Merge rust-lang/libm#78
78: force PR code to be formatted r=japaric a=japaric Co-authored-by: Jorge Aparicio <jorge@japaric.io>
This commit is contained in:
commit
664a00f9ac
13 changed files with 155 additions and 149 deletions
|
|
@ -15,6 +15,8 @@ matrix:
|
|||
- env: TARGET=powerpc64-unknown-linux-gnu
|
||||
- env: TARGET=powerpc64le-unknown-linux-gnu
|
||||
- env: TARGET=x86_64-unknown-linux-gnu
|
||||
- env: TARGET=cargo-fmt
|
||||
rust: beta
|
||||
|
||||
before_install: set -e
|
||||
|
||||
|
|
|
|||
|
|
@ -1,6 +1,11 @@
|
|||
set -euxo pipefail
|
||||
|
||||
main() {
|
||||
if [ $TARGET = cargo-fmt ]; then
|
||||
rustup component add rustfmt-preview
|
||||
return
|
||||
fi
|
||||
|
||||
if ! hash cross >/dev/null 2>&1; then
|
||||
cargo install cross
|
||||
fi
|
||||
|
|
|
|||
|
|
@ -1,6 +1,11 @@
|
|||
set -euxo pipefail
|
||||
|
||||
main() {
|
||||
if [ $TARGET = cargo-fmt ]; then
|
||||
cargo fmt -- --check
|
||||
return
|
||||
fi
|
||||
|
||||
# quick check
|
||||
cargo check
|
||||
|
||||
|
|
|
|||
|
|
@ -1,60 +1,67 @@
|
|||
use super::scalbnf;
|
||||
|
||||
const HALF : [f32; 2] = [0.5,-0.5];
|
||||
const LN2_HI : f32 = 6.9314575195e-01; /* 0x3f317200 */
|
||||
const LN2_LO : f32 = 1.4286067653e-06; /* 0x35bfbe8e */
|
||||
const INV_LN2 : f32 = 1.4426950216e+00; /* 0x3fb8aa3b */
|
||||
const HALF: [f32; 2] = [0.5, -0.5];
|
||||
const LN2_HI: f32 = 6.9314575195e-01; /* 0x3f317200 */
|
||||
const LN2_LO: f32 = 1.4286067653e-06; /* 0x35bfbe8e */
|
||||
const INV_LN2: f32 = 1.4426950216e+00; /* 0x3fb8aa3b */
|
||||
/*
|
||||
* Domain [-0.34568, 0.34568], range ~[-4.278e-9, 4.447e-9]:
|
||||
* |x*(exp(x)+1)/(exp(x)-1) - p(x)| < 2**-27.74
|
||||
*/
|
||||
const P1 : f32 = 1.6666625440e-1; /* 0xaaaa8f.0p-26 */
|
||||
const P2 : f32 = -2.7667332906e-3; /* -0xb55215.0p-32 */
|
||||
const P1: f32 = 1.6666625440e-1; /* 0xaaaa8f.0p-26 */
|
||||
const P2: f32 = -2.7667332906e-3; /* -0xb55215.0p-32 */
|
||||
|
||||
#[inline]
|
||||
pub fn expf(mut x: f32) -> f32 {
|
||||
let x1p127 = f32::from_bits(0x7f000000); // 0x1p127f === 2 ^ 127
|
||||
let x1p_126 = f32::from_bits(0x800000); // 0x1p-126f === 2 ^ -126 /*original 0x1p-149f ??????????? */
|
||||
|
||||
let x1p127 = f32::from_bits(0x7f000000); // 0x1p127f === 2 ^ 127
|
||||
let x1p_126 = f32::from_bits(0x800000); // 0x1p-126f === 2 ^ -126 /*original 0x1p-149f ??????????? */
|
||||
|
||||
let mut hx = x.to_bits();
|
||||
let sign = (hx >> 31) as i32; /* sign bit of x */
|
||||
let signb : bool = sign != 0;
|
||||
hx &= 0x7fffffff; /* high word of |x| */
|
||||
|
||||
let sign = (hx >> 31) as i32; /* sign bit of x */
|
||||
let signb: bool = sign != 0;
|
||||
hx &= 0x7fffffff; /* high word of |x| */
|
||||
|
||||
/* special cases */
|
||||
if hx >= 0x42aeac50 { /* if |x| >= -87.33655f or NaN */
|
||||
if hx > 0x7f800000 {/* NaN */
|
||||
if hx >= 0x42aeac50 {
|
||||
/* if |x| >= -87.33655f or NaN */
|
||||
if hx > 0x7f800000 {
|
||||
/* NaN */
|
||||
return x;
|
||||
}
|
||||
if (hx >= 0x42b17218) && (!signb) { /* x >= 88.722839f */
|
||||
if (hx >= 0x42b17218) && (!signb) {
|
||||
/* x >= 88.722839f */
|
||||
/* overflow */
|
||||
x *= x1p127;
|
||||
return x;
|
||||
}
|
||||
if signb {
|
||||
/* underflow */
|
||||
force_eval!(-x1p_126/x);
|
||||
if hx >= 0x42cff1b5 { /* x <= -103.972084f */
|
||||
return 0.
|
||||
force_eval!(-x1p_126 / x);
|
||||
if hx >= 0x42cff1b5 {
|
||||
/* x <= -103.972084f */
|
||||
return 0.;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
/* argument reduction */
|
||||
let k : i32;
|
||||
let hi : f32;
|
||||
let lo : f32;
|
||||
if hx > 0x3eb17218 { /* if |x| > 0.5 ln2 */
|
||||
if hx > 0x3f851592 { /* if |x| > 1.5 ln2 */
|
||||
k = (INV_LN2*x + HALF[sign as usize]) as i32;
|
||||
let k: i32;
|
||||
let hi: f32;
|
||||
let lo: f32;
|
||||
if hx > 0x3eb17218 {
|
||||
/* if |x| > 0.5 ln2 */
|
||||
if hx > 0x3f851592 {
|
||||
/* if |x| > 1.5 ln2 */
|
||||
k = (INV_LN2 * x + HALF[sign as usize]) as i32;
|
||||
} else {
|
||||
k = 1 - sign - sign;
|
||||
}
|
||||
let kf = k as f32;
|
||||
hi = x - kf*LN2_HI; /* k*ln2hi is exact here */
|
||||
lo = kf*LN2_LO;
|
||||
hi = x - kf * LN2_HI; /* k*ln2hi is exact here */
|
||||
lo = kf * LN2_LO;
|
||||
x = hi - lo;
|
||||
} else if hx > 0x39000000 { /* |x| > 2**-14 */
|
||||
} else if hx > 0x39000000 {
|
||||
/* |x| > 2**-14 */
|
||||
k = 0;
|
||||
hi = x;
|
||||
lo = 0.;
|
||||
|
|
@ -63,11 +70,11 @@ pub fn expf(mut x: f32) -> f32 {
|
|||
force_eval!(x1p127 + x);
|
||||
return 1. + x;
|
||||
}
|
||||
|
||||
|
||||
/* x is now in primary range */
|
||||
let xx = x*x;
|
||||
let c = x - xx*(P1+xx*P2);
|
||||
let y = 1. + (x*c/(2.-c) - lo + hi);
|
||||
let xx = x * x;
|
||||
let c = x - xx * (P1 + xx * P2);
|
||||
let y = 1. + (x * c / (2. - c) - lo + hi);
|
||||
if k == 0 {
|
||||
y
|
||||
} else {
|
||||
|
|
|
|||
|
|
@ -1,27 +1,27 @@
|
|||
use core::f64;
|
||||
|
||||
const TOINT : f64 = 1. / f64::EPSILON;
|
||||
const TOINT: f64 = 1. / f64::EPSILON;
|
||||
|
||||
#[inline]
|
||||
pub fn floor(x : f64) -> f64 {
|
||||
pub fn floor(x: f64) -> f64 {
|
||||
let ui = x.to_bits();
|
||||
let e = ((ui >> 52) & 0x7ff) as i32;
|
||||
let e = ((ui >> 52) & 0x7ff) as i32;
|
||||
|
||||
if (e >= 0x3ff+52) || (x == 0.) {
|
||||
return x;
|
||||
if (e >= 0x3ff + 52) || (x == 0.) {
|
||||
return x;
|
||||
}
|
||||
/* y = int(x) - x, where int(x) is an integer neighbor of x */
|
||||
let y = if (ui >> 63) != 0 {
|
||||
x - TOINT + TOINT - x
|
||||
} else {
|
||||
x + TOINT - TOINT - x
|
||||
/* y = int(x) - x, where int(x) is an integer neighbor of x */
|
||||
let y = if (ui >> 63) != 0 {
|
||||
x - TOINT + TOINT - x
|
||||
} else {
|
||||
x + TOINT - TOINT - x
|
||||
};
|
||||
/* special case because of non-nearest rounding modes */
|
||||
if e <= 0x3ff-1 {
|
||||
force_eval!(y);
|
||||
return if (ui >> 63) != 0 { -1. } else { 0. };
|
||||
}
|
||||
if y > 0. {
|
||||
/* special case because of non-nearest rounding modes */
|
||||
if e <= 0x3ff - 1 {
|
||||
force_eval!(y);
|
||||
return if (ui >> 63) != 0 { -1. } else { 0. };
|
||||
}
|
||||
if y > 0. {
|
||||
x + y - 1.
|
||||
} else {
|
||||
x + y
|
||||
|
|
|
|||
|
|
@ -12,8 +12,8 @@ fn sq(x: f64) -> (f64, 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;
|
||||
let hi = x * x;
|
||||
let lo = xh * xh - hi + 2. * xh * xl + xl * xl;
|
||||
(hi, lo)
|
||||
}
|
||||
|
||||
|
|
@ -39,8 +39,8 @@ pub fn hypot(mut x: f64, mut y: f64) -> f64 {
|
|||
}
|
||||
|
||||
/* special cases */
|
||||
ex = (uxi>>52) as i64;
|
||||
ey = (uyi>>52) as i64;
|
||||
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 */
|
||||
|
|
@ -59,16 +59,16 @@ pub fn hypot(mut x: f64, mut y: f64) -> f64 {
|
|||
/* 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 {
|
||||
if ex > 0x3ff + 510 {
|
||||
z = x1p700;
|
||||
x *= x1p_700;
|
||||
y *= x1p_700;
|
||||
} else if ey < 0x3ff-450 {
|
||||
} 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);
|
||||
return z * sqrt(ly + lx + hy + hx);
|
||||
}
|
||||
|
|
|
|||
|
|
@ -22,22 +22,22 @@ pub fn hypotf(mut x: f32, mut y: f32) -> f32 {
|
|||
|
||||
x = f32::from_bits(uxi);
|
||||
y = f32::from_bits(uyi);
|
||||
if uyi == 0xff<<23 {
|
||||
if uyi == 0xff << 23 {
|
||||
return y;
|
||||
}
|
||||
if uxi >= 0xff<<23 || uyi == 0 || uxi - uyi >= 25<<23 {
|
||||
if uxi >= 0xff << 23 || uyi == 0 || uxi - uyi >= 25 << 23 {
|
||||
return x + y;
|
||||
}
|
||||
|
||||
z = 1.;
|
||||
if uxi >= (0x7f+60)<<23 {
|
||||
if uxi >= (0x7f + 60) << 23 {
|
||||
z = x1p90;
|
||||
x *= x1p_90;
|
||||
y *= x1p_90;
|
||||
} else if uyi < (0x7f-60)<<23 {
|
||||
} else if uyi < (0x7f - 60) << 23 {
|
||||
z = x1p_90;
|
||||
x *= x1p90;
|
||||
y *= x1p90;
|
||||
}
|
||||
z*sqrtf((x as f64 * x as f64 + y as f64 * y as f64) as f32)
|
||||
z * sqrtf((x as f64 * x as f64 + y as f64 * y as f64) as f32)
|
||||
}
|
||||
|
|
|
|||
|
|
@ -1,24 +1,25 @@
|
|||
const LN2_HI : f32 = 6.9313812256e-01; /* 0x3f317180 */
|
||||
const LN2_LO : f32 = 9.0580006145e-06; /* 0x3717f7d1 */
|
||||
const LN2_HI: f32 = 6.9313812256e-01; /* 0x3f317180 */
|
||||
const LN2_LO: f32 = 9.0580006145e-06; /* 0x3717f7d1 */
|
||||
/* |(log(1+s)-log(1-s))/s - Lg(s)| < 2**-34.24 (~[-4.95e-11, 4.97e-11]). */
|
||||
const LG1 : f32 = 0.66666662693; /* 0xaaaaaa.0p-24*/
|
||||
const LG2 : f32 = 0.40000972152; /* 0xccce13.0p-25 */
|
||||
const LG3 : f32 = 0.28498786688; /* 0x91e9ee.0p-25 */
|
||||
const LG4 : f32 = 0.24279078841; /* 0xf89e26.0p-26 */
|
||||
const LG1: f32 = 0.66666662693; /* 0xaaaaaa.0p-24*/
|
||||
const LG2: f32 = 0.40000972152; /* 0xccce13.0p-25 */
|
||||
const LG3: f32 = 0.28498786688; /* 0x91e9ee.0p-25 */
|
||||
const LG4: f32 = 0.24279078841; /* 0xf89e26.0p-26 */
|
||||
|
||||
#[inline]
|
||||
pub fn logf(mut x: f32) -> f32 {
|
||||
let x1p25 = f32::from_bits(0x4c000000); // 0x1p25f === 2 ^ 25
|
||||
|
||||
|
||||
let mut ix = x.to_bits();
|
||||
let mut k = 0i32;
|
||||
|
||||
if (ix < 0x00800000) || ((ix>>31) != 0) { /* x < 2**-126 */
|
||||
if ix<<1 == 0 {
|
||||
return -1./(x*x); /* log(+-0)=-inf */
|
||||
|
||||
if (ix < 0x00800000) || ((ix >> 31) != 0) {
|
||||
/* x < 2**-126 */
|
||||
if ix << 1 == 0 {
|
||||
return -1. / (x * x); /* log(+-0)=-inf */
|
||||
}
|
||||
if (ix>>31) != 0 {
|
||||
return (x-x)/0.; /* log(-#) = NaN */
|
||||
if (ix >> 31) != 0 {
|
||||
return (x - x) / 0.; /* log(-#) = NaN */
|
||||
}
|
||||
/* subnormal number, scale up x */
|
||||
k -= 25;
|
||||
|
|
@ -32,18 +33,18 @@ pub fn logf(mut x: f32) -> f32 {
|
|||
|
||||
/* reduce x into [sqrt(2)/2, sqrt(2)] */
|
||||
ix += 0x3f800000 - 0x3f3504f3;
|
||||
k += ((ix>>23) as i32) - 0x7f;
|
||||
k += ((ix >> 23) as i32) - 0x7f;
|
||||
ix = (ix & 0x007fffff) + 0x3f3504f3;
|
||||
x = f32::from_bits(ix);
|
||||
x = f32::from_bits(ix);
|
||||
|
||||
let f = x - 1.;
|
||||
let s = f/(2. + f);
|
||||
let z = s*s;
|
||||
let w = z*z;
|
||||
let t1 = w*(LG2+w*LG4);
|
||||
let t2 = z*(LG1+w*LG3);
|
||||
let s = f / (2. + f);
|
||||
let z = s * s;
|
||||
let w = z * z;
|
||||
let t1 = w * (LG2 + w * LG4);
|
||||
let t2 = z * (LG1 + w * LG3);
|
||||
let r = t2 + t1;
|
||||
let hfsq = 0.5*f*f;
|
||||
let hfsq = 0.5 * f * f;
|
||||
let dk = k as f32;
|
||||
s*(hfsq+r) + dk*LN2_LO - hfsq + f + dk*LN2_HI
|
||||
s * (hfsq + r) + dk * LN2_LO - hfsq + f + dk * LN2_HI
|
||||
}
|
||||
|
|
|
|||
|
|
@ -7,10 +7,15 @@ macro_rules! force_eval {
|
|||
}
|
||||
|
||||
mod ceilf;
|
||||
mod expf;
|
||||
mod fabs;
|
||||
mod fabsf;
|
||||
mod floor;
|
||||
mod floorf;
|
||||
mod fmodf;
|
||||
mod hypot;
|
||||
mod hypotf;
|
||||
mod logf;
|
||||
mod powf;
|
||||
mod round;
|
||||
mod roundf;
|
||||
|
|
@ -18,36 +23,15 @@ mod scalbn;
|
|||
mod scalbnf;
|
||||
mod sqrt;
|
||||
mod sqrtf;
|
||||
mod logf;
|
||||
mod expf;
|
||||
mod floor;
|
||||
mod trunc;
|
||||
mod truncf;
|
||||
mod hypot;
|
||||
mod hypotf;
|
||||
|
||||
//mod service;
|
||||
|
||||
pub use self::{
|
||||
ceilf::ceilf,
|
||||
fabs::fabs,
|
||||
fabsf::fabsf,
|
||||
floorf::floorf,
|
||||
fmodf::fmodf,
|
||||
powf::powf,
|
||||
round::round,
|
||||
roundf::roundf,
|
||||
scalbn::scalbn,
|
||||
scalbnf::scalbnf,
|
||||
sqrt::sqrt,
|
||||
sqrtf::sqrtf,
|
||||
logf::logf,
|
||||
expf::expf,
|
||||
floor::floor,
|
||||
trunc::trunc,
|
||||
truncf::truncf,
|
||||
hypot::hypot,
|
||||
hypotf::hypotf,
|
||||
ceilf::ceilf, expf::expf, fabs::fabs, fabsf::fabsf, floor::floor, floorf::floorf, fmodf::fmodf,
|
||||
hypot::hypot, hypotf::hypotf, logf::logf, powf::powf, round::round, roundf::roundf,
|
||||
scalbn::scalbn, scalbnf::scalbnf, sqrt::sqrt, sqrtf::sqrtf, trunc::trunc, truncf::truncf,
|
||||
};
|
||||
|
||||
fn isnanf(x: f32) -> bool {
|
||||
|
|
|
|||
|
|
@ -26,7 +26,7 @@ pub fn round(mut x: f64) -> f64 {
|
|||
} else {
|
||||
y = y + x;
|
||||
}
|
||||
|
||||
|
||||
if i >> 63 != 0 {
|
||||
-y
|
||||
} else {
|
||||
|
|
|
|||
|
|
@ -1,11 +1,11 @@
|
|||
#[inline]
|
||||
pub fn scalbn(x : f64, mut n: i32) -> f64 {
|
||||
let x1p1023 = f64::from_bits(0x7fe0000000000000); // 0x1p1023 === 2 ^ 1023
|
||||
let x1p53 = f64::from_bits(0x4340000000000000); // 0x1p53 === 2 ^ 53
|
||||
pub fn scalbn(x: f64, mut n: i32) -> f64 {
|
||||
let x1p1023 = f64::from_bits(0x7fe0000000000000); // 0x1p1023 === 2 ^ 1023
|
||||
let x1p53 = f64::from_bits(0x4340000000000000); // 0x1p53 === 2 ^ 53
|
||||
let x1p_1022 = f64::from_bits(0x0010000000000000); // 0x1p-1022 === 2 ^ (-1022)
|
||||
|
||||
|
||||
let mut y = x;
|
||||
|
||||
|
||||
if n > 1023 {
|
||||
y *= x1p1023;
|
||||
n -= 1023;
|
||||
|
|
@ -29,5 +29,5 @@ pub fn scalbn(x : f64, mut n: i32) -> f64 {
|
|||
}
|
||||
}
|
||||
}
|
||||
y*f64::from_bits(((0x3ff+n) as u64)<<52)
|
||||
y * f64::from_bits(((0x3ff + n) as u64) << 52)
|
||||
}
|
||||
|
|
|
|||
|
|
@ -1,9 +1,9 @@
|
|||
#[inline]
|
||||
pub fn scalbnf(mut x: f32, mut n : i32) -> f32 {
|
||||
let x1p127 = f32::from_bits(0x7f000000); // 0x1p127f === 2 ^ 127
|
||||
let x1p_126 = f32::from_bits(0x800000); // 0x1p-126f === 2 ^ -126
|
||||
let x1p24 = f32::from_bits(0x4b800000); // 0x1p24f === 2 ^ 24
|
||||
|
||||
pub fn scalbnf(mut x: f32, mut n: i32) -> f32 {
|
||||
let x1p127 = f32::from_bits(0x7f000000); // 0x1p127f === 2 ^ 127
|
||||
let x1p_126 = f32::from_bits(0x800000); // 0x1p-126f === 2 ^ -126
|
||||
let x1p24 = f32::from_bits(0x4b800000); // 0x1p24f === 2 ^ 24
|
||||
|
||||
if n > 127 {
|
||||
x *= x1p127;
|
||||
n -= 127;
|
||||
|
|
@ -25,5 +25,5 @@ pub fn scalbnf(mut x: f32, mut n : i32) -> f32 {
|
|||
}
|
||||
}
|
||||
}
|
||||
x * f32::from_bits(((0x7f+n) as u32)<<23)
|
||||
x * f32::from_bits(((0x7f + n) as u32) << 23)
|
||||
}
|
||||
|
|
|
|||
|
|
@ -22,60 +22,62 @@ pub fn sqrt(x: f64) -> f64 {
|
|||
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 */
|
||||
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 & !(sign as i32)) | ix1 as i32) == 0 {
|
||||
return x; /* sqrt(+-0) = +-0 */
|
||||
}
|
||||
if ix0 < 0 {
|
||||
return (x - x) / (x - x); /* sqrt(-ve) = sNaN */
|
||||
return (x - x) / (x - x); /* sqrt(-ve) = sNaN */
|
||||
}
|
||||
}
|
||||
/* normalize x */
|
||||
m = ix0>>20;
|
||||
if m == 0 { /* subnormal x */
|
||||
m = ix0 >> 20;
|
||||
if m == 0 {
|
||||
/* subnormal x */
|
||||
while ix0 == 0 {
|
||||
m -= 21;
|
||||
ix0 |= (ix1>>11) as i32;
|
||||
ix0 |= (ix1 >> 11) as i32;
|
||||
ix1 <<= 21;
|
||||
}
|
||||
i=0;
|
||||
while (ix0&0x00100000) == 0 {
|
||||
i = 0;
|
||||
while (ix0 & 0x00100000) == 0 {
|
||||
i += 1;
|
||||
ix0 <<= 1;
|
||||
}
|
||||
m -= i - 1;
|
||||
ix0 |= (ix1>>(32-i)) as i32;
|
||||
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;
|
||||
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] */
|
||||
m >>= 1; /* m = [m/2] */
|
||||
|
||||
/* generate sqrt(x) bit by bit */
|
||||
ix0 += ix0 + ((ix1&sign)>>31) as i32;
|
||||
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 */
|
||||
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;
|
||||
s0 = t + r as i32;
|
||||
ix0 -= t;
|
||||
q += r as i32;
|
||||
q += r as i32;
|
||||
}
|
||||
ix0 += ix0 + ((ix1&sign)>>31) as i32;
|
||||
ix0 += ix0 + ((ix1 & sign) >> 31) as i32;
|
||||
ix1 += ix1;
|
||||
r >>= 1;
|
||||
}
|
||||
|
|
@ -83,10 +85,10 @@ pub fn sqrt(x: f64) -> f64 {
|
|||
r = sign;
|
||||
while r != 0 {
|
||||
t1 = s1 + r;
|
||||
t = s0;
|
||||
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) == 0 {
|
||||
s0 += 1;
|
||||
}
|
||||
ix0 -= t;
|
||||
|
|
@ -96,19 +98,19 @@ pub fn sqrt(x: f64) -> f64 {
|
|||
ix1 -= t1;
|
||||
q1 += r;
|
||||
}
|
||||
ix0 += ix0 + ((ix1&sign)>>31) as i32;
|
||||
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 {
|
||||
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;
|
||||
q += 1;
|
||||
} else if z > 1.0 {
|
||||
if q1 == 0xfffffffe {
|
||||
q += 1;
|
||||
|
|
@ -119,9 +121,9 @@ pub fn sqrt(x: f64) -> f64 {
|
|||
}
|
||||
}
|
||||
}
|
||||
ix0 = (q>>1) + 0x3fe00000;
|
||||
ix1 = q1>>1;
|
||||
if (q&1) == 1 {
|
||||
ix0 = (q >> 1) + 0x3fe00000;
|
||||
ix1 = q1 >> 1;
|
||||
if (q & 1) == 1 {
|
||||
ix1 |= sign;
|
||||
}
|
||||
ix0 += m << 20;
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue