Merge rust-lang#64

64: add exp and log/ln functions for f32 r=japaric a=burrbull

closes rust-lang#17 
closes rust-lang#29 
closes rust-lang#55 

Co-authored-by: Andrey Zgarbul <[email protected]>
Co-authored-by: Jorge Aparicio <[email protected]>

Author: 3 people

src/lib.rs

@@ -71,16 +71,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)]
@@ -220,7 +217,6 @@ impl F32Ext for f32 {
         sqrtf(self)
     }
 
-    #[cfg(todo)]
     #[inline]
     fn exp(self) -> Self {
         expf(self)
@@ -232,13 +228,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()
@@ -371,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)]
@@ -485,7 +478,6 @@ pub trait F64Ext {
 }
 
 impl F64Ext for f64 {
-    #[cfg(todo)]
     #[inline]
     fn floor(self) -> Self {
         floor(self)

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

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

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

src/math/mod.rs

@@ -9,16 +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;
+//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

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

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

test-generator/src/main.rs

@@ -570,10 +570,11 @@ f32_f32! {
     // cosf,
     // coshf,
     // exp2f,
-    // expf,
+    expf,
     // fdimf,
     // log10f,
     // log2f,
+    logf,
     // roundf,
     // sinf,
     // sinhf,
@@ -613,7 +614,7 @@ f64_f64! {
     // exp,
     // exp2,
     // expm1,
-    // floor,
+    floor,
     // log,
     // log10,
     // log1p,
@@ -644,5 +645,5 @@ f64f64f64_f64! {
 
 // With signature `fn(f64, i32) -> f64`
 f64i32_f64! {
-    // scalbn,
+    scalbn,
 }