search

RobinHankin / gsl

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

Bug in hyperg_2F1 #13

Closed michaelgordy closed 3 years ago

michaelgordy commented 3 years ago

This bug is probably in GSL rather than your wrapper, but I wonder what you think. Attached zip provides Rmd example and HTML output.

bug2F1.zip

RobinHankin commented 3 years ago

thanks for this. The hypergeometric functions are difficult to deal with near x=1: there is a branch point on the complex plane there, and GSL is restricted to the real case, which causes problems. The hypergeo package has to jump through hoops to deal the branch cuts. But yes, as you say, it is almost certainly an infelicity [I hesistate to say "bug"] in GSL rather than gsl.

best wishes

Robin

@.***>

On Thu, Mar 18, 2021 at 8:46 AM michaelgordy @.***> wrote:

This bug is probably in GSL rather than your wrapper, but I wonder what you think. Attached zip provides Rmd example and HTML output.

bug2F1.zip https://github.com/RobinHankin/gsl/files/6159321/bug2F1.zip

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/RobinHankin/gsl/issues/13, or unsubscribe https://github.com/notifications/unsubscribe-auth/ADFFZUWXNJNKZJYNFR25653TEEBLFANCNFSM4ZLHWLGA .

michaelgordy commented 3 years ago

Thanks, Robin. I can confirm that the infelicity arises within GSL, as I tested yesterday using the Julia wrapper to the library. It is odd that the issue is so local to specific parameter values, but I can appreciate that things get ugly as x->1. I will take the liberty of closing the Issue.

RobinHankin commented 3 years ago

Thanks for this Michael. MWE follows together with mathematica values for comparison:

> hyperg_2F1(1.3,1,2.3,c(0.994,0.995))
[1] 6.170794e+00 5.854680e+15
>

Mathematica gives

6.1707936168
6.4009067299

for Hypergeometric2F1[1.3,1,2.3,0.994] and Hypergeometric2F1[1.3,1,2.3,0.994] respectively. Happy to leave closed as it is not an error with the gsl R package but one day I should report it to the GSL maintainers.