gsl_sf_coupling_3j, BUG #7

["G.W.M. Vissers" <vissers@theochem.kun.nl>]

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