This is the mail archive of the
gsl-discuss@sources.redhat.com
mailing list for the GSL project.
gsl_sf_coupling_3j, BUG #7
- From: "G.W.M. Vissers" <vissers at theochem dot kun dot nl>
- To: gsl-discuss at sources dot redhat dot com
- Date: Wed, 23 Apr 2003 14:41:01 +0200 (CEST)
- Subject: gsl_sf_coupling_3j, BUG #7
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