Re: [Bug-gsl] gsl_sf_coupling_3j bug report

[Alexey A. Illarionov]

Hi Grigory,

Good job, but I think it should be UNDERFLOW_ERROR instead of
OVERFLOW_ERROR. I did not look closely at this part of the code
previously. I usually avoid underflow by computation of ln(norm) instead
of norm, for example here is how I usually compute this part. (I'm not
quite sure what method is better).

=== modified file 'specfunc/coupling.c'
--- specfunc/coupling.c 2011-09-20 13:52:43 +0000
+++ specfunc/coupling.c 2011-10-02 00:08:27 +0000
@@ -146,21 +146,22 @@
         k, sign = GSL_IS_ODD (kmin - jpma + jmmb) ? -1 : 1,
         status = 0;
     double sum_pos = 0.0, sum_neg = 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);
+    gsl_sf_result bc1, bc2, bc3, bcn1, bcn2, bcd1, bcd2, bcd3, bcd4,
elnorm;
+
+    status += gsl_sf_lnchoose_e (two_ja, jcc , &bcn1);
+    status += gsl_sf_lnchoose_e (two_jb, jcc , &bcn2);
+    status += gsl_sf_lnchoose_e (jsum+1, jcc , &bcd1);
+    status += gsl_sf_lnchoose_e (two_ja, jmma, &bcd2);
+    status += gsl_sf_lnchoose_e (two_jb, jmmb, &bcd3);
+    status += gsl_sf_lnchoose_e (two_jc, jpmc, &bcd4);
+
+    status += gsl_sf_exp_e( 0.5 * (bcn1.val + bcn2.val - bcd1.val -
bcd2.val - bcd3.val - bcd4.val),
+      &elnorm);
+    norm = elnorm.val / sqrt((double) two_jc + 1.);
+
+    if (norm == 0.) {
+      UNDERFLOW_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);


P.S. I don't think long double will help you, since as far as I know
it's support is highly dependant on compiler and architecture. If you
really need large j's with high precision I think it's better to write
something with GMP or similar library.

On 01/10/11 06:07 PM, Grigory I. Rubtsov wrote:
> Hi, Alexey,
> 
> Thank you for the hint. In fact there is OVERFLOW_ERROR for larger l
> (>500). I found the origin of my bug: in my case norm (coupling.c line
> 153) is equal to 0 because of large denominator. I suggest the
> following patch which extends the range of the function from ~250 to
> ~500.
> 
> --- coupling_prev.c     2011-10-02 01:42:52.000000000 +0400
> +++ coupling.c  2011-10-02 01:43:32.000000000 +0400
> @@ -151,7 +151,10 @@
>      }
> 
>      norm = sqrt (bcn1.val * bcn2.val)
> -           / sqrt (bcd1.val * bcd2.val * bcd3.val * bcd4.val *
> ((double) two_jc + 1.0));
> +      / sqrt 
> (bcd1.val)/sqrt(bcd2.val)/sqrt(bcd3.val)/sqrt(bcd4.val)/sqrt((double)
> two_jc + 1.0);
> +    if(norm==0) {
> +      OVERFLOW_ERROR (result);
> +    }
> 
>      for (k = kmin; k <= kmax; k++) {
>        status += gsl_sf_choose_e (jcc, k, &bc1);
> 
> Is there an easy way to switch GSL from double to long double?
> 
> Best wishes,
> Grigory Rubtsov
> 
> 
> 2011/10/2 Alexey A Illarionov <address@hidden>:
>> Hi,
>>
>> As far as I know it is a feature of the imlemented algorithm. Due to
>> cancellation in the sum at coupling.c:164-184, it is bad for large
>> arguments. Although I'm curious myself, why you did not get OVERFLOW_ERROR.
>>
>> And you can easily write subroutine yourself by using
>> 1. The same algorithm (formula Racah -- eq 2. in [Wei1999]) but using
>> BigInt instead of int (see GMP library)
>> 2. Use some exotic algorithms, for example [Wei1999].
>>
>> --
>> Best regards, Alexey A. Illarionov
>>
>> References:
>> Wei, Computer Physics Communications 120 (1999) 222-230.
>>
>> On 29/09/11 09:09 AM, Grigory I. Rubtsov wrote:
>>> Dear GSL developers,
>>>
>>> Thank you for the great effort in developing easy to use and fast
>>> scientific library. I am using GSL in most of scientific calculations.
>>>
>>> Recently I encounter a bug in calculation of 3j symbol for relatively large 
>>> l.
>>> In particular:
>>> gsl_sf_coupling_3j_e(2*l1,2*l2,2*l3,2*m1,2*m2,2*m3,&r)
>>>   for l1=249, l2=248, l3=2, m1=5, m2=-6, m3=1
>>>   gives 0 and error 0 (correct answer should be -0.022878)
>>> while
>>>   for l1=248, l2=247, l3=2, m1=5, m2=-6, m3=1
>>>   gives -0.0229267 and error 2.98369e-17 (the result is correct)
>>>
>>> The calculation of 3j symbols for l up to 1000 is important for the
>>> analysis of cosmic microwave background anisotropy data, especially
>>> WMAP and Planck missions.
>>>
>>> Bug details:
>>>     * GSL version: 1.15 on SL 5.7 (64-bit) at Pentium E5400  @ 2.70GHz
>>>     * The compiler: gcc version 4.1.2 20080704 (Red Hat 4.1.2-51)
>>>     * Compiler options: g++  -lm -lgsl -lgslcblas  gsl_bug.cpp   -o gsl_bug
>>>     * A short program which exercises the bug
>>>
>>> #include <gsl/gsl_sf_coupling.h>
>>> #include <stdio.h>
>>> #include <math.h>
>>>
>>> int main() {
>>>     gsl_sf_result r;
>>>     int j=248,m=5;
>>>     int l1=j+1, l2=j, l3=2;
>>>     int m1=m, m2=-m-1, m3=1;
>>>     double jj = double(j), mm=double(m);
>>>     double ans=-2*(jj+2*mm+2)*sqrt(
>>> (jj-mm+1)*(jj-mm)/2/jj/(2*jj+1)/(2*jj+2)/(2*jj+3)/(2*jj+4) );
>>>     gsl_sf_coupling_3j_e(2*l1,2*l2,2*l3,2*m1,2*m2,2*m3,&r);
>>>     printf("(%i %i %i) \n(%i %i %i) = %g %g\n", l1,l2,l3,m1,m2,m3,r.val, 
>>> r.err);
>>>     printf("Expected: %g\n", ans);
>>> }
>>>
>>> Sincerely yours,
>>> Grigory Rubtsov,
>>> Institute for Nuclear Research of the
>>> Russian Academy of Sciences
>>>
>>> _______________________________________________
>>> Bug-gsl mailing list
>>> address@hidden
>>> https://lists.gnu.org/mailman/listinfo/bug-gsl
>>
>>
>>
> 
> _______________________________________________
> Bug-gsl mailing list
> address@hidden
> https://lists.gnu.org/mailman/listinfo/bug-gsl