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
> [email protected]
> https://lists.gnu.org/mailman/listinfo/bug-gsl
_______________________________________________
Bug-gsl mailing list
[email protected]
https://lists.gnu.org/mailman/listinfo/bug-gsl
