I think the main issue that makes this tricky is that trying to
evaluate the power series for the hypergeometric function is
numerically unstable due to catastrophic cancellation for
one negative parameter and function value $\ll1$. Otherwise, having real parameters, and a real argument in $[0,1]$ should be fine.
If you look at the DLMF entry for hypergeometric functions, you will
find a relation between contiguous
functions, $F(a,b;c;z)$ and $F(a\pm1,b;c;z)$:
$$ f_a = \frac{F(a+1,b;c;z)}{F(a,b;c;z)} $$
$$ (c-a) \frac1{f_{a-1}} + (2a-c+(b-a)z) + a(z-1)f_a = 0. $$
So for large negative $a$, you can find $\bar a=a-\lfloor a\rfloor\geq0$
and then relate $f_a$ to $f_{\bar a}$, which are separated in the above recurrence relation by integer distance $-\lfloor a\rfloor$. This algorithm takes time
linear in $|a|$, so it would work well for $|a|<1000$, but poorly for
$|a|$ much larger than $10^6$, say.
Since $\bar a$ is positive, evaluating $f_{\bar a}$ is a lot simpler since it involves
hypergeometric function of real variable with all positive parameters
(so no cancellation). In particular, you can do with power series (or
GSL for that matter).
Here is some python code that I wrote to test this idea. (It doesn't work for some values of $a, c$ that cause division by zero, but that's straightforward to fix; you can use limits for example.)
from mpmath import *
def hsum(a, b, c, x):
ans, ansAbs, t, i = mpf(1), mpf(1), mpf(1), 0
while abs(t) / abs(ans) > mp.eps():
t *= (a + i) * (b + i) / ((c + i) * (i+1)) * x
i += 1
ans += t
ansAbs += abs(t)
print("hsum(a=%f): %f (%f)" % (a, ans, ansAbs / abs(ans)))
return ans
def h3(a, b, c, x):
a, b, c, x = mpf(a), mpf(b), mpf(c), mpf(x)
ans0 = hyp2f1(a+1, b, c, x) / hyp2f1(a, b, c, x)
if a >= 0:
return ans0
def f(a):
if a >= 0: return hsum(a+1, b, c, x) / hsum(a, b, c, x)
return - (c-a-1) / ((2*(a+1)-c+(b-a-1)*x) + (a+1)*(x-1)*f(a+1))
def f2(a):
n = int(floor(a))
a0 = a - n
ans = hsum(a0+1, b, c, x) / hsum(a0, b, c, x)
def alpha(a): return -(c-a-1)/((a+1) * (x-1))
def beta(a): return (2*(a+1)-c+(b-a-1)*x) / ((a+1)*(x-1))
aa = a0
for i in range(0, -n):
aa -= 1
ans = alpha(aa) / (beta(aa) + ans)
return ans
ans1 = f2(a)
if (abs(a) <= 200):
assert(abs(f2(a)/f(a)-1) < eps() * 1000)
print("a %s: %s" % (a, float(abs(ans1/ans0-1))))
return ans1
P.S. You said in a comment "I am computing the numerator and the
denominator, which are very small or very large, and thus their ratio
is not very accurate". This isn't really right. The numerical stability of
evaluating the expression $a/b$ can be found by computing the
condition number of the function $(a,b)\mapsto a/b$. It is
$$ \left|\frac{d\log(a/b)}{d\log a}\right| = 1, \qquad
\left|\frac{d\log(a/b)}{d\log b}\right| = 1. $$
Thus the computation of $a/b$ should be well-conditioned, regardless of their values (assuming no overflow/underflow, which is not what's going on for $|a|,|b|,|c|\approx 100$ anyway). You might be
thinking of $a-b$, which does have a singular condition number at $a\approx
b$:
$$ \left|\frac{d\log(a-b)}{d\log a}\right| = \frac{|a|}{|a-b|}. $$
This doesn't happen to $a/b$. The numerical problems here are of a
different sort.
GSL_EUNDRFLW
? Please post your code, otherwise it is hard to guess whether this is a problem with gsl. $\endgroup$