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 22734c1e84b4..56ab82eace0f 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 { @@ -347,13 +354,14 @@ impl F32Ext for f32 { pub trait F64Ext: private::Sealed { fn floor(self) -> Self; - #[cfg(todo)] fn ceil(self) -> Self; fn round(self) -> Self; fn trunc(self) -> Self; + fn fdim(self, rhs: Self) -> Self; + #[cfg(todo)] fn fract(self) -> Self; @@ -379,7 +387,6 @@ pub trait F64Ext: private::Sealed { fn sqrt(self) -> Self; - #[cfg(todo)] fn exp(self) -> Self; #[cfg(todo)] @@ -455,7 +462,6 @@ impl F64Ext for f64 { floor(self) } - #[cfg(todo)] #[inline] fn ceil(self) -> Self { ceil(self) @@ -471,6 +477,10 @@ impl F64Ext for f64 { trunc(self) } + #[inline] + fn fdim(self, rhs: Self) -> Self { + fdim(self, rhs) + } #[cfg(todo)] #[inline] fn fract(self) -> Self { @@ -520,7 +530,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/ceil.rs b/library/compiler-builtins/libm/src/math/ceil.rs new file mode 100644 index 000000000000..4db2ca840368 --- /dev/null +++ b/library/compiler-builtins/libm/src/math/ceil.rs @@ -0,0 +1,39 @@ +use core::f64; + +const TOINT: f64 = 1. / f64::EPSILON; + +#[inline] +pub fn ceil(x: f64) -> f64 { + let u: u64 = x.to_bits(); + let e: i64 = (u >> 52 & 0x7ff) as i64; + let y: f64; + + if e >= 0x3ff + 52 || x == 0. { + return x; + } + // y = int(x) - x, where int(x) is an integer neighbor of x + y = if (u >> 63) != 0 { + x - TOINT + TOINT - x + } else { + x + TOINT - TOINT - x + }; + // special case because of non-nearest rounding modes + if e <= 0x3ff - 1 { + force_eval!(y); + return if (u >> 63) != 0 { -0. } else { 1. }; + } + if y < 0. { + x + y + 1. + } else { + x + y + } +} + +#[cfg(test)] +mod tests { + #[test] + fn sanity_check() { + assert_eq!(super::ceil(1.1), 2.0); + assert_eq!(super::ceil(2.9), 3.0); + } +} 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 8d0c2d5c82c1..d84cfeb01b83 100644 --- a/library/compiler-builtins/libm/src/math/fmodf.rs +++ b/library/compiler-builtins/libm/src/math/fmodf.rs @@ -1,3 +1,4 @@ +use core::f32; use core::u32; #[inline] diff --git a/library/compiler-builtins/libm/src/math/mod.rs b/library/compiler-builtins/libm/src/math/mod.rs index d3746e432cdf..a5debc701ccc 100644 --- a/library/compiler-builtins/libm/src/math/mod.rs +++ b/library/compiler-builtins/libm/src/math/mod.rs @@ -6,11 +6,15 @@ macro_rules! force_eval { }; } +mod ceil; mod ceilf; mod cosf; +mod exp; mod expf; mod fabs; mod fabsf; +mod fdim; +mod fdimf; mod floor; mod floorf; mod fmod; @@ -36,11 +40,15 @@ mod trunc; mod truncf; // Use separated imports instead of {}-grouped imports for easier merging. +pub use self::ceil::ceil; pub use self::ceilf::ceilf; pub use self::cosf::cosf; +pub use self::exp::exp; pub use self::expf::expf; 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::fmod::fmod; diff --git a/library/compiler-builtins/libm/test-generator/src/main.rs b/library/compiler-builtins/libm/test-generator/src/main.rs index a5d713afba84..719d596a72f5 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, @@ -700,10 +700,10 @@ f64_f64! { // asin, // atan, // cbrt, - // ceil, + 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,