This file is the GSL bug tracking system.  The CVS version of this
file should be kept up-to-date.

----------------------------------------------------------------------
BUG#1 -- gsl_sf_hyperg_2F1_e fails for some arguments 

From: keith.briggs@bt.com
Subject: gsl_sf_hyperg_2F1 bug report
Date: Thu, 31 Jan 2002 12:30:04 -0000

gsl_sf_hyperg_2F1_e fails with arguments (1,13,14,0.999227196008978,&r).
It should return 53.4645... .

#include <gsl/gsl_sf.h>
#include <stdio.h>

int main (void)
{
  gsl_sf_result r;
  gsl_sf_hyperg_2F1_e (1,13,14,0.999227196008978,&r);
  printf("r = %g %g\n", r.val, r.err);
}

NOTES: The program overflows the maximum number of iterations in
gsl_sf_hyperg_2F1, due to the presence of a nearby singularity at
(c=a+b,x=1) so the sum is slowly convergent.

The exact result is 53.46451441879150950530608621 as calculated by
gp-pari using sumpos(k=0,gamma(a+k)*gamma(b+k)*gamma(c)*gamma(1)/
(gamma(c+k)*gamma(1+k)*gamma(a)*gamma(b))*x^k)

The code needs to be extended to handle the case c=a+b. This is the
main problem. The case c=a+b is special and needs to be computed
differently.  There is a special formula given for it in Abramowitz &
Stegun 15.3.10

As reported by Lee Warren <warren@atom.chem.utk.edu> another set of
arguments which fail are:

#include <gsl/gsl_sf.h>
#include <stdio.h>

int main (void)
{
  gsl_sf_result r;
  gsl_sf_hyperg_2F1_e (-1, -1, -0.5, 1.5, &r);
  printf("r = %g %g\n", r.val, r.err);
}

The correct value is -2.

See also, 

From: Olaf Wucknitz <wucknitz@jive.nl>
To: bug-gsl@gnu.org
Subject: [Bug-gsl] gsl_sf_hyperg_2F1

Hi,

I am having a problem with gsl_sf_hyperg_2F1.
With the parameters (-0.5, 0.5, 1, x) the returned values show a jump at 
x=0.5. For x<0.5 the results seem to be correct, while for x>0.5 they 
aren't.
The function gsl_sf_hyperg_2F1_e calls hyperg_2F1_series for x<0.5, but
hyperg_2F1_reflect for x>0.5. The latter function checks for c-a-b being 
an integer (which it is in my case). If I change one of the parameters 
a,b,c by a small amount, the other branch of the function is taken and the 
results are correct again.
Unfortunately I know too little about the numerics of 2F1 to suggest a 
patch at the moment.

Regards,
Olaf Wucknitz
-- 
Joint Institute for VLBI in Europe                wucknitz@jive.nl

----------------------------------------------------------------------
BUG#44 -- gamma_inc_P and gamma_inc_Q only satisfy P+Q=1 within errors

The sum of gamma_inc_P and gamma_inc_Q doesn't always satisfy the
identity P+Q=1 exactly (although it is correct within errors), due the
slightly different branch conditions for the series and continued
fraction expansions.  These could be made identical so that P+Q=1 exactly.

#include <stdio.h>
#include <gsl/gsl_sf_gamma.h>

int
main (void)
{
  gsl_sf_result r1, r2;
  double a = 0.3, x = 1.0;
  gsl_sf_gamma_inc_P_e (a, x, &r1);
  gsl_sf_gamma_inc_Q_e (a, x, &r2);
  printf("%.18e\n", r1.val);
  printf("%.18e\n", r2.val);
  printf("%.18e\n", r1.val + r2.val);
}

$ ./a.out
9.156741562411074842e-01
8.432584375889111417e-02
9.999999999999985567e-01


---------------------------------------------------------------------
BUG#50 - gsl_linalg_solve_symm_tridiag requires positive definite matrix

A zero on the diagonal will cause NaNs even though a reasonable
solution could be computed in principle. 

#include <gsl/gsl_linalg.h>

int main (void)
{
  double d[] = { 0.00, 1.21, 0.80, 1.55, 0.76 } ;
  double e[] = { 0.82, 0.39, 0.09, 0.68 } ;
  double b[] = { 0.07, 0.62, 0.81, 0.11, 0.65} ;
  double x[] = { 0.00, 0.00, 0.00, 0.00, 0.00} ;

  gsl_vector_view dv = gsl_vector_view_array(d, 5);
  gsl_vector_view ev = gsl_vector_view_array(e, 4);
  gsl_vector_view bv = gsl_vector_view_array(b, 5);
  gsl_vector_view xv = gsl_vector_view_array(x, 5);

  gsl_linalg_solve_symm_tridiag(&dv.vector, &ev.vector, &bv.vector, &xv.vector);
  gsl_vector_fprintf(stdout, &xv.vector, "% .5f");

  d[0] += 1e-5;
  gsl_linalg_solve_symm_tridiag(&dv.vector, &ev.vector, &bv.vector, &xv.vector);
  gsl_vector_fprintf(stdout, &xv.vector, "% .5f");
}

$ ./a.out
 nan
 nan
 nan
 nan
 nan
 0.13626
 0.08536
 1.03840
-0.60009
 1.39219

AUG 2007: We now return an error code for this case.  To return a solution
we would need to do a permutation, see slatec/dgtsl.f

----------------------------------------------------------------------
Last assigned bug number = 62