Permalink
Browse files

Merge pull request #489 from pv/mathieu-bug

BUG: special/specfun: fix bug in MTU12
  • Loading branch information...
2 parents 8da73d9 + 8bcdd37 commit 4fe03a73ec5546091e9435f3d0d34b5566bc1c7a @rgommers rgommers committed Mar 29, 2013
Showing with 12 additions and 3 deletions.
  1. +3 −3 scipy/special/specfun/specfun.f
  2. +9 −0 scipy/special/tests/test_basic.py
@@ -12491,7 +12491,7 @@ SUBROUTINE MTU12(KF,KC,M,Q,X,F1R,D1R,F2R,D2R)
QM=17.0+3.1*SQRT(Q)-.126*Q+.0037*SQRT(Q)*Q
ENDIF
KM=INT(QM+0.5*M)
- IF(KM.GT.251) THEN
+ IF(KM.GE.251) THEN
F1R=DNAN()
D1R=DNAN()
F2R=DNAN()
@@ -12505,8 +12505,8 @@ SUBROUTINE MTU12(KF,KC,M,Q,X,F1R,D1R,F2R,D2R)
C2=DEXP(X)
U1=DSQRT(Q)*C1
U2=DSQRT(Q)*C2
- CALL JYNB(KM,U1,NM,BJ1,DJ1,BY1,DY1)
- CALL JYNB(KM,U2,NM,BJ2,DJ2,BY2,DY2)
+ CALL JYNB(KM+1,U1,NM,BJ1,DJ1,BY1,DY1)
+ CALL JYNB(KM+1,U2,NM,BJ2,DJ2,BY2,DY2)
W1=0.0D0
W2=0.0D0
IF (KC.EQ.2) GO TO 50
@@ -455,6 +455,15 @@ def test_mathieu_overflow(self):
assert_equal(cephes.mathieu_modcem2(10000, 1.5, 1.3), (np.nan, np.nan))
assert_equal(cephes.mathieu_modsem2(10000, 1.5, 1.3), (np.nan, np.nan))
+ def test_mathieu_ticket_1847(self):
+ # Regression test --- this call had some out-of-bounds access
+ # and could return nan occasionally
+ for k in range(60):
+ v = cephes.mathieu_modsem2(2, 100, -1)
+ # Values from ACM TOMS 804 (derivate by numerical differentiation)
+ assert_allclose(v[0], 0.1431742913063671074347, rtol=1e-10)
+ assert_allclose(v[1], 0.9017807375832909144719, rtol=1e-4)
+
def test_modfresnelm(self):
cephes.modfresnelm(0)
def test_modfresnelp(self):

0 comments on commit 4fe03a7

Please sign in to comment.