Implement part of sin/cos with quadrant selection unimplemented
This commit is contained in:
C Jones
parent
802f6a2b46
commit
440a1b7273
5 changed files with 452 additions and 8 deletions
|
|
@ -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)
|
||||
|
|
|
|||
73
library/compiler-builtins/libm/src/math/cos.rs
Normal file
73
library/compiler-builtins/libm/src/math/cos.rs
Normal file
|
|
@ -0,0 +1,73 @@
|
|||
/* 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};
|
||||
|
||||
// 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
|
||||
//
|
||||
// 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 _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 => _cos(y0, y1),
|
||||
1 => -_sin(y0, y1, 1),
|
||||
2 => -_cos(y0, y1),
|
||||
_ => _sin(y0, y1, 1),
|
||||
}
|
||||
}
|
||||
|
|
@ -15,6 +15,7 @@ mod cbrt;
|
|||
mod cbrtf;
|
||||
mod ceil;
|
||||
mod ceilf;
|
||||
mod cos;
|
||||
mod cosf;
|
||||
mod exp;
|
||||
mod exp2;
|
||||
|
|
@ -46,6 +47,7 @@ mod round;
|
|||
mod roundf;
|
||||
mod scalbn;
|
||||
mod scalbnf;
|
||||
mod sin;
|
||||
mod sinf;
|
||||
mod sqrt;
|
||||
mod sqrtf;
|
||||
|
|
@ -63,6 +65,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 +97,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;
|
||||
|
|
@ -106,11 +110,13 @@ mod k_sinf;
|
|||
mod k_tanf;
|
||||
mod rem_pio2_large;
|
||||
mod rem_pio2f;
|
||||
mod trig_common;
|
||||
|
||||
use self::{
|
||||
k_cosf::k_cosf, k_sinf::k_sinf, k_tanf::k_tanf, rem_pio2_large::rem_pio2_large,
|
||||
rem_pio2f::rem_pio2f,
|
||||
};
|
||||
use self::k_cosf::k_cosf;
|
||||
use self::k_sinf::k_sinf;
|
||||
use self::k_tanf::k_tanf;
|
||||
use self::rem_pio2_large::rem_pio2_large;
|
||||
use self::rem_pio2f::rem_pio2f;
|
||||
|
||||
#[inline]
|
||||
pub fn get_high_word(x: f64) -> u32 {
|
||||
|
|
|
|||
80
library/compiler-builtins/libm/src/math/sin.rs
Normal file
80
library/compiler-builtins/libm/src/math/sin.rs
Normal file
|
|
@ -0,0 +1,80 @@
|
|||
/* 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};
|
||||
|
||||
// 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
|
||||
//
|
||||
// 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 _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 => _sin(y0, y1, 1),
|
||||
1 => _cos(y0, y1),
|
||||
2 => -_sin(y0, y1, 1),
|
||||
_ => -_cos(y0, y1),
|
||||
}
|
||||
}
|
||||
289
library/compiler-builtins/libm/src/math/trig_common.rs
Normal file
289
library/compiler-builtins/libm/src/math/trig_common.rs
Normal file
|
|
@ -0,0 +1,289 @@
|
|||
/*
|
||||
* ====================================================
|
||||
* 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!()
|
||||
}
|
||||
Loading…
Add table
Add a link
Reference in a new issue