[Bug-gsl] Bug in the GSL 1.5 Legendre Functions

[MAGNEVILLE Christophe]

Bug report: Using GSL 1.5

       Hello,

 I have found a bug in the functions:

    "Legendre Functions and Spherical Harmonics"

    gsl_sf_legendre_sphPlm_array(Lmax,m,cx,Xlm);
    gsl_sf_legendre_sphPlm(l,m,cx);

 When I check the precision of various routines
distribued on the WEB, I saw that the GSL Legendre Functions
gave big differences with the standard reference (HealPix, Sophya)
for some region in cos(teta),m,l:

 Here is the result of the test I have made (lmax=2000):

The test is done for cos(teta) ranging from -1 to +1 (step 0.04)
As you may see, there are differences around cos(teta)=+/- 0.92
That happens for big value of l and m

 I think that the GSL answer is wrong because it returns 0
and the Sophya library returns non zero:

 I reproduce that bug on:
Alpha/OSF1 V5.1 2650 Compiler Driver V6.3-008 (cxx) cxx Driver
Linux 2.4.18-14 gcc 3.2 20020903 (Red Hat Linux 8.0 3.2-7)
  (but it may exist on other systems)

 In the list below, when the diff is greater than 0.1,
I have printed Xlm= (GSL value) and Xlm1= (Sophya value)

  ==> the GSL return Xlm=0 when Sophya return a non zero reasonable value.

  (Also look at my second test about continuity after that output)

------------------------------------------------------
Lmax=2000
>>>>>>>> Try: 0 cx=-1 (th=3.14159)  n_echec=0 diffmax=2.40945e-11 (l=1913,m=0)
>>>>>>>> Try: 1 cx=-0.96 (th=2.8578)  n_echec=1030 diffmax=0.00229663 
>>>>>>>> (l=2000,m=586)
m=796 l=1997 cx=-0.92 Xlm=-0 Xlm1=-0.100142 diff=0.100142 approx=-0.14598       
      ***
m=796 l=1998 cx=-0.92 Xlm=0 Xlm1=0.108967 diff=0.108967 approx=0.185656         
      ***
m=796 l=1999 cx=-0.92 Xlm=0 Xlm1=-0.118442 diff=0.118442 approx=-0.191241       
      ***
m=796 l=2000 cx=-0.92 Xlm=-0 Xlm1=0.128604 diff=0.128604 approx=0.161718        
      ***
m=797 l=2000 cx=-0.92 Xlm=-0 Xlm1=0.104166 diff=0.104166 approx=-0.0813683      
      ***
>>>>>>>> Try: 2 cx=-0.92 (th=2.73888)  n_echec=2638 diffmax=0.128604 
>>>>>>>> (l=2000,m=796)   ***
>>>>>>>> Try: 3 cx=-0.88 (th=2.64666)  n_echec=0 diffmax=8.86466e-07 
>>>>>>>> (l=2000,m=1002)
>>>>>>>> Try: 4 cx=-0.84 (th=2.56808)  n_echec=0 diffmax=9.35252e-13 
>>>>>>>> (l=1373,m=739)
>>>>>>>> Try: 5 cx=-0.8 (th=2.49809)  n_echec=0 diffmax=1.02962e-12 
>>>>>>>> (l=1588,m=946)
>>>>>>>> Try: 6 cx=-0.76 (th=2.43411)  n_echec=0 diffmax=1.18283e-12 
>>>>>>>> (l=1805,m=1166)
>>>>>>>> Try: 7 cx=-0.72 (th=2.3746)  n_echec=0 diffmax=1.67688e-12 
>>>>>>>> (l=1942,m=1341)
>>>>>>>> Try: 8 cx=-0.68 (th=2.31856)  n_echec=0 diffmax=1.6005e-12 
>>>>>>>> (l=1995,m=1456)
>>>>>>>> Try: 9 cx=-0.64 (th=2.26529)  n_echec=0 diffmax=2.13363e-12 
>>>>>>>> (l=1937,m=1482)
>>>>>>>> Try: 10 cx=-0.6 (th=2.2143)  n_echec=0 diffmax=2.085e-12 
>>>>>>>> (l=1860,m=1482)
>>>>>>>> Try: 11 cx=-0.56 (th=2.16518)  n_echec=0 diffmax=2.26485e-12 
>>>>>>>> (l=1796,m=1482)
>>>>>>>> Try: 12 cx=-0.52 (th=2.11765)  n_echec=0 diffmax=2.1847e-12 
>>>>>>>> (l=1742,m=1482)
>>>>>>>> Try: 13 cx=-0.48 (th=2.07145)  n_echec=0 diffmax=2.19358e-12 
>>>>>>>> (l=1695,m=1482)
>>>>>>>> Try: 14 cx=-0.44 (th=2.0264)  n_echec=0 diffmax=2.07234e-12 
>>>>>>>> (l=1656,m=1482)
>>>>>>>> Try: 15 cx=-0.4 (th=1.98231)  n_echec=0 diffmax=2.1676e-12 
>>>>>>>> (l=1622,m=1482)
>>>>>>>> Try: 16 cx=-0.36 (th=1.93906)  n_echec=0 diffmax=2.21356e-12 
>>>>>>>> (l=1593,m=1482)
>>>>>>>> Try: 17 cx=-0.32 (th=1.89653)  n_echec=0 diffmax=2.15872e-12 
>>>>>>>> (l=1568,m=1482)
>>>>>>>> Try: 18 cx=-0.28 (th=1.85459)  n_echec=0 diffmax=2.22244e-12 
>>>>>>>> (l=1547,m=1482)
>>>>>>>> Try: 19 cx=-0.24 (th=1.81316)  n_echec=0 diffmax=2.29949e-12 
>>>>>>>> (l=1530,m=1482)
>>>>>>>> Try: 20 cx=-0.2 (th=1.77215)  n_echec=0 diffmax=2.39142e-12 
>>>>>>>> (l=1515,m=1482)
>>>>>>>> Try: 21 cx=-0.16 (th=1.73149)  n_echec=0 diffmax=2.4325e-12 
>>>>>>>> (l=1504,m=1482)
>>>>>>>> Try: 22 cx=-0.12 (th=1.69109)  n_echec=0 diffmax=2.40941e-12 
>>>>>>>> (l=1495,m=1482)
>>>>>>>> Try: 23 cx=-0.08 (th=1.65088)  n_echec=0 diffmax=2.63967e-12 
>>>>>>>> (l=1488,m=1482)
>>>>>>>> Try: 24 cx=-0.04 (th=1.61081)  n_echec=0 diffmax=2.98939e-12 
>>>>>>>> (l=1484,m=1482)
>>>>>>>> Try: 25 cx=0 (th=1.5708)  n_echec=0 diffmax=3.61888e-12 (l=1482,m=1482)
>>>>>>>> Try: 26 cx=0.04 (th=1.53079)  n_echec=0 diffmax=2.97806e-12 
>>>>>>>> (l=1484,m=1482)
>>>>>>>> Try: 27 cx=0.08 (th=1.49071)  n_echec=0 diffmax=2.61791e-12 
>>>>>>>> (l=1488,m=1482)
>>>>>>>> Try: 28 cx=0.12 (th=1.45051)  n_echec=0 diffmax=2.38121e-12 
>>>>>>>> (l=1495,m=1482)
>>>>>>>> Try: 29 cx=0.16 (th=1.41011)  n_echec=0 diffmax=2.45759e-12 
>>>>>>>> (l=1504,m=1482)
>>>>>>>> Try: 30 cx=0.2 (th=1.36944)  n_echec=0 diffmax=2.43205e-12 
>>>>>>>> (l=1515,m=1482)
>>>>>>>> Try: 31 cx=0.24 (th=1.32843)  n_echec=0 diffmax=2.25198e-12 
>>>>>>>> (l=1530,m=1482)
>>>>>>>> Try: 32 cx=0.28 (th=1.287)  n_echec=0 diffmax=2.27773e-12 
>>>>>>>> (l=1547,m=1482)
>>>>>>>> Try: 33 cx=0.32 (th=1.24507)  n_echec=0 diffmax=2.2462e-12 
>>>>>>>> (l=1568,m=1482)
>>>>>>>> Try: 34 cx=0.36 (th=1.20253)  n_echec=0 diffmax=2.12541e-12 
>>>>>>>> (l=1593,m=1482)
>>>>>>>> Try: 35 cx=0.4 (th=1.15928)  n_echec=0 diffmax=2.07634e-12 
>>>>>>>> (l=1622,m=1482)
>>>>>>>> Try: 36 cx=0.44 (th=1.1152)  n_echec=0 diffmax=2.01839e-12 
>>>>>>>> (l=1656,m=1482)
>>>>>>>> Try: 37 cx=0.48 (th=1.07014)  n_echec=0 diffmax=2.08433e-12 
>>>>>>>> (l=1695,m=1482)
>>>>>>>> Try: 38 cx=0.52 (th=1.02395)  n_echec=0 diffmax=2.1847e-12 
>>>>>>>> (l=1742,m=1482)
>>>>>>>> Try: 39 cx=0.56 (th=0.976411)  n_echec=0 diffmax=2.31903e-12 
>>>>>>>> (l=1796,m=1482)
>>>>>>>> Try: 40 cx=0.6 (th=0.927295)  n_echec=0 diffmax=2.1867e-12 
>>>>>>>> (l=1860,m=1482)
>>>>>>>> Try: 41 cx=0.64 (th=0.876298)  n_echec=0 diffmax=2.28262e-12 
>>>>>>>> (l=1937,m=1482)
>>>>>>>> Try: 42 cx=0.68 (th=0.823034)  n_echec=0 diffmax=1.82254e-12 
>>>>>>>> (l=1995,m=1456)
>>>>>>>> Try: 43 cx=0.72 (th=0.766994)  n_echec=0 diffmax=1.69109e-12 
>>>>>>>> (l=1942,m=1341)
>>>>>>>> Try: 44 cx=0.76 (th=0.707483)  n_echec=0 diffmax=1.18283e-12 
>>>>>>>> (l=1805,m=1166)
>>>>>>>> Try: 45 cx=0.8 (th=0.643501)  n_echec=0 diffmax=1.15041e-12 
>>>>>>>> (l=1955,m=1166)
>>>>>>>> Try: 46 cx=0.84 (th=0.573513)  n_echec=0 diffmax=1.0596e-12 
>>>>>>>> (l=1658,m=893)
>>>>>>>> Try: 47 cx=0.88 (th=0.494934)  n_echec=0 diffmax=8.86466e-07 
>>>>>>>> (l=2000,m=1002)
m=796 l=1997 cx=0.92 Xlm=0 Xlm1=0.100142 diff=0.100142 approx=1.28935           
     ***
m=796 l=1998 cx=0.92 Xlm=0 Xlm1=0.108967 diff=0.108967 approx=1.28019           
     ***
m=796 l=1999 cx=0.92 Xlm=0 Xlm1=0.118442 diff=0.118442 approx=1.26662           
     ***
m=796 l=2000 cx=0.92 Xlm=0 Xlm1=0.128604 diff=0.128604 approx=1.24846           
     ***
m=797 l=2000 cx=0.92 Xlm=0 Xlm1=-0.104166 diff=0.104166 approx=1.2855           
     ***
>>>>>>>> Try: 48 cx=0.92 (th=0.402716)  n_echec=2638 diffmax=0.128604 
>>>>>>>> (l=2000,m=796) ***
>>>>>>>> Try: 49 cx=0.96 (th=0.283794)  n_echec=1030 diffmax=0.00229663 
>>>>>>>> (l=2000,m=586)
>>>>>>>> Try: 50 cx=1 (th=0)  n_echec=0 diffmax=2.40945e-11 (l=1913,m=0)
Nombre de differences>1e-06 = 7336, diff_maxi=0.128604
---------------------------------------------------------------------------------------


---------------------------------------------------------------------------------------
=========> Another test, I have made, may help you:
That was done to test the continuity of the legendre functions:

int main (void) {
 double cx; int lmax=2000, l=1997, m=796;
 double Xlm[2005];
 for(cx=-1;cx<-0.91;cx+=0.001) {
   gsl_sf_legendre_sphPlm_array(lmax,m,cx,Xlm);
   double v = Xlm[l-m];
   double lv = 0.;
   if(v>0.) lv = log10(Xlm[l-m]);
   if(v<0.) lv = log10(-Xlm[l-m]);
   cout<<"cx="<<cx<<" l="<<l<<" m="<<m<<" Plm="<<v<<" log10(Plm)="<<lv<<endl;
 }

Answer: See the data discontinuity at cos(teta)=-0.92

---------------------------------------------------------------------------------------
cx=-1 l=1997 m=796 Plm=0 log10(Plm)=0
cx=-0.999 l=1997 m=796 Plm=-0 log10(Plm)=0
cx=-0.998 l=1997 m=796 Plm=-0 log10(Plm)=0
cx=-0.997 l=1997 m=796 Plm=-0 log10(Plm)=0
cx=-0.996 l=1997 m=796 Plm=-0 log10(Plm)=0
cx=-0.995 l=1997 m=796 Plm=-0 log10(Plm)=0
cx=-0.994 l=1997 m=796 Plm=-0 log10(Plm)=0
cx=-0.993 l=1997 m=796 Plm=-0 log10(Plm)=0
cx=-0.992 l=1997 m=796 Plm=-0 log10(Plm)=0
cx=-0.991 l=1997 m=796 Plm=-0 log10(Plm)=0
cx=-0.99 l=1997 m=796 Plm=-0 log10(Plm)=0
cx=-0.989 l=1997 m=796 Plm=-0 log10(Plm)=0
cx=-0.988 l=1997 m=796 Plm=-0 log10(Plm)=0
cx=-0.987 l=1997 m=796 Plm=-0 log10(Plm)=0
cx=-0.986 l=1997 m=796 Plm=-0 log10(Plm)=0
cx=-0.985 l=1997 m=796 Plm=-0 log10(Plm)=0
cx=-0.984 l=1997 m=796 Plm=-0 log10(Plm)=0
cx=-0.983 l=1997 m=796 Plm=-0 log10(Plm)=0
cx=-0.982 l=1997 m=796 Plm=-0 log10(Plm)=0
cx=-0.981 l=1997 m=796 Plm=-0 log10(Plm)=0
cx=-0.98 l=1997 m=796 Plm=-0 log10(Plm)=0
cx=-0.979 l=1997 m=796 Plm=-0 log10(Plm)=0
cx=-0.978 l=1997 m=796 Plm=-0 log10(Plm)=0
cx=-0.977 l=1997 m=796 Plm=-0 log10(Plm)=0
cx=-0.976 l=1997 m=796 Plm=-0 log10(Plm)=0
cx=-0.975 l=1997 m=796 Plm=-0 log10(Plm)=0
cx=-0.974 l=1997 m=796 Plm=-0 log10(Plm)=0
cx=-0.973 l=1997 m=796 Plm=-0 log10(Plm)=0
cx=-0.972 l=1997 m=796 Plm=-0 log10(Plm)=0
cx=-0.971 l=1997 m=796 Plm=-0 log10(Plm)=0
cx=-0.97 l=1997 m=796 Plm=-0 log10(Plm)=0
cx=-0.969 l=1997 m=796 Plm=-0 log10(Plm)=0
cx=-0.968 l=1997 m=796 Plm=-0 log10(Plm)=0
cx=-0.967 l=1997 m=796 Plm=-0 log10(Plm)=0
cx=-0.966 l=1997 m=796 Plm=-0 log10(Plm)=0
cx=-0.965 l=1997 m=796 Plm=-0 log10(Plm)=0
cx=-0.964 l=1997 m=796 Plm=-0 log10(Plm)=0
cx=-0.963 l=1997 m=796 Plm=-0 log10(Plm)=0
cx=-0.962 l=1997 m=796 Plm=-0 log10(Plm)=0
cx=-0.961 l=1997 m=796 Plm=-0 log10(Plm)=0
cx=-0.96 l=1997 m=796 Plm=-0 log10(Plm)=0
cx=-0.959 l=1997 m=796 Plm=-0 log10(Plm)=0
cx=-0.958 l=1997 m=796 Plm=-0 log10(Plm)=0
cx=-0.957 l=1997 m=796 Plm=-0 log10(Plm)=0
cx=-0.956 l=1997 m=796 Plm=-0 log10(Plm)=0
cx=-0.955 l=1997 m=796 Plm=-0 log10(Plm)=0
cx=-0.954 l=1997 m=796 Plm=-0 log10(Plm)=0
cx=-0.953 l=1997 m=796 Plm=-0 log10(Plm)=0
cx=-0.952 l=1997 m=796 Plm=-0 log10(Plm)=0
cx=-0.951 l=1997 m=796 Plm=-0 log10(Plm)=0
cx=-0.95 l=1997 m=796 Plm=-0 log10(Plm)=0
cx=-0.949 l=1997 m=796 Plm=-0 log10(Plm)=0
cx=-0.948 l=1997 m=796 Plm=-0 log10(Plm)=0
cx=-0.947 l=1997 m=796 Plm=-0 log10(Plm)=0
cx=-0.946 l=1997 m=796 Plm=-0 log10(Plm)=0
cx=-0.945 l=1997 m=796 Plm=-0 log10(Plm)=0
cx=-0.944 l=1997 m=796 Plm=-0 log10(Plm)=0
cx=-0.943 l=1997 m=796 Plm=-0 log10(Plm)=0
cx=-0.942 l=1997 m=796 Plm=-0 log10(Plm)=0
cx=-0.941 l=1997 m=796 Plm=-0 log10(Plm)=0
cx=-0.94 l=1997 m=796 Plm=-0 log10(Plm)=0
cx=-0.939 l=1997 m=796 Plm=-0 log10(Plm)=0
cx=-0.938 l=1997 m=796 Plm=-0 log10(Plm)=0
cx=-0.937 l=1997 m=796 Plm=-0 log10(Plm)=0
cx=-0.936 l=1997 m=796 Plm=-0 log10(Plm)=0
cx=-0.935 l=1997 m=796 Plm=-0 log10(Plm)=0
cx=-0.934 l=1997 m=796 Plm=-0 log10(Plm)=0
cx=-0.933 l=1997 m=796 Plm=-0 log10(Plm)=0
cx=-0.932 l=1997 m=796 Plm=-0 log10(Plm)=0
cx=-0.931 l=1997 m=796 Plm=-0 log10(Plm)=0
cx=-0.93 l=1997 m=796 Plm=-0 log10(Plm)=0
cx=-0.929 l=1997 m=796 Plm=-0 log10(Plm)=0
cx=-0.928 l=1997 m=796 Plm=-0 log10(Plm)=0
cx=-0.927 l=1997 m=796 Plm=-0 log10(Plm)=0
cx=-0.926 l=1997 m=796 Plm=-0 log10(Plm)=0
cx=-0.925 l=1997 m=796 Plm=-0 log10(Plm)=0
cx=-0.924 l=1997 m=796 Plm=-0 log10(Plm)=0
cx=-0.923 l=1997 m=796 Plm=-0 log10(Plm)=0
cx=-0.922 l=1997 m=796 Plm=-0 log10(Plm)=0
cx=-0.921 l=1997 m=796 Plm=-0 log10(Plm)=0
cx=-0.92 l=1997 m=796 Plm=-0 log10(Plm)=0
cx=-0.919 l=1997 m=796 Plm=-0.253614 log10(Plm)=-0.595827
cx=-0.918 l=1997 m=796 Plm=-0.54171 log10(Plm)=-0.266233
cx=-0.917 l=1997 m=796 Plm=-0.943909 log10(Plm)=-0.0250699
cx=-0.916 l=1997 m=796 Plm=-1.27768 log10(Plm)=0.106422
cx=-0.915 l=1997 m=796 Plm=-1.19331 log10(Plm)=0.0767529
cx=-0.914 l=1997 m=796 Plm=-0.452665 log10(Plm)=-0.344223
cx=-0.913 l=1997 m=796 Plm=0.598563 log10(Plm)=-0.22289
cx=-0.912 l=1997 m=796 Plm=1.01032 log10(Plm)=0.00445825
cx=-0.911 l=1997 m=796 Plm=0.228302 log10(Plm)=-0.641489
---------------------------------------------------------------------------------------


 At present, Legendre Functions are very important
in astronomy and the studies of the early universe
with the Cosmic Background Fluctuations (CMB) measurement.
 New high-resolution CMB experiments require the computation
of Legendre fonctions up to Lmax as high as 3000 !
 To my knowledge, GSL Legendre Functions are the only C code
freeware. Some other exits (HealPix in Fortran 90 !!!)
and Sophya in C++. Sophya seems to work OK but it is not
easy to call Sophya C++ routines from C or f77 code.

          Could you do something about that bug ?

               With my thanks,

                      Best Regards,

                         Christophe Magneville


-- 
-----------------------------------------------------
Christophe Magneville    Email: address@hidden
CEA Saclay                     Tel: 33 (0)1 6908 3344
DAPNIA/SPP Bat 141             Fax: 33 (0)1 6908 6428
91191 Gif-sur-Yvette Cedex (France)
-----------------------------------------------------

gsl.txt
Description: Text document