Re: [Bug-gsl] incorrect values by bessel_steed_array

[Brian Gough]

Hu Zhan writes:
 > There is one line in gsl_sf_bessel_jl_steed_array:
 > 
 >     W = x_inv / sqrt(FP*FP + F*F);
 > 
 > When LMAX is large, FP and F can be larger than 1e160, and cause W = 0.
 > The option -mfpmath=387 (90-bit arithmetics internally) can make it work
 > on the opteron, but not pentium.

Thanks for the explanation, I see the problem.  Could you try the
following patch (to replace sqrt(a^2+b^2) by hypot) and let me know if
it works.

-- 
best regards

Brian Gough


Index: bessel_j.c
===================================================================
RCS file: /home/gsl-cvs/gsl/specfunc/bessel_j.c,v
retrieving revision 1.36
diff -u -r1.36 bessel_j.c
--- bessel_j.c  26 Jul 2003 13:44:38 -0000      1.36
+++ bessel_j.c  11 Feb 2005 15:42:01 -0000
@@ -357,7 +357,7 @@
     }
     
     /* normalization */
-    W = x_inv / sqrt(FP*FP + F*F);
+    W = x_inv / hypot(FP, F);
     jl_x[0] = W*F;
     if(lmax > 0) {
       int L;
Index: test_bessel.c
===================================================================
RCS file: /home/gsl-cvs/gsl/specfunc/test_bessel.c,v
retrieving revision 1.29
diff -u -r1.29 test_bessel.c
--- test_bessel.c       25 Jun 2003 21:09:31 -0000      1.29
+++ test_bessel.c       11 Feb 2005 15:34:54 -0000
@@ -380,11 +380,12 @@
   s += sa;
 
   sa = 0;
-  gsl_sf_bessel_jl_steed_array(50, 1.0, J);
+  gsl_sf_bessel_jl_steed_array(99, 1.0, J);
   sa += ( test_sf_frac_diff(J[0],  0.84147098480789650670   ) > TEST_TOL0 );
   sa += ( test_sf_frac_diff(J[1],  0.30116867893975678925   ) > TEST_TOL0 );
   sa += ( test_sf_frac_diff(J[10], 7.116552640047313024e-11 ) > TEST_TOL0 );
   sa += ( test_sf_frac_diff(J[50], 3.615274717489787311e-81 ) > TEST_TOL0 );
+  sa += ( test_sf_frac_diff(J[80], 1.136352423414503264e-144 ) > TEST_TOL0 );
   gsl_test(sa, "  gsl_sf_bessel_jl_steed_array");
   s += sa;