[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Bug-gsl] convergence of hyperg_2F1
From: |
Andrew Benson |
Subject: |
[Bug-gsl] convergence of hyperg_2F1 |
Date: |
Sat, 31 Jan 2015 17:52:08 -0800 |
User-agent: |
KMail/4.14.3 (Linux/3.18.3-201.fc21.x86_64; KDE/4.14.3; x86_64; ; ) |
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