diff --git a/src/math/sqrt.rs b/src/math/sqrt.rs index a05a521fb..14404d4eb 100644 --- a/src/math/sqrt.rs +++ b/src/math/sqrt.rs @@ -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 = 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; + let mut t1: Wrapping; + let mut s1: Wrapping; + let mut ix1: Wrapping; + let mut q1: Wrapping; 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) }