diff --git a/library/compiler-builtins/libm/src/lib.rs b/library/compiler-builtins/libm/src/lib.rs index 0d13590c7ecc..bb00149dfd24 100644 --- a/library/compiler-builtins/libm/src/lib.rs +++ b/library/compiler-builtins/libm/src/lib.rs @@ -71,16 +71,13 @@ pub trait F32Ext { fn sqrt(self) -> Self; - #[cfg(todo)] fn exp(self) -> Self; #[cfg(todo)] fn exp2(self) -> Self; - #[cfg(todo)] fn ln(self) -> Self; - #[cfg(todo)] fn log(self, base: Self) -> Self; #[cfg(todo)] @@ -220,7 +217,6 @@ impl F32Ext for f32 { sqrtf(self) } - #[cfg(todo)] #[inline] fn exp(self) -> Self { expf(self) @@ -232,13 +228,11 @@ impl F32Ext for f32 { exp2f(self) } - #[cfg(todo)] #[inline] fn ln(self) -> Self { logf(self) } - #[cfg(todo)] #[inline] fn log(self, base: Self) -> Self { self.ln() / base.ln() diff --git a/library/compiler-builtins/libm/src/math/expf.rs b/library/compiler-builtins/libm/src/math/expf.rs new file mode 100644 index 000000000000..c26ff8e64124 --- /dev/null +++ b/library/compiler-builtins/libm/src/math/expf.rs @@ -0,0 +1,76 @@ +use super::scalbnf; + +const HALF : [f32; 2] = [0.5,-0.5]; +const LN2_HI : f32 = 6.9314575195e-01; /* 0x3f317200 */ +const LN2_LO : f32 = 1.4286067653e-06; /* 0x35bfbe8e */ +const INV_LN2 : f32 = 1.4426950216e+00; /* 0x3fb8aa3b */ +/* + * Domain [-0.34568, 0.34568], range ~[-4.278e-9, 4.447e-9]: + * |x*(exp(x)+1)/(exp(x)-1) - p(x)| < 2**-27.74 + */ +const P1 : f32 = 1.6666625440e-1; /* 0xaaaa8f.0p-26 */ +const P2 : f32 = -2.7667332906e-3; /* -0xb55215.0p-32 */ + +#[inline] +pub fn expf(mut x: f32) -> f32 { + let x1p127 = f32::from_bits(0x7f000000); // 0x1p127f === 2 ^ 127 + let x1p_126 = f32::from_bits(0x800000); // 0x1p-126f === 2 ^ -126 /*original 0x1p-149f ??????????? */ + + let mut hx = x.to_bits() as i32; + let sign = (hx >> 31) as i32; /* sign bit of x */ + let signb : bool = sign != 0; + hx &= 0x7fffffff; /* high word of |x| */ + + /* special cases */ + if hx >= 0x42aeac50 { /* if |x| >= -87.33655f or NaN */ + if hx > 0x7f800000 {/* NaN */ + return x; + } + if (hx >= 0x42b17218) && (!signb) { /* x >= 88.722839f */ + /* overflow */ + x *= x1p127; + return x; + } + if signb { + /* underflow */ + force_eval!(-x1p_126/x); + if hx >= 0x42cff1b5 { /* x <= -103.972084f */ + return 0. + } + } + } + + /* argument reduction */ + let k : i32; + let hi : f32; + let lo : f32; + if hx > 0x3eb17218 { /* if |x| > 0.5 ln2 */ + if hx > 0x3f851592 { /* if |x| > 1.5 ln2 */ + k = (INV_LN2*x + HALF[sign as usize]) as i32; + } else { + k = 1 - sign - sign; + } + let kf = k as f32; + hi = x - kf*LN2_HI; /* k*ln2hi is exact here */ + lo = kf*LN2_LO; + x = hi - lo; + } else if hx > 0x39000000 { /* |x| > 2**-14 */ + k = 0; + hi = x; + lo = 0.; + } else { + /* raise inexact */ + force_eval!(x1p127 + x); + return 1. + x; + } + + /* x is now in primary range */ + let xx = x*x; + let c = x - xx*(P1+xx*P2); + let y = 1. + (x*c/(2.-c) - lo + hi); + if k == 0 { + y + } else { + scalbnf(y, k) + } +} diff --git a/library/compiler-builtins/libm/src/math/logf.rs b/library/compiler-builtins/libm/src/math/logf.rs new file mode 100644 index 000000000000..f8bbe733aa85 --- /dev/null +++ b/library/compiler-builtins/libm/src/math/logf.rs @@ -0,0 +1,49 @@ +const LN2_HI : f32 = 6.9313812256e-01; /* 0x3f317180 */ +const LN2_LO : f32 = 9.0580006145e-06; /* 0x3717f7d1 */ +/* |(log(1+s)-log(1-s))/s - Lg(s)| < 2**-34.24 (~[-4.95e-11, 4.97e-11]). */ +const LG1 : f32 = 0.66666662693; /* 0xaaaaaa.0p-24*/ +const LG2 : f32 = 0.40000972152; /* 0xccce13.0p-25 */ +const LG3 : f32 = 0.28498786688; /* 0x91e9ee.0p-25 */ +const LG4 : f32 = 0.24279078841; /* 0xf89e26.0p-26 */ + +#[inline] +pub fn logf(mut x: f32) -> f32 { + let x1p25 = f32::from_bits(0x4c000000); // 0x1p24f === 2 ^ 25 + + let mut ix = x.to_bits(); + let mut k = 0i32; + + if (ix < 0x00800000) || ((ix>>31) != 0) { /* x < 2**-126 */ + if ix<<1 == 0 { + return -1./(x*x); /* log(+-0)=-inf */ + } + if (ix>>31) != 0 { + return (x-x)/0.; /* log(-#) = NaN */ + } + /* subnormal number, scale up x */ + k -= 25; + x *= x1p25; + ix = x.to_bits(); + } else if ix >= 0x7f800000 { + return x; + } else if ix == 0x3f800000 { + return 0.; + } + + /* reduce x into [sqrt(2)/2, sqrt(2)] */ + ix += 0x3f800000 - 0x3f3504f3; + k += ((ix>>23) as i32) - 0x7f; + ix = (ix & 0x007fffff) + 0x3f3504f3; + x = f32::from_bits(ix); + + let f = x - 1.; + let s = f/(2. + f); + let z = s*s; + let w = z*z; + let t1 = w*(LG2+w*LG4); + let t2 = z*(LG1+w*LG3); + let r = t2 + t1; + let hfsq = 0.5*f*f; + let dk = k as f32; + s*(hfsq+r) + dk*LN2_LO - hfsq + f + dk*LN2_HI +} diff --git a/library/compiler-builtins/libm/src/math/mod.rs b/library/compiler-builtins/libm/src/math/mod.rs index 41359bf8c059..5edf6fd740df 100644 --- a/library/compiler-builtins/libm/src/math/mod.rs +++ b/library/compiler-builtins/libm/src/math/mod.rs @@ -11,6 +11,8 @@ mod powf; mod round; mod scalbnf; mod sqrtf; +mod logf; +mod expf; pub use self::fabs::fabs; pub use self::fabsf::fabsf; @@ -19,6 +21,8 @@ pub use self::powf::powf; pub use self::round::round; pub use self::scalbnf::scalbnf; pub use self::sqrtf::sqrtf; +pub use self::logf::logf; +pub use self::expf::expf; fn isnanf(x: f32) -> bool { x.to_bits() & 0x7fffffff > 0x7f800000