BUG #1  gsl_sf_hyperg_2F1_e fails for some arguments
From: email is unavailable
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
gppari 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: email is unavailable
Subject: [Buggsl] 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 cab 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 email is unavailable
