[do not merge] Add test case for #170

src/math/jn.rs

@@ -341,3 +341,18 @@ pub fn yn(n: i32, x: f64) -> f64 {
         b
     }
 }
+
+#[cfg(test)]
+mod test {
+
+    #[test]
+    pub fn issue_170() {
+        let expected = f64::from_bits(1617394868955);
+        let result = super::jn(114, f64::from_bits(4594974205335009568));
+        assert_eq!(result, expected);
+
+        let expected = f64::from_bits(15924833);
+        let result = super::jn(136, f64::from_bits(4602429132083530282));
+        assert_eq!(result, expected);
+    }
+}

src/math/jnf.rs

@@ -169,10 +169,10 @@ pub fn jnf(n: i32, mut x: f32) -> f32 {
                     b = 2.0 * (i as f32) * b / x - a;
                     a = temp;
                     /* scale b to avoid spurious overflow */
-                    let x1p60 = f32::from_bits(0x5d800000); // 0x1p60 == 2^60
-                    if b > x1p60 {
+                    // let x1p60 = f32::from_bits(0x5d800000); // 0x1p60 == 2^60
+                    if b > 1.152922e+18 {
                         a /= b;
-                        t /= b;
+                        t = t / b;
                         b = 1.0;
                     }
                     i -= 1;
@@ -257,3 +257,14 @@ pub fn ynf(n: i32, x: f32) -> f32 {
         b
     }
 }
+
+#[cfg(test)]
+mod test {
+
+    #[test]
+    pub fn issue_170() {
+        let expected = f32::from_bits(86);
+        let result = super::jnf(29, f32::from_bits(1061546867));
+        assert_eq!(result, expected);
+    }
+}