Implement cosh, exp and expm1
This commit is contained in:
Lukas Wirth
parent
965c736cc1
commit
3a15b30d16
5 changed files with 273 additions and 7 deletions
54
library/compiler-builtins/libm/src/math/cosh.rs
Normal file
54
library/compiler-builtins/libm/src/math/cosh.rs
Normal file
|
|
@ -0,0 +1,54 @@
|
|||
use core::f64;
|
||||
|
||||
use super::exp;
|
||||
use super::expm1;
|
||||
|
||||
pub fn cosh(mut x: f64) -> f64 {
|
||||
let t: f64;
|
||||
/* |x| */
|
||||
let mut ui = x.to_bits();
|
||||
ui &= !0u64;
|
||||
x = f64::from_bits(ui);
|
||||
let w = (ui >> 32) as u32;
|
||||
|
||||
/* |x| < log(2) */
|
||||
if w < 0x3fe62e42 {
|
||||
if w < 0x3ff00000 - (26 << 20) {
|
||||
/* raise inexact if x!=0 */
|
||||
force_eval!(x + f64::from_bits(0x4770000000000000));
|
||||
return 1.0;
|
||||
}
|
||||
let t = expm1(x);
|
||||
return 1.0 + t * t / (2.0 * (1.0 + t));
|
||||
}
|
||||
|
||||
/* |x| < log(DBL_MAX) */
|
||||
if w < 0x40862e42 {
|
||||
t = exp(x);
|
||||
/* note: if x>log(0x1p26) then the 1/t is not needed */
|
||||
return 0.5 * (t + 1.0 / t);
|
||||
}
|
||||
|
||||
/* |x| > log(DBL_MAX) or nan */
|
||||
/* note: the result is stored to handle overflow */
|
||||
t = __expo2(x);
|
||||
return t;
|
||||
}
|
||||
|
||||
const K: u32 = 2043;
|
||||
|
||||
pub fn __expo2(x: f64) -> f64 {
|
||||
let kln2 = f64::from_bits(0x40962066151add8b);
|
||||
/* note that k is odd and scale*scale overflows */
|
||||
let scale = f64::from_bits(((0x3ff + K / 2) << 20) as u64);
|
||||
/* exp(x - k ln2) * 2**(k-1) */
|
||||
return exp(x - kln2) * scale * scale;
|
||||
}
|
||||
|
||||
#[cfg(test)]
|
||||
mod tests {
|
||||
#[test]
|
||||
fn sanity_check() {
|
||||
assert_eq!(super::cosh(1.1), 1.6685185538222564);
|
||||
}
|
||||
}
|
||||
84
library/compiler-builtins/libm/src/math/exp.rs
Normal file
84
library/compiler-builtins/libm/src/math/exp.rs
Normal file
|
|
@ -0,0 +1,84 @@
|
|||
use super::scalbn;
|
||||
|
||||
const HALF: [f64; 2] = [0.5, -0.5];
|
||||
const LN2_HI: f64 = 6.93147180369123816490e-01; /* 0x3fe62e42, 0xfee00000 */
|
||||
const LN2_LO: f64 = 1.90821492927058770002e-10; /* 0x3dea39ef, 0x35793c76 */
|
||||
const INV_LN2: f64 = 1.44269504088896338700e+00; /* 0x3ff71547, 0x652b82fe */
|
||||
const P1: f64 = 1.66666666666666019037e-01; /* 0x3FC55555, 0x5555553E */
|
||||
const P2: f64 = -2.77777777770155933842e-03; /* 0xBF66C16C, 0x16BEBD93 */
|
||||
const P3: f64 = 6.61375632143793436117e-05; /* 0x3F11566A, 0xAF25DE2C */
|
||||
const P4: f64 = -1.65339022054652515390e-06; /* 0xBEBBBD41, 0xC5D26BF1 */
|
||||
const P5: f64 = 4.13813679705723846039e-08; /* 0x3E663769, 0x72BEA4D0 */
|
||||
|
||||
#[inline]
|
||||
pub fn exp(mut x: f64) -> f64 {
|
||||
let mut hx: u32 = (x.to_bits() >> 32) as u32;
|
||||
let sign = (hx >> 31) as i32; /* sign bit of x */
|
||||
hx &= 0x7fffffff; /* high word of |x| */
|
||||
|
||||
/* special cases */
|
||||
if hx >= 0x4086232b {
|
||||
/* if |x| >= 708.39... */
|
||||
if x.is_nan() {
|
||||
return x;
|
||||
}
|
||||
if x > 709.782712893383973096 {
|
||||
/* overflow if x!=inf */
|
||||
x *= f64::from_bits(0x7fe0000000000000);
|
||||
return x;
|
||||
}
|
||||
if x < -708.39641853226410622 {
|
||||
/* underflow if x!=-inf */
|
||||
force_eval!((f64::from_bits(0xb6a0000000000000) / x) as f32);
|
||||
if x < -745.13321910194110842 {
|
||||
return 0.0;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/* argument reduction */
|
||||
let k: i32;
|
||||
let hi: f64;
|
||||
let lo: f64;
|
||||
if hx > 0x3fd62e42 {
|
||||
/* if |x| > 0.5 ln2 */
|
||||
/* if |x| > 0.5 ln2 */
|
||||
if hx > 0x3ff0a2b2 {
|
||||
/* if |x| > 1.5 ln2 */
|
||||
k = (INV_LN2 * x + HALF[sign as usize]) as i32;
|
||||
} else {
|
||||
k = 1 - sign - sign;
|
||||
}
|
||||
let kf = k as f64;
|
||||
hi = x - kf * LN2_HI; /* k*ln2hi is exact here */
|
||||
lo = kf * LN2_LO;
|
||||
x = hi - lo;
|
||||
} else if hx > 0x3e300000 {
|
||||
/* |x| > 2**-14 */
|
||||
k = 0;
|
||||
hi = x;
|
||||
lo = 0.0;
|
||||
} else {
|
||||
/* raise inexact */
|
||||
force_eval!(f64::from_bits(0x7fe0000000000000) + x);
|
||||
return 1.0 + x;
|
||||
}
|
||||
|
||||
/* x is now in primary range */
|
||||
let xx = x * x;
|
||||
let c = x - xx * (P1 + xx * (P2 + xx * (P3 + xx * (P4 + xx * P5))));
|
||||
let y = 1.0 + (x * c / (2.0 - c) - lo + hi);
|
||||
if k == 0 {
|
||||
y
|
||||
} else {
|
||||
scalbn(y, k)
|
||||
}
|
||||
}
|
||||
|
||||
#[cfg(test)]
|
||||
mod tests {
|
||||
#[test]
|
||||
fn sanity_check() {
|
||||
assert_eq!(super::exp(1.1), 3.0041660239464334);
|
||||
}
|
||||
}
|
||||
124
library/compiler-builtins/libm/src/math/expm1.rs
Normal file
124
library/compiler-builtins/libm/src/math/expm1.rs
Normal file
|
|
@ -0,0 +1,124 @@
|
|||
use core::f64;
|
||||
|
||||
const O_THRESHOLD: f64 = 7.09782712893383973096e+02; /* 0x40862E42, 0xFEFA39EF */
|
||||
const LN2_HI: f64 = 6.93147180369123816490e-01; /* 0x3fe62e42, 0xfee00000 */
|
||||
const LN2_LO: f64 = 1.90821492927058770002e-10; /* 0x3dea39ef, 0x35793c76 */
|
||||
const INVLN2: f64 = 1.44269504088896338700e+00; /* 0x3ff71547, 0x652b82fe */
|
||||
/* Scaled Q's: Qn_here = 2**n * Qn_above, for R(2*z) where z = hxs = x*x/2: */
|
||||
const Q1: f64 = -3.33333333333331316428e-02; /* BFA11111 111110F4 */
|
||||
const Q2: f64 = 1.58730158725481460165e-03; /* 3F5A01A0 19FE5585 */
|
||||
const Q3: f64 = -7.93650757867487942473e-05; /* BF14CE19 9EAADBB7 */
|
||||
const Q4: f64 = 4.00821782732936239552e-06; /* 3ED0CFCA 86E65239 */
|
||||
const Q5: f64 = -2.01099218183624371326e-07; /* BE8AFDB7 6E09C32D */
|
||||
|
||||
pub fn expm1(mut x: f64) -> f64 {
|
||||
let hi: f64;
|
||||
let lo: f64;
|
||||
let k: i32;
|
||||
let c: f64;
|
||||
let mut t: f64;
|
||||
let mut y: f64;
|
||||
|
||||
let mut ui = x.to_bits() >> 32;
|
||||
let hx = ui & 0x7fffffff;
|
||||
let sign = (ui >> 63) as i32;
|
||||
|
||||
/* filter out huge and non-finite argument */
|
||||
if hx >= 0x4043687A {
|
||||
/* if |x|>=56*ln2 */
|
||||
if x.is_nan() {
|
||||
return x;
|
||||
}
|
||||
if sign != 0 {
|
||||
return -1.0;
|
||||
}
|
||||
if x > O_THRESHOLD {
|
||||
x *= f64::from_bits(0x7fe0000000000000);
|
||||
return x;
|
||||
}
|
||||
}
|
||||
|
||||
/* argument reduction */
|
||||
if hx > 0x3fd62e42 {
|
||||
/* if |x| > 0.5 ln2 */
|
||||
if hx < 0x3FF0A2B2 {
|
||||
/* and |x| < 1.5 ln2 */
|
||||
if sign == 0 {
|
||||
hi = x - LN2_HI;
|
||||
lo = LN2_LO;
|
||||
k = 1;
|
||||
} else {
|
||||
hi = x + LN2_HI;
|
||||
lo = -LN2_LO;
|
||||
k = -1;
|
||||
}
|
||||
} else {
|
||||
k = (INVLN2 * x + if sign != 0 { -0.5 } else { 0.5 }) as i32;
|
||||
t = k as f64;
|
||||
hi = x - t * LN2_HI; /* t*ln2_hi is exact here */
|
||||
lo = t * LN2_LO;
|
||||
}
|
||||
x = hi - lo;
|
||||
c = (hi - x) - lo;
|
||||
} else if hx < 0x3c900000 {
|
||||
/* |x| < 2**-54, return x */
|
||||
if hx < 0x00100000 {
|
||||
force_eval!(x as f32);
|
||||
}
|
||||
return x;
|
||||
} else {
|
||||
c = 0.0;
|
||||
k = 0;
|
||||
}
|
||||
|
||||
/* x is now in primary range */
|
||||
let hfx = 0.5 * x;
|
||||
let hxs = x * hfx;
|
||||
let r1 = 1.0 + hxs * (Q1 + hxs * (Q2 + hxs * (Q3 + hxs * (Q4 + hxs * Q5))));
|
||||
t = 3.0 - r1 * hfx;
|
||||
let mut e = hxs * ((r1 - t) / (6.0 - x * t));
|
||||
if k == 0 {
|
||||
/* c is 0 */
|
||||
return x - (x * e - hxs);
|
||||
}
|
||||
e = x * (e - c) - c;
|
||||
e -= hxs;
|
||||
/* exp(x) ~ 2^k (x_reduced - e + 1) */
|
||||
if k == -1 {
|
||||
return 0.5 * (x - e) - 0.5;
|
||||
}
|
||||
if k == 1 {
|
||||
if x < -0.25 {
|
||||
return -2.0 * (e - (x + 0.5));
|
||||
}
|
||||
return 1.0 + 2.0 * (x - e);
|
||||
}
|
||||
ui = ((0x3ff + k) as u64) << 52; /* 2^k */
|
||||
let twopk = f64::from_bits(ui);
|
||||
if k < 0 || k > 56 {
|
||||
/* suffice to return exp(x)-1 */
|
||||
y = x - e + 1.0;
|
||||
if k == 1024 {
|
||||
y = y * 2.0 * f64::from_bits(0x7fe0000000000000);
|
||||
} else {
|
||||
y = y * twopk;
|
||||
}
|
||||
return y - 1.0;
|
||||
}
|
||||
ui = ((0x3ff - k) as u64) << 52; /* 2^-k */
|
||||
let uf = f64::from_bits(ui);
|
||||
if k < 20 {
|
||||
y = (x - e + (1.0 - uf)) * twopk;
|
||||
} else {
|
||||
y = (x - (e + uf) + 1.0) * twopk;
|
||||
}
|
||||
y
|
||||
}
|
||||
|
||||
#[cfg(test)]
|
||||
mod tests {
|
||||
#[test]
|
||||
fn sanity_check() {
|
||||
assert_eq!(super::expm1(1.1), 2.0041660239464334);
|
||||
}
|
||||
}
|
||||
|
|
@ -7,7 +7,10 @@ macro_rules! force_eval {
|
|||
}
|
||||
|
||||
mod ceilf;
|
||||
mod cosh;
|
||||
mod exp;
|
||||
mod expf;
|
||||
mod expm1;
|
||||
mod fabs;
|
||||
mod fabsf;
|
||||
mod floor;
|
||||
|
|
@ -36,10 +39,11 @@ mod truncf;
|
|||
//mod service;
|
||||
|
||||
pub use self::{
|
||||
ceilf::ceilf, expf::expf, fabs::fabs, fabsf::fabsf, floor::floor, floorf::floorf, fmodf::fmodf,
|
||||
hypot::hypot, hypotf::hypotf, log::log, log10::log10, log10f::log10f, log1p::log1p,
|
||||
log1pf::log1pf, log2::log2, log2f::log2f, logf::logf, powf::powf, round::round, roundf::roundf,
|
||||
scalbn::scalbn, scalbnf::scalbnf, sqrt::sqrt, sqrtf::sqrtf, trunc::trunc, truncf::truncf,
|
||||
ceilf::ceilf, cosh::cosh, exp::exp, expf::expf, expm1::expm1, fabs::fabs, fabsf::fabsf,
|
||||
floor::floor, floorf::floorf, fmodf::fmodf, hypot::hypot, hypotf::hypotf, log::log,
|
||||
log10::log10, log10f::log10f, log1p::log1p, log1pf::log1pf, log2::log2, log2f::log2f,
|
||||
logf::logf, powf::powf, round::round, roundf::roundf, scalbn::scalbn, scalbnf::scalbnf,
|
||||
sqrt::sqrt, sqrtf::sqrtf, trunc::trunc, truncf::truncf,
|
||||
};
|
||||
|
||||
fn isnanf(x: f32) -> bool {
|
||||
|
|
|
|||
|
|
@ -702,10 +702,10 @@ f64_f64! {
|
|||
// cbrt,
|
||||
// ceil,
|
||||
// cos,
|
||||
// cosh,
|
||||
// exp,
|
||||
cosh,
|
||||
exp,
|
||||
// exp2,
|
||||
// expm1,
|
||||
expm1,
|
||||
floor,
|
||||
log,
|
||||
log10,
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue