user0/rust
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Merge rust-lang/libm#80

Browse source 
80: Add sin and cos r=japaric a=porglezomp

Closes rust-lang/libm#11 
Closes rust-lang/libm#33

Co-authored-by: C Jones <code@calebjones.net>
This commit is contained in:
bors[bot] 2018-07-14 23:08:33 +00:00
commit 486a0066e5
9 changed files with 476 additions and 11 deletions
Download patch file Download diff file Expand all files Collapse all files

4
library/compiler-builtins/libm/src/lib.rs
View file

@ -386,10 +386,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)]
@ -548,13 +546,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)

72
library/compiler-builtins/libm/src/math/cos.rs Normal file
View file

@ -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),
}
}

2
library/compiler-builtins/libm/src/math/exp2.rs
View file

@ -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;

62
library/compiler-builtins/libm/src/math/k_cos.rs Normal file
View file

@ -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))
}

57
library/compiler-builtins/libm/src/math/k_sin.rs Normal file
View file

@ -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)
}
}

22
library/compiler-builtins/libm/src/math/mod.rs
View file

@ -6,6 +6,7 @@ macro_rules! force_eval {
};
}
// Public modules
mod acos;
mod acosf;
mod asin;
@ -15,6 +16,7 @@ mod cbrt;
mod cbrtf;
mod ceil;
mod ceilf;
mod cos;
mod cosf;
mod exp;
mod exp2;
@ -46,6 +48,7 @@ mod round;
mod roundf;
mod scalbn;
mod scalbnf;
mod sin;
mod sinf;
mod sqrt;
mod sqrtf;
@ -63,6 +66,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::exp::exp;
pub use self::exp2::exp2;
@ -94,6 +98,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;
@ -101,16 +106,25 @@ pub use self::tanf::tanf;
pub use self::trunc::trunc;
pub use self::truncf::truncf;
// Private modules
mod k_cos;
mod k_cosf;
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_sinf::k_sinf, k_tanf::k_tanf, rem_pio2_large::rem_pio2_large,
rem_pio2f::rem_pio2f,
};
// Private re-imports
use self::k_cos::k_cos;
use self::k_cosf::k_cosf;
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 {

187
library/compiler-builtins/libm/src/math/rem_pio2.rs Normal file
View file

@ -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]);
}

77
library/compiler-builtins/libm/src/math/sin.rs Normal file
View file

@ -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),
}
}

4
library/compiler-builtins/libm/test-generator/src/main.rs
View file

@ -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,