Add a generic version of scalbn

Author: tgross35

crates/libm-test/src/mpfloat.rs

@@ -158,7 +158,8 @@ libm_macros::for_each_function! {
         ilogbf,
         jn,
         jnf,
-        ldexp,ldexpf,
+        ldexp,
+        ldexpf,
         lgamma_r,
         lgammaf_r,
         modf,

etc/function-definitions.json

@@ -698,12 +698,14 @@
     "scalbn": {
         "sources": [
             "src/libm_helper.rs",
+            "src/math/generic/scalbn.rs",
             "src/math/scalbn.rs"
         ],
         "type": "f64"
     },
     "scalbnf": {
         "sources": [
+            "src/math/generic/scalbn.rs",
             "src/math/scalbnf.rs"
         ],
         "type": "f32"

src/math/generic/mod.rs

@@ -4,6 +4,7 @@ mod fabs;
 mod fdim;
 mod floor;
 mod rint;
+mod scalbn;
 mod sqrt;
 mod trunc;
 
@@ -13,5 +14,6 @@ pub use fabs::fabs;
 pub use fdim::fdim;
 pub use floor::floor;
 pub use rint::rint;
+pub use scalbn::scalbn;
 pub use sqrt::sqrt;
 pub use trunc::trunc;

src/math/generic/scalbn.rs

@@ -0,0 +1,83 @@
+use super::super::{CastFrom, CastInto, Float, IntTy, MinInt};
+
+/// Scale the exponent.
+///
+/// From N3220:
+///
+/// > The scalbn and scalbln functions compute `x * b^n`, where `b = FLT_RADIX` if the return type
+/// > of the function is a standard floating type, or `b = 10` if the return type of the function
+/// > is a decimal floating type. A range error occurs for some finite x, depending on n.
+/// >
+/// > [...]
+/// >
+/// > * `scalbn(±0, n)` returns `±0`.
+/// > * `scalbn(x, 0)` returns `x`.
+/// > * `scalbn(±∞, n)` returns `±∞`.
+/// >
+/// > If the calculation does not overflow or underflow, the returned value is exact and
+/// > independent of the current rounding direction mode.
+pub fn scalbn<F: Float>(mut x: F, mut n: i32) -> F
+where
+    u32: CastInto<F::Int>,
+    F::Int: CastFrom<i32>,
+    F::Int: CastFrom<u32>,
+{
+    let zero = IntTy::<F>::ZERO;
+
+    // if n == 0 || x == F::ZERO || x.is_nan() || x.is_infinite() {
+    //     return x;
+    // }
+
+    // Bits including the implicit bit
+    let sig_total_bits = F::SIG_BITS + 1;
+
+    // Maximum and minimum values when biased
+    let exp_max: i32 = F::EXP_BIAS as i32;
+    let exp_min = -(exp_max - 1);
+    // let exp_min_with_subnorm = -((F::EXP_BIAS + F::SIG_BITS + 1) as i32);
+
+    // // let x_exp = x.exp();
+    // // let x_sig = x.frac();
+
+    // if n > exp_max {
+    //     return F::INFINITY * x.signum();
+    // }
+
+    // if n < exp_min_with_subnorm {
+    //     return F::ZERO * x.signum();
+    // }
+
+    // 2 ^ Emax, where Emax is the maximum biased exponent value (1023 for f64)
+    let f_exp_max = F::from_bits(F::Int::cast_from(F::EXP_BIAS << 1) << F::SIG_BITS);
+    // 2 ^ Emin, where Emin is the minimum biased exponent value (-1022 for f64)
+    let f_exp_min = F::from_bits(IntTy::<F>::ONE << F::SIG_BITS);
+    // 2 ^ sig_total_bits, representation of what can be accounted for with subnormals
+    let f_exp_subnorm = F::from_parts(false, sig_total_bits + F::EXP_BIAS, zero);
+
+    if n > exp_max {
+        x *= f_exp_max;
+        n -= exp_max;
+        if n > exp_max {
+            x *= f_exp_max;
+            n -= exp_max;
+            if n > exp_max {
+                n = exp_max;
+            }
+        }
+    } else if n < exp_min {
+        let mul = f_exp_min * f_exp_subnorm;
+        let add = (exp_max - 1) - sig_total_bits as i32;
+
+        x *= mul;
+        n += add;
+        if n < exp_min {
+            x *= mul;
+            n += add;
+            if n < exp_min {
+                n = exp_min;
+            }
+        }
+    }
+
+    x * F::from_bits(F::Int::cast_from(F::EXP_BIAS as i32 + n) << F::SIG_BITS)
+}

src/math/scalbn.rs

@@ -1,33 +1,4 @@
 #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
-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)
+pub fn scalbn(x: f64, n: i32) -> f64 {
+    super::generic::scalbn(x, n)
 }

src/math/scalbnf.rs

@@ -1,29 +1,4 @@
 #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
-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;
-        if n > 127 {
-            x *= x1p127;
-            n -= 127;
-            if n > 127 {
-                n = 127;
-            }
-        }
-    } else if n < -126 {
-        x *= x1p_126 * x1p24;
-        n += 126 - 24;
-        if n < -126 {
-            x *= x1p_126 * x1p24;
-            n += 126 - 24;
-            if n < -126 {
-                n = -126;
-            }
-        }
-    }
-    x * f32::from_bits(((0x7f + n) as u32) << 23)
+pub fn scalbnf(x: f32, n: i32) -> f32 {
+    super::generic::scalbn(x, n)
 }