Numerically evaluate Gauss' hypergeometric function ${}_{2}F_{1}(a,b;c;x) $ for large $|a|$ or $|b|$ and $x\ll 0$ or $ x \approx 1$?

Question

I need to compute Gauss' hypergeometric function $${}_{2}F_{1}(a,b;c;x)$$ for the case where one of $|a|$ or $|b|$ is large and $x\ll 0$ or $ x \approx 1$. By employing some linear transformations, I can choose whether I want $x\ll 0$ or $ x \approx 1$, but the upper parameters remain large in absolute value.

For example, consider $a=1000, b=0.006, c=0.01, x = -9000$. Using the transformations or not, these parameter values result in errors with hypergeom in MATLAB, hyperg_2F1 from the gsl package in R and hypergeo from the hypergeo package in R. It works when $a$ is around 100. However, I need to compute ${}_{2}F_{1}$ with $a$ upwards of $10^6$.

Are there any known tricks for this problem?

Wolfram Mathematica has no trouble with that input, although for $a = 10^6$ the computing time is about seconds. — uranix, Jul 30 '15 at 14:35
@uranix That's interesting - then at least I can be sure there are ways to do this. How long was the computing time for $a=10^6$? — vgnils, Jul 30 '15 at 15:24
2 seconds. I found an interesting recurrence here on page 31 https://people.maths.ox.ac.uk/porterm/research/pearson_final.pdf . I'm developing an answer now — uranix, Jul 30 '15 at 15:26

score 2 · Accepted Answer · answered Jul 30 '15 at 18:26

From Abramowitz and Stegun, 15.2.10 $$ (c-a){}_2F_1(a-1, b; c; z) + (2a-c+(b-a)z){}_2F_1(a, b; c; z) +a(z-1){}_2F_1(a+1, b; c; z) = 0 $$ Let $G(a) = {}_2F_1(a, b; c; z)$. Then we can use $$ G(a+1) = \frac{2a-c+(b-a)z}{a(1-z)}G(a)+\frac{c-a}{a(1-z)}G(a-1) $$ recurrence to compute $G(a)$ for big values of $a$, starting from a pair $G(a - \lfloor a \rfloor + 1), G(a - \lfloor a \rfloor + 2)$.

Here's the MATLAB code implementing this idea (relies on hypergeom for small $a$)

function g = Hyp2F1(a, b, c, z)
    steps = max(floor(a) - 2, 0);
    t = a - steps;
    g = hypergeom([t, b], c, z);
    gprev = hypergeom([t-1, b], c, z);

    for j = 1:steps
        gsave = g;
        alpha = (2 * t - c + (b - t) * z) / t / (1 - z);
        beta = (c - t) / t / (1 - z);
        g = alpha * g + beta * gprev;
        t = t + 1;
        gprev = gsave;
    end
end

A very clean solution that seems to be fast and reliable - it's a bit strange that such techniques haven't been implemented in the 2F1-functions available for MATLAB and R. Anyway, many thanks. — vgnils, Jul 31 '15 at 08:43
@vgnils Thanks, I'm thinking of some optimization that should improve the precision. I think this case is just quite rare and there are too many possible cases for parameters, thus covering every case is impossible. Also, take a look at this question, http://math.stackexchange.com/questions/478052/fast-matlab-code-for-hypergeometric-function-2f-1 . There is a link to MATLAB code that is using this idea among others. — uranix, Jul 31 '15 at 09:00

Numerically evaluate Gauss' hypergeometric function ${}_{2}F_{1}(a,b;c;x) $ for large $|a|$ or $|b|$ and $x\ll 0$ or $ x \approx 1$?

1 Answers1