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

Submitter:  Anders Søndergaard <andersas>
Submitted:  Fri 19 Sep 2014 12:45:52 PM UTC
   
 
Category:  Runtime error Severity:  3 - Normal
Operating System:  Ubuntu linux Status:  None
Assigned to:  None Open/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>

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

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
  • -email is unavailable- added by andersas (Submitted the item)
  •  

    Follows 1 latest change.

    Date Changed by Updated Field Previous Value => Replaced by
    2014-09-22 andersas Attached File- Added sixjsymbols.c, #32162

    Back to the top

    Powered by Savane 3.13-f8d8.
    Corresponding source code