bugGNU Scientific Library - Bugs: bug #43256, gsl_sf_coupling_6j overflows

 
 

You are not allowed to post comments on this tracker with your current authentication level.

bug #43256: gsl_sf_coupling_6j overflows

Submitted by:  Anders Søndergaard <andersas>
Submitted on:  Fri 19 Sep 2014 12:45:52 PM UTC  
 
Category: Runtime errorSeverity: 3 - Normal
Operating System: Ubuntu linuxStatus: None
Assigned to: NoneOpen/Closed: Open
Release: 1.16 (ubuntu latest)

Mon 22 Sep 2014 11:54:15 AM UTC, comment #2:

Added a file containing the added and modified function because
the markup failed in my previous comment.

Anders

(file #32162)

Anders Søndergaard <andersas>
Mon 22 Sep 2014 11:52:56 AM UTC, comment #1:

I implemented the suggested changes of working with the logarithm of the factorials instead of the factorials themselves.

This seems to work fine.

I still need to work out the error calculation and to verify that the changes actually work (by e.g. comparing with sympy).

Also, perhaps it is faster not to work with the logarithm for small j. Who knows? Maybe the sixj function should
jump to the logarithm version in case of overflow.

You can have the code that works so far here:

In the specfunc/coupling.c file I changed the delta function to read:

static
int
lndelta(int ta, int tb, int tc, gsl_sf_result * d)
{
gsl_sf_result f1, f2, f3, f4;
int status = 0;
status += gsl_sf_lnfact_e((ta + tb - tc)/2, &f1);
status += gsl_sf_lnfact_e((ta + tc - tb)/2, &f2);
status += gsl_sf_lnfact_e((tb + tc - ta)/2, &f3);
status += gsl_sf_lnfact_e((ta + tb + tc)/2 + 1, &f4);
if(status != 0) {
OVERFLOW_ERROR(d);
}
d->val = f1.val + f2.val + f3.val - f4.val;
d->err = 4.0 * GSL_DBL_EPSILON * fabs(d->val); // Is this correct?
return GSL_SUCCESS;
}

And the gsl_sf_coupling_6j_e function:

int
gsl_sf_coupling_6j_e(int two_ja, int two_jb, int two_jc,
int two_jd, int two_je, int two_jf,
gsl_sf_result * result)
{
/* CHECK_POINTER(result) */

if( two_ja < 0 || two_jb < 0 || two_jc < 0
|| two_jd < 0 || two_je < 0 || two_jf < 0
) {
DOMAIN_ERROR(result);
}
else if( triangle_selection_fails(two_ja, two_jb, two_jc)
|| triangle_selection_fails(two_ja, two_je, two_jf)
|| triangle_selection_fails(two_jb, two_jd, two_jf)
|| triangle_selection_fails(two_je, two_jd, two_jc)
) {
result->val = 0.0;
result->err = 0.0;
return GSL_SUCCESS;
}
else {
gsl_sf_result n1;
gsl_sf_result d1, d2, d3, d4, d5, d6;
double lnnorm;
int tk, tkmin, tkmax;
double phase;
double sum_pos = 0.0;
double sum_neg = 0.0;
double sumsq_err = 0.0;
int status = 0;
status += lndelta(two_ja, two_jb, two_jc, &d1);
status += lndelta(two_ja, two_je, two_jf, &d2);
status += lndelta(two_jb, two_jd, two_jf, &d3);
status += lndelta(two_je, two_jd, two_jc, &d4);
if(status != GSL_SUCCESS) {
OVERFLOW_ERROR(result);
}
// Dividing by two inside exp() is the same as sqrt() outside
lnnorm = (d1.val + d2.val + d3.val + d4.val)/2;

tkmin = locMax3(0,
two_ja + two_jd - two_jc - two_jf,
two_jb + two_je - two_jc - two_jf);

tkmax = locMin5(two_ja + two_jb + two_je + two_jd + 2,
two_ja + two_jb - two_jc,
two_je + two_jd - two_jc,
two_ja + two_je - two_jf,
two_jb + two_jd - two_jf);

phase = GSL_IS_ODD((two_ja + two_jb + two_je + two_jd + tkmin)/2)
? -1.0
: 1.0;

for(tk=tkmin; tk<=tkmax; tk += 2) {
double term;
double term_err;
gsl_sf_result den_1, den_2;
gsl_sf_result d1_a, d1_b;
status = 0;

status += gsl_sf_lnfact_e((two_ja + two_jb + two_je + two_jd - tk)/2 + 1, &n1);
status += gsl_sf_lnfact_e(tk/2, &d1_a);
status += gsl_sf_lnfact_e((two_jc + two_jf - two_ja - two_jd + tk)/2, &d1_b);
status += gsl_sf_lnfact_e((two_jc + two_jf - two_jb - two_je + tk)/2, &d2);
status += gsl_sf_lnfact_e((two_ja + two_jb - two_jc - tk)/2, &d3);
status += gsl_sf_lnfact_e((two_je + two_jd - two_jc - tk)/2, &d4);
status += gsl_sf_lnfact_e((two_ja + two_je - two_jf - tk)/2, &d5);
status += gsl_sf_lnfact_e((two_jb + two_jd - two_jf - tk)/2, &d6);

if(status != GSL_SUCCESS) {
OVERFLOW_ERROR(result);
}

d1.val = d1_a.val + d1_b.val;
d1.err = d1_a.err * fabs(d1_b.val) + fabs(d1_a.val) * d1_b.err; // ???

den_1.val = d1.val+d2.val+d3.val;
den_1.err = d1.err * fabs(d2.val+d3.val); // ???
den_1.err += d2.err * fabs(d1.val+d3.val); // ???
den_1.err += d3.err * fabs(d1.val+d2.val); // ???

den_2.val = d4.val+d5.val+d6.val;
den_2.err = d4.err * fabs(d5.val+d6.val); // ???
den_2.err += d5.err * fabs(d4.val+d6.val); // ???
den_2.err += d6.err * fabs(d4.val+d5.val); // ???

term = n1.val - den_1.val - den_2.val;
term_err = n1.err / fabs(den_1.val) / fabs(den_2.val); // ???
term_err += fabs(term / den_1.val) * den_1.err; // ???
term_err += fabs(term / den_2.val) * den_2.err; // ???

if(phase >= 0.0) {
sum_pos += phase*exp(lnnorm+term);
}
else {
sum_neg -= phase*exp(lnnorm+term);
}

phase = -phase;
sumsq_err += exp(lnnorm+lnnorm) + term_err*term_err; // ???
}

result->val = sum_pos - sum_neg;
result->err = 2.0 * GSL_DBL_EPSILON * (sum_pos + sum_neg); // ???
result->err += sqrt(sumsq_err / (0.5*(tkmax-tkmin)+1.0)); // ???
result->err += 2.0 * GSL_DBL_EPSILON * (tkmax - tkmin + 2.0) * fabs(result->val); // ???

return GSL_SUCCESS;
}
}

Anders Søndergaard <andersas>
Fri 19 Sep 2014 12:45:52 PM UTC, original submission:

I need to calculate 6-j symbols for somewhat large J.
For example the symbol:

(14, 16, 16
78, 62, 76) = 6/1739*sqrt(16971005238954/1556382731177197) = 0.000360286

but this fails because gsl_sf_fact overflows.
It is used internally in gsl_sf_coupling_6j.

I have looked a bit at the code, and I think in these large
cases it is possible to work with the logarithm of the factorials
and in the end exponentiate.

Other libraries I have tested also gets these symbols wrong (but silently returns the wrong value),
except pythons sympy, which calculates the symbols as the square root of a rational number.

I would rather not interface my Fortran program with Python.

Anders Søndergaard <andersas>

 

Attached Files
file #32162:  sixjsymbols.c added by andersas (4KiB - text/x-csrc)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -unavailable- added by andersas (Submitted the item)
  •  

    Please enter the title of George Orwell's famous dystopian book (it's a date):

     

     

    Follows 1 latest change.

    Date Changed By Updated Field Previous Value => Replaced By
    Mon 22 Sep 2014 11:54:15 AM UTCandersasAttached File-=>Added sixjsymbols.c, #32162

    Back to the top


    Powered by Savane 3.1-cleanup1