Compute hypergeometric function ratio: $\frac{_{2}F_{1}(a+1,b;c;x)}{_{2}F_{1}(a,b;c;x)}$?

Question

I need a numerically stable way to compute the following ratio:

$$\frac{_{2}F_{1}(a+1,b;c;x)}{_{2}F_{1}(a,b;c;x)}$$

All the parameters are real numbers, with $a< 0$,$\ $ $b,c > 0$ and $0<x<1$.

Right now I am using GSL's implementation of the hypergeometric function, but I keep getting under/over-flows.

Is there a simplification that I can use? Or an aproximation scheme that remains accurate for a wide range of the parameters (say $-1000 < a,b,c < 1000$)?

As a particular case, I am also interested in positive values of the parameters, in the range $0<a,b,c<1000$. What approximation can I use in that case?

Answer 1

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.