Merge #71

71: implement sqrt and hypot r=japaric a=erikdesjardins

closes #57, closes #22

Co-authored-by: Erik <erikdesjardins@users.noreply.github.com>

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)

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);
+}

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,
 };
 

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)
+}

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,
 }