[Bug-gsl] [bug #43256] gsl_sf_coupling_6j overflows

[Anders Søndergaard]

Follow-up Comment #1, bug #43256 (project gsl):

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;
  }
}




    _______________________________________________________

Reply to this item at:

  <http://savannah.gnu.org/bugs/?43256>

_______________________________________________
  Message sent via/by Savannah
  http://savannah.gnu.org/