bug-gsl
[Top][All Lists]

[Bug-gsl] Re: hyperg_U(a, b, x) Questions about x<0 and values of a,

 From: Raymond Rogers Subject: [Bug-gsl] Re: hyperg_U(a, b, x) Questions about x<0 and values of a, Date: Sun, 11 Jul 2010 14:43:43 -0500 User-agent: Thunderbird 2.0.0.24 (X11/20100411)

```hyperg_U basically fails with b=1, a non-integer; because
gsl_sf_poch_e(1+a-b,-a,&r1); is throwing a domain error when given
gamma(0)/gamma(a).
Checking on and using b=1 after a-integer is checked is illustrated
below in Octave.   I also put in recursion to evaluate b>=2.
I checked the b=1 expression against Maple; for a few values x<0,a<0,b=1
and x<0,a<0,b>=2 integer.
--------------
Unfortunately the routine in Octave to call hyperg_U  is only set up for
real returns, which was okay for versions <1.14 .  Sad to say I am the
one who implemented the hyperg_U interface, and will probably have to go
back :-( .  Integrating these functions into Octave was not pleasant;
but perhaps somebody made it easier.  I did translate the active parts
of hyperg_U into octave though; so it can be used in that way.

Ray

#
#
# Test function to evaluate b=1 for gsl 1.14 hyperg_U x<0
#
function anss=hyperg_U_negx_1(a,b,x)
int_a=(floor(a)==a);
int_b=(floor(b)==b);
#neg, int, a is already taken care of so use it
if (int_a && a<=0)
anss=hyperg_U(a,b,x);
elseif (int_b && (b==1))
#from the new NIST DLMF 13.2.41

anss=gamma(1-a)*exp(x)*(hyperg_U(1-a,1,-x)/gamma(a)-hyperg_1F1(1-a,1,-x)*exp((1-a)*pi*I));
elseif (b>=2)
#DLMF 13.3.10
anss=((b-1-a)*hyperg_U_negx_1(a,b-1,x) +
hyperg_U_negx_1(a-1,b-1,x))/x;
else
anss=hyperg_U(a,b,x);
endif
#
endfunction

```