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

[scipy-svn at scipy.org]

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.

Additional tests included.

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 
-       the wrong answer.   
-       
-       Simple fix for now:
-     */
-    if ((n==3) && (x == 4)) {
-       return 0.43017147387562193;
-    }
-
     nint = 0;			/* Flag for integer n */
     sign = 1;			/* Flag for sign inversion */
     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