# GNU Scientific Library - Bugs: bug #30510, problems with hyperg_U(a,b,x) for...

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

*bug #30510*: problems with hyperg_U(a,b,x) for x<0

Submitter: | -Deleted Account- <bjg> | ||

Submitted: | Wed 21 Jul 2010 08:17:45 PM UTC | ||

Category: | Runtime error | Severity: | 3 - Normal |

Operating System: | Status: | None | |

Assigned to: | None | Open/Closed: | Open |

Release: | 1.14 |

No files currently attached

Depends on the following items: None found

Items that depend on this one: None found

Carbon-Copy List

Follows 1 latest change.

Date | Changed by | Updated Field | Previous Value | => | Replaced by |
---|---|---|---|---|---|

2010-07-21 | bjg | Carbon-Copy | - | Added -email is unavailable- |

Reply-To: -email is unavailable-

From: Raymond Rogers <raymond.rogers72@gmail.com>

To: -email is unavailable-

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

Date: Thu, 08 Jul 2010 11:49:15 -0500

Brian Gough <bjg@gnu.org>

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

of a

To: -email is unavailable-

Cc: -email is unavailable-

Message-ID: <43wrt61fkt.wl%bjg@gnu.org>

Content-Type: text/plain; charset=US-ASCII

At Wed, 07 Jul 2010 10:14:34 -0500,

Raymond Rogers wrote:

| Brian Gough

| Thanks for the email. There are some comments about the domain for

| the hyperg_U_negx function in specfunc/hyperg_U.c -- do they help?

They explain some things, but I believe the section

if (b_int && b >= 2 && !(a_int && a <= (b - 2))){}

else {}

is implemented incorrectly; and probably the preceding section as well. Some restructuring of the code would make things clea\

rer; but things like that should probably done in a different forum: email, blog, etc...

I think the switches might be wrong. In any case it seems that b=1 has a hole. Is there a source for this code?

Note: the new NIST Mathematical handbook might have better algorithms. I am certainly no expert on implementing mathematical \

functions (except for finding ways to make them fail).

Ray

Reply-To: -email is unavailable-

From: Raymond Rogers <raymond.rogers72@gmail.com>

To: -email is unavailable-

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

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