# [Scipy-svn] r5515 - in trunk/scipy/special: cephes tests

scipy-svn@scip... scipy-svn@scip...
Fri Jan 23 17:14:40 CST 2009

```Author: ptvirtan
Date: 2009-01-23 17:14:26 -0600 (Fri, 23 Jan 2009)
New Revision: 5515

Modified:
trunk/scipy/special/cephes/jv.c
trunk/scipy/special/tests/test_basic.py
Log:
Fix #623: better Cephes jv.c:recur continued fraction termination conditions

There were two problems in the old termination conditions for the
continued fraction for Jn/Jn-1:

1) It could terminate spuriously on the first iteration, because the
initial value for `ans` corresponded to a possible value of the fraction
but not one obtained from the iteration.

2) `jv` could call `recur` with values of `n` and `x` such that the
continued fraction iteration would not converge and the maximum
iteration limit was reached.

The code is now more careful with 1) and enforces a minimum number of
iterations, and 2) is worked around by bumping the maximum iteration
limit so that the CF should converge for all calls that come from `jv`.
The iteration limits can be estimated as the number of iterations
needed is known for this continued fraction.

Modified: trunk/scipy/special/cephes/jv.c
===================================================================
--- trunk/scipy/special/cephes/jv.c	2009-01-23 23:13:54 UTC (rev 5514)
+++ trunk/scipy/special/cephes/jv.c	2009-01-23 23:14:26 UTC (rev 5515)
@@ -98,15 +98,6 @@
double k, q, t, y, an;
int i, sign, nint;

-    /* The recursion used when n = 3 and x = 4 in recur gives
-
-       Simple fix for now:
-     */
-    if ((n==3) && (x == 4)) {
-       return 0.43017147387562193;
-    }
-
nint = 0;			/* Flag for integer n */
an = fabs(n);
@@ -265,8 +256,29 @@
double k, ans, qk, xk, yk, r, t, kf;
static double big = BIG;
int nflag, ctr;
+    int miniter, maxiter;

-/* continued fraction for Jn(x)/Jn-1(x)  */
+/* Continued fraction for Jn(x)/Jn-1(x)
+ * AMS 9.1.73
+ *
+ *    x       -x^2      -x^2
+ * ------  ---------  ---------   ...
+ * 2 n +   2(n+1) +   2(n+2) +
+ *
+ * Compute it with the simplest possible algorithm.
+ *
+ * This continued fraction starts to converge when (|n| + m) > |x|.
+ * Hence, at least |x|-|n| iterations are necessary before convergence is
+ * achieved. There is a hard limit set below, m <= 30000, which is chosen
+ * so that no branch in `jv` requires more iterations to converge.
+ * The exact maximum number is (500/3.6)^2 - 500 ~ 19000
+ */
+
+    maxiter = 22000;
+    miniter = fabs(x) - fabs(*n);
+    if (miniter < 1)
+        miniter = 1;
+
if (*n < 0.0)
nflag = 1;
else
@@ -284,7 +296,7 @@
qkm1 = *n + *n;
xk = -x * x;
yk = qkm1;
-    ans = 1.0;
+    ans = 0.0; /* ans=0.0 ensures that t=1.0 in the first iteration */
ctr = 0;
do {
yk += 2.0;
@@ -294,23 +306,28 @@
pkm1 = pk;
qkm2 = qkm1;
qkm1 = qk;
-	if (qk != 0)
+
+	/* check convergence */
+	if (qk != 0 && ctr > miniter)
r = pk / qk;
else
r = 0.0;
+
if (r != 0) {
t = fabs((ans - r) / r);
ans = r;
-	} else
+	} else {
t = 1.0;
+	}

-	if (++ctr > 1000) {
+	if (++ctr > maxiter) {
mtherr("jv", UNDERFLOW);
goto done;
}
if (t < MACHEP)
goto done;

+	/* renormalize coefficients */
if (fabs(pk) > big) {
pkm2 /= big;
pkm1 /= big;
@@ -321,6 +338,8 @@
while (t > MACHEP);

done:
+    if (ans == 0)
+        ans = 1.0;

#if CEPHES_DEBUG
printf("%.6e\n", ans);

Modified: trunk/scipy/special/tests/test_basic.py
===================================================================
--- trunk/scipy/special/tests/test_basic.py	2009-01-23 23:13:54 UTC (rev 5514)
+++ trunk/scipy/special/tests/test_basic.py	2009-01-23 23:14:26 UTC (rev 5515)
@@ -1629,17 +1629,29 @@
yvp1 = yvp(2,.2)
assert_array_almost_equal(yvp1,yvpr,10)

+class TestBesselJ(object):
+
+    def test_cephes_vs_specfun(self):
+        for v in [-120, -20., -10., -1., 0., 1., 12.49, 120., 301]:
+            for z in [1., 10., 200.5, 400., 600.5, 700.6, 1300, 10000]:
+                c1, c2, c3 = jv(v, z), jv(v,z+0j), jn(int(v), z)
+                if np.isinf(c1):
+                    assert np.abs(c2) >= 1e150
+                else:
+                    assert_tol_equal(c1, c2, err_msg=(v, z), rtol=1e-11)
+                    if v == int(v):
+                        assert_tol_equal(c2, c3, err_msg=(v, z), rtol=1e-11)
+
+    def test_ticket_623(self):
+        assert_tol_equal(jv(3, 4), 0.43017147387562193)
+        assert_tol_equal(jv(301, 1300), 0.0183487151115275)
+        assert_tol_equal(jv(301, 1296.0682), -0.0224174325312048)
+
class TestBesselI(object):

-    def _lggamma(self, q):
-        res = zeros_like(q)
-        res[q>=2] = gammaln(q[q>=2])
-        res[q<2] = log(gamma(q[q<2]))
-        return res
-
def _series(self, v, z, n=200):
k = arange(0, n).astype(float_)
-        r = (v+2*k)*log(.5*z) - self._lggamma(k+1) - self._lggamma(v+k+1)
+        r = (v+2*k)*log(.5*z) - gammaln(k+1) - gammaln(v+k+1)
r[isnan(r)] = inf
r = exp(r)
err = abs(r).max() * finfo(float_).eps * n + abs(r[-1])*10

```