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?

  • $\begingroup$ Wolfram Mathematica has no trouble with that input, although for $a = 10^6$ the computing time is about seconds. $\endgroup$ – uranix Jul 30 '15 at 14:35
  • $\begingroup$ @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$? $\endgroup$ – vgnils Jul 30 '15 at 15:24
  • $\begingroup$ 2 seconds. I found an interesting recurrence here on page 31 people.maths.ox.ac.uk/porterm/research/pearson_final.pdf . I'm developing an answer now $\endgroup$ – uranix Jul 30 '15 at 15:26
  • $\begingroup$ But there's a mistake, will try to fix it $\endgroup$ – uranix Jul 30 '15 at 15:39

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;
  • $\begingroup$ 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. $\endgroup$ – vgnils Jul 31 '15 at 8:43
  • 1
    $\begingroup$ @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, math.stackexchange.com/questions/478052/… . There is a link to MATLAB code that is using this idea among others. $\endgroup$ – uranix Jul 31 '15 at 9:00

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.