Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Already on GitHub? Sign in to your account

BUG: special/specfun: fix bug in MTU12 #489

Merged
merged 1 commit into from Mar 29, 2013
Jump to file or symbol
Failed to load files and symbols.
+12 −3
Split
@@ -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):