From a749e176b170c4527813a5d3191b6f203458870d Mon Sep 17 00:00:00 2001 From: Joseph Ryan Date: Fri, 13 Jul 2018 23:15:36 -0500 Subject: [PATCH 1/5] Implement fmod and tweak fmodf --- .../compiler-builtins/libm/src/math/fmod.rs | 80 +++++++++++++++++++ .../compiler-builtins/libm/src/math/fmodf.rs | 4 +- .../compiler-builtins/libm/src/math/mod.rs | 36 +++------ .../libm/test-generator/src/main.rs | 2 +- 4 files changed, 91 insertions(+), 31 deletions(-) create mode 100644 library/compiler-builtins/libm/src/math/fmod.rs diff --git a/library/compiler-builtins/libm/src/math/fmod.rs b/library/compiler-builtins/libm/src/math/fmod.rs new file mode 100644 index 000000000000..23f0c4846c4b --- /dev/null +++ b/library/compiler-builtins/libm/src/math/fmod.rs @@ -0,0 +1,80 @@ +use core::u64; + +#[inline] +pub fn fmod(x: f64, y: f64) -> f64 { + let mut uxi = x.to_bits(); + let mut uyi = y.to_bits(); + let mut ex = (uxi >> 52 & 0x7ff) as i64; + let mut ey = (uyi >> 52 & 0x7ff) as i64; + let sx = uxi >> 63; + let mut i; + + if uyi << 1 == 0 || y.is_nan() || ex == 0x7ff { + return (x * y) / (x * y); + } + if uxi << 1 <= uyi << 1 { + if uxi << 1 == uyi << 1 { + return 0.0 * x; + } + return x; + } + + /* normalize x and y */ + if ex == 0 { + i = uxi << 12; + while i >> 63 == 0 { + ex -= 1; + i <<= 1; + } + uxi <<= -ex + 1; + } else { + uxi &= u64::MAX >> 12; + uxi |= 1 << 52; + } + if ey == 0 { + i = uyi << 12; + while i >> 63 == 0 { + ey -= 1; + i <<= 1; + } + uyi <<= -ey + 1; + } else { + uyi &= u64::MAX >> 12; + uyi |= 1 << 52; + } + + /* x mod y */ + while ex > ey { + i = uxi - uyi; + if i >> 63 == 0 { + if i == 0 { + return 0.0 * x; + } + uxi = i; + } + uxi <<= 1; + ex -= 1; + } + i = uxi - uyi; + if i >> 63 == 0 { + if i == 0 { + return 0.0 * x; + } + uxi = i; + } + while uxi >> 52 == 0 { + uxi <<= 1; + ex -= 1; + } + + /* scale result */ + if ex > 0 { + uxi -= 1 << 52; + uxi |= (ex as u64) << 52; + } else { + uxi >>= -ex + 1; + } + uxi |= (sx as u64) << 63; + + f64::from_bits(uxi) +} diff --git a/library/compiler-builtins/libm/src/math/fmodf.rs b/library/compiler-builtins/libm/src/math/fmodf.rs index 909775249085..8d0c2d5c82c1 100644 --- a/library/compiler-builtins/libm/src/math/fmodf.rs +++ b/library/compiler-builtins/libm/src/math/fmodf.rs @@ -1,7 +1,5 @@ use core::u32; -use super::isnanf; - #[inline] pub fn fmodf(x: f32, y: f32) -> f32 { let mut uxi = x.to_bits(); @@ -11,7 +9,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 d0121048d96b..a81fc174bcb3 100644 --- a/library/compiler-builtins/libm/src/math/mod.rs +++ b/library/compiler-builtins/libm/src/math/mod.rs @@ -7,10 +7,16 @@ macro_rules! force_eval { } mod ceilf; +mod expf; mod fabs; mod fabsf; +mod floor; mod floorf; +mod fmod; mod fmodf; +mod hypot; +mod hypotf; +mod logf; mod powf; mod round; mod roundf; @@ -18,38 +24,14 @@ mod scalbn; mod scalbnf; mod sqrt; mod sqrtf; -mod logf; -mod expf; -mod floor; mod trunc; mod truncf; -mod hypot; -mod hypotf; //mod service; pub use self::{ - ceilf::ceilf, - fabs::fabs, - fabsf::fabsf, - floorf::floorf, - fmodf::fmodf, - powf::powf, - round::round, - roundf::roundf, - scalbn::scalbn, - scalbnf::scalbnf, - sqrt::sqrt, - sqrtf::sqrtf, - logf::logf, - expf::expf, - floor::floor, - trunc::trunc, + ceilf::ceilf, expf::expf, fabs::fabs, fabsf::fabsf, floor::floor, floorf::floorf, fmod::fmod, + fmodf::fmodf, hypot::hypot, hypotf::hypotf, logf::logf, powf::powf, round::round, + roundf::roundf, scalbn::scalbn, scalbnf::scalbnf, sqrt::sqrt, sqrtf::sqrtf, trunc::trunc, truncf::truncf, - hypot::hypot, - hypotf::hypotf, }; - -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 257232ca694a..ee291027d21b 100644 --- a/library/compiler-builtins/libm/test-generator/src/main.rs +++ b/library/compiler-builtins/libm/test-generator/src/main.rs @@ -724,7 +724,7 @@ f64_f64! { f64f64_f64! { // atan2, // fdim, - // fmod, + fmod, hypot, // pow, } From e9575059748b90eba751e31883b4b0b733e6253b Mon Sep 17 00:00:00 2001 From: Joseph Ryan Date: Fri, 13 Jul 2018 23:56:00 -0500 Subject: [PATCH 2/5] Fix log2 --- library/compiler-builtins/libm/src/lib.rs | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/library/compiler-builtins/libm/src/lib.rs b/library/compiler-builtins/libm/src/lib.rs index 0194722af66c..735f8162e5a2 100644 --- a/library/compiler-builtins/libm/src/lib.rs +++ b/library/compiler-builtins/libm/src/lib.rs @@ -229,7 +229,7 @@ impl F32Ext for f32 { #[inline] fn log2(self) -> Self { - log2f(self) + self.log2f(self) } #[cfg(todo)] @@ -556,7 +556,7 @@ impl F64Ext for f64 { #[inline] fn log2(self) -> Self { - log2(self) + self.log2(self) } #[cfg(todo)] From 119bb9bd248aca933b27d8dc50892fcf2cca5402 Mon Sep 17 00:00:00 2001 From: Joseph Ryan Date: Sat, 14 Jul 2018 00:15:24 -0500 Subject: [PATCH 3/5] Revert log2 breakage --- library/compiler-builtins/libm/src/lib.rs | 4 ++-- library/compiler-builtins/libm/src/math/mod.rs | 6 +++--- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/library/compiler-builtins/libm/src/lib.rs b/library/compiler-builtins/libm/src/lib.rs index 735f8162e5a2..0194722af66c 100644 --- a/library/compiler-builtins/libm/src/lib.rs +++ b/library/compiler-builtins/libm/src/lib.rs @@ -229,7 +229,7 @@ impl F32Ext for f32 { #[inline] fn log2(self) -> Self { - self.log2f(self) + log2f(self) } #[cfg(todo)] @@ -556,7 +556,7 @@ impl F64Ext for f64 { #[inline] fn log2(self) -> Self { - self.log2(self) + log2(self) } #[cfg(todo)] diff --git a/library/compiler-builtins/libm/src/math/mod.rs b/library/compiler-builtins/libm/src/math/mod.rs index e75429df59c7..56488c43806b 100644 --- a/library/compiler-builtins/libm/src/math/mod.rs +++ b/library/compiler-builtins/libm/src/math/mod.rs @@ -33,7 +33,7 @@ mod truncf; pub use self::{ ceilf::ceilf, expf::expf, fabs::fabs, fabsf::fabsf, floor::floor, floorf::floorf, fmod::fmod, - fmodf::fmodf, hypot::hypot, hypotf::hypotf, logf::logf, powf::powf, round::round, - roundf::roundf, scalbn::scalbn, scalbnf::scalbnf, sqrt::sqrt, sqrtf::sqrtf, trunc::trunc, - truncf::truncf, + 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, }; From 88faf19ca55c94eebff3a113e40158f4413dbb16 Mon Sep 17 00:00:00 2001 From: Joseph Ryan Date: Sat, 14 Jul 2018 00:44:36 -0500 Subject: [PATCH 4/5] Run rustfmt --- library/compiler-builtins/libm/src/math/mod.rs | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/library/compiler-builtins/libm/src/math/mod.rs b/library/compiler-builtins/libm/src/math/mod.rs index 855929081470..afa8cd8ef233 100644 --- a/library/compiler-builtins/libm/src/math/mod.rs +++ b/library/compiler-builtins/libm/src/math/mod.rs @@ -35,7 +35,7 @@ mod truncf; pub use self::{ ceilf::ceilf, expf::expf, fabs::fabs, fabsf::fabsf, floor::floor, floorf::floorf, fmod::fmod, - fmodf::fmodf, hypot::hypot, hypotf::hypotf, log10::log10, log10f::log10f, log2::log2, - log2f::log2f, logf::logf, powf::powf, round::round, roundf::roundf, scalbn::scalbn, + fmodf::fmodf, 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, }; From e22803605379ea10989050707a6a30d02ab8a141 Mon Sep 17 00:00:00 2001 From: Erik Date: Sat, 14 Jul 2018 14:07:14 -0400 Subject: [PATCH 5/5] implement cbrt and cbrtf --- library/compiler-builtins/libm/src/lib.rs | 4 - .../compiler-builtins/libm/src/math/cbrt.rs | 110 ++++++++++++++++++ .../compiler-builtins/libm/src/math/cbrtf.rs | 72 ++++++++++++ .../compiler-builtins/libm/src/math/mod.rs | 4 + .../libm/test-generator/src/main.rs | 4 +- 5 files changed, 188 insertions(+), 6 deletions(-) create mode 100644 library/compiler-builtins/libm/src/math/cbrt.rs create mode 100644 library/compiler-builtins/libm/src/math/cbrtf.rs diff --git a/library/compiler-builtins/libm/src/lib.rs b/library/compiler-builtins/libm/src/lib.rs index f112aaaca6b1..b175f6864117 100644 --- a/library/compiler-builtins/libm/src/lib.rs +++ b/library/compiler-builtins/libm/src/lib.rs @@ -79,7 +79,6 @@ pub trait F32Ext: private::Sealed { fn log10(self) -> Self; - #[cfg(todo)] fn cbrt(self) -> Self; fn hypot(self, other: Self) -> Self; @@ -234,7 +233,6 @@ impl F32Ext for f32 { log10f(self) } - #[cfg(todo)] #[inline] fn cbrt(self) -> Self { cbrtf(self) @@ -391,7 +389,6 @@ pub trait F64Ext: private::Sealed { fn log10(self) -> Self; - #[cfg(todo)] fn cbrt(self) -> Self; fn hypot(self, other: Self) -> Self; @@ -548,7 +545,6 @@ impl F64Ext for f64 { log10(self) } - #[cfg(todo)] #[inline] fn cbrt(self) -> Self { cbrt(self) diff --git a/library/compiler-builtins/libm/src/math/cbrt.rs b/library/compiler-builtins/libm/src/math/cbrt.rs new file mode 100644 index 000000000000..8c37f0b266b8 --- /dev/null +++ b/library/compiler-builtins/libm/src/math/cbrt.rs @@ -0,0 +1,110 @@ +/* origin: FreeBSD /usr/src/lib/msun/src/s_cbrt.c */ +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + * + * Optimized by Bruce D. Evans. + */ +/* cbrt(x) + * Return cube root of x + */ + +use core::f64; + +const B1: u32 = 715094163; /* B1 = (1023-1023/3-0.03306235651)*2**20 */ +const B2: u32 = 696219795; /* B2 = (1023-1023/3-54/3-0.03306235651)*2**20 */ + +/* |1/cbrt(x) - p(x)| < 2**-23.5 (~[-7.93e-8, 7.929e-8]). */ +const P0: f64 = 1.87595182427177009643; /* 0x3ffe03e6, 0x0f61e692 */ +const P1: f64 = -1.88497979543377169875; /* 0xbffe28e0, 0x92f02420 */ +const P2: f64 = 1.621429720105354466140; /* 0x3ff9f160, 0x4a49d6c2 */ +const P3: f64 = -0.758397934778766047437; /* 0xbfe844cb, 0xbee751d9 */ +const P4: f64 = 0.145996192886612446982; /* 0x3fc2b000, 0xd4e4edd7 */ + +#[inline] +pub fn cbrt(x: f64) -> f64 { + let x1p54 = f64::from_bits(0x4350000000000000); // 0x1p54 === 2 ^ 54 + + let mut ui: u64 = x.to_bits(); + let mut r: f64; + let s: f64; + let mut t: f64; + let w: f64; + let mut hx: u32 = (ui >> 32) as u32 & 0x7fffffff; + + if hx >= 0x7ff00000 { + /* cbrt(NaN,INF) is itself */ + return x + x; + } + + /* + * Rough cbrt to 5 bits: + * cbrt(2**e*(1+m) ~= 2**(e/3)*(1+(e%3+m)/3) + * where e is integral and >= 0, m is real and in [0, 1), and "/" and + * "%" are integer division and modulus with rounding towards minus + * infinity. The RHS is always >= the LHS and has a maximum relative + * error of about 1 in 16. Adding a bias of -0.03306235651 to the + * (e%3+m)/3 term reduces the error to about 1 in 32. With the IEEE + * floating point representation, for finite positive normal values, + * ordinary integer divison of the value in bits magically gives + * almost exactly the RHS of the above provided we first subtract the + * exponent bias (1023 for doubles) and later add it back. We do the + * subtraction virtually to keep e >= 0 so that ordinary integer + * division rounds towards minus infinity; this is also efficient. + */ + if hx < 0x00100000 { + /* zero or subnormal? */ + ui = (x * x1p54).to_bits(); + hx = (ui >> 32) as u32 & 0x7fffffff; + if hx == 0 { + return x; /* cbrt(0) is itself */ + } + hx = hx / 3 + B2; + } else { + hx = hx / 3 + B1; + } + ui &= 1 << 63; + ui |= (hx as u64) << 32; + t = f64::from_bits(ui); + + /* + * New cbrt to 23 bits: + * cbrt(x) = t*cbrt(x/t**3) ~= t*P(t**3/x) + * where P(r) is a polynomial of degree 4 that approximates 1/cbrt(r) + * to within 2**-23.5 when |r - 1| < 1/10. The rough approximation + * has produced t such than |t/cbrt(x) - 1| ~< 1/32, and cubing this + * gives us bounds for r = t**3/x. + * + * Try to optimize for parallel evaluation as in __tanf.c. + */ + r = (t * t) * (t / x); + t = t * ((P0 + r * (P1 + r * P2)) + ((r * r) * r) * (P3 + r * P4)); + + /* + * Round t away from zero to 23 bits (sloppily except for ensuring that + * the result is larger in magnitude than cbrt(x) but not much more than + * 2 23-bit ulps larger). With rounding towards zero, the error bound + * would be ~5/6 instead of ~4/6. With a maximum error of 2 23-bit ulps + * in the rounded t, the infinite-precision error in the Newton + * approximation barely affects third digit in the final error + * 0.667; the error in the rounded t can be up to about 3 23-bit ulps + * before the final error is larger than 0.667 ulps. + */ + ui = t.to_bits(); + ui = (ui + 0x80000000) & 0xffffffffc0000000; + t = f64::from_bits(ui); + + /* one step Newton iteration to 53 bits with error < 0.667 ulps */ + s = t * t; /* t*t is exact */ + r = x / s; /* error <= 0.5 ulps; |r| < |t| */ + w = t + t; /* t+t is exact */ + r = (r - t) / (w + r); /* r-t is exact; w+r ~= 3*t */ + t = t + t * r; /* error <= 0.5 + 0.5/3 + epsilon */ + t +} diff --git a/library/compiler-builtins/libm/src/math/cbrtf.rs b/library/compiler-builtins/libm/src/math/cbrtf.rs new file mode 100644 index 000000000000..878372eefbf8 --- /dev/null +++ b/library/compiler-builtins/libm/src/math/cbrtf.rs @@ -0,0 +1,72 @@ +/* origin: FreeBSD /usr/src/lib/msun/src/s_cbrtf.c */ +/* + * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. + * Debugged and optimized by Bruce D. Evans. + */ +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ +/* cbrtf(x) + * Return cube root of x + */ + +use core::f32; + +const B1: u32 = 709958130; /* B1 = (127-127.0/3-0.03306235651)*2**23 */ +const B2: u32 = 642849266; /* B2 = (127-127.0/3-24/3-0.03306235651)*2**23 */ + +#[inline] +pub fn cbrtf(x: f32) -> f32 { + let x1p24 = f32::from_bits(0x4b800000); // 0x1p24f === 2 ^ 24 + + let mut r: f64; + let mut t: f64; + let mut ui: u32 = x.to_bits(); + let mut hx: u32 = ui & 0x7fffffff; + + if hx >= 0x7f800000 { + /* cbrt(NaN,INF) is itself */ + return x + x; + } + + /* rough cbrt to 5 bits */ + if hx < 0x00800000 { + /* zero or subnormal? */ + if hx == 0 { + return x; /* cbrt(+-0) is itself */ + } + ui = (x * x1p24).to_bits(); + hx = ui & 0x7fffffff; + hx = hx / 3 + B2; + } else { + hx = hx / 3 + B1; + } + ui &= 0x80000000; + ui |= hx; + + /* + * First step Newton iteration (solving t*t-x/t == 0) to 16 bits. In + * double precision so that its terms can be arranged for efficiency + * without causing overflow or underflow. + */ + t = f32::from_bits(ui) as f64; + r = t * t * t; + t = t * (x as f64 + x as f64 + r) / (x as f64 + r + r); + + /* + * Second step Newton iteration to 47 bits. In double precision for + * efficiency and accuracy. + */ + r = t * t * t; + t = t * (x as f64 + x as f64 + r) / (x as f64 + r + r); + + /* rounding to 24 bits is perfect in round-to-nearest mode */ + t as f32 +} diff --git a/library/compiler-builtins/libm/src/math/mod.rs b/library/compiler-builtins/libm/src/math/mod.rs index 2a84b463d9fd..0f112a0cbbee 100644 --- a/library/compiler-builtins/libm/src/math/mod.rs +++ b/library/compiler-builtins/libm/src/math/mod.rs @@ -6,6 +6,8 @@ macro_rules! force_eval { }; } +mod cbrt; +mod cbrtf; mod ceil; mod ceilf; mod cosf; @@ -37,6 +39,8 @@ mod trunc; mod truncf; // Use separated imports instead of {}-grouped imports for easier merging. +pub use self::cbrt::cbrt; +pub use self::cbrtf::cbrtf; pub use self::ceil::ceil; pub use self::ceilf::ceilf; pub use self::cosf::cosf; diff --git a/library/compiler-builtins/libm/test-generator/src/main.rs b/library/compiler-builtins/libm/test-generator/src/main.rs index a54d8b2710f5..0521538deff4 100644 --- a/library/compiler-builtins/libm/test-generator/src/main.rs +++ b/library/compiler-builtins/libm/test-generator/src/main.rs @@ -656,7 +656,7 @@ f32_f32! { truncf, // asinf, // atanf, - // cbrtf, + cbrtf, cosf, ceilf, // coshf, @@ -699,7 +699,7 @@ f64_f64! { // acos, // asin, // atan, - // cbrt, + cbrt, ceil, // cos, // cosh,