search

RobinHankin / gsl

an R wrapper for the Gnu Scientific Library
15 stars 9 forks source link

GSL: hyperg_2F1 returns Inf #19

Closed karchjd closed 2 years ago

karchjd commented 2 years ago

Not really an issue, but I did not know a better platform to ask for help. I need to evaluate hyperg_2F1 at the input hyperg_2F1(a,a,b,z) with a=N/2, and b=p/2, N integer from 1 to unlimited, p integer from 1 to unlimited, and N>p. This typically works but, for example, fails for N=500, p=2, z=0.9. In this case, gsl::hyperg_2F1(250, 250, 1, 0.9) returns Inf. I am assuming that this is caused by overflow.

Any suggestion on how to solve this? My main idea was to find a transformation (I tried Euler's transformation but failed) that avoids inputs that cause the overflow. A brute force idea might be to increase the number of bits used for floating-point numbers in gsl, but I am not sure whether I, as a user, have control over this.

RobinHankin commented 2 years ago

not really my bag, as the R package is nothing but a brain-dead wrapper for GSL, but I would observe that mathematica gives about 7e641 (that is, 7 times 10^641) so returning Inf is reasonable?

karchjd commented 2 years ago

Thanks for the swift reply. I was expecting this to be the answer. Yes, Inf is kind of reasonable. In the code I use it, however, it is multiplied with a very small number, and the result should be between 0 and 1, that's, of course, not happening if one of the factors is inf. I will look for help elsewhere.

RobinHankin commented 2 years ago

The hypergeo package is a little shabby too in this case (I'll raise an issue and link back here):

> library(hypergeo)
> hypergeo(250,250,1,0.9)
[1] -1.266637e+20+0i
>