bugGNU Scientific Library - Bugs: bug #21835, gsl_sf_hyperg_2F1 problematic...

 
 

You are not allowed to post comments on this tracker with your current authentication level.

bug #21835: gsl_sf_hyperg_2F1 problematic arguments

Submitter:  -Deleted Account- <bjg>
Submitted:  Tue 18 Dec 2007 06:49:13 PM UTC
   
 
Category:  Accuracy problem Severity:  3 - Normal
Operating System:  Status:  None
Assigned to:  None Open/Closed:  Open
Release:  1.10

Tue 18 Dec 2007 06:49:13 PM UTC, original submission:  

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
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: -email is unavailable-
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                -email is unavailable-

-Deleted Account- <bjg>

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

No files currently attached

 

Depends on the following items: None found

Items that depend on this one: None found

 

CC list is empty

 

No changes have been made to this item

Back to the top

Powered by Savane 3.13-3230.
Corresponding source code