[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Bug-gsl] gsl_sf_coupling_3j bug report
From: |
Grigory I. Rubtsov |
Subject: |
[Bug-gsl] gsl_sf_coupling_3j bug report |
Date: |
Thu, 29 Sep 2011 17:09:15 +0400 |
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] gsl_sf_coupling_3j bug report,
Grigory I. Rubtsov <=
- 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, 2011/10/01
- 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