Merge remote-tracking branch 'upstream/master'

library/compiler-builtins/libm/src/lib.rs

@@ -32,9 +32,8 @@ pub fn _eq(a: u64, b: u64) -> bool {
 
 /// Math support 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 F32Ext {
+/// This trait is sealed and cannot be implemented outside of `libm`.
+pub trait F32Ext: private::Sealed {
     #[cfg(todo)]
     fn floor(self) -> Self;
 
@@ -43,7 +42,6 @@ pub trait F32Ext {
 
     fn round(self) -> Self;
 
-    #[cfg(todo)]
     fn trunc(self) -> Self;
 
     #[cfg(todo)]
@@ -70,16 +68,13 @@ pub trait F32Ext {
 
     fn sqrt(self) -> Self;
 
-    #[cfg(todo)]
     fn exp(self) -> Self;
 
     #[cfg(todo)]
     fn exp2(self) -> Self;
 
-    #[cfg(todo)]
     fn ln(self) -> Self;
 
-    #[cfg(todo)]
     fn log(self, base: Self) -> Self;
 
     #[cfg(todo)]
@@ -91,7 +86,6 @@ pub trait F32Ext {
     #[cfg(todo)]
     fn cbrt(self) -> Self;
 
-    #[cfg(todo)]
     fn hypot(self, other: Self) -> Self;
 
     #[cfg(todo)]
@@ -164,7 +158,6 @@ impl F32Ext for f32 {
         roundf(self)
     }
 
-    #[cfg(todo)]
     #[inline]
     fn trunc(self) -> Self {
         truncf(self)
@@ -218,7 +211,6 @@ impl F32Ext for f32 {
         sqrtf(self)
     }
 
-    #[cfg(todo)]
     #[inline]
     fn exp(self) -> Self {
         expf(self)
@@ -230,13 +222,11 @@ impl F32Ext for f32 {
         exp2f(self)
     }
 
-    #[cfg(todo)]
     #[inline]
     fn ln(self) -> Self {
         logf(self)
     }
 
-    #[cfg(todo)]
     #[inline]
     fn log(self, base: Self) -> Self {
         self.ln() / base.ln()
@@ -260,7 +250,6 @@ impl F32Ext for f32 {
         cbrtf(self)
     }
 
-    #[cfg(todo)]
     #[inline]
     fn hypot(self, other: Self) -> Self {
         hypotf(self, other)
@@ -364,12 +353,10 @@ impl F32Ext for f32 {
     }
 }
 
-/// Math support for `f32`
+/// Math support for `f64`
 ///
-/// 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)]
+/// This trait is sealed and cannot be implemented outside of `libm`.
+pub trait F64Ext: private::Sealed {
     fn floor(self) -> Self;
 
     #[cfg(todo)]
@@ -377,7 +364,6 @@ pub trait F64Ext {
 
     fn round(self) -> Self;
 
-    #[cfg(todo)]
     fn trunc(self) -> Self;
 
     #[cfg(todo)]
@@ -403,7 +389,6 @@ pub trait F64Ext {
     #[cfg(todo)]
     fn powf(self, n: Self) -> Self;
 
-    #[cfg(todo)]
     fn sqrt(self) -> Self;
 
     #[cfg(todo)]
@@ -427,7 +412,6 @@ pub trait F64Ext {
     #[cfg(todo)]
     fn cbrt(self) -> Self;
 
-    #[cfg(todo)]
     fn hypot(self, other: Self) -> Self;
 
     #[cfg(todo)]
@@ -483,7 +467,6 @@ pub trait F64Ext {
 }
 
 impl F64Ext for f64 {
-    #[cfg(todo)]
     #[inline]
     fn floor(self) -> Self {
         floor(self)
@@ -500,7 +483,6 @@ impl F64Ext for f64 {
         round(self)
     }
 
-    #[cfg(todo)]
     #[inline]
     fn trunc(self) -> Self {
         trunc(self)
@@ -550,7 +532,6 @@ impl F64Ext for f64 {
         pow(self, n)
     }
 
-    #[cfg(todo)]
     #[inline]
     fn sqrt(self) -> Self {
         sqrt(self)
@@ -598,7 +579,6 @@ impl F64Ext for f64 {
         cbrt(self)
     }
 
-    #[cfg(todo)]
     #[inline]
     fn hypot(self, other: Self) -> Self {
         hypot(self, other)
@@ -701,3 +681,10 @@ impl F64Ext for f64 {
         0.5 * ((2.0 * self) / (1.0 - self)).ln_1p()
     }
 }
+
+mod private {
+    pub trait Sealed {}
+
+    impl Sealed for f32 {}
+    impl Sealed for f64 {}
+}

library/compiler-builtins/libm/src/math/expf.rs

@@ -0,0 +1,76 @@
+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 */
+/*
+ * 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 */
+
+#[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 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 > 0x7f800000 {/* NaN */
+            return x;
+        }
+        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.
+            }
+        }
+    }
+    
+    /* 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;
+        } else {
+            k = 1 - sign - sign;
+        }
+        let kf = k as f32;
+        hi = x - kf*LN2_HI;  /* k*ln2hi is exact here */
+        lo = kf*LN2_LO;
+        x = hi - lo;
+    } else if hx > 0x39000000 {  /* |x| > 2**-14 */
+        k = 0;
+        hi = x;
+        lo = 0.;
+    } else {
+        /* raise inexact */
+        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);
+    if k == 0 {
+        y
+    } else {
+        scalbnf(y, k)
+    }
+}

library/compiler-builtins/libm/src/math/floor.rs

@@ -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
+    }
+}

library/compiler-builtins/libm/src/math/hypot.rs

@@ -0,0 +1,74 @@
+use core::f64;
+
+use super::sqrt;
+
+const SPLIT: f64 = 134217728. + 1.; // 0x1p27 + 1 === (2 ^ 27) + 1
+
+fn sq(x: f64) -> (f64, f64) {
+    let xh: f64;
+    let xl: f64;
+    let xc: 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;
+    (hi, lo)
+}
+
+#[inline]
+pub fn hypot(mut x: f64, mut y: f64) -> f64 {
+    let x1p700 = f64::from_bits(0x6bb0000000000000); // 0x1p700 === 2 ^ 700
+    let x1p_700 = f64::from_bits(0x1430000000000000); // 0x1p-700 === 2 ^ -700
+
+    let mut uxi = x.to_bits();
+    let mut uyi = y.to_bits();
+    let uti;
+    let ex: i64;
+    let ey: i64;
+    let mut z: f64;
+
+    /* arrange |x| >= |y| */
+    uxi &= -1i64 as u64 >> 1;
+    uyi &= -1i64 as u64 >> 1;
+    if uxi < uyi {
+        uti = uxi;
+        uxi = uyi;
+        uyi = uti;
+    }
+
+    /* special cases */
+    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 */
+    if ey == 0x7ff {
+        return y;
+    }
+    if ex == 0x7ff || uyi == 0 {
+        return x;
+    }
+    /* note: hypot(x,y) ~= x + y*y/x/2 with inexact for small y/x */
+    /* 64 difference is enough for ld80 double_t */
+    if ex - ey > 64 {
+        return x + y;
+    }
+
+    /* 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 {
+        z = x1p700;
+        x *= x1p_700;
+        y *= x1p_700;
+    } 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);
+}

library/compiler-builtins/libm/src/math/hypotf.rs

@@ -0,0 +1,43 @@
+use core::f32;
+
+use super::sqrtf;
+
+#[inline]
+pub fn hypotf(mut x: f32, mut y: f32) -> f32 {
+    let x1p90 = f32::from_bits(0x6c800000); // 0x1p90f === 2 ^ 90
+    let x1p_90 = f32::from_bits(0x12800000); // 0x1p-90f === 2 ^ -90
+
+    let mut uxi = x.to_bits();
+    let mut uyi = y.to_bits();
+    let uti;
+    let mut z: f32;
+
+    uxi &= -1i32 as u32 >> 1;
+    uyi &= -1i32 as u32 >> 1;
+    if uxi < uyi {
+        uti = uxi;
+        uxi = uyi;
+        uyi = uti;
+    }
+
+    x = f32::from_bits(uxi);
+    y = f32::from_bits(uyi);
+    if uyi == 0xff<<23 {
+        return y;
+    }
+    if uxi >= 0xff<<23 || uyi == 0 || uxi - uyi >= 25<<23 {
+        return x + y;
+    }
+
+    z = 1.;
+    if uxi >= (0x7f+60)<<23 {
+        z = x1p90;
+        x *= x1p_90;
+        y *= x1p_90;
+    } 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)
+}

library/compiler-builtins/libm/src/math/logf.rs

@@ -0,0 +1,49 @@
+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 */
+
+#[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>>31) != 0 {
+            return (x-x)/0.;   /* log(-#) = NaN */
+        }
+        /* subnormal number, scale up x */
+        k -= 25;
+        x *= x1p25;
+        ix = x.to_bits();
+    } else if ix >= 0x7f800000 {
+        return x;
+    } else if ix == 0x3f800000 {
+        return 0.;
+    }
+
+    /* reduce x into [sqrt(2)/2, sqrt(2)] */
+    ix += 0x3f800000 - 0x3f3504f3;
+    k += ((ix>>23) as i32) - 0x7f;
+    ix = (ix & 0x007fffff) + 0x3f3504f3;
+    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 r = t2 + t1;
+    let hfsq = 0.5*f*f;
+    let dk = k as f32;
+    s*(hfsq+r) + dk*LN2_LO - hfsq + f + dk*LN2_HI
+}

library/compiler-builtins/libm/src/math/mod.rs

@@ -12,17 +12,39 @@ mod fmodf;
 mod powf;
 mod round;
 mod roundf;
+mod scalbn;
 mod scalbnf;
+mod sqrt;
 mod sqrtf;
+mod logf;
+mod expf;
+mod floor;
+mod trunc;
+mod truncf;
+mod hypot;
+mod hypotf;
 
-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::roundf::roundf;
-pub use self::scalbnf::scalbnf;
-pub use self::sqrtf::sqrtf;
+//mod service;
+
+pub use self::{
+    fabs::fabs,
+    fabsf::fabsf,
+    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,
+};
 
 fn isnanf(x: f32) -> bool {
     x.to_bits() & 0x7fffffff > 0x7f800000

library/compiler-builtins/libm/src/math/scalbn.rs

@@ -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)
+}

library/compiler-builtins/libm/src/math/scalbnf.rs

@@ -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)
 }

library/compiler-builtins/libm/src/math/sqrt.rs

@@ -0,0 +1,129 @@
+use core::f64;
+
+const TINY: f64 = 1.0e-300;
+
+#[inline]
+pub fn sqrt(x: f64) -> f64 {
+    let mut z: f64;
+    let sign: u32 = 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;
+
+    ix0 = (x.to_bits() >> 32) as i32;
+    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 */
+    }
+    /* take care of zero */
+    if ix0 <= 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 */
+        }
+    }
+    /* normalize x */
+    m = ix0>>20;
+    if m == 0 {  /* subnormal x */
+        while ix0 == 0 {
+            m -= 21;
+            ix0 |= (ix1>>11) as i32;
+            ix1 <<= 21;
+        }
+        i=0;
+        while (ix0&0x00100000) == 0 {
+            i += 1;
+            ix0 <<= 1;
+        }
+        m -= i - 1;
+        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;
+        ix1 += ix1;
+    }
+    m >>= 1;      /* m = [m/2] */
+
+    /* generate sqrt(x) bit by bit */
+    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 */
+
+    while r != 0 {
+        t = s0 + r as i32;
+        if t <= ix0 {
+            s0   = t + r as i32;
+            ix0 -= t;
+            q   += r as i32;
+        }
+        ix0 += ix0 + ((ix1&sign)>>31) as i32;
+        ix1 += ix1;
+        r >>= 1;
+    }
+
+    r = sign;
+    while r != 0 {
+        t1 = s1 + r;
+        t  = s0;
+        if t < ix0 || (t == ix0 && t1 <= ix1) {
+            s1 = t1 + r;
+            if (t1&sign) == sign && (s1&sign) == 0 {
+                s0 += 1;
+            }
+            ix0 -= t;
+            if ix1 < t1 {
+                ix0 -= 1;
+            }
+            ix1 -= t1;
+            q1 += r;
+        }
+        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 {
+        z = 1.0 - TINY; /* raise inexact flag */
+        if z >= 1.0 {
+            z = 1.0 + TINY;
+            if q1 == 0xffffffff {
+                q1 = 0;
+                q+=1;
+            } else if z > 1.0 {
+                if q1 == 0xfffffffe {
+                    q += 1;
+                }
+                q1 += 2;
+            } else {
+                q1 += q1 & 1;
+            }
+        }
+    }
+    ix0 = (q>>1) + 0x3fe00000;
+    ix1 = q1>>1;
+    if (q&1) == 1 {
+        ix1 |= sign;
+    }
+    ix0 += m << 20;
+    f64::from_bits((ix0 as u64) << 32 | ix1 as u64)
+}

library/compiler-builtins/libm/src/math/trunc.rs

@@ -0,0 +1,32 @@
+use core::f64;
+
+#[inline]
+pub fn trunc(x: f64) -> f64 {
+    let x1p120 = f64::from_bits(0x4770000000000000); // 0x1p120f === 2 ^ 120
+
+    let mut i: u64 = x.to_bits();
+    let mut e: i64 = (i >> 52 & 0x7ff) as i64 - 0x3ff + 12;
+    let m: u64;
+
+    if e >= 52 + 12 {
+        return x;
+    }
+    if e < 12 {
+        e = 1;
+    }
+    m = -1i64 as u64 >> e;
+    if (i & m) == 0 {
+        return x;
+    }
+    force_eval!(x + x1p120);
+    i &= !m;
+    f64::from_bits(i)
+}
+
+#[cfg(test)]
+mod tests {
+    #[test]
+    fn sanity_check() {
+        assert_eq!(super::trunc(1.1), 1.0);
+    }
+}

library/compiler-builtins/libm/src/math/truncf.rs

@@ -0,0 +1,32 @@
+use core::f32;
+
+#[inline]
+pub fn truncf(x: f32) -> f32 {
+    let x1p120 = f32::from_bits(0x7b800000); // 0x1p120f === 2 ^ 120
+
+    let mut i: u32 = x.to_bits();
+    let mut e: i32 = (i >> 23 & 0xff) as i32 - 0x7f + 9;
+    let m: u32;
+
+    if e >= 23 + 9 {
+        return x;
+    }
+    if e < 9 {
+        e = 1;
+    }
+    m = -1i32 as u32 >> e;
+    if (i & m) == 0 {
+        return x;
+    }
+    force_eval!(x + x1p120);
+    i &= !m;
+    f32::from_bits(i)
+}
+
+#[cfg(test)]
+mod tests {
+    #[test]
+    fn sanity_check() {
+        assert_eq!(super::truncf(1.1), 1.0);
+    }
+}