Add forced rounding to storage format for x87 to rem_pio2.rs as well.

Author: Peter Michael Green

src/math/rem_pio2.rs

@@ -50,7 +50,11 @@ 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")))]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 */

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
+            );
+        }
+    }
+}