Incorrect values in hypergeometric_2F1

[amellnik]

My coworker @nicolasdidier noticed this strange issue with sf_hyperg_2F1 on half of the domain:

(I wasn't able to test it on 1.0 due to unsatisfiable package requirements.) Is there any way this could be an issue with this package or is this definitely due to GSL?

[amellnik]

I'm reasonably sure that it's the same thing as this bug in GSL itself: http://lists.gnu.org/archive/html/bug-gsl/2015-09/msg00003.html . I'll keep this issue open to track when it's fixed upstream.

[simonbyrne]

It looks like it has been fixed upstream. This package is still using a quite old version of GSL, so should be upgraded.

[ludvigak]

This still appears to be an issue, even after switching to GSL 2.5. Interestingly, the error estimate provided by GSL is also completely off

using GSL, PyPlot, PyCall

x = 0:0.01:1
a, b, c = -1/4, 1/4, 1
out = sf_hyperg_2F1_e.(a, b, c, x)
gslval = [r.val for r in out]
gslerr = [r.err for r in out]
pyval = pyimport("scipy.special")[:hyp2f1](a, b, c, x)

subplot(2,1,1);
plot(x, gslval, ".", label="GSL value")
plot(x, pyval, "--", label="SciPy value")
legend(); grid(); title("sf_hyperg_2F1_e(-1/4, 1/4, 1, x)")
subplot(2,1,2);
semilogy(x, gslerr, ".", label="GSL error est")
legend(); grid()

[simonbyrne]

It seems to be the second bug mentioned in #21835:

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.

[ludvigak]

That's too bad. There actually seems to be quite a few bug reports on gsl_sf_hyperg_2F1. Let's just hope they get fixed in a future version.