Re: [Bug-gsl] gsl_sf_coupling_3j bug report

[Grigory I. Rubtsov]

Hi Alexey,

You gave me an idea of further improvement. If one mupltiply the norm
with every term, there will be no large numbers at all. Please
consider the following patch for inclusion into GSL. It works with
practically arbitrary large l.

Best wishes,
Grigory Rubtsov

--- gsl-1.15_orig/specfunc/coupling.c   2010-12-26 20:57:08.000000000 +0300
+++ gsl-1.15/specfunc/coupling.c        2011-10-02 13:19:39.000000000 +0400
@@ -136,38 +136,33 @@
         kmax = locMin3 (jcc, jmma, jpmb),
         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;
+    double sum_pos = 0.0, sum_neg = 0.0, lnorm;
+    gsl_sf_result bc1, bc2, bc3, bcn1, bcn2, bcd1, bcd2, bcd3, bcd4, term;

-    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));
+    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);
+
+    lnorm = 0.5 * (bcn1.val + bcn2.val - bcd1.val - bcd2.val - bcd3.val
+                  - bcd4.val - log((double) two_jc + 1));

     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);
+      status += gsl_sf_lnchoose_e (jcc, k, &bc1);
+      status += gsl_sf_lnchoose_e (jcb, jmma - k, &bc2);
+      status += gsl_sf_lnchoose_e (jca, jpmb - k, &bc3);
+      status += gsl_sf_exp_e(bc1.val + bc2.val + bc3.val + lnorm, &term);

       if (status != 0) {
         OVERFLOW_ERROR (result);
-      }
-
-      term = bc1.val * bc2.val * bc3.val;
+      }

       if (sign < 0) {
-        sum_neg += norm * term;
+        sum_neg += term.val;
       } else {
-        sum_pos += norm * term;
+        sum_pos += term.val;
       }

       sign = -sign;



2011/10/2 Alexey A. Illarionov <address@hidden>:
> 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
>