add floor, scalbn; fixes in expf, scalbnf

This commit is contained in:
Andrey Zgarbul 2018-07-13 21:02:11 +03:00
parent 505e4ce81c
commit 7e419befa1
7 changed files with 94 additions and 32 deletions

View file

@ -365,7 +365,6 @@ impl F32Ext for f32 {
/// NOTE this meant to be a closed extension trait. The only stable way to use this trait is to
/// import it to access its methods.
pub trait F64Ext {
#[cfg(todo)]
fn floor(self) -> Self;
#[cfg(todo)]
@ -479,7 +478,6 @@ pub trait F64Ext {
}
impl F64Ext for f64 {
#[cfg(todo)]
#[inline]
fn floor(self) -> Self {
floor(self)

View file

@ -13,16 +13,16 @@ 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() as i32;
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| */
/* special cases */
if hx >= 0x42aeac50 { /* if |x| >= -87.33655f or NaN */
if hx >= 0x42aeac50 { /* if |x| >= -87.33655f or NaN */
if hx > 0x7f800000 {/* NaN */
return x;
}

View file

@ -0,0 +1,29 @@
use core::f64;
const TOINT : f64 = 1. / f64::EPSILON;
#[inline]
pub fn floor(x : f64) -> f64 {
let ui = x.to_bits();
let e = ((ui >> 52) & 0x7ff) as i32;
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
};
/* 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
}
}

View file

@ -8,7 +8,7 @@ const LG4 : f32 = 0.24279078841; /* 0xf89e26.0p-26 */
#[inline]
pub fn logf(mut x: f32) -> f32 {
let x1p25 = f32::from_bits(0x4c000000); // 0x1p24f === 2 ^ 25
let x1p25 = f32::from_bits(0x4c000000); // 0x1p25f === 2 ^ 25
let mut ix = x.to_bits();
let mut k = 0i32;

View file

@ -9,20 +9,28 @@ mod fabsf;
mod fmodf;
mod powf;
mod round;
mod scalbn;
mod scalbnf;
mod sqrtf;
mod logf;
mod expf;
mod floor;
pub use self::fabs::fabs;
pub use self::fabsf::fabsf;
pub use self::fmodf::fmodf;
pub use self::powf::powf;
pub use self::round::round;
pub use self::scalbnf::scalbnf;
pub use self::sqrtf::sqrtf;
pub use self::logf::logf;
pub use self::expf::expf;
//mod service;
pub use self::{
fabs::fabs,
fabsf::fabsf,
fmodf::fmodf,
powf::powf,
round::round,
scalbn::scalbn,
scalbnf::scalbnf,
sqrtf::sqrtf,
logf::logf,
expf::expf,
floor::floor,
};
fn isnanf(x: f32) -> bool {
x.to_bits() & 0x7fffffff > 0x7f800000

View file

@ -0,0 +1,33 @@
#[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
let x1p_1022 = f64::from_bits(0x0010000000000000); // 0x1p-1022 === 2 ^ (-1022)
let mut y = x;
if n > 1023 {
y *= x1p1023;
n -= 1023;
if n > 1023 {
y *= x1p1023;
n -= 1023;
if n > 1023 {
n = 1023;
}
}
} else if n < -1022 {
/* make sure final n < -53 to avoid double
rounding in the subnormal range */
y *= x1p_1022 * x1p53;
n += 1022 - 53;
if n < -1022 {
y *= x1p_1022 * x1p53;
n += 1022 - 53;
if n < -1022 {
n = -1022;
}
}
}
y*f64::from_bits(((0x3ff+n) as u64)<<52)
}

View file

@ -1,35 +1,29 @@
#[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
let mut y: f32 = x;
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 {
y *= x1p127;
x *= x1p127;
n -= 127;
if n > 127 {
y *= x1p127;
x *= x1p127;
n -= 127;
if n > 127 {
n = 127;
}
}
} else if n < -126 {
y *= x1p_126;
y *= x1p24;
x *= x1p_126 * x1p24;
n += 126 - 24;
if n < -126 {
y *= x1p_126;
y *= x1p24;
x *= x1p_126 * x1p24;
n += 126 - 24;
if n < -126 {
n = -126;
}
}
}
x = y * f32::from_bits((0x7f + n as u32) << 23);
x
x * f32::from_bits(((0x7f+n) as u32)<<23)
}