Should be doing regularised _2F1
_2F1 sucks compared to regularised _2F1 (bold F in https://dlmf.nist.gov/15.2 or Hypergeometric2F1Regularized in mathematica), as the regularised variant avoids the artificial polls at negative integers for c. We really should compute the regularised version and then define
_₂F₁(a, b, c, z) = gamma(c)*regularized_₂F₁(a, b, c, z)
I think the reason I chose this way is because if you are near the origin you don't need to call gamma(c) as the Maclaurin series will do.
If we want to implement the regularized version of the Michiels & Stoitsov algorithm, we'd just need to scale these initial conditions by 1/gamma(c)
https://github.com/JuliaMath/HypergeometricFunctions.jl/blob/ee236778e08b0795e2b3e24f02f919423fcfabd0/src/specialfunctions.jl#L309-L312
for example.
Possibly can be done for the rational algorithms too.