I recently tried implementing Gauss Hypergeometric function with c++ in two different ways, but found that they each had some problems in certain parameter regions.
The first way uses the naive series approach and is extremely simple to write:
http://www.cplusplus.com/forum/general/255896/
The problem is that this series diverges, e.g., when $x = 1$.
The second way simply calls the function from GNU standard library (GSL), but also encountered (expected) failure mode for the following set of parameters, discovered by my colleagues:
$a = 1, b = -62.5, c = 2, x \in [0.27^2, 1].$ Curiously, the series approach had no problem with this case.
But mathematica seems able to handle both cases well:
My suspicion is that mathematica uses quadrature method for better numeric stability, though I have no way to verify that. Also is quadrature going to be significantly slower than series methods?
As a bonus, what's a recommended c++ portable implementation to use for 2F1 with the largest possible parameter range? In fact my range is only going to be say $a \in [0, 2]$, $b \in [-100, 0]$, $c \in [0, 5]$ and $x \in [-1, 1]$. My current plan is to simply use the naive series implementation as a backup for GSL, but I am worried there may be points where both fail.