Skip to content
This repository

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

Merged
merged 1 commit into from about 1 year ago

2 participants

Pauli Virtanen Ralf Gommers
Pauli Virtanen
Owner
pv commented March 28, 2013

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.

Closes Trac #1847

Pauli Virtanen 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)
8bcdd37
Ralf Gommers
Owner

Thanks @pv. Tests pass, so merging.

Ralf Gommers rgommers merged commit 4fe03a7 into from March 29, 2013
Ralf Gommers rgommers closed this March 29, 2013
Ralf Gommers
Owner

backported to 0.12.x in 0476cae

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Showing 1 unique commit by 1 author.

Mar 29, 2013
Pauli Virtanen 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)
8bcdd37
This page is out of date. Refresh to see the latest.
6  scipy/special/specfun/specfun.f
@@ -12491,7 +12491,7 @@ SUBROUTINE MTU12(KF,KC,M,Q,X,F1R,D1R,F2R,D2R)
12491 12491
            QM=17.0+3.1*SQRT(Q)-.126*Q+.0037*SQRT(Q)*Q
12492 12492
         ENDIF
12493 12493
         KM=INT(QM+0.5*M)
12494  
-        IF(KM.GT.251) THEN
  12494
+        IF(KM.GE.251) THEN
12495 12495
            F1R=DNAN()
12496 12496
            D1R=DNAN()
12497 12497
            F2R=DNAN()
@@ -12505,8 +12505,8 @@ SUBROUTINE MTU12(KF,KC,M,Q,X,F1R,D1R,F2R,D2R)
12505 12505
         C2=DEXP(X)
12506 12506
         U1=DSQRT(Q)*C1
12507 12507
         U2=DSQRT(Q)*C2
12508  
-        CALL JYNB(KM,U1,NM,BJ1,DJ1,BY1,DY1)
12509  
-        CALL JYNB(KM,U2,NM,BJ2,DJ2,BY2,DY2)
  12508
+        CALL JYNB(KM+1,U1,NM,BJ1,DJ1,BY1,DY1)
  12509
+        CALL JYNB(KM+1,U2,NM,BJ2,DJ2,BY2,DY2)
12510 12510
         W1=0.0D0
12511 12511
         W2=0.0D0
12512 12512
         IF (KC.EQ.2) GO TO 50
9  scipy/special/tests/test_basic.py
@@ -455,6 +455,15 @@ def test_mathieu_overflow(self):
455 455
         assert_equal(cephes.mathieu_modcem2(10000, 1.5, 1.3), (np.nan, np.nan))
456 456
         assert_equal(cephes.mathieu_modsem2(10000, 1.5, 1.3), (np.nan, np.nan))
457 457
 
  458
+    def test_mathieu_ticket_1847(self):
  459
+        # Regression test --- this call had some out-of-bounds access
  460
+        # and could return nan occasionally
  461
+        for k in range(60):
  462
+            v = cephes.mathieu_modsem2(2, 100, -1)
  463
+            # Values from ACM TOMS 804 (derivate by numerical differentiation)
  464
+            assert_allclose(v[0], 0.1431742913063671074347, rtol=1e-10)
  465
+            assert_allclose(v[1], 0.9017807375832909144719, rtol=1e-4)
  466
+
458 467
     def test_modfresnelm(self):
459 468
         cephes.modfresnelm(0)
460 469
     def test_modfresnelp(self):
Commit_comment_tip

Tip: You can add notes to lines in a file. Hover to the left of a line to make a note

Something went wrong with that request. Please try again.