gsl_sf_coupling_3j, BUG #7

gsl_sf_coupling_3j, BUG #7

[G.W.M. Vissers]

Hi,

I saw at documents.wolfram.com/v4/RefGuide/ThreeJSymbol.html that 
mathematica's 3j symbol is built upon binomial coefficients. I implemented 
their expression in the following routine:

int
gsl_sf_coupling_3j_e_new (int two_ja, int two_jb, int two_jc,
                          int two_ma, int two_mb, int two_mc,
                          gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */

  if(two_ja < 0 || two_jb < 0 || two_jc < 0) {
    DOMAIN_ERROR(result);
  }
  else if (   triangle_selection_fails(two_ja, two_jb, two_jc)
           || m_selection_fails(two_ja, two_jb, two_jc, two_ma, two_mb, two_mc)
     ) {
    result->val = 0.0;
    result->err = 0.0;
    return GSL_SUCCESS;
  }
  else {
    int jca  = (-two_ja + two_jb + two_jc) / 2,
        jcb  = ( two_ja - two_jb + two_jc) / 2,
        jcc  = ( two_ja + two_jb - two_jc) / 2,
        jmma = ( two_ja - two_ma) / 2,
        jmmb = ( two_jb - two_mb) / 2,
        jmmc = ( two_jc - two_mc) / 2,
        jpma = ( two_ja + two_ma) / 2,
        jpmb = ( two_jb + two_mb) / 2,
        jpmc = ( two_jc + two_mc) / 2,
        jsum = ( two_ja + two_jb + two_jc) / 2,
        kmin = locMax3 (0, jpmb - jmmc, jmma - jpmc),
        kmax = locMin3 (jcc, jmma, jpmb),
        k, sign = GSL_IS_ODD (kmin - jpma + jmmb) ? -1 : 1,
        status = 0;
    double ksump = 0.0, ksumn = 0.0, norm, term;
    gsl_sf_result bc1, bc2, bc3, bcn1, bcn2, bcd1, bcd2, bcd3, bcd4;

    status += gsl_sf_choose_e (two_ja, jcc , &bcn1);
    status += gsl_sf_choose_e (two_jb, jcc , &bcn2);
    status += gsl_sf_choose_e (jsum+1, jcc , &bcd1);
    status += gsl_sf_choose_e (two_ja, jmma, &bcd2);
    status += gsl_sf_choose_e (two_jb, jmmb, &bcd3);
    status += gsl_sf_choose_e (two_jc, jpmc, &bcd4);
    
    if (status != 0) {
      OVERFLOW_ERROR (result);
    }
    
    norm = sqrt (bcn1.val * bcn2.val)
           / sqrt (bcd1.val * bcd2.val * bcd3.val * bcd4.val * ((double) two_jc + 1.0));

    for (k = kmin; k <= kmax; k++) {
      status += gsl_sf_choose_e (jcc, k, &bc1);
      status += gsl_sf_choose_e (jcb, jmma - k, &bc2);
      status += gsl_sf_choose_e (jca, jpmb - k, &bc3);
      
      if (status != 0) {
        OVERFLOW_ERROR (result);
      }
      
      term = bc1.val * bc2.val * bc3.val;
      
      if (sign < 0) {
        ksumn += norm * term;
      } else {
        ksump += norm * term;
      }
      
      sign = -sign;
    }
    
    result->val  = ksump - ksumn;
    result->err  = 2.0 * GSL_DBL_EPSILON * (ksump + ksumn);
    result->err += 2.0 * GSL_DBL_EPSILON * (kmax - kmin) * fabs(result->val);

    return GSL_SUCCESS;
  }
}

This seems to be a lot more stable than the previous routine (to use the 
example in the BUGS file:
 
gsl_sf_coupling_3j_new (80,80,80,0,0,0) 

now produces 1.4968524483265e-02 with an estimated error of 
1.3934261781184e-10. The actual value is 1.4968524489706e-02 which is 
well within the error bars.)

The routine passes all 6 tests in test_sf.c. Please let me know if 
anything is wrong with it, 

Cheers,

Gé Vissers

Re: gsl_sf_coupling_3j, BUG #7

[Brian Gough]

G.W.M. Vissers writes:
 > I saw at documents.wolfram.com/v4/RefGuide/ThreeJSymbol.html that 
 > mathematica's 3j symbol is built upon binomial coefficients. I implemented 
 > their expression in the following routine:
 > 
 > int
 > gsl_sf_coupling_3j_e_new (int two_ja, int two_jb, int two_jc,
 >                           int two_ma, int two_mb, int two_mc,
 >                           gsl_sf_result * result)

Thanks.  

It looks good to me but maybe the original author of the gsl
3j functions has some comments.

Brian

p.s. Did you look at the papers mentioned in the BUGS file?

http://www.ks.uiuc.edu/Publications/Papers/abstract.cgi?tbcode=SCHU76B
http://www.ks.uiuc.edu/Publications/Papers/abstract.cgi?tbcode=SCHU75A
http://www.ks.uiuc.edu/Publications/Papers/abstract.cgi?tbcode=SCHU75

Re: gsl_sf_coupling_3j, BUG #7

[G.W.M. Vissers]

>  > I saw at documents.wolfram.com/v4/RefGuide/ThreeJSymbol.html that 
>  > mathematica's 3j symbol is built upon binomial coefficients. I implemented 
>  > their expression in the following routine:
>  > 
>  > int
>  > gsl_sf_coupling_3j_e_new (int two_ja, int two_jb, int two_jc,
>  >                           int two_ma, int two_mb, int two_mc,
>  >                           gsl_sf_result * result)
> 
> Thanks.  
> 
> It looks good to me but maybe the original author of the gsl
> 3j functions has some comments.
> 
> Brian
> 
> p.s. Did you look at the papers mentioned in the BUGS file?
> 
> http://www.ks.uiuc.edu/Publications/Papers/abstract.cgi?tbcode=SCHU76B
> http://www.ks.uiuc.edu/Publications/Papers/abstract.cgi?tbcode=SCHU75A
> http://www.ks.uiuc.edu/Publications/Papers/abstract.cgi?tbcode=SCHU75

I have to admit that although this function does better than the original 
one, it is still not as stable as one would like it to be. I have printed 
the papers mentioned, and will take a better look at them tonight. I 
expect that for large j the recurrence relations in there may be worthwile 
to implement. Maybe later this week ;-) 

Gé

Re: gsl_sf_coupling_3j, BUG #7

[Brian Gough]

G.W.M. Vissers writes:
 > 
 > >  > I saw at documents.wolfram.com/v4/RefGuide/ThreeJSymbol.html that 
 > >  > mathematica's 3j symbol is built upon binomial coefficients. I implemented 
 > >  > their expression in the following routine:
 > >  > 
 > >  > int
 > >  > gsl_sf_coupling_3j_e_new (int two_ja, int two_jb, int two_jc,
 > >  >                           int two_ma, int two_mb, int two_mc,
 > >  >                           gsl_sf_result * result)

Thanks, I've added the new 3j routine to the next release.

-- 
Brian