From e22803605379ea10989050707a6a30d02ab8a141 Mon Sep 17 00:00:00 2001 From: Erik Date: Sat, 14 Jul 2018 14:07:14 -0400 Subject: [PATCH] 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,