[Bug-gsl] convergence of hyperg_2F1

[Andrew Benson]

In some cases, the convergence criterion logic in hyperg_2F1_series(a,b,c,x) 
seems to be broken. The Gauss series, summed over k, evaluates:

    del *= (a+k)*(b+k) * x / ((c+k) * (k+1.0)); 

and accumulates negative and positive terms separately. The convergence 
criterion is:

fabs((del_pos + del_neg)/(sum_pos-sum_neg)) <= GSL_DBL_EPSILON

where del_pos and del_neg are the most recent (in k) positive and negative 
values of del respectively. 

If x>0, then for sufficiently large k (such that a+k>0, b+k>0, and c+k>0) del 
no 
longer switches signs for successive values of k. This can break the 
convergence criterion, since, for example, if del remains negative for all 
subsequent values of k, del_neg continues to decrease, but del_pos remains 
fixed at some (possibly large) value. In this case, the series is truncated 
when del undeflows to zero (which is caught by the test for a or b being a 
negative integer). 

For example, in the case:

a=1.5000000000000000      
b=-1.8637325556200852      
c=-0.86373255562008522        
x=0.023977964684202376

this means that ~500 terms in the series are evaluated, rather than just the 
10 terms needed to evaluate hyperg_2F1 to machine precision. For some 
applications (e.g. my own!) this results in a huge speed-up.

I can supply a patch that detects these situations and adjusts the convergence 
criterion appropriately if there's interest.

-Andrew

-- 

* Andrew Benson: http://users.obs.carnegiescience.edu/abenson/contact.html

* Galacticus: http://sites.google.com/site/galacticusmodel