# 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$?

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 people.maths.ox.ac.uk/porterm/research/pearson_final.pdf . I'm developing an answer now – uranix Jul 30 '15 at 15:26
• But there's a mistake, will try to fix it – uranix Jul 30 '15 at 15:39

## 1 Answer

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