diff --git a/src/math/k_cos.rs b/src/math/k_cos.rs index 49b2fc64..3ab516bd 100644 --- 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 index 9dd96c94..e1d06106 100644 --- 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 index 1329b41a..889ae046 100644 --- 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); }