From 892cd4f5b83761f2abf96e27c109ad3a6d87bceb Mon Sep 17 00:00:00 2001 From: Erik Date: Fri, 13 Jul 2018 18:31:01 -0400 Subject: [PATCH] implement sqrt and hypot --- library/compiler-builtins/libm/src/lib.rs | 4 - .../compiler-builtins/libm/src/math/hypot.rs | 74 ++++++++++ .../compiler-builtins/libm/src/math/mod.rs | 4 + .../compiler-builtins/libm/src/math/sqrt.rs | 129 ++++++++++++++++++ .../libm/test-generator/src/main.rs | 4 +- 5 files changed, 209 insertions(+), 6 deletions(-) create mode 100644 library/compiler-builtins/libm/src/math/hypot.rs create mode 100644 library/compiler-builtins/libm/src/math/sqrt.rs diff --git a/library/compiler-builtins/libm/src/lib.rs b/library/compiler-builtins/libm/src/lib.rs index 0fb555a7b499..3e25c16db29d 100644 --- a/library/compiler-builtins/libm/src/lib.rs +++ b/library/compiler-builtins/libm/src/lib.rs @@ -391,7 +391,6 @@ pub trait F64Ext: private::Sealed { #[cfg(todo)] fn powf(self, n: Self) -> Self; - #[cfg(todo)] fn sqrt(self) -> Self; #[cfg(todo)] @@ -415,7 +414,6 @@ pub trait F64Ext: private::Sealed { #[cfg(todo)] fn cbrt(self) -> Self; - #[cfg(todo)] fn hypot(self, other: Self) -> Self; #[cfg(todo)] @@ -536,7 +534,6 @@ impl F64Ext for f64 { pow(self, n) } - #[cfg(todo)] #[inline] fn sqrt(self) -> Self { sqrt(self) @@ -584,7 +581,6 @@ impl F64Ext for f64 { cbrt(self) } - #[cfg(todo)] #[inline] fn hypot(self, other: Self) -> Self { hypot(self, other) diff --git a/library/compiler-builtins/libm/src/math/hypot.rs b/library/compiler-builtins/libm/src/math/hypot.rs new file mode 100644 index 000000000000..dcc17d914fb9 --- /dev/null +++ b/library/compiler-builtins/libm/src/math/hypot.rs @@ -0,0 +1,74 @@ +use core::f64; + +use super::sqrt; + +const SPLIT: f64 = 134217728. + 1.; // 0x1p27 + 1 === (2 ^ 27) + 1 + +fn sq(x: f64) -> (f64, f64) { + let xh: f64; + let xl: f64; + let xc: f64; + + xc = x * SPLIT; + xh = x - xc + xc; + xl = x - xh; + let hi = x*x; + let lo = xh*xh - hi + 2.*xh*xl + xl*xl; + (hi, lo) +} + +#[inline] +pub fn hypot(mut x: f64, mut y: f64) -> f64 { + let x1p700 = f64::from_bits(0x6bb0000000000000); // 0x1p700 === 2 ^ 700 + let x1p_700 = f64::from_bits(0x1430000000000000); // 0x1p-700 === 2 ^ -700 + + let mut uxi = x.to_bits(); + let mut uyi = y.to_bits(); + let uti; + let ex: i64; + let ey: i64; + let mut z: f64; + + /* arrange |x| >= |y| */ + uxi &= -1i64 as u64 >> 1; + uyi &= -1i64 as u64 >> 1; + if uxi < uyi { + uti = uxi; + uxi = uyi; + uyi = uti; + } + + /* special cases */ + ex = (uxi>>52) as i64; + ey = (uyi>>52) as i64; + x = f64::from_bits(uxi); + y = f64::from_bits(uyi); + /* note: hypot(inf,nan) == inf */ + if ey == 0x7ff { + return y; + } + if ex == 0x7ff || uyi == 0 { + return x; + } + /* note: hypot(x,y) ~= x + y*y/x/2 with inexact for small y/x */ + /* 64 difference is enough for ld80 double_t */ + if ex - ey > 64 { + return x + y; + } + + /* precise sqrt argument in nearest rounding mode without overflow */ + /* xh*xh must not overflow and xl*xl must not underflow in sq */ + z = 1.; + if ex > 0x3ff+510 { + z = x1p700; + x *= x1p_700; + y *= x1p_700; + } else if ey < 0x3ff-450 { + z = x1p_700; + x *= x1p700; + y *= x1p700; + } + let (hx, lx) = sq(x); + let (hy, ly) = sq(y); + return z*sqrt(ly+lx+hy+hx); +} diff --git a/library/compiler-builtins/libm/src/math/mod.rs b/library/compiler-builtins/libm/src/math/mod.rs index a34ebd767e54..444652e81eff 100644 --- a/library/compiler-builtins/libm/src/math/mod.rs +++ b/library/compiler-builtins/libm/src/math/mod.rs @@ -11,12 +11,14 @@ mod powf; mod round; mod scalbn; mod scalbnf; +mod sqrt; mod sqrtf; mod logf; mod expf; mod floor; mod trunc; mod truncf; +mod hypot; mod hypotf; //mod service; @@ -29,12 +31,14 @@ pub use self::{ round::round, scalbn::scalbn, scalbnf::scalbnf, + sqrt::sqrt, sqrtf::sqrtf, logf::logf, expf::expf, floor::floor, trunc::trunc, truncf::truncf, + hypot::hypot, hypotf::hypotf, }; diff --git a/library/compiler-builtins/libm/src/math/sqrt.rs b/library/compiler-builtins/libm/src/math/sqrt.rs new file mode 100644 index 000000000000..49fc58fff03c --- /dev/null +++ b/library/compiler-builtins/libm/src/math/sqrt.rs @@ -0,0 +1,129 @@ +use core::f64; + +const TINY: f64 = 1.0e-300; + +#[inline] +pub fn sqrt(x: f64) -> f64 { + let mut z: f64; + let sign: u32 = 0x80000000; + let mut ix0: i32; + let mut s0: i32; + let mut q: i32; + let mut m: i32; + let mut t: i32; + let mut i: i32; + let mut r: u32; + let mut t1: u32; + let mut s1: u32; + let mut ix1: u32; + let mut q1: u32; + + ix0 = (x.to_bits() >> 32) as i32; + ix1 = x.to_bits() as u32; + + /* take care of Inf and NaN */ + if (ix0&0x7ff00000) == 0x7ff00000 { + return x*x + x; /* sqrt(NaN)=NaN, sqrt(+inf)=+inf, sqrt(-inf)=sNaN */ + } + /* take care of zero */ + if ix0 <= 0 { + if ((ix0&!(sign as i32))|ix1 as i32) == 0 { + return x; /* sqrt(+-0) = +-0 */ + } + if ix0 < 0 { + return (x - x) / (x - x); /* sqrt(-ve) = sNaN */ + } + } + /* normalize x */ + m = ix0>>20; + if m == 0 { /* subnormal x */ + while ix0 == 0 { + m -= 21; + ix0 |= (ix1>>11) as i32; + ix1 <<= 21; + } + i=0; + while (ix0&0x00100000) == 0 { + i += 1; + ix0 <<= 1; + } + m -= i - 1; + ix0 |= (ix1>>(32-i)) as i32; + ix1 <<= i; + } + m -= 1023; /* unbias exponent */ + ix0 = (ix0&0x000fffff)|0x00100000; + if (m & 1) == 1 { /* odd m, double x to make it even */ + ix0 += ix0 + ((ix1&sign)>>31) as i32; + ix1 += ix1; + } + m >>= 1; /* m = [m/2] */ + + /* generate sqrt(x) bit by bit */ + ix0 += ix0 + ((ix1&sign)>>31) as i32; + ix1 += ix1; + q = 0; /* [q,q1] = sqrt(x) */ + q1 = 0; + s0 = 0; + s1 = 0; + r = 0x00200000; /* r = moving bit from right to left */ + + while r != 0 { + t = s0 + r as i32; + if t <= ix0 { + s0 = t + r as i32; + ix0 -= t; + q += r as i32; + } + ix0 += ix0 + ((ix1&sign)>>31) as i32; + ix1 += ix1; + r >>= 1; + } + + r = sign; + while r != 0 { + t1 = s1 + r; + t = s0; + if t < ix0 || (t == ix0 && t1 <= ix1) { + s1 = t1 + r; + if (t1&sign) == sign && (s1&sign) == 0 { + s0 += 1; + } + ix0 -= t; + if ix1 < t1 { + ix0 -= 1; + } + ix1 -= t1; + q1 += r; + } + ix0 += ix0 + ((ix1&sign)>>31) as i32; + ix1 += ix1; + r >>= 1; + } + + /* use floating add to find out rounding direction */ + if (ix0 as u32|ix1) != 0 { + z = 1.0 - TINY; /* raise inexact flag */ + if z >= 1.0 { + z = 1.0 + TINY; + if q1 == 0xffffffff { + q1 = 0; + q+=1; + } else if z > 1.0 { + if q1 == 0xfffffffe { + q += 1; + } + q1 += 2; + } else { + q1 += q1 & 1; + } + } + } + ix0 = (q>>1) + 0x3fe00000; + ix1 = q1>>1; + if (q&1) == 1 { + ix1 |= sign; + } + ix0 += m << 20; + f64::from_bits((ix0 as u64) << 32 | ix1 as u64) +} diff --git a/library/compiler-builtins/libm/test-generator/src/main.rs b/library/compiler-builtins/libm/test-generator/src/main.rs index 8da1a920eef2..ce20e80d1952 100644 --- a/library/compiler-builtins/libm/test-generator/src/main.rs +++ b/library/compiler-builtins/libm/test-generator/src/main.rs @@ -713,7 +713,7 @@ f64_f64! { round, // sin, // sinh, - // sqrt, + sqrt, // tan, // tanh, trunc, @@ -725,7 +725,7 @@ f64f64_f64! { // atan2, // fdim, // fmod, - // hypot, + hypot, // pow, }