The present implementation of the function in scipy, involves one Euler
Transform followed by the application of Power Series (as a and b turn negative after Euler Transform is applied once). This most probably blows up the values of function as |c|>>>|a| or |b|.
In one of the attempts to debug this the function hyp2f1 (when c<0; a,b>0; |c|>>a,b) was made to call hys2f1 (power series) without Euler Transformation. The result was the same
as what mpmath gives when the precision is less (~100), as the error tolerance
for hys2f1 is high. (MACHEP of order 10^-17). On increasing the sensitivity by
changing the local tolerance to the order of 10^-200, it works perfect.
What are the implications of adopting this method in the cases where hyp2f1
fails to give accurate results? Except of course the fact that this implementation
would be heavy.
Our current hyp2f1 implementation does use recurrences (hyp2f1ra), but
perhaps they are not invoked for this case. The problem here can be the
accurate determination of the convergence region for each parameter value.
The straight-forwardness of power series and its resemblance to mpmath tempted
me to first try with hys2f1, However I have a strong feeling that owing to strong sensitivity
of recurrence function in general, the implementation will be much faster. However
direct implementation without a change in sensitivity fails to give the required answer.
