Fix impl of 'sin()' to be equal to newlib libm

JOE1994 · Aug 28, 2020 · a55ad16 · a55ad16
1 parent 98c21c7
commit a55ad16
Show file tree

Hide file tree

Showing 3 changed files with 28 additions and 17 deletions.
diff --git a/src/math/k_cos.rs b/src/math/k_cos.rs
@@ -53,10 +53,28 @@ const C6: f64 = -1.13596475577881948265e-11; /* 0xBDA8FAE9, 0xBE8838D4 */
 //         any extra precision in w.
 #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
 pub(crate) fn k_cos(x: f64, y: f64) -> f64 {
+    let ix = (f64::to_bits(x) >> 32) as u32 & 0x7fffffff;
+    if ix < 0x3e400000			/* |x| < 2**-27 */
+	    { if (x as i32) == 0 { return x; }}		/* generate inexact */
+
     let z = x * x;
-    let w = z * z;
-    let r = z * (C1 + z * (C2 + z * C3)) + w * w * (C4 + z * (C5 + z * C6));
+    let r  = z * (C1 + z * (C2 + z * (C3 + z * (C4 + z * (C5 + z * C6)))));
+    if ix < 0x3FD33333 { 			/* if |x| < 0.3 */ 
+	    return 1.0 - (0.5*z - (z*r - x*y));
+	} else {
+        let qx = if ix > 0x3fe90000 {		/* x > 0.78125 */
+		    0.28125
+	    } else {
+            // INSERT_WORDS(qx,ix-0x00200000,0);	/* x/4 */
+            f64::from_bits(((ix - 0x00200000) as u64) << 32 + 0)
+	    };
+	    let hz = 0.5*z - qx;
+	    let a  = 1.0 - qx;
+	    return a - (hz - (z*r-x*y));
+	}
+
+    /*
     let hz = 0.5 * z;
     let w = 1.0 - hz;
-    w + (((1.0 - w) - hz) + (z * r - x * y))
+    w + (((1.0 - w) - hz) + (z * r - x * y)) */
 }
diff --git a/src/math/k_sin.rs b/src/math/k_sin.rs
@@ -45,10 +45,15 @@ const S6: f64 = 1.58969099521155010221e-10; /* 0x3DE5D93A, 0x5ACFD57C */
 //              sin(x) = x + (S1*x + (x *(r-y/2)+y))
 #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
 pub(crate) fn k_sin(x: f64, y: f64, iy: i32) -> f64 {
+    /* High word of x. */
+    let ix = (f64::to_bits(x) >> 32) as u32 & 0x7fffffff;
+    if ix < 0x3e400000			/* |x| < 2**-27 */
+	    { if (x as i32) == 0 { return x; }}		/* generate inexact */
+
     let z = x * x;
-    let w = z * z;
-    let r = S2 + z * (S3 + z * S4) + z * w * (S5 + z * S6);
     let v = z * x;
+    let r =  S2 + z * (S3 + z * (S4 + z * (S5 + z * S6)));
+
     if iy == 0 {
         x + v * (S1 + z * r)
     } else {

diff --git a/src/math/sin.rs b/src/math/sin.rs
@@ -42,23 +42,11 @@ use super::{k_cos, k_sin, rem_pio2};
 //      TRIG(x) returns trig(x) nearly rounded
 #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
 pub fn sin(x: f64) -> f64 {
-    let x1p120 = f64::from_bits(0x4770000000000000); // 0x1p120f === 2 ^ 120
-
     /* High word of x. */
     let ix = (f64::to_bits(x) >> 32) as u32 & 0x7fffffff;
 
     /* |x| ~< pi/4 */
     if ix <= 0x3fe921fb {
-        if ix < 0x3e500000 {
-            /* |x| < 2**-26 */
-            /* raise inexact if x != 0 and underflow if subnormal*/
-            if ix < 0x00100000 {
-                force_eval!(x / x1p120);
-            } else {
-                force_eval!(x + x1p120);
-            }
-            return x;
-        }
         return k_sin(x, 0.0, 0);
     }