diff --git a/library/compiler-builtins/libm/src/lib.rs b/library/compiler-builtins/libm/src/lib.rs index 969e7c45b5e4..a06aeffe2480 100644 --- a/library/compiler-builtins/libm/src/lib.rs +++ b/library/compiler-builtins/libm/src/lib.rs @@ -380,10 +380,8 @@ pub trait F64Ext: private::Sealed { fn hypot(self, other: Self) -> Self; - #[cfg(todo)] fn sin(self) -> Self; - #[cfg(todo)] fn cos(self) -> Self; #[cfg(todo)] @@ -542,13 +540,11 @@ impl F64Ext for f64 { hypot(self, other) } - #[cfg(todo)] #[inline] fn sin(self) -> Self { sin(self) } - #[cfg(todo)] #[inline] fn cos(self) -> Self { cos(self) diff --git a/library/compiler-builtins/libm/src/math/cos.rs b/library/compiler-builtins/libm/src/math/cos.rs new file mode 100644 index 000000000000..e6e9b373674c --- /dev/null +++ b/library/compiler-builtins/libm/src/math/cos.rs @@ -0,0 +1,72 @@ +// origin: FreeBSD /usr/src/lib/msun/src/s_cos.c */ +// +// ==================================================== +// Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. +// +// Developed at SunPro, a Sun Microsystems, Inc. business. +// Permission to use, copy, modify, and distribute this +// software is freely granted, provided that this notice +// is preserved. +// ==================================================== + +use super::{k_cos, k_sin, rem_pio2}; + +// cos(x) +// Return cosine function of x. +// +// kernel function: +// k_sin ... sine function on [-pi/4,pi/4] +// k_cos ... cosine function on [-pi/4,pi/4] +// rem_pio2 ... argument reduction routine +// +// Method. +// Let S,C and T denote the sin, cos and tan respectively on +// [-PI/4, +PI/4]. Reduce the argument x to y1+y2 = x-k*pi/2 +// in [-pi/4 , +pi/4], and let n = k mod 4. +// We have +// +// n sin(x) cos(x) tan(x) +// ---------------------------------------------------------- +// 0 S C T +// 1 C -S -1/T +// 2 -S -C T +// 3 -C S -1/T +// ---------------------------------------------------------- +// +// Special cases: +// Let trig be any of sin, cos, or tan. +// trig(+-INF) is NaN, with signals; +// trig(NaN) is that NaN; +// +// Accuracy: +// TRIG(x) returns trig(x) nearly rounded +// +pub fn cos(x: f64) -> f64 { + let ix = (f64::to_bits(x) >> 32) as u32 & 0x7fffffff; + + /* |x| ~< pi/4 */ + if ix <= 0x3fe921fb { + if ix < 0x3e46a09e { + /* if x < 2**-27 * sqrt(2) */ + /* raise inexact if x != 0 */ + if x as i32 == 0 { + return 1.0; + } + } + return k_cos(x, 0.0); + } + + /* cos(Inf or NaN) is NaN */ + if ix >= 0x7ff00000 { + return x - x; + } + + /* argument reduction needed */ + let (n, y0, y1) = rem_pio2(x); + match n & 3 { + 0 => k_cos(y0, y1), + 1 => -k_sin(y0, y1, 1), + 2 => -k_cos(y0, y1), + _ => k_sin(y0, y1, 1), + } +} diff --git a/library/compiler-builtins/libm/src/math/exp2.rs b/library/compiler-builtins/libm/src/math/exp2.rs index 0ab11989679b..3952e93007e1 100644 --- a/library/compiler-builtins/libm/src/math/exp2.rs +++ b/library/compiler-builtins/libm/src/math/exp2.rs @@ -24,7 +24,7 @@ // OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF // SUCH DAMAGE. -use super::scalbn::scalbn; +use super::scalbn; const TBLSIZE: usize = 256; diff --git a/library/compiler-builtins/libm/src/math/k_cos.rs b/library/compiler-builtins/libm/src/math/k_cos.rs new file mode 100644 index 000000000000..693950d1d418 --- /dev/null +++ b/library/compiler-builtins/libm/src/math/k_cos.rs @@ -0,0 +1,62 @@ +// origin: FreeBSD /usr/src/lib/msun/src/k_cos.c +// +// ==================================================== +// Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. +// +// Developed at SunSoft, a Sun Microsystems, Inc. business. +// Permission to use, copy, modify, and distribute this +// software is freely granted, provided that this notice +// is preserved. +// ==================================================== + +const C1: f64 = 4.16666666666666019037e-02; /* 0x3FA55555, 0x5555554C */ +const C2: f64 = -1.38888888888741095749e-03; /* 0xBF56C16C, 0x16C15177 */ +const C3: f64 = 2.48015872894767294178e-05; /* 0x3EFA01A0, 0x19CB1590 */ +const C4: f64 = -2.75573143513906633035e-07; /* 0xBE927E4F, 0x809C52AD */ +const C5: f64 = 2.08757232129817482790e-09; /* 0x3E21EE9E, 0xBDB4B1C4 */ +const C6: f64 = -1.13596475577881948265e-11; /* 0xBDA8FAE9, 0xBE8838D4 */ + +// kernel cos function on [-pi/4, pi/4], pi/4 ~ 0.785398164 +// Input x is assumed to be bounded by ~pi/4 in magnitude. +// Input y is the tail of x. +// +// Algorithm +// 1. Since cos(-x) = cos(x), we need only to consider positive x. +// 2. if x < 2^-27 (hx<0x3e400000 0), return 1 with inexact if x!=0. +// 3. cos(x) is approximated by a polynomial of degree 14 on +// [0,pi/4] +// 4 14 +// cos(x) ~ 1 - x*x/2 + C1*x + ... + C6*x +// where the remez error is +// +// | 2 4 6 8 10 12 14 | -58 +// |cos(x)-(1-.5*x +C1*x +C2*x +C3*x +C4*x +C5*x +C6*x )| <= 2 +// | | +// +// 4 6 8 10 12 14 +// 4. let r = C1*x +C2*x +C3*x +C4*x +C5*x +C6*x , then +// cos(x) ~ 1 - x*x/2 + r +// since cos(x+y) ~ cos(x) - sin(x)*y +// ~ cos(x) - x*y, +// a correction term is necessary in cos(x) and hence +// cos(x+y) = 1 - (x*x/2 - (r - x*y)) +// For better accuracy, rearrange to +// cos(x+y) ~ w + (tmp + (r-x*y)) +// where w = 1 - x*x/2 and tmp is a tiny correction term +// (1 - x*x/2 == w + tmp exactly in infinite precision). +// The exactness of w + tmp in infinite precision depends on w +// and tmp having the same precision as x. If they have extra +// precision due to compiler bugs, then the extra precision is +// only good provided it is retained in all terms of the final +// expression for cos(). Retention happens in all cases tested +// under FreeBSD, so don't pessimize things by forcibly clipping +// any extra precision in w. +#[inline] +pub fn k_cos(x: f64, y: f64) -> f64 { + let z = x * x; + let w = z * z; + let r = z * (C1 + z * (C2 + z * C3)) + w * w * (C4 + z * (C5 + z * C6)); + let hz = 0.5 * z; + let w = 1.0 - hz; + w + (((1.0 - w) - hz) + (z * r - x * y)) +} diff --git a/library/compiler-builtins/libm/src/math/k_sin.rs b/library/compiler-builtins/libm/src/math/k_sin.rs new file mode 100644 index 000000000000..3e07c3594475 --- /dev/null +++ b/library/compiler-builtins/libm/src/math/k_sin.rs @@ -0,0 +1,57 @@ +// origin: FreeBSD /usr/src/lib/msun/src/k_sin.c +// +// ==================================================== +// Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. +// +// Developed at SunSoft, a Sun Microsystems, Inc. business. +// Permission to use, copy, modify, and distribute this +// software is freely granted, provided that this notice +// is preserved. +// ==================================================== + +const S1: f64 = -1.66666666666666324348e-01; /* 0xBFC55555, 0x55555549 */ +const S2: f64 = 8.33333333332248946124e-03; /* 0x3F811111, 0x1110F8A6 */ +const S3: f64 = -1.98412698298579493134e-04; /* 0xBF2A01A0, 0x19C161D5 */ +const S4: f64 = 2.75573137070700676789e-06; /* 0x3EC71DE3, 0x57B1FE7D */ +const S5: f64 = -2.50507602534068634195e-08; /* 0xBE5AE5E6, 0x8A2B9CEB */ +const S6: f64 = 1.58969099521155010221e-10; /* 0x3DE5D93A, 0x5ACFD57C */ + +// kernel sin function on ~[-pi/4, pi/4] (except on -0), pi/4 ~ 0.7854 +// Input x is assumed to be bounded by ~pi/4 in magnitude. +// Input y is the tail of x. +// Input iy indicates whether y is 0. (if iy=0, y assume to be 0). +// +// Algorithm +// 1. Since sin(-x) = -sin(x), we need only to consider positive x. +// 2. Callers must return sin(-0) = -0 without calling here since our +// odd polynomial is not evaluated in a way that preserves -0. +// Callers may do the optimization sin(x) ~ x for tiny x. +// 3. sin(x) is approximated by a polynomial of degree 13 on +// [0,pi/4] +// 3 13 +// sin(x) ~ x + S1*x + ... + S6*x +// where +// +// |sin(x) 2 4 6 8 10 12 | -58 +// |----- - (1+S1*x +S2*x +S3*x +S4*x +S5*x +S6*x )| <= 2 +// | x | +// +// 4. sin(x+y) = sin(x) + sin'(x')*y +// ~ sin(x) + (1-x*x/2)*y +// For better accuracy, let +// 3 2 2 2 2 +// r = x *(S2+x *(S3+x *(S4+x *(S5+x *S6)))) +// then 3 2 +// sin(x) = x + (S1*x + (x *(r-y/2)+y)) +#[inline] +pub fn k_sin(x: f64, y: f64, iy: i32) -> f64 { + let z = x * x; + let w = z * z; + let r = S2 + z * (S3 + z * S4) + z * w * (S5 + z * S6); + let v = z * x; + if iy == 0 { + x + v * (S1 + z * r) + } else { + x - ((z * (0.5 * y - v * r) - y) - v * S1) + } +} diff --git a/library/compiler-builtins/libm/src/math/mod.rs b/library/compiler-builtins/libm/src/math/mod.rs index 1ac5472576b2..3def31082a72 100644 --- a/library/compiler-builtins/libm/src/math/mod.rs +++ b/library/compiler-builtins/libm/src/math/mod.rs @@ -6,6 +6,7 @@ macro_rules! force_eval { }; } +// Public modules mod acos; mod acosf; mod asin; @@ -16,6 +17,7 @@ mod cbrt; mod cbrtf; mod ceil; mod ceilf; +mod cos; mod cosf; mod coshf; mod exp; @@ -48,6 +50,7 @@ mod round; mod roundf; mod scalbn; mod scalbnf; +mod sin; mod sinf; mod sqrt; mod sqrtf; @@ -67,6 +70,7 @@ pub use self::cbrt::cbrt; pub use self::cbrtf::cbrtf; pub use self::ceil::ceil; pub use self::ceilf::ceilf; +pub use self::cos::cos; pub use self::cosf::cosf; pub use self::coshf::coshf; pub use self::exp::exp; @@ -99,6 +103,7 @@ pub use self::round::round; pub use self::roundf::roundf; pub use self::scalbn::scalbn; pub use self::scalbnf::scalbnf; +pub use self::sin::sin; pub use self::sinf::sinf; pub use self::sqrt::sqrt; pub use self::sqrtf::sqrtf; @@ -107,17 +112,26 @@ pub use self::tanhf::tanhf; pub use self::trunc::trunc; pub use self::truncf::truncf; +// Private modules +mod k_cos; mod k_cosf; mod k_expo2f; +mod k_sin; mod k_sinf; mod k_tanf; +mod rem_pio2; mod rem_pio2_large; mod rem_pio2f; -use self::{ - k_cosf::k_cosf, k_expo2f::k_expo2f, k_sinf::k_sinf, k_tanf::k_tanf, - rem_pio2_large::rem_pio2_large, rem_pio2f::rem_pio2f, -}; +use self::k_cos::k_cos; +use self::k_cosf::k_cosf; +use self::k_expo2f::k_expo2f; +use self::k_sin::k_sin; +use self::k_sinf::k_sinf; +use self::k_tanf::k_tanf; +use self::rem_pio2::rem_pio2; +use self::rem_pio2_large::rem_pio2_large; +use self::rem_pio2f::rem_pio2f; #[inline] pub fn get_high_word(x: f64) -> u32 { diff --git a/library/compiler-builtins/libm/src/math/rem_pio2.rs b/library/compiler-builtins/libm/src/math/rem_pio2.rs new file mode 100644 index 000000000000..47bab6e65f19 --- /dev/null +++ b/library/compiler-builtins/libm/src/math/rem_pio2.rs @@ -0,0 +1,187 @@ +// origin: FreeBSD /usr/src/lib/msun/src/e_rem_pio2.c +// +// ==================================================== +// Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. +// +// Developed at SunPro, a Sun Microsystems, Inc. business. +// Permission to use, copy, modify, and distribute this +// software is freely granted, provided that this notice +// is preserved. +// ==================================================== +// +// Optimized by Bruce D. Evans. */ + +use super::rem_pio2_large; + +// #if FLT_EVAL_METHOD==0 || FLT_EVAL_METHOD==1 +// #define EPS DBL_EPSILON +const EPS: f64 = 2.2204460492503131e-16; +// #elif FLT_EVAL_METHOD==2 +// #define EPS LDBL_EPSILON +// #endif + +// TODO: Support FLT_EVAL_METHOD? + +const TO_INT: f64 = 1.5 / EPS; +/// 53 bits of 2/pi +const INV_PIO2: f64 = 6.36619772367581382433e-01; /* 0x3FE45F30, 0x6DC9C883 */ +/// first 33 bits of pi/2 +const PIO2_1: f64 = 1.57079632673412561417e+00; /* 0x3FF921FB, 0x54400000 */ +/// pi/2 - PIO2_1 +const PIO2_1T: f64 = 6.07710050650619224932e-11; /* 0x3DD0B461, 0x1A626331 */ +/// second 33 bits of pi/2 +const PIO2_2: f64 = 6.07710050630396597660e-11; /* 0x3DD0B461, 0x1A600000 */ +/// pi/2 - (PIO2_1+PIO2_2) +const PIO2_2T: f64 = 2.02226624879595063154e-21; /* 0x3BA3198A, 0x2E037073 */ +/// third 33 bits of pi/2 +const PIO2_3: f64 = 2.02226624871116645580e-21; /* 0x3BA3198A, 0x2E000000 */ +/// pi/2 - (PIO2_1+PIO2_2+PIO2_3) +const PIO2_3T: f64 = 8.47842766036889956997e-32; /* 0x397B839A, 0x252049C1 */ + +// return the remainder of x rem pi/2 in y[0]+y[1] +// use rem_pio2_large() for large x +// +// caller must handle the case when reduction is not needed: |x| ~<= pi/4 */ +pub fn rem_pio2(x: f64) -> (i32, f64, f64) { + let x1p24 = f64::from_bits(0x4170000000000000); + + let sign = (f64::to_bits(x) >> 63) as i32; + let ix = (f64::to_bits(x) >> 32) as u32 & 0x7fffffff; + + 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 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 */ + let mut y0 = r - w; + let ui = f64::to_bits(y0); + let ey = (ui >> 52) as i32 & 0x7ff; + let ex = (ix >> 20) as i32; + if ex - ey > 16 { + /* 2nd round, good to 118 bits */ + let t = r; + w = f_n * PIO2_2; + r = t - w; + w = f_n * PIO2_2T - ((t - r) - w); + y0 = r - w; + let ey = (f64::to_bits(y0) >> 52) as i32 & 0x7ff; + if ex - ey > 49 { + /* 3rd round, good to 151 bits, covers all cases */ + let t = r; + w = f_n * PIO2_3; + r = t - w; + w = f_n * PIO2_3T - ((t - r) - w); + y0 = r - w; + } + } + let y1 = (r - y0) - w; + return (n, y0, y1); + } + + if ix <= 0x400f6a7a { + /* |x| ~<= 5pi/4 */ + if (ix & 0xfffff) == 0x921fb { + /* |x| ~= pi/2 or 2pi/2 */ + medium(x, ix); /* cancellation -- use medium case */ + } + if ix <= 0x4002d97c { + /* |x| ~<= 3pi/4 */ + if sign == 0 { + let z = x - PIO2_1; /* one round good to 85 bits */ + let y0 = z - PIO2_1T; + let y1 = (z - y0) - PIO2_1T; + return (1, y0, y1); + } else { + let z = x + PIO2_1; + let y0 = z + PIO2_1T; + let y1 = (z - y0) + PIO2_1T; + return (-1, y0, y1); + } + } else { + if sign == 0 { + let z = x - 2.0 * PIO2_1; + let y0 = z - 2.0 * PIO2_1T; + let y1 = (z - y0) - 2.0 * PIO2_1T; + return (2, y0, y1); + } else { + let z = x + 2.0 * PIO2_1; + let y0 = z + 2.0 * PIO2_1T; + let y1 = (z - y0) + 2.0 * PIO2_1T; + return (-2, y0, y1); + } + } + } + if ix <= 0x401c463b { + /* |x| ~<= 9pi/4 */ + if ix <= 0x4015fdbc { + /* |x| ~<= 7pi/4 */ + if ix == 0x4012d97c { + /* |x| ~= 3pi/2 */ + return medium(x, ix); + } + if sign == 0 { + let z = x - 3.0 * PIO2_1; + let y0 = z - 3.0 * PIO2_1T; + let y1 = (z - y0) - 3.0 * PIO2_1T; + return (3, y0, y1); + } else { + let z = x + 3.0 * PIO2_1; + let y0 = z + 3.0 * PIO2_1T; + let y1 = (z - y0) + 3.0 * PIO2_1T; + return (-3, y0, y1); + } + } else { + if ix == 0x401921fb { + /* |x| ~= 4pi/2 */ + return medium(x, ix); + } + if sign == 0 { + let z = x - 4.0 * PIO2_1; + let y0 = z - 4.0 * PIO2_1T; + let y1 = (z - y0) - 4.0 * PIO2_1T; + return (4, y0, y1); + } else { + let z = x + 4.0 * PIO2_1; + let y0 = z + 4.0 * PIO2_1T; + let y1 = (z - y0) + 4.0 * PIO2_1T; + return (-4, y0, y1); + } + } + } + if ix < 0x413921fb { + /* |x| ~< 2^20*(pi/2), medium size */ + return medium(x, ix); + } + /* + * all other (large) arguments + */ + if ix >= 0x7ff00000 { + /* x is inf or NaN */ + let y0 = x - x; + let y1 = y0; + return (0, y0, y1); + } + /* set z = scalbn(|x|,-ilogb(x)+23) */ + let mut ui = f64::to_bits(x); + ui &= (!1) >> 12; + ui |= (0x3ff + 23) << 52; + let mut z = f64::from_bits(ui); + let mut tx = [0.0; 3]; + for i in 0..2 { + tx[i] = z as i32 as f64; + z = (z - tx[i]) * x1p24; + } + tx[2] = z; + /* skip zero terms, first term is non-zero */ + let mut i = 2; + while tx[i] == 0.0 { + i -= 1; + } + let mut ty = [0.0; 3]; + let n = rem_pio2_large(&tx[..=i], &mut ty, ((ix >> 20) - (0x3ff + 23)) as i32, 1); + if sign != 0 { + return (-n, -ty[0], -ty[1]); + } + return (n, ty[0], ty[1]); +} diff --git a/library/compiler-builtins/libm/src/math/sin.rs b/library/compiler-builtins/libm/src/math/sin.rs new file mode 100644 index 000000000000..13eb30248ab7 --- /dev/null +++ b/library/compiler-builtins/libm/src/math/sin.rs @@ -0,0 +1,77 @@ +// origin: FreeBSD /usr/src/lib/msun/src/s_sin.c */ +// +// ==================================================== +// Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. +// +// Developed at SunPro, a Sun Microsystems, Inc. business. +// Permission to use, copy, modify, and distribute this +// software is freely granted, provided that this notice +// is preserved. +// ==================================================== + +use super::{k_cos, k_sin, rem_pio2}; + +// sin(x) +// Return sine function of x. +// +// kernel function: +// k_sin ... sine function on [-pi/4,pi/4] +// k_cos ... cose function on [-pi/4,pi/4] +// rem_pio2 ... argument reduction routine +// +// Method. +// Let S,C and T denote the sin, cos and tan respectively on +// [-PI/4, +PI/4]. Reduce the argument x to y1+y2 = x-k*pi/2 +// in [-pi/4 , +pi/4], and let n = k mod 4. +// We have +// +// n sin(x) cos(x) tan(x) +// ---------------------------------------------------------- +// 0 S C T +// 1 C -S -1/T +// 2 -S -C T +// 3 -C S -1/T +// ---------------------------------------------------------- +// +// Special cases: +// Let trig be any of sin, cos, or tan. +// trig(+-INF) is NaN, with signals; +// trig(NaN) is that NaN; +// +// Accuracy: +// TRIG(x) returns trig(x) nearly rounded +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); + } + + /* sin(Inf or NaN) is NaN */ + if ix >= 0x7ff00000 { + return x - x; + } + + /* argument reduction needed */ + let (n, y0, y1) = rem_pio2(x); + match n & 3 { + 0 => k_sin(y0, y1, 1), + 1 => k_cos(y0, y1), + 2 => -k_sin(y0, y1, 1), + _ => -k_cos(y0, y1), + } +} diff --git a/library/compiler-builtins/libm/test-generator/src/main.rs b/library/compiler-builtins/libm/test-generator/src/main.rs index 1c261b2ebafe..7cc299079a47 100644 --- a/library/compiler-builtins/libm/test-generator/src/main.rs +++ b/library/compiler-builtins/libm/test-generator/src/main.rs @@ -702,7 +702,7 @@ f64_f64! { // atan, cbrt, ceil, - // cos, + cos, // cosh, exp, exp2, @@ -713,7 +713,7 @@ f64_f64! { log1p, log2, round, - // sin, + sin, // sinh, sqrt, // tan,