Add scalbnf16 and scalbnf128

The `scalbn` functions are similar enough that they can easily be made
generic. Do so and `f16` and `f128` versions.

Author: tgross35

crates/libm-macros/src/shared.rs

@@ -134,6 +134,13 @@ const ALL_OPERATIONS_NESTED: &[(FloatTy, Signature, Option<Signature>, &[&str])]
         None,
         &["jn", "yn"],
     ),
+    (
+        // `(f16, i32) -> f16`
+        FloatTy::F16,
+        Signature { args: &[Ty::F16, Ty::I32], returns: &[Ty::F16] },
+        None,
+        &["scalbnf16"],
+    ),
     (
         // `(f32, i32) -> f32`
         FloatTy::F32,
@@ -148,6 +155,13 @@ const ALL_OPERATIONS_NESTED: &[(FloatTy, Signature, Option<Signature>, &[&str])]
         None,
         &["scalbn", "ldexp"],
     ),
+    (
+        // `(f128, i32) -> f128`
+        FloatTy::F128,
+        Signature { args: &[Ty::F128, Ty::I32], returns: &[Ty::F128] },
+        None,
+        &["scalbnf128"],
+    ),
     (
         // `(f32, &mut f32) -> f32` as `(f32) -> (f32, f32)`
         FloatTy::F32,

crates/libm-test/benches/random.rs

@@ -117,7 +117,7 @@ libm_macros::for_each_function! {
         exp10 | exp10f | exp2 | exp2f => (true, Some(musl_math_sys::MACRO_FN_NAME)),
 
         // Musl does not provide `f16` and `f128` functions
-        copysignf16 | copysignf128 | fabsf16 | fabsf128 => (false, None),
+        copysignf16 | copysignf128 | fabsf16 | fabsf128 | scalbnf16 | scalbnf128 => (false, None),
 
         // By default we never skip (false) and always have a musl function available
         _ => (false, Some(musl_math_sys::MACRO_FN_NAME))

crates/libm-test/src/mpfloat.rs

@@ -137,7 +137,7 @@ libm_macros::for_each_function! {
         fmod, fmodf, frexp, frexpf, ilogb, ilogbf, jn, jnf, ldexp, ldexpf,
         lgamma_r, lgammaf_r, modf, modff, nextafter, nextafterf, pow,powf,
         remquo, remquof, scalbn, scalbnf, sincos, sincosf, yn, ynf,
-        copysignf16, copysignf128, fabsf16, fabsf128,
+        copysignf16, copysignf128, fabsf16, fabsf128, scalbnf16, scalbnf128,
     ],
     fn_extra: match MACRO_FN_NAME {
         // Remap function names that are different between mpfr and libm

crates/libm-test/src/precision.rs

@@ -403,5 +403,9 @@ fn bessel_prec_dropoff<F: Float>(
 
 impl MaybeOverride<(f32, f32, f32)> for SpecialCase {}
 impl MaybeOverride<(f64, f64, f64)> for SpecialCase {}
+#[cfg(f16_enabled)]
+impl MaybeOverride<(f16, i32)> for SpecialCase {}
 impl MaybeOverride<(f32, i32)> for SpecialCase {}
 impl MaybeOverride<(f64, i32)> for SpecialCase {}
+#[cfg(f128_enabled)]
+impl MaybeOverride<(f128, i32)> for SpecialCase {}

crates/libm-test/tests/compare_built_musl.rs

@@ -47,7 +47,7 @@ where
 libm_macros::for_each_function! {
     callback: musl_rand_tests,
     // Musl does not support `f16` and `f128` on all platforms.
-    skip: [copysignf16, copysignf128, fabsf16, fabsf128],
+    skip: [copysignf16, copysignf128, fabsf16, fabsf128, scalbnf16, scalbnf128],
     attributes: [
         #[cfg_attr(x86_no_sse, ignore)] // FIXME(correctness): wrong result on i586
         [exp10, exp10f, exp2, exp2f, rint]

crates/libm-test/tests/multiprecision.rs

@@ -154,6 +154,8 @@ libm_macros::for_each_function! {
         remquof,
         scalbn,
         scalbnf,
+        scalbnf16,
+        scalbnf128,
         yn,
         ynf,
 

crates/libm-test/tests/z_extensive/run.rs

@@ -61,6 +61,8 @@ fn register_all_tests() -> Vec<Trial> {
             remquof,
             scalbn,
             scalbnf,
+            scalbnf16,
+            scalbnf128,
 
             // FIXME: test needed, see
             // https://github.com/rust-lang/libm/pull/311#discussion_r1818273392

etc/function-definitions.json

@@ -630,16 +630,32 @@
     "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"
     },
+    "scalbnf128": {
+        "sources": [
+            "src/math/generic/scalbn.rs",
+            "src/math/scalbnf128.rs"
+        ],
+        "type": "f128"
+    },
+    "scalbnf16": {
+        "sources": [
+            "src/math/generic/scalbn.rs",
+            "src/math/scalbnf16.rs"
+        ],
+        "type": "f16"
+    },
     "sin": {
         "sources": [
             "src/libm_helper.rs",

etc/function-list.txt

@@ -95,6 +95,8 @@ round
 roundf
 scalbn
 scalbnf
+scalbnf128
+scalbnf16
 sin
 sincos
 sincosf

src/math/generic/mod.rs

@@ -1,5 +1,7 @@
 mod copysign;
 mod fabs;
+mod scalbn;
 
 pub use copysign::copysign;
 pub use fabs::fabs;
+pub use scalbn::scalbn;

src/math/generic/scalbn.rs

@@ -0,0 +1,49 @@
+use super::super::{CastFrom, CastInto, Float, IntTy, MinInt};
+
+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>,
+{
+    // 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);
+
+    // 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_bits((F::EXP_BIAS + sig_total_bits).cast() << F::SIG_BITS);
+
+    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/mod.rs

@@ -114,7 +114,7 @@ use self::rem_pio2::rem_pio2;
 use self::rem_pio2_large::rem_pio2_large;
 use self::rem_pio2f::rem_pio2f;
 #[allow(unused_imports)]
-use self::support::{CastFrom, CastInto, DInt, Float, HInt, Int, MinInt};
+use self::support::{CastFrom, CastInto, DInt, Float, HInt, Int, IntTy, MinInt};
 
 // Public modules
 mod acos;
@@ -336,19 +336,23 @@ cfg_if! {
     if #[cfg(f16_enabled)] {
         mod copysignf16;
         mod fabsf16;
+        mod scalbnf16;
 
         pub use self::copysignf16::copysignf16;
         pub use self::fabsf16::fabsf16;
+        pub use self::scalbnf16::scalbnf16;
     }
 }
 
 cfg_if! {
     if #[cfg(f128_enabled)] {
         mod copysignf128;
         mod fabsf128;
+        mod scalbnf128;
 
         pub use self::copysignf128::copysignf128;
         pub use self::fabsf128::fabsf128;
+        pub use self::scalbnf128::scalbnf128;
     }
 }
 

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

src/math/scalbnf128.rs

@@ -0,0 +1,4 @@
+#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
+pub fn scalbnf128(x: f128, n: i32) -> f128 {
+    super::generic::scalbn(x, n)
+}

src/math/scalbnf16.rs

@@ -0,0 +1,4 @@
+#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
+pub fn scalbnf16(x: f16, n: i32) -> f16 {
+    super::generic::scalbn(x, n)
+}