[Scipy-svn] r6970 - in trunk/scipy/special: . specfun tests

scipy-svn@scip... scipy-svn@scip...
Sun Nov 28 16:55:42 CST 2010


Author: ptvirtan
Date: 2010-11-28 16:55:42 -0600 (Sun, 28 Nov 2010)
New Revision: 6970

Modified:
   trunk/scipy/special/cephes_doc.h
   trunk/scipy/special/specfun/specfun.f
   trunk/scipy/special/specfun_wrappers.c
   trunk/scipy/special/tests/test_basic.py
   trunk/scipy/special/tests/test_data.py
   trunk/scipy/special/tests/test_mpmath.py
Log:
ENH: special/lpmv: implement degree down-recursion in associated legendre polynomial (fixes #1296)

Also implement support for negative orders and degrees, except in the
cases V+M+1 <= 0, M < 0.

Modified: trunk/scipy/special/cephes_doc.h
===================================================================
--- trunk/scipy/special/cephes_doc.h	2010-11-28 17:46:31 UTC (rev 6969)
+++ trunk/scipy/special/cephes_doc.h	2010-11-28 22:55:42 UTC (rev 6970)
@@ -104,7 +104,7 @@
 #define kv_doc "y=kv(v,z) returns the modified Bessel function of the second kind (sometimes called the third kind) for\nreal order v at complex z."
 #define kve_doc "y=kve(v,z) returns the exponentially scaled, modified Bessel function\nof the second kind (sometimes called the third kind) for real order v at complex z: kve(v,z) = kv(v,z) * exp(z)"
 #define log1p_doc "y=log1p(x) calculates log(1+x) for use when x is near zero."
-#define lpmv_doc "y=lpmv(m,v,x) returns the associated legendre function of integer order\nm and nonnegative degree v: |x|<=1."
+#define lpmv_doc "y=lpmv(m,v,x) returns the associated legendre function of integer order\nm and real degree v (s.t. v>-m-1 or v<m): |x|<=1."
 #define mathieu_a_doc "lmbda=mathieu_a(m,q) returns the characteristic value for the even solution, \nce_m(z,q), of Mathieu's equation"
 #define mathieu_b_doc "lmbda=mathieu_b(m,q) returns the characteristic value for the odd solution, \nse_m(z,q), of Mathieu's equation"
 #define mathieu_cem_doc "(y,yp)=mathieu_cem(m,q,x) returns the even Mathieu function, ce_m(x,q), \nof order m and parameter q evaluated at x (given in degrees).\nAlso returns the derivative with respect to x of ce_m(x,q)"

Modified: trunk/scipy/special/specfun/specfun.f
===================================================================
--- trunk/scipy/special/specfun/specfun.f	2010-11-28 17:46:31 UTC (rev 6969)
+++ trunk/scipy/special/specfun/specfun.f	2010-11-28 22:55:42 UTC (rev 6970)
@@ -7771,7 +7771,7 @@
 
 C       **********************************
 
-        SUBROUTINE LPMV(V,M,X,PMV)
+        SUBROUTINE LPMV0(V,M,X,PMV)
 C
 C       =======================================================
 C       Purpose: Compute the associated Legendre function
@@ -7866,8 +7866,63 @@
         RETURN
         END
 
+C       **********************************
 
+        SUBROUTINE LPMV(V,M,X,PMV)
+C
+C       =======================================================
+C       Purpose: Compute the associated Legendre function
+C                Pmv(x) with an integer order and an arbitrary
+C                degree v, using down-recursion for large degrees
+C       Input :  x   --- Argument of Pm(x)  ( -1 ≤ x ≤ 1 )
+C                m   --- Order of Pmv(x)
+C                v   --- Degree of Pmv(x)
+C       Output:  PMV --- Pmv(x)
+C       Routine called:  LPMV0
+C       =======================================================
+C
+        IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+        IF (X.EQ.-1.0D0.AND.V.NE.INT(V)) THEN
+           IF (M.EQ.0) PMV=-1.0D+300
+           IF (M.NE.0) PMV=1.0D+300
+           RETURN
+        ENDIF
+        VX=V
+        MX=M
+        IF (V.LT.0) THEN
+           VX=-VX-1
+        ENDIF
+        NEG_M=0
+        IF (M.LT.0.AND.(VX+M+1).GE.0D0) THEN
+C          XXX: does not handle the cases where AMS 8.2.5
+C               does not help
+           NEG_M=1
+           MX=-M
+        ENDIF
+        NV=INT(VX)
+        V0=VX-NV
+        IF (NV.GT.2.AND.NV.GT.MX) THEN
+C          Up-recursion on degree, AMS 8.5.3
+           CALL LPMV0(V0+MX, MX, X, P0)
+           CALL LPMV0(V0+MX+1, MX, X, P1)
+           PMV = P1
+           DO 10 J=MX+2,NV
+              PMV = ((2*(V0+J)-1)*X*P1 - (V0+J-1+MX)*P0) / (V0+J-MX)
+              P0 = P1
+              P1 = PMV
+10         CONTINUE
+        ELSE
+           CALL LPMV0(VX, MX, X, PMV)
+        ENDIF
+        IF (NEG_M.NE.0.AND.ABS(PMV).LT.1.0D+300) THEN
+C          AMS 8.2.5, for integer order
+           CALL GAMMA2(VX-MX+1, G1)
+           CALL GAMMA2(VX+MX+1, G2)
+           PMV = PMV*G1/G2 * (-1)**MX
+        ENDIF
+        END
 
+
 C       **********************************
 
         SUBROUTINE CGAMA(X,Y,KF,GR,GI)

Modified: trunk/scipy/special/specfun_wrappers.c
===================================================================
--- trunk/scipy/special/specfun_wrappers.c	2010-11-28 17:46:31 UTC (rev 6969)
+++ trunk/scipy/special/specfun_wrappers.c	2010-11-28 22:55:42 UTC (rev 6970)
@@ -608,6 +608,7 @@
   if (m != floor(m)) return NPY_NAN;
   int_m = (int ) m;
   F_FUNC(lpmv,LPMV)(&v, &int_m, &x, &out);
+  CONVINF(out);
   return out;
 }
 

Modified: trunk/scipy/special/tests/test_basic.py
===================================================================
--- trunk/scipy/special/tests/test_basic.py	2010-11-28 17:46:31 UTC (rev 6969)
+++ trunk/scipy/special/tests/test_basic.py	2010-11-28 22:55:42 UTC (rev 6970)
@@ -1884,7 +1884,9 @@
 
     def test_lpmv(self):
         lp = special.lpmv(0,2,.5)
-        assert_almost_equal(lp,-0.125,3)
+        assert_almost_equal(lp,-0.125,7)
+        lp = special.lpmv(0,40,.001)
+        assert_almost_equal(lp,0.1252678976534484,7)
 
     def test_lqmn(self):
         lqmnf = special.lqmn(0,2,.5)

Modified: trunk/scipy/special/tests/test_data.py
===================================================================
--- trunk/scipy/special/tests/test_data.py	2010-11-28 17:46:31 UTC (rev 6969)
+++ trunk/scipy/special/tests/test_data.py	2010-11-28 22:55:42 UTC (rev 6970)
@@ -5,7 +5,7 @@
     arccosh, arcsinh, arctanh, erf, erfc, log1p, expm1,
     jn, jv, yn, yv, iv, kv, kn, gamma, gammaln, digamma, beta, cbrt,
     ellipe, ellipeinc, ellipk, ellipj, erfinv, erfcinv, exp1, expi, expn,
-    zeta, gammaincinv,
+    zeta, gammaincinv, lpmv
 )
 
 from testutils import FuncData
@@ -27,6 +27,11 @@
     return ellipj(k*k)
 def zeta_(x):
     return zeta(x, 1.)
+def assoc_legendre_p_boost_(nu, mu, x):
+    # the boost test data is for integer orders only
+    return lpmv(mu, nu.astype(int), x)
+def legendre_p_via_assoc_(nu, x):
+    return lpmv(0, nu, x)
 
 def test_boost():
     TESTS = [
@@ -39,6 +44,11 @@
         data(arctanh, 'atanh_data_ipp-atanh_data', 0, 1, rtol=1e-11),
         data(arctanh, 'atanh_data_ipp-atanh_data', 0j, 1, rtol=1e-11),
 
+        data(assoc_legendre_p_boost_, 'assoc_legendre_p_ipp-assoc_legendre_p',
+             (0,1,2), 3, rtol=1e-11),
+        data(legendre_p_via_assoc_, 'legendre_p_ipp-legendre_p',
+             (0,1), 2, rtol=1e-11),
+
         data(beta, 'beta_exp_data_ipp-beta_exp_data', (0,1), 2, rtol=1e-13),
         data(beta, 'beta_exp_data_ipp-beta_exp_data', (0,1), 2, rtol=1e-13),
         data(beta, 'beta_small_data_ipp-beta_small_data', (0,1), 2),

Modified: trunk/scipy/special/tests/test_mpmath.py
===================================================================
--- trunk/scipy/special/tests/test_mpmath.py	2010-11-28 17:46:31 UTC (rev 6969)
+++ trunk/scipy/special/tests/test_mpmath.py	2010-11-28 22:55:42 UTC (rev 6970)
@@ -176,3 +176,48 @@
                           vectorized=False, rtol=2e-8)
     finally:
         mpmath.mp.dps, mpmath.mp.prec = old_dps, old_prec
+
+
+
+#------------------------------------------------------------------------------
+# lpmv
+#------------------------------------------------------------------------------
+
+@mpmath_check('0.15')
+def test_lpmv():
+    pts = []
+    for x in [-0.99, -0.557, 1e-6, 0.132, 1]:
+        pts.extend([
+            (1, 1, x),
+            (1, -1, x),
+            (-1, 1, x),
+            (-1, -2, x),
+            (1, 1.7, x),
+            (1, -1.7, x),
+            (-1, 1.7, x),
+            (-1, -2.7, x),
+            (1, 10, x),
+            (1, 11, x),
+            (3, 8, x),
+            (5, 11, x),
+            (-3, 8, x),
+            (-5, 11, x),
+            (3, -8, x),
+            (5, -11, x),
+            (-3, -8, x),
+            (-5, -11, x),
+            (3, 8.3, x),
+            (5, 11.3, x),
+            (-3, 8.3, x),
+            (-5, 11.3, x),
+            (3, -8.3, x),
+            (5, -11.3, x),
+            (-3, -8.3, x),
+            (-5, -11.3, x),
+        ])
+
+    dataset = [p + (mpmath.legenp(p[1], p[0], p[2]),) for p in pts]
+    dataset = np.array(dataset, dtype=np.float_)
+
+    evf = lambda mu,nu,x: sc.lpmv(mu.astype(int), nu, x)
+    FuncData(evf, dataset, (0,1,2), 3, rtol=1e-10, atol=1e-14).check()



More information about the Scipy-svn mailing list