Merge pull request #249 from plugwash/master

rust-lang · Jan 4, 2022 · 2f3fc96 · 2f3fc96
2 parents 10fdf12 + 5e68d37
commit 2f3fc96
Show file tree

Hide file tree

Showing 11 changed files with 189 additions and 16 deletions.
diff --git a/src/math/ceil.rs b/src/math/ceil.rs
@@ -1,3 +1,4 @@
+#![allow(unreachable_code)]
 use core::f64;
 
 const TOINT: f64 = 1. / f64::EPSILON;
@@ -15,6 +16,24 @@ pub fn ceil(x: f64) -> f64 {
             return unsafe { ::core::intrinsics::ceilf64(x) }
         }
     }
+    #[cfg(all(target_arch = "x86", not(target_feature = "sse2")))]
+    {
+        //use an alternative implementation on x86, because the
+        //main implementation fails with the x87 FPU used by
+        //debian i386, probablly due to excess precision issues.
+        //basic implementation taken from https://github.com/rust-lang/libm/issues/219
+        use super::fabs;
+        if fabs(x).to_bits() < 4503599627370496.0_f64.to_bits() {
+            let truncated = x as i64 as f64;
+            if truncated < x {
+                return truncated + 1.0;
+            } else {
+                return truncated;
+            }
+        } else {
+            return x;
+        }
+    }
     let u: u64 = x.to_bits();
     let e: i64 = (u >> 52 & 0x7ff) as i64;
     let y: f64;

diff --git a/src/math/floor.rs b/src/math/floor.rs
@@ -1,3 +1,4 @@
+#![allow(unreachable_code)]
 use core::f64;
 
 const TOINT: f64 = 1. / f64::EPSILON;
@@ -15,6 +16,24 @@ pub fn floor(x: f64) -> f64 {
             return unsafe { ::core::intrinsics::floorf64(x) }
         }
     }
+    #[cfg(all(target_arch = "x86", not(target_feature = "sse2")))]
+    {
+        //use an alternative implementation on x86, because the
+        //main implementation fails with the x87 FPU used by
+        //debian i386, probablly due to excess precision issues.
+        //basic implementation taken from https://github.com/rust-lang/libm/issues/219
+        use super::fabs;
+        if fabs(x).to_bits() < 4503599627370496.0_f64.to_bits() {
+            let truncated = x as i64 as f64;
+            if truncated > x {
+                return truncated - 1.0;
+            } else {
+                return truncated;
+            }
+        } else {
+            return x;
+        }
+    }
     let ui = x.to_bits();
     let e = ((ui >> 52) & 0x7ff) as i32;
 

diff --git a/src/math/fma.rs b/src/math/fma.rs
@@ -218,7 +218,11 @@ mod tests {
             -0.00000000000000022204460492503126,
         );
 
-        assert_eq!(fma(-0.992, -0.992, -0.992), -0.007936000000000007,);
+        let result = fma(-0.992, -0.992, -0.992);
+        //force rounding to storage format on x87 to prevent superious errors.
+        #[cfg(all(target_arch = "x86", not(target_feature = "sse2")))]
+        let result = force_eval!(result);
+        assert_eq!(result, -0.007936000000000007,);
     }
 
     #[test]

diff --git a/src/math/j1f.rs b/src/math/j1f.rs
@@ -369,6 +369,12 @@ mod tests {
     }
     #[test]
     fn test_y1f_2002() {
-        assert_eq!(y1f(2.0000002_f32), -0.10703229_f32);
+        //allow slightly different result on x87
+        let res = y1f(2.0000002_f32);
+        if cfg!(all(target_arch = "x86", not(target_feature = "sse2"))) && (res == -0.10703231_f32)
+        {
+            return;
+        }
+        assert_eq!(res, -0.10703229_f32);
     }
 }
diff --git a/src/math/mod.rs b/src/math/mod.rs
@@ -1,8 +1,6 @@
 macro_rules! force_eval {
     ($e:expr) => {
-        unsafe {
-            ::core::ptr::read_volatile(&$e);
-        }
+        unsafe { ::core::ptr::read_volatile(&$e) }
     };
 }
 

diff --git a/src/math/pow.rs b/src/math/pow.rs
@@ -484,6 +484,10 @@ mod tests {
                 let exp = expected(*val);
                 let res = computed(*val);
 
+                #[cfg(all(target_arch = "x86", not(target_feature = "sse2")))]
+                let exp = force_eval!(exp);
+                #[cfg(all(target_arch = "x86", not(target_feature = "sse2")))]
+                let res = force_eval!(res);
                 assert!(
                     if exp.is_nan() {
                         res.is_nan()

diff --git a/src/math/rem_pio2.rs b/src/math/rem_pio2.rs
@@ -50,7 +50,12 @@ pub(crate) fn rem_pio2(x: f64) -> (i32, f64, f64) {
 
     fn medium(x: f64, ix: u32) -> (i32, f64, f64) {
         /* rint(x/(pi/2)), Assume round-to-nearest. */
-        let f_n = x as f64 * INV_PIO2 + TO_INT - TO_INT;
+        let tmp = x as f64 * INV_PIO2 + TO_INT;
+        // force rounding of tmp to it's storage format on x87 to avoid
+        // excess precision issues.
+        #[cfg(all(target_arch = "x86", not(target_feature = "sse2")))]
+        let tmp = force_eval!(tmp);
+        let f_n = tmp - TO_INT;
         let n = f_n as i32;
         let mut r = x - f_n * PIO2_1;
         let mut w = f_n * PIO2_1T; /* 1st round, good to 85 bits */
@@ -190,20 +195,28 @@ mod tests {
 
     #[test]
     fn test_near_pi() {
+        let arg = 3.141592025756836;
+        let arg = force_eval!(arg);
         assert_eq!(
-            rem_pio2(3.141592025756836),
+            rem_pio2(arg),
             (2, -6.278329573009626e-7, -2.1125998133974653e-23)
         );
+        let arg = 3.141592033207416;
+        let arg = force_eval!(arg);
         assert_eq!(
-            rem_pio2(3.141592033207416),
+            rem_pio2(arg),
             (2, -6.20382377148128e-7, -2.1125998133974653e-23)
         );
+        let arg = 3.141592144966125;
+        let arg = force_eval!(arg);
         assert_eq!(
-            rem_pio2(3.141592144966125),
+            rem_pio2(arg),
             (2, -5.086236681942706e-7, -2.1125998133974653e-23)
         );
+        let arg = 3.141592979431152;
+        let arg = force_eval!(arg);
         assert_eq!(
-            rem_pio2(3.141592979431152),
+            rem_pio2(arg),
             (2, 3.2584135866119817e-7, -2.1125998133974653e-23)
         );
     }

diff --git a/src/math/rem_pio2f.rs b/src/math/rem_pio2f.rs
@@ -43,7 +43,12 @@ pub(crate) fn rem_pio2f(x: f32) -> (i32, f64) {
     if ix < 0x4dc90fdb {
         /* |x| ~< 2^28*(pi/2), medium size */
         /* Use a specialized rint() to get fn.  Assume round-to-nearest. */
-        let f_n = x64 * INV_PIO2 + TOINT - TOINT;
+        let tmp = x64 * INV_PIO2 + TOINT;
+        // force rounding of tmp to it's storage format on x87 to avoid
+        // excess precision issues.
+        #[cfg(all(target_arch = "x86", not(target_feature = "sse2")))]
+        let tmp = force_eval!(tmp);
+        let f_n = tmp - TOINT;
         return (f_n as i32, x64 - f_n * PIO2_1 - f_n * PIO2_1T);
     }
     if ix >= 0x7f800000 {

diff --git a/src/math/sin.rs b/src/math/sin.rs
@@ -81,5 +81,8 @@ pub fn sin(x: f64) -> f64 {
 fn test_near_pi() {
     let x = f64::from_bits(0x400921fb000FD5DD); // 3.141592026217707
     let sx = f64::from_bits(0x3ea50d15ced1a4a2); // 6.273720864039205e-7
-    assert_eq!(sin(x), sx);
+    let result = sin(x);
+    #[cfg(all(target_arch = "x86", not(target_feature = "sse2")))]
+    let result = force_eval!(result);
+    assert_eq!(result, sx);
 }
diff --git a/src/math/sincos.rs b/src/math/sincos.rs
@@ -57,3 +57,77 @@ pub fn sincos(x: f64) -> (f64, f64) {
         _ => (0.0, 1.0),
     }
 }
+
+// These tests are based on those from sincosf.rs
+#[cfg(test)]
+mod tests {
+    use super::sincos;
+
+    const TOLERANCE: f64 = 1e-6;
+
+    #[test]
+    fn with_pi() {
+        let (s, c) = sincos(core::f64::consts::PI);
+        assert!(
+            (s - 0.0).abs() < TOLERANCE,
+            "|{} - {}| = {} >= {}",
+            s,
+            0.0,
+            (s - 0.0).abs(),
+            TOLERANCE
+        );
+        assert!(
+            (c + 1.0).abs() < TOLERANCE,
+            "|{} + {}| = {} >= {}",
+            c,
+            1.0,
+            (s + 1.0).abs(),
+            TOLERANCE
+        );
+    }
+
+    #[test]
+    fn rotational_symmetry() {
+        use core::f64::consts::PI;
+        const N: usize = 24;
+        for n in 0..N {
+            let theta = 2. * PI * (n as f64) / (N as f64);
+            let (s, c) = sincos(theta);
+            let (s_plus, c_plus) = sincos(theta + 2. * PI);
+            let (s_minus, c_minus) = sincos(theta - 2. * PI);
+
+            assert!(
+                (s - s_plus).abs() < TOLERANCE,
+                "|{} - {}| = {} >= {}",
+                s,
+                s_plus,
+                (s - s_plus).abs(),
+                TOLERANCE
+            );
+            assert!(
+                (s - s_minus).abs() < TOLERANCE,
+                "|{} - {}| = {} >= {}",
+                s,
+                s_minus,
+                (s - s_minus).abs(),
+                TOLERANCE
+            );
+            assert!(
+                (c - c_plus).abs() < TOLERANCE,
+                "|{} - {}| = {} >= {}",
+                c,
+                c_plus,
+                (c - c_plus).abs(),
+                TOLERANCE
+            );
+            assert!(
+                (c - c_minus).abs() < TOLERANCE,
+                "|{} - {}| = {} >= {}",
+                c,
+                c_minus,
+                (c - c_minus).abs(),
+                TOLERANCE
+            );
+        }
+    }
+}
diff --git a/src/math/sincosf.rs b/src/math/sincosf.rs
@@ -147,10 +147,38 @@ mod tests {
             let (s_minus, c_minus) = sincosf(theta - 2. * PI);
 
             const TOLERANCE: f32 = 1e-6;
-            assert!((s - s_plus).abs() < TOLERANCE);
-            assert!((s - s_minus).abs() < TOLERANCE);
-            assert!((c - c_plus).abs() < TOLERANCE);
-            assert!((c - c_minus).abs() < TOLERANCE);
+            assert!(
+                (s - s_plus).abs() < TOLERANCE,
+                "|{} - {}| = {} >= {}",
+                s,
+                s_plus,
+                (s - s_plus).abs(),
+                TOLERANCE
+            );
+            assert!(
+                (s - s_minus).abs() < TOLERANCE,
+                "|{} - {}| = {} >= {}",
+                s,
+                s_minus,
+                (s - s_minus).abs(),
+                TOLERANCE
+            );
+            assert!(
+                (c - c_plus).abs() < TOLERANCE,
+                "|{} - {}| = {} >= {}",
+                c,
+                c_plus,
+                (c - c_plus).abs(),
+                TOLERANCE
+            );
+            assert!(
+                (c - c_minus).abs() < TOLERANCE,
+                "|{} - {}| = {} >= {}",
+                c,
+                c_minus,
+                (c - c_minus).abs(),
+                TOLERANCE
+            );
         }
     }
 }