Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Change sqrt to use wrapping newtypes #142

Merged
merged 1 commit into from
May 2, 2019
Merged

Change sqrt to use wrapping newtypes #142

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
63 changes: 32 additions & 31 deletions src/math/sqrt.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -77,6 +77,7 @@
*/

use core::f64;
use core::num::Wrapping;

const TINY: f64 = 1.0e-300;

Expand All @@ -96,29 +97,29 @@ 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 {
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 {
if ((ix0 & !(sign.0 as i32)) | ix1.0 as i32) == 0 {
return x; /* sqrt(+-0) = +-0 */
}
if ix0 < 0 {
Expand All @@ -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;
Expand All @@ -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;
Expand All @@ -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);
}
}
}
Expand All @@ -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)
}