From 65368071063a2b4a7d117301d67b6f3b66c70d1a Mon Sep 17 00:00:00 2001 From: Erik Date: Sat, 14 Jul 2018 00:24:22 -0400 Subject: [PATCH] implement log10 and log10f --- library/compiler-builtins/libm/src/lib.rs | 4 - .../compiler-builtins/libm/src/math/log10.rs | 98 +++++++++++++++++++ .../compiler-builtins/libm/src/math/log10f.rs | 76 ++++++++++++++ .../compiler-builtins/libm/src/math/mod.rs | 8 +- .../libm/test-generator/src/main.rs | 4 +- 5 files changed, 181 insertions(+), 9 deletions(-) create mode 100644 library/compiler-builtins/libm/src/math/log10.rs create mode 100644 library/compiler-builtins/libm/src/math/log10f.rs diff --git a/library/compiler-builtins/libm/src/lib.rs b/library/compiler-builtins/libm/src/lib.rs index 0194722af66c..3ed89f729507 100644 --- a/library/compiler-builtins/libm/src/lib.rs +++ b/library/compiler-builtins/libm/src/lib.rs @@ -77,7 +77,6 @@ pub trait F32Ext: private::Sealed { fn log2(self) -> Self; - #[cfg(todo)] fn log10(self) -> Self; #[cfg(todo)] @@ -232,7 +231,6 @@ impl F32Ext for f32 { log2f(self) } - #[cfg(todo)] #[inline] fn log10(self) -> Self { log10f(self) @@ -399,7 +397,6 @@ pub trait F64Ext: private::Sealed { fn log2(self) -> Self; - #[cfg(todo)] fn log10(self) -> Self; #[cfg(todo)] @@ -559,7 +556,6 @@ impl F64Ext for f64 { log2(self) } - #[cfg(todo)] #[inline] fn log10(self) -> Self { log10(self) diff --git a/library/compiler-builtins/libm/src/math/log10.rs b/library/compiler-builtins/libm/src/math/log10.rs new file mode 100644 index 000000000000..137d8966f521 --- /dev/null +++ b/library/compiler-builtins/libm/src/math/log10.rs @@ -0,0 +1,98 @@ +use core::f64; + +const IVLN10HI: f64 = 4.34294481878168880939e-01; /* 0x3fdbcb7b, 0x15200000 */ +const IVLN10LO: f64 = 2.50829467116452752298e-11; /* 0x3dbb9438, 0xca9aadd5 */ +const LOG10_2HI: f64 = 3.01029995663611771306e-01; /* 0x3FD34413, 0x509F6000 */ +const LOG10_2LO: f64 = 3.69423907715893078616e-13; /* 0x3D59FEF3, 0x11F12B36 */ +const LG1: f64 = 6.666666666666735130e-01; /* 3FE55555 55555593 */ +const LG2: f64 = 3.999999999940941908e-01; /* 3FD99999 9997FA04 */ +const LG3: f64 = 2.857142874366239149e-01; /* 3FD24924 94229359 */ +const LG4: f64 = 2.222219843214978396e-01; /* 3FCC71C5 1D8E78AF */ +const LG5: f64 = 1.818357216161805012e-01; /* 3FC74664 96CB03DE */ +const LG6: f64 = 1.531383769920937332e-01; /* 3FC39A09 D078C69F */ +const LG7: f64 = 1.479819860511658591e-01; /* 3FC2F112 DF3E5244 */ + +#[inline] +pub fn log10(mut x: f64) -> f64 { + let x1p54 = f64::from_bits(0x4350000000000000); // 0x1p54 === 2 ^ 54 + + let mut ui: u64 = x.to_bits(); + let hfsq: f64; + let f: f64; + let s: f64; + let z: f64; + let r: f64; + let mut w: f64; + let t1: f64; + let t2: f64; + let dk: f64; + let y: f64; + let mut hi: f64; + let lo: f64; + let mut val_hi: f64; + let mut val_lo: f64; + let mut hx: u32; + let mut k: i32; + + hx = (ui >> 32) as u32; + k = 0; + if hx < 0x00100000 || (hx >> 31) > 0 { + if ui << 1 == 0 { + return -1. / (x * x); /* log(+-0)=-inf */ + } + if (hx >> 31) > 0 { + return (x - x) / 0.0; /* log(-#) = NaN */ + } + /* subnormal number, scale x up */ + k -= 54; + x *= x1p54; + ui = x.to_bits(); + hx = (ui >> 32) as u32; + } else if hx >= 0x7ff00000 { + return x; + } else if hx == 0x3ff00000 && ui << 32 == 0 { + return 0.; + } + + /* reduce x into [sqrt(2)/2, sqrt(2)] */ + hx += 0x3ff00000 - 0x3fe6a09e; + k += (hx >> 20) as i32 - 0x3ff; + hx = (hx & 0x000fffff) + 0x3fe6a09e; + ui = (hx as u64) << 32 | (ui & 0xffffffff); + x = f64::from_bits(ui); + + f = x - 1.0; + hfsq = 0.5 * f * f; + s = f / (2.0 + f); + z = s * s; + w = z * z; + t1 = w * (LG2 + w * (LG4 + w * LG6)); + t2 = z * (LG1 + w * (LG3 + w * (LG5 + w * LG7))); + r = t2 + t1; + + /* See log2.c for details. */ + /* hi+lo = f - hfsq + s*(hfsq+R) ~ log(1+f) */ + hi = f - hfsq; + ui = hi.to_bits(); + ui &= (-1i64 as u64) << 32; + hi = f64::from_bits(ui); + lo = f - hi - hfsq + s * (hfsq + r); + + /* val_hi+val_lo ~ log10(1+f) + k*log10(2) */ + val_hi = hi * IVLN10HI; + dk = k as f64; + y = dk * LOG10_2HI; + val_lo = dk * LOG10_2LO + (lo + hi) * IVLN10LO + lo * IVLN10HI; + + /* + * Extra precision in for adding y is not strictly needed + * since there is no very large cancellation near x = sqrt(2) or + * x = 1/sqrt(2), but we do it anyway since it costs little on CPUs + * with some parallelism and it reduces the error for many args. + */ + w = y + val_hi; + val_lo += (y - w) + val_hi; + val_hi = w; + + return val_lo + val_hi; +} diff --git a/library/compiler-builtins/libm/src/math/log10f.rs b/library/compiler-builtins/libm/src/math/log10f.rs new file mode 100644 index 000000000000..58db093441f3 --- /dev/null +++ b/library/compiler-builtins/libm/src/math/log10f.rs @@ -0,0 +1,76 @@ +use core::f32; + +const IVLN10HI: f32 = 4.3432617188e-01; /* 0x3ede6000 */ +const IVLN10LO: f32 = -3.1689971365e-05; /* 0xb804ead9 */ +const LOG10_2HI: f32 = 3.0102920532e-01; /* 0x3e9a2080 */ +const LOG10_2LO: f32 = 7.9034151668e-07; /* 0x355427db */ +/* |(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 log10f(mut x: f32) -> f32 { + let x1p25f = f32::from_bits(0x4c000000); // 0x1p25f === 2 ^ 25 + + let mut ui: u32 = x.to_bits(); + let hfsq: f32; + let f: f32; + let s: f32; + let z: f32; + let r: f32; + let w: f32; + let t1: f32; + let t2: f32; + let dk: f32; + let mut hi: f32; + let lo: f32; + let mut ix: u32; + let mut k: i32; + + ix = ui; + k = 0; + 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.0; /* log(-#) = NaN */ + } + /* subnormal number, scale up x */ + k -= 25; + x *= x1p25f; + ui = x.to_bits(); + ix = ui; + } 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; + ui = ix; + x = f32::from_bits(ui); + + f = x - 1.0; + s = f / (2.0 + f); + z = s * s; + w = z * z; + t1 = w * (LG2 + w * LG4); + t2 = z * (LG1 + w * LG3); + r = t2 + t1; + hfsq = 0.5 * f * f; + + hi = f - hfsq; + ui = hi.to_bits(); + ui &= 0xfffff000; + hi = f32::from_bits(ui); + lo = f - hi - hfsq + s * (hfsq + r); + dk = k as f32; + return dk * LOG10_2LO + (lo + hi) * IVLN10LO + lo * IVLN10HI + hi * IVLN10HI + dk * LOG10_2HI; +} diff --git a/library/compiler-builtins/libm/src/math/mod.rs b/library/compiler-builtins/libm/src/math/mod.rs index a82c1ba7e592..828eef859085 100644 --- a/library/compiler-builtins/libm/src/math/mod.rs +++ b/library/compiler-builtins/libm/src/math/mod.rs @@ -15,6 +15,8 @@ mod floorf; mod fmodf; mod hypot; mod hypotf; +mod log10; +mod log10f; mod log2; mod log2f; mod logf; @@ -32,9 +34,9 @@ mod truncf; pub use self::{ ceilf::ceilf, expf::expf, fabs::fabs, fabsf::fabsf, floor::floor, floorf::floorf, fmodf::fmodf, - hypot::hypot, hypotf::hypotf, log2::log2, log2f::log2f, logf::logf, powf::powf, round::round, - roundf::roundf, scalbn::scalbn, scalbnf::scalbnf, sqrt::sqrt, sqrtf::sqrtf, trunc::trunc, - truncf::truncf, + hypot::hypot, hypotf::hypotf, log10::log10, log10f::log10f, 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 { diff --git a/library/compiler-builtins/libm/test-generator/src/main.rs b/library/compiler-builtins/libm/test-generator/src/main.rs index 7ec2292c3498..799a354a90af 100644 --- a/library/compiler-builtins/libm/test-generator/src/main.rs +++ b/library/compiler-builtins/libm/test-generator/src/main.rs @@ -663,7 +663,7 @@ f32_f32! { // exp2f, expf, // fdimf, - // log10f, + log10f, log2f, logf, roundf, @@ -707,7 +707,7 @@ f64_f64! { // expm1, floor, // log, - // log10, + log10, // log1p, log2, round,