Merge branch 'master' into master
This commit is contained in:
Zgarbul Andrey
GitHub
commit
6b7bd7fc3d
9 changed files with 476 additions and 11 deletions
|
|
@ -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)
|
||||
|
|
|
|||
72
library/compiler-builtins/libm/src/math/cos.rs
Normal file
72
library/compiler-builtins/libm/src/math/cos.rs
Normal 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),
|
||||
}
|
||||
}
|
||||
|
|
@ -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
62
library/compiler-builtins/libm/src/math/k_cos.rs
Normal 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
57
library/compiler-builtins/libm/src/math/k_sin.rs
Normal 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)
|
||||
}
|
||||
}
|
||||
|
|
@ -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 {
|
||||
|
|
|
|||
187
library/compiler-builtins/libm/src/math/rem_pio2.rs
Normal file
187
library/compiler-builtins/libm/src/math/rem_pio2.rs
Normal 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
77
library/compiler-builtins/libm/src/math/sin.rs
Normal 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),
|
||||
}
|
||||
}
|
||||
|
|
@ -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,
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue