Add a generic version of frexp and ilogb

This is used to replace `frexp`, `frexpf`, `ilogb`, and `ilogbf`.

Author: tgross35

etc/function-definitions.json

@@ -450,13 +450,15 @@
     "frexp": {
         "sources": [
             "src/libm_helper.rs",
-            "src/math/frexp.rs"
+            "src/math/frexp.rs",
+            "src/math/generic/frexp.rs"
         ],
         "type": "f64"
     },
     "frexpf": {
         "sources": [
-            "src/math/frexpf.rs"
+            "src/math/frexpf.rs",
+            "src/math/generic/frexp.rs"
         ],
         "type": "f32"
     },
@@ -476,12 +478,14 @@
     "ilogb": {
         "sources": [
             "src/libm_helper.rs",
+            "src/math/generic/ilogb.rs",
             "src/math/ilogb.rs"
         ],
         "type": "f64"
     },
     "ilogbf": {
         "sources": [
+            "src/math/generic/ilogb.rs",
             "src/math/ilogbf.rs"
         ],
         "type": "f32"

src/math/frexp.rs

@@ -1,20 +1,7 @@
+/// Decompose a float into a normalized value within the range `[0.5, 1)`, and a power of 2.
+///
+/// That is, `x * 2^p` will represent the input value.
+#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
 pub fn frexp(x: f64) -> (f64, i32) {
-    let mut y = x.to_bits();
-    let ee = ((y >> 52) & 0x7ff) as i32;
-
-    if ee == 0 {
-        if x != 0.0 {
-            let x1p64 = f64::from_bits(0x43f0000000000000);
-            let (x, e) = frexp(x * x1p64);
-            return (x, e - 64);
-        }
-        return (x, 0);
-    } else if ee == 0x7ff {
-        return (x, 0);
-    }
-
-    let e = ee - 0x3fe;
-    y &= 0x800fffffffffffff;
-    y |= 0x3fe0000000000000;
-    return (f64::from_bits(y), e);
+    super::generic::frexp(x)
 }

src/math/frexpf.rs

@@ -1,21 +1,7 @@
+/// Decompose a float into a normalized value within the range `[0.5, 1)`, and a power of 2.
+///
+/// That is, `x * 2^p` will represent the input value.
+#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
 pub fn frexpf(x: f32) -> (f32, i32) {
-    let mut y = x.to_bits();
-    let ee: i32 = ((y >> 23) & 0xff) as i32;
-
-    if ee == 0 {
-        if x != 0.0 {
-            let x1p64 = f32::from_bits(0x5f800000);
-            let (x, e) = frexpf(x * x1p64);
-            return (x, e - 64);
-        } else {
-            return (x, 0);
-        }
-    } else if ee == 0xff {
-        return (x, 0);
-    }
-
-    let e = ee - 0x7e;
-    y &= 0x807fffff;
-    y |= 0x3f000000;
-    (f32::from_bits(y), e)
+    super::generic::frexp(x)
 }

src/math/generic/frexp.rs

@@ -0,0 +1,25 @@
+use super::super::{CastFrom, Float, MinInt};
+
+pub fn frexp<F: Float>(x: F) -> (F, i32) {
+    let mut ix = x.to_bits();
+    let ee = x.exp();
+
+    if ee == 0 {
+        if x != F::ZERO {
+            // normalize via multiplication; 1p64 for `f64`
+            let magic = F::from_parts(false, F::EXP_BIAS + F::BITS, F::Int::ZERO);
+            magic.to_bits();
+
+            let (x, e) = frexp(x * magic);
+            return (x, e - F::BITS as i32);
+        }
+        return (x, 0);
+    } else if ee == F::EXP_MAX {
+        return (x, 0);
+    }
+
+    let e = ee as i32 - (F::EXP_BIAS as i32 - 1);
+    ix &= F::SIGN_MASK | F::SIG_MASK;
+    ix |= F::Int::cast_from(F::EXP_BIAS - 1) << F::SIG_BITS;
+    (F::from_bits(ix), e)
+}

src/math/generic/ilogb.rs

@@ -0,0 +1,30 @@
+use super::super::{Float, MinInt};
+
+const FP_ILOGBNAN: i32 = i32::MIN;
+const FP_ILOGB0: i32 = FP_ILOGBNAN;
+
+pub fn ilogb<F: Float>(x: F) -> i32 {
+    let zero = F::Int::ZERO;
+    let mut i = x.to_bits();
+    let e = x.exp() as i32;
+
+    if e == 0 {
+        i <<= F::EXP_BITS + 1;
+        if i == F::Int::ZERO {
+            force_eval!(0.0 / 0.0);
+            return FP_ILOGB0;
+        }
+        /* subnormal x */
+        let mut e = -(F::EXP_BIAS as i32);
+        while i >> (F::BITS - 1) == zero {
+            e -= 1;
+            i <<= 1;
+        }
+        e
+    } else if e == F::EXP_MAX as i32 {
+        force_eval!(0.0 / 0.0);
+        if i << (F::EXP_BITS + 1) != zero { FP_ILOGBNAN } else { i32::MAX }
+    } else {
+        e - F::EXP_BIAS as i32
+    }
+}

src/math/generic/mod.rs

@@ -5,6 +5,8 @@ mod fdim;
 mod floor;
 mod fmax;
 mod fmin;
+mod frexp;
+mod ilogb;
 mod rint;
 mod round;
 mod scalbn;
@@ -18,6 +20,8 @@ pub use fdim::fdim;
 pub use floor::floor;
 pub use fmax::fmax;
 pub use fmin::fmin;
+pub use frexp::frexp;
+pub use ilogb::ilogb;
 pub use rint::rint;
 pub use round::round;
 pub use scalbn::scalbn;

src/math/ilogb.rs

@@ -1,28 +1,5 @@
-const FP_ILOGBNAN: i32 = -1 - 0x7fffffff;
-const FP_ILOGB0: i32 = FP_ILOGBNAN;
-
+/// Extract the binary exponent of `x`.
 #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
 pub fn ilogb(x: f64) -> i32 {
-    let mut i: u64 = x.to_bits();
-    let e = ((i >> 52) & 0x7ff) as i32;
-
-    if e == 0 {
-        i <<= 12;
-        if i == 0 {
-            force_eval!(0.0 / 0.0);
-            return FP_ILOGB0;
-        }
-        /* subnormal x */
-        let mut e = -0x3ff;
-        while (i >> 63) == 0 {
-            e -= 1;
-            i <<= 1;
-        }
-        e
-    } else if e == 0x7ff {
-        force_eval!(0.0 / 0.0);
-        if (i << 12) != 0 { FP_ILOGBNAN } else { i32::MAX }
-    } else {
-        e - 0x3ff
-    }
+    super::generic::ilogb(x)
 }

src/math/ilogbf.rs

@@ -1,28 +1,5 @@
-const FP_ILOGBNAN: i32 = -1 - 0x7fffffff;
-const FP_ILOGB0: i32 = FP_ILOGBNAN;
-
+/// Extract the binary exponent of `x`.
 #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
 pub fn ilogbf(x: f32) -> i32 {
-    let mut i = x.to_bits();
-    let e = ((i >> 23) & 0xff) as i32;
-
-    if e == 0 {
-        i <<= 9;
-        if i == 0 {
-            force_eval!(0.0 / 0.0);
-            return FP_ILOGB0;
-        }
-        /* subnormal x */
-        let mut e = -0x7f;
-        while (i >> 31) == 0 {
-            e -= 1;
-            i <<= 1;
-        }
-        e
-    } else if e == 0xff {
-        force_eval!(0.0 / 0.0);
-        if (i << 9) != 0 { FP_ILOGBNAN } else { i32::MAX }
-    } else {
-        e - 0x7f
-    }
+    super::generic::ilogb(x)
 }

src/math/support/int_traits.rs

@@ -45,6 +45,8 @@ pub trait Int:
     + ops::BitOrAssign
     + ops::BitXorAssign
     + ops::ShlAssign<i32>
+    + ops::ShlAssign<u32>
+    + ops::ShrAssign<i32>
     + ops::ShrAssign<u32>
     + ops::Add<Output = Self>
     + ops::Sub<Output = Self>