diff --git a/library/compiler-builtins/libm/.travis.yml b/library/compiler-builtins/libm/.travis.yml index 2bbb60fc64d1..223b586ea6da 100644 --- a/library/compiler-builtins/libm/.travis.yml +++ b/library/compiler-builtins/libm/.travis.yml @@ -35,6 +35,5 @@ before_cache: branches: only: - - master - staging - trying diff --git a/library/compiler-builtins/libm/src/lib.rs b/library/compiler-builtins/libm/src/lib.rs index c53fefb2f297..8091079a6241 100644 --- a/library/compiler-builtins/libm/src/lib.rs +++ b/library/compiler-builtins/libm/src/lib.rs @@ -42,6 +42,8 @@ pub trait F32Ext: private::Sealed { fn trunc(self) -> Self; + fn fdim(self, rhs: Self) -> Self; + #[cfg(todo)] fn fract(self) -> Self; @@ -155,6 +157,11 @@ impl F32Ext for f32 { truncf(self) } + #[inline] + fn fdim(self, rhs: Self) -> Self { + fdimf(self, rhs) + } + #[cfg(todo)] #[inline] fn fract(self) -> Self { @@ -353,6 +360,8 @@ pub trait F64Ext: private::Sealed { fn trunc(self) -> Self; + fn fdim(self, rhs: Self) -> Self; + #[cfg(todo)] fn fract(self) -> Self; @@ -378,7 +387,6 @@ pub trait F64Ext: private::Sealed { fn sqrt(self) -> Self; - #[cfg(todo)] fn exp(self) -> Self; #[cfg(todo)] @@ -467,6 +475,10 @@ impl F64Ext for f64 { trunc(self) } + #[inline] + fn fdim(self, rhs: Self) -> Self { + fdim(self, rhs) + } #[cfg(todo)] #[inline] fn fract(self) -> Self { @@ -516,7 +528,6 @@ impl F64Ext for f64 { sqrt(self) } - #[cfg(todo)] #[inline] fn exp(self) -> Self { exp(self) diff --git a/library/compiler-builtins/libm/src/math/exp.rs b/library/compiler-builtins/libm/src/math/exp.rs new file mode 100644 index 000000000000..cd63b8fb3f15 --- /dev/null +++ b/library/compiler-builtins/libm/src/math/exp.rs @@ -0,0 +1,150 @@ +/* origin: FreeBSD /usr/src/lib/msun/src/e_exp.c */ +/* + * ==================================================== + * Copyright (C) 2004 by Sun Microsystems, Inc. All rights reserved. + * + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ +/* exp(x) + * Returns the exponential of x. + * + * Method + * 1. Argument reduction: + * Reduce x to an r so that |r| <= 0.5*ln2 ~ 0.34658. + * Given x, find r and integer k such that + * + * x = k*ln2 + r, |r| <= 0.5*ln2. + * + * Here r will be represented as r = hi-lo for better + * accuracy. + * + * 2. Approximation of exp(r) by a special rational function on + * the interval [0,0.34658]: + * Write + * R(r**2) = r*(exp(r)+1)/(exp(r)-1) = 2 + r*r/6 - r**4/360 + ... + * We use a special Remez algorithm on [0,0.34658] to generate + * a polynomial of degree 5 to approximate R. The maximum error + * of this polynomial approximation is bounded by 2**-59. In + * other words, + * R(z) ~ 2.0 + P1*z + P2*z**2 + P3*z**3 + P4*z**4 + P5*z**5 + * (where z=r*r, and the values of P1 to P5 are listed below) + * and + * | 5 | -59 + * | 2.0+P1*z+...+P5*z - R(z) | <= 2 + * | | + * The computation of exp(r) thus becomes + * 2*r + * exp(r) = 1 + ---------- + * R(r) - r + * r*c(r) + * = 1 + r + ----------- (for better accuracy) + * 2 - c(r) + * where + * 2 4 10 + * c(r) = r - (P1*r + P2*r + ... + P5*r ). + * + * 3. Scale back to obtain exp(x): + * From step 1, we have + * exp(x) = 2^k * exp(r) + * + * Special cases: + * exp(INF) is INF, exp(NaN) is NaN; + * exp(-INF) is 0, and + * for finite argument, only exp(0)=1 is exact. + * + * Accuracy: + * according to an error analysis, the error is always less than + * 1 ulp (unit in the last place). + * + * Misc. info. + * For IEEE double + * if x > 709.782712893383973096 then exp(x) overflows + * if x < -745.133219101941108420 then exp(x) underflows + */ + +use super::scalbn; + +const HALF: [f64; 2] = [0.5, -0.5]; +const LN2HI: f64 = 6.93147180369123816490e-01; /* 0x3fe62e42, 0xfee00000 */ +const LN2LO: f64 = 1.90821492927058770002e-10; /* 0x3dea39ef, 0x35793c76 */ +const INVLN2: 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 x1p1023 = f64::from_bits(0x7fe0000000000000); // 0x1p1023 === 2 ^ 1023 + let x1p_149 = f64::from_bits(0x36a0000000000000); // 0x1p-149 === 2 ^ -149 + + let hi: f64; + let lo: f64; + let c: f64; + let xx: f64; + let y: f64; + let k: i32; + let sign: i32; + let mut hx: u32; + + hx = (x.to_bits() >> 32) as u32; + sign = (hx >> 31) as i32; + 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 *= x1p1023; + return x; + } + if x < -708.39641853226410622 { + /* underflow if x!=-inf */ + force_eval!((-x1p_149 / x) as f32); + if x < -745.13321910194110842 { + return 0.; + } + } + } + + /* argument reduction */ + if hx > 0x3fd62e42 { + /* if |x| > 0.5 ln2 */ + if hx >= 0x3ff0a2b2 { + /* if |x| >= 1.5 ln2 */ + k = (INVLN2 * x + HALF[sign as usize]) as i32; + } else { + k = 1 - sign - sign; + } + hi = x - k as f64 * LN2HI; /* k*ln2hi is exact here */ + lo = k as f64 * LN2LO; + x = hi - lo; + } else if hx > 0x3e300000 { + /* if |x| > 2**-28 */ + k = 0; + hi = x; + lo = 0.; + } else { + /* inexact if x!=0 */ + force_eval!(x1p1023 + x); + return 1. + x; + } + + /* x is now in primary range */ + xx = x * x; + c = x - xx * (P1 + xx * (P2 + xx * (P3 + xx * (P4 + xx * P5)))); + y = 1. + (x * c / (2. - c) - lo + hi); + if k == 0 { + y + } else { + scalbn(y, k) + } +} diff --git a/library/compiler-builtins/libm/src/math/fdim.rs b/library/compiler-builtins/libm/src/math/fdim.rs new file mode 100644 index 000000000000..2b277eab03e1 --- /dev/null +++ b/library/compiler-builtins/libm/src/math/fdim.rs @@ -0,0 +1,15 @@ +use core::f64; + +pub fn fdim(x: f64, y: f64) -> f64 { + if x.is_nan() { + x + } else if y.is_nan() { + y + } else { + if x > y { + x - y + } else { + 0.0 + } + } +} diff --git a/library/compiler-builtins/libm/src/math/fdimf.rs b/library/compiler-builtins/libm/src/math/fdimf.rs new file mode 100644 index 000000000000..44bf2d6806ef --- /dev/null +++ b/library/compiler-builtins/libm/src/math/fdimf.rs @@ -0,0 +1,15 @@ +use core::f32; + +pub fn fdimf(x: f32, y: f32) -> f32 { + if x.is_nan() { + x + } else if y.is_nan() { + y + } else { + if x > y { + x - y + } else { + 0.0 + } + } +} diff --git a/library/compiler-builtins/libm/src/math/fmodf.rs b/library/compiler-builtins/libm/src/math/fmodf.rs index 909775249085..d84cfeb01b83 100644 --- a/library/compiler-builtins/libm/src/math/fmodf.rs +++ b/library/compiler-builtins/libm/src/math/fmodf.rs @@ -1,7 +1,6 @@ +use core::f32; use core::u32; -use super::isnanf; - #[inline] pub fn fmodf(x: f32, y: f32) -> f32 { let mut uxi = x.to_bits(); @@ -11,7 +10,7 @@ pub fn fmodf(x: f32, y: f32) -> f32 { let sx = uxi & 0x80000000; let mut i; - if uyi << 1 == 0 || isnanf(y) || ex == 0xff { + if uyi << 1 == 0 || y.is_nan() || ex == 0xff { return (x * y) / (x * y); } diff --git a/library/compiler-builtins/libm/src/math/mod.rs b/library/compiler-builtins/libm/src/math/mod.rs index be17295659a8..86c181d13236 100644 --- a/library/compiler-builtins/libm/src/math/mod.rs +++ b/library/compiler-builtins/libm/src/math/mod.rs @@ -10,10 +10,13 @@ mod ceil; mod ceilf; mod cosf; mod cosh; +mod exp; mod expf; mod expm1; mod fabs; mod fabsf; +mod fdim; +mod fdimf; mod floor; mod floorf; mod fmodf; @@ -42,10 +45,13 @@ pub use self::ceil::ceil; pub use self::ceilf::ceilf; pub use self::cosf::cosf; pub use self::cosh::cosh; +pub use self::exp::exp; pub use self::expf::expf; pub use self::expm1::expm1; pub use self::fabs::fabs; pub use self::fabsf::fabsf; +pub use self::fdim::fdim; +pub use self::fdimf::fdimf; pub use self::floor::floor; pub use self::floorf::floorf; pub use self::fmodf::fmodf; @@ -79,7 +85,3 @@ use self::{ k_cosf::k_cosf, k_expo2::k_expo2, k_sinf::k_sinf, rem_pio2_large::rem_pio2_large, rem_pio2f::rem_pio2f, }; - -fn isnanf(x: f32) -> bool { - x.to_bits() & 0x7fffffff > 0x7f800000 -} diff --git a/library/compiler-builtins/libm/test-generator/src/main.rs b/library/compiler-builtins/libm/test-generator/src/main.rs index fc2842b0b81f..dff949459a28 100644 --- a/library/compiler-builtins/libm/test-generator/src/main.rs +++ b/library/compiler-builtins/libm/test-generator/src/main.rs @@ -662,7 +662,6 @@ f32_f32! { // coshf, // exp2f, expf, - // fdimf, log10f, log1pf, log2f, @@ -679,6 +678,7 @@ f32_f32! { // With signature `fn(f32, f32) -> f32` f32f32_f32! { // atan2f, + fdimf, hypotf, fmodf, powf, @@ -703,7 +703,7 @@ f64_f64! { ceil, // cos, cosh, - // exp, + exp, // exp2, expm1, floor, @@ -724,7 +724,7 @@ f64_f64! { // With signature `fn(f64, f64) -> f64` f64f64_f64! { // atan2, - // fdim, + fdim, // fmod, hypot, // pow,