yggdrasil/libm
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Merge pull request #142 from mb64/wrapping-sqrt

Browse Source 
Change `sqrt` to use wrapping newtypes
This commit is contained in:
Alex Crichton
2019-05-02 14:30:49 -05:00
committed by GitHub
parent 45417e5e4d 339502c69e
commit 77a4032353
1 changed files with 32 additions and 31 deletions
Download Patch File Download Diff File Expand all files Collapse all files
src/math/sqrt.rs
+32 -31
View File
@@ -77,6 +77,7 @@
*/
use core::f64;
use core::num::Wrapping;
const TINY: f64 = 1.0e-300;
@@ -96,21 +97,21 @@ pub fn sqrt(x: f64) -> f64 {
}
}
let mut z: f64;
let sign: u32 = 0x80000000;
let sign: Wrapping<u32> = Wrapping(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;
let mut r: Wrapping<u32>;
let mut t1: Wrapping<u32>;
let mut s1: Wrapping<u32>;
let mut ix1: Wrapping<u32>;
let mut q1: Wrapping<u32>;
ix0 = (x.to_bits() >> 32) as i32;
ix1 = x.to_bits() as u32;
ix1 = Wrapping(x.to_bits() as u32);
/* take care of Inf and NaN */
if (ix0 & 0x7ff00000) == 0x7ff00000 {
@@ -118,7 +119,7 @@ pub fn sqrt(x: f64) -> f64 {
}
/* take care of zero */
if ix0 <= 0 {
if ((ix0 & !(sign as i32)) | ix1 as i32) == 0 {
if ((ix0 & !(sign.0 as i32)) | ix1.0 as i32) == 0 {
return x; /* sqrt(+-0) = +-0 */
}
if ix0 < 0 {
@@ -131,7 +132,7 @@ pub fn sqrt(x: f64) -> f64 {
/* subnormal x */
while ix0 == 0 {
m -= 21;
ix0 |= (ix1 >> 11) as i32;
ix0 |= (ix1 >> 11).0 as i32;
ix1 <<= 21;
}
i = 0;
@@ -140,46 +141,46 @@ pub fn sqrt(x: f64) -> f64 {
ix0 <<= 1;
}
m -= i - 1;
ix0 |= (ix1 >> (32 - i)) as i32;
ix1 <<= i;
ix0 |= (ix1 >> (32 - i) as usize).0 as i32;
ix1 = ix1 << i as usize;
}
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;
ix0 += ix0 + ((ix1 & sign) >> 31).0 as i32;
ix1 += ix1;
}
m >>= 1; /* m = [m/2] */
/* generate sqrt(x) bit by bit */
ix0 += ix0 + ((ix1 & sign) >> 31) as i32;
ix0 += ix0 + ((ix1 & sign) >> 31).0 as i32;
ix1 += ix1;
q = 0; /* [q,q1] = sqrt(x) */
q1 = 0;
q1 = Wrapping(0);
s0 = 0;
s1 = 0;
r = 0x00200000; /* r = moving bit from right to left */
s1 = Wrapping(0);
r = Wrapping(0x00200000); /* r = moving bit from right to left */
while r != 0 {
t = s0 + r as i32;
while r != Wrapping(0) {
t = s0 + r.0 as i32;
if t <= ix0 {
s0 = t + r as i32;
s0 = t + r.0 as i32;
ix0 -= t;
q += r as i32;
q += r.0 as i32;
}
ix0 += ix0 + ((ix1 & sign) >> 31) as i32;
ix0 += ix0 + ((ix1 & sign) >> 31).0 as i32;
ix1 += ix1;
r >>= 1;
}
r = sign;
while r != 0 {
while r != Wrapping(0) {
t1 = s1 + r;
t = s0;
if t < ix0 || (t == ix0 && t1 <= ix1) {
s1 = t1 + r;
if (t1 & sign) == sign && (s1 & sign) == 0 {
if (t1 & sign) == sign && (s1 & sign) == Wrapping(0) {
s0 += 1;
}
ix0 -= t;
@@ -189,26 +190,26 @@ pub fn sqrt(x: f64) -> f64 {
ix1 -= t1;
q1 += r;
}
ix0 += ix0 + ((ix1 & sign) >> 31) as i32;
ix0 += ix0 + ((ix1 & sign) >> 31).0 as i32;
ix1 += ix1;
r >>= 1;
}
/* use floating add to find out rounding direction */
if (ix0 as u32 | ix1) != 0 {
if (ix0 as u32 | ix1.0) != 0 {
z = 1.0 - TINY; /* raise inexact flag */
if z >= 1.0 {
z = 1.0 + TINY;
if q1 == 0xffffffff {
q1 = 0;
if q1.0 == 0xffffffff {
q1 = Wrapping(0);
q += 1;
} else if z > 1.0 {
if q1 == 0xfffffffe {
if q1.0 == 0xfffffffe {
q += 1;
}
q1 += 2;
q1 += Wrapping(2);
} else {
q1 += q1 & 1;
q1 += q1 & Wrapping(1);
}
}
}
@@ -218,5 +219,5 @@ pub fn sqrt(x: f64) -> f64 {
ix1 |= sign;
}
ix0 += m << 20;
f64::from_bits((ix0 as u64) << 32 | ix1 as u64)
f64::from_bits((ix0 as u64) << 32 | ix1.0 as u64)
}