[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: [Bug-gsl] gsl_sf_coupling_3j bug report
From: |
Grigory I. Rubtsov |
Subject: |
Re: [Bug-gsl] gsl_sf_coupling_3j bug report |
Date: |
Sun, 2 Oct 2011 02:07:07 +0400 |
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] gsl_sf_coupling_3j bug report, Grigory I. Rubtsov, 2011/10/01
- Re: [Bug-gsl] gsl_sf_coupling_3j bug report, Alexey A Illarionov, 2011/10/02
- Re: [Bug-gsl] gsl_sf_coupling_3j bug report,
Grigory I. Rubtsov <=
- Re: [Bug-gsl] gsl_sf_coupling_3j bug report, Alexey A. Illarionov, 2011/10/01
- Re: [Bug-gsl] gsl_sf_coupling_3j bug report, Grigory I. Rubtsov, 2011/10/02
- Re: [Bug-gsl] gsl_sf_coupling_3j bug report, Alexey A. Illarionov, 2011/10/02
- Re: [Bug-gsl] gsl_sf_coupling_3j bug report, Brian Gough, 2011/10/10
- Re: [Bug-gsl] gsl_sf_coupling_3j bug report, Alexey A. Illarionov, 2011/10/12
- Re: [Bug-gsl] gsl_sf_coupling_3j bug report, Brian Gough, 2011/10/12
Re: [Bug-gsl] gsl_sf_coupling_3j bug report, Brian Gough, 2011/10/07