Permalink
Browse files

BUG: special/specfun: fix bug in MTU12

Off-by-one error in the order up to which the Bessel functions were
computed (the arrays are accessed up to KM+1). Doesn't seem to affect
results much as the series seems mostly to converge before reaching the
limit.

Covered by test_mathieu_modsem2 (see Trac #1847)

(backport of 8bcdd37)
  • Loading branch information...
1 parent e6cf3d3 commit 0476caeda2966dc451e421511ad6281ca7c719a7 @pv pv committed with rgommers Mar 28, 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
View
6 scipy/special/specfun/specfun.f
@@ -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
View
9 scipy/special/tests/test_basic.py
@@ -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 0476cae

Please sign in to comment.