[Scipy-svn] r6970 - in trunk/scipy/special: . specfun tests
scipy-svn at scipy.org
scipy-svn at scipy.org
Sun Nov 28 17:55:42 EST 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
+#------------------------------------------------------------------------------
+
+ at 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