Re: [Bug-gsl] gsl_sf_coupling_3j bug report

[Alexey A. Illarionov]

Hi Grigory,

Though it is possible to avoid overflow during calculation this way the
algorithm is still not working for large l-s. This is just a feature of
floating points numbers, and there is nothing can be done unless you
switch to a different method or BigIntegers. You can see this by
calculating
(200 200 200)
(-10 60 -50)

Best regards.

On 02/10/11 05:25 AM, Grigory I. Rubtsov wrote:
> 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
>>