Re: [Bug-gsl] bug in gsl_sf_coulomb_wave_FG_e (gsl 1.14 compiled with gcc 4.5.0)

[Brian Gough]

At Wed, 25 Aug 2010 23:41:14 -0400,
Alexey A. Illarionov wrote:
> It seems that I find a bug in gsl_sf_coulomb_wave_FG_e. For large values 
> of lambda and small values of x, with eta = 0 it produces nan values 
> without raising of the flag.
> 
> Here the sample cases (see attachment with example)
>   gsl_sf_coulomb_wave_FG_e(0.,1.2693881947287221e-07, 37, 1, ...)
>   gsl_sf_coulomb_wave_FG_e(0.,5.911135077240691e-07, 39, 1, ...)
>   gsl_sf_coulomb_wave_FG_e(0.,6.118185507743884e-07, 40, 1, ...)
> and lots of others.

Thanks for the bug report. I'm adding a test case for this.  Can you
confirm if the following values are correct (I'm not working with
these functions often):

  lam_F = 37.0;
  eta = 0.0;
  x = 1.2693881947287221e-07;
  k_G = 1;
  gsl_sf_coulomb_wave_FG_e(eta, x, lam_F, k_G, &F, &Fp, &G, &Gp, &Fe, &Ge);
  s = 0;
  message_buff[0] = 0;
  s += test_sf_check_result(message_buff,  F,  
-2.020327525317343313380459069e299, TEST_TOL3);
  s += test_sf_check_result(message_buff, Fp, 
5.729672862364037942289798904e307, TEST_TOL3);
  s += test_sf_check_result(message_buff,  G, 
4.397355593129492554472984282e299, TEST_TOL3);
  s += test_sf_check_result(message_buff, Gp,  - 
1.247095270068213211088454516e308, TEST_TOL3);

I've computed them using Pari as

  c(L,n) = { 2^L * exp(-Pi*n/2) * abs(gamma(L+1+I*n)) / gamma(2*L+2)}
  cofac(L,n,r) = { I * exp(-I*r)*r^(-L)/(gamma(2*L+2)*c(L,n)) }
  FplusIG(L,n,r) = { cofac(L,n,r) * 
intnum(t=0,[[1],1],exp(-t)*t^(L-I*n)*(t+2*I*r)^(L+I*n)) }

  FplusIG(37,0,1.2693881947287221e-07+x) # for F and F' (real part)
  FplusIG(36,0,1.2693881947287221e-07+x) # for G and G' (imaginary part)

-- 
Brian Gough