Convert rem_pio2 code, split up modules
This commit is contained in:
C Jones
parent
440a1b7273
commit
b2751897d2
9 changed files with 356 additions and 335 deletions
|
|
@ -1,24 +1,23 @@
|
|||
/* 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.
|
||||
* ====================================================
|
||||
*/
|
||||
// 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 math::trig_common::{_cos, _sin, rem_pio2};
|
||||
use super::{k_cos, k_sin, rem_pio2};
|
||||
|
||||
// cos(x)
|
||||
// Return cosine function of x.
|
||||
//
|
||||
// kernel function:
|
||||
// __sin ... sine function on [-pi/4,pi/4]
|
||||
// __cos ... cosine function on [-pi/4,pi/4]
|
||||
// __rem_pio2 ... argument reduction routine
|
||||
// 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
|
||||
|
|
@ -54,7 +53,7 @@ pub fn cos(x: f64) -> f64 {
|
|||
return 1.0;
|
||||
}
|
||||
}
|
||||
return _cos(x, 0.0);
|
||||
return k_cos(x, 0.0);
|
||||
}
|
||||
|
||||
/* cos(Inf or NaN) is NaN */
|
||||
|
|
@ -65,9 +64,9 @@ pub fn cos(x: f64) -> f64 {
|
|||
/* argument reduction needed */
|
||||
let (n, y0, y1) = rem_pio2(x);
|
||||
match n & 3 {
|
||||
0 => _cos(y0, y1),
|
||||
1 => -_sin(y0, y1, 1),
|
||||
2 => -_cos(y0, y1),
|
||||
_ => _sin(y0, y1, 1),
|
||||
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;
|
||||
|
|
@ -105,16 +106,23 @@ 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;
|
||||
mod trig_common;
|
||||
|
||||
// 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;
|
||||
|
||||
|
|
|
|||
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(0x7041);
|
||||
|
||||
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]);
|
||||
}
|
||||
|
|
@ -1,26 +1,23 @@
|
|||
/* 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.
|
||||
* ====================================================
|
||||
*/
|
||||
// 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 core::f64;
|
||||
|
||||
use math::trig_common::{_cos, _sin, rem_pio2};
|
||||
use super::{k_cos, k_sin, rem_pio2};
|
||||
|
||||
// sin(x)
|
||||
// Return sine function of x.
|
||||
//
|
||||
// kernel function:
|
||||
// __sin ... sine function on [-pi/4,pi/4]
|
||||
// __cos ... cose function on [-pi/4,pi/4]
|
||||
// __rem_pio2 ... argument reduction routine
|
||||
// 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
|
||||
|
|
@ -61,7 +58,7 @@ pub fn sin(x: f64) -> f64 {
|
|||
}
|
||||
return x;
|
||||
}
|
||||
return _sin(x, 0.0, 0);
|
||||
return k_sin(x, 0.0, 0);
|
||||
}
|
||||
|
||||
/* sin(Inf or NaN) is NaN */
|
||||
|
|
@ -72,9 +69,9 @@ pub fn sin(x: f64) -> f64 {
|
|||
/* argument reduction needed */
|
||||
let (n, y0, y1) = rem_pio2(x);
|
||||
match n & 3 {
|
||||
0 => _sin(y0, y1, 1),
|
||||
1 => _cos(y0, y1),
|
||||
2 => -_sin(y0, y1, 1),
|
||||
_ => -_cos(y0, y1),
|
||||
0 => k_sin(y0, y1, 1),
|
||||
1 => k_cos(y0, y1),
|
||||
2 => -k_sin(y0, y1, 1),
|
||||
_ => -k_cos(y0, y1),
|
||||
}
|
||||
}
|
||||
|
|
|
|||
|
|
@ -1,289 +0,0 @@
|
|||
/*
|
||||
* ====================================================
|
||||
* 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.
|
||||
* ====================================================
|
||||
*/
|
||||
|
||||
/* origin: FreeBSD /usr/src/lib/msun/src/k_sin.c */
|
||||
|
||||
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))
|
||||
pub fn _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)
|
||||
}
|
||||
}
|
||||
|
||||
/* origin: FreeBSD /usr/src/lib/msun/src/k_cos.c */
|
||||
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.
|
||||
pub fn _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))
|
||||
}
|
||||
|
||||
// origin: FreeBSD /usr/src/lib/msun/src/e_rem_pio2.c
|
||||
// Optimized by Bruce D. Evans. */
|
||||
|
||||
// #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?
|
||||
|
||||
#[allow(unused, non_upper_case_globals)]
|
||||
const toint: f64 = 1.5 / EPS;
|
||||
/// 53 bits of 2/pi
|
||||
#[allow(unused, non_upper_case_globals)]
|
||||
const invpio2: f64 = 6.36619772367581382433e-01; /* 0x3FE45F30, 0x6DC9C883 */
|
||||
/// first 33 bits of pi/2
|
||||
#[allow(unused, non_upper_case_globals)]
|
||||
const pio2_1: f64 = 1.57079632673412561417e+00; /* 0x3FF921FB, 0x54400000 */
|
||||
/// pi/2 - pio2_1
|
||||
#[allow(unused, non_upper_case_globals)]
|
||||
const pio2_1t: f64 = 6.07710050650619224932e-11; /* 0x3DD0B461, 0x1A626331 */
|
||||
/// second 33 bits of pi/2
|
||||
#[allow(unused, non_upper_case_globals)]
|
||||
const pio2_2: f64 = 6.07710050630396597660e-11; /* 0x3DD0B461, 0x1A600000 */
|
||||
/// pi/2 - (pio2_1+pio2_2)
|
||||
#[allow(unused, non_upper_case_globals)]
|
||||
const pio2_2t: f64 = 2.02226624879595063154e-21; /* 0x3BA3198A, 0x2E037073 */
|
||||
/// third 33 bits of pi/2
|
||||
#[allow(unused, non_upper_case_globals)]
|
||||
const pio2_3: f64 = 2.02226624871116645580e-21; /* 0x3BA3198A, 0x2E000000 */
|
||||
/// pi/2 - (pio2_1+pio2_2+pio2_3)
|
||||
#[allow(unused, non_upper_case_globals)]
|
||||
const pio2_3t: f64 = 8.47842766036889956997e-32; /* 0x397B839A, 0x252049C1 */
|
||||
|
||||
/* __rem_pio2(x,y)
|
||||
*
|
||||
* 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 */
|
||||
/*
|
||||
{
|
||||
union {double f; uint64_t i;} u = {x};
|
||||
double_t z,w,t,r,fn;
|
||||
double tx[3],ty[2];
|
||||
uint32_t ix;
|
||||
int sign, n, ex, ey, i;
|
||||
|
||||
sign = u.i>>63;
|
||||
ix = u.i>>32 & 0x7fffffff;
|
||||
if (ix <= 0x400f6a7a) { /* |x| ~<= 5pi/4 */
|
||||
if ((ix & 0xfffff) == 0x921fb) /* |x| ~= pi/2 or 2pi/2 */
|
||||
goto medium; /* cancellation -- use medium case */
|
||||
if (ix <= 0x4002d97c) { /* |x| ~<= 3pi/4 */
|
||||
if (!sign) {
|
||||
z = x - pio2_1; /* one round good to 85 bits */
|
||||
y[0] = z - pio2_1t;
|
||||
y[1] = (z-y[0]) - pio2_1t;
|
||||
return 1;
|
||||
} else {
|
||||
z = x + pio2_1;
|
||||
y[0] = z + pio2_1t;
|
||||
y[1] = (z-y[0]) + pio2_1t;
|
||||
return -1;
|
||||
}
|
||||
} else {
|
||||
if (!sign) {
|
||||
z = x - 2*pio2_1;
|
||||
y[0] = z - 2*pio2_1t;
|
||||
y[1] = (z-y[0]) - 2*pio2_1t;
|
||||
return 2;
|
||||
} else {
|
||||
z = x + 2*pio2_1;
|
||||
y[0] = z + 2*pio2_1t;
|
||||
y[1] = (z-y[0]) + 2*pio2_1t;
|
||||
return -2;
|
||||
}
|
||||
}
|
||||
}
|
||||
if (ix <= 0x401c463b) { /* |x| ~<= 9pi/4 */
|
||||
if (ix <= 0x4015fdbc) { /* |x| ~<= 7pi/4 */
|
||||
if (ix == 0x4012d97c) /* |x| ~= 3pi/2 */
|
||||
goto medium;
|
||||
if (!sign) {
|
||||
z = x - 3*pio2_1;
|
||||
y[0] = z - 3*pio2_1t;
|
||||
y[1] = (z-y[0]) - 3*pio2_1t;
|
||||
return 3;
|
||||
} else {
|
||||
z = x + 3*pio2_1;
|
||||
y[0] = z + 3*pio2_1t;
|
||||
y[1] = (z-y[0]) + 3*pio2_1t;
|
||||
return -3;
|
||||
}
|
||||
} else {
|
||||
if (ix == 0x401921fb) /* |x| ~= 4pi/2 */
|
||||
goto medium;
|
||||
if (!sign) {
|
||||
z = x - 4*pio2_1;
|
||||
y[0] = z - 4*pio2_1t;
|
||||
y[1] = (z-y[0]) - 4*pio2_1t;
|
||||
return 4;
|
||||
} else {
|
||||
z = x + 4*pio2_1;
|
||||
y[0] = z + 4*pio2_1t;
|
||||
y[1] = (z-y[0]) + 4*pio2_1t;
|
||||
return -4;
|
||||
}
|
||||
}
|
||||
}
|
||||
if (ix < 0x413921fb) { /* |x| ~< 2^20*(pi/2), medium size */
|
||||
medium:
|
||||
/* rint(x/(pi/2)), Assume round-to-nearest. */
|
||||
fn = (double_t)x*invpio2 + toint - toint;
|
||||
n = (int32_t)fn;
|
||||
r = x - fn*pio2_1;
|
||||
w = fn*pio2_1t; /* 1st round, good to 85 bits */
|
||||
y[0] = r - w;
|
||||
u.f = y[0];
|
||||
ey = u.i>>52 & 0x7ff;
|
||||
ex = ix>>20;
|
||||
if (ex - ey > 16) { /* 2nd round, good to 118 bits */
|
||||
t = r;
|
||||
w = fn*pio2_2;
|
||||
r = t - w;
|
||||
w = fn*pio2_2t - ((t-r)-w);
|
||||
y[0] = r - w;
|
||||
u.f = y[0];
|
||||
ey = u.i>>52 & 0x7ff;
|
||||
if (ex - ey > 49) { /* 3rd round, good to 151 bits, covers all cases */
|
||||
t = r;
|
||||
w = fn*pio2_3;
|
||||
r = t - w;
|
||||
w = fn*pio2_3t - ((t-r)-w);
|
||||
y[0] = r - w;
|
||||
}
|
||||
}
|
||||
y[1] = (r - y[0]) - w;
|
||||
return n;
|
||||
}
|
||||
/*
|
||||
* all other (large) arguments
|
||||
*/
|
||||
if (ix >= 0x7ff00000) { /* x is inf or NaN */
|
||||
y[0] = y[1] = x - x;
|
||||
return 0;
|
||||
}
|
||||
/* set z = scalbn(|x|,-ilogb(x)+23) */
|
||||
u.f = x;
|
||||
u.i &= (uint64_t)-1>>12;
|
||||
u.i |= (uint64_t)(0x3ff + 23)<<52;
|
||||
z = u.f;
|
||||
for (i=0; i < 2; i++) {
|
||||
tx[i] = (double)(int32_t)z;
|
||||
z = (z-tx[i])*0x1p24;
|
||||
}
|
||||
tx[i] = z;
|
||||
/* skip zero terms, first term is non-zero */
|
||||
while (tx[i] == 0.0)
|
||||
i--;
|
||||
n = __rem_pio2_large(tx,ty,(int)(ix>>20)-(0x3ff+23),i+1,1);
|
||||
if (sign) {
|
||||
y[0] = -ty[0];
|
||||
y[1] = -ty[1];
|
||||
return -n;
|
||||
}
|
||||
y[0] = ty[0];
|
||||
y[1] = ty[1];
|
||||
return n;
|
||||
}
|
||||
*/
|
||||
|
||||
pub fn rem_pio2(_x: f64) -> (i32, f64, f64) {
|
||||
unimplemented!()
|
||||
}
|
||||
|
|
@ -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