Add scalbnf16 and scalbnf128

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

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", "ldexpf16"],
+    ),
     (
         // `(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", "ldexpf128"],
+    ),
     (
         // `(f32, &mut f32) -> f32` as `(f32) -> (f32, f32)`
         FloatTy::F32,

crates/libm-test/benches/random.rs

@@ -127,8 +127,12 @@ libm_macros::for_each_function! {
         | fdimf16
         | floorf128
         | floorf16
+        | ldexpf128
+        | ldexpf16
         | rintf128
         | rintf16
+        | scalbnf128
+        | scalbnf16
         | sqrtf128
         | sqrtf16
         | truncf128

crates/libm-test/src/mpfloat.rs

@@ -158,7 +158,10 @@ libm_macros::for_each_function! {
         ilogbf,
         jn,
         jnf,
-        ldexp,ldexpf,
+        ldexp,
+        ldexpf,
+        ldexpf128,
+        ldexpf16,
         lgamma_r,
         lgammaf_r,
         modf,
@@ -176,6 +179,8 @@ libm_macros::for_each_function! {
         roundf,
         scalbn,
         scalbnf,
+        scalbnf128,
+        scalbnf16,
         sincos,sincosf,
         trunc,
         truncf,
@@ -367,34 +372,6 @@ macro_rules! impl_op_for_ty {
                 }
             }
 
-            // `ldexp` and `scalbn` are the same for binary floating point, so just forward all
-            // methods.
-            impl MpOp for crate::op::[<ldexp $suffix>]::Routine {
-                type MpTy = <crate::op::[<scalbn $suffix>]::Routine as MpOp>::MpTy;
-
-                fn new_mp() -> Self::MpTy {
-                    <crate::op::[<scalbn $suffix>]::Routine as MpOp>::new_mp()
-                }
-
-                fn run(this: &mut Self::MpTy, input: Self::RustArgs) -> Self::RustRet {
-                    <crate::op::[<scalbn $suffix>]::Routine as MpOp>::run(this, input)
-                }
-            }
-
-            impl MpOp for crate::op::[<scalbn $suffix>]::Routine {
-                type MpTy = MpFloat;
-
-                fn new_mp() -> Self::MpTy {
-                    new_mpfloat::<Self::FTy>()
-                }
-
-                fn run(this: &mut Self::MpTy, input: Self::RustArgs) -> Self::RustRet {
-                    this.assign(input.0);
-                    *this <<= input.1;
-                    prep_retval::<Self::FTy>(this, Ordering::Equal)
-                }
-            }
-
             impl MpOp for crate::op::[<sincos $suffix>]::Routine {
                 type MpTy = (MpFloat, MpFloat);
 
@@ -475,6 +452,34 @@ macro_rules! impl_op_for_ty_all {
                     prep_retval::<Self::RustRet>(&mut this.0, Ordering::Equal)
                 }
             }
+
+            // `ldexp` and `scalbn` are the same for binary floating point, so just forward all
+            // methods.
+            impl MpOp for crate::op::[<ldexp $suffix>]::Routine {
+                type MpTy = <crate::op::[<scalbn $suffix>]::Routine as MpOp>::MpTy;
+
+                fn new_mp() -> Self::MpTy {
+                    <crate::op::[<scalbn $suffix>]::Routine as MpOp>::new_mp()
+                }
+
+                fn run(this: &mut Self::MpTy, input: Self::RustArgs) -> Self::RustRet {
+                    <crate::op::[<scalbn $suffix>]::Routine as MpOp>::run(this, input)
+                }
+            }
+
+            impl MpOp for crate::op::[<scalbn $suffix>]::Routine {
+                type MpTy = MpFloat;
+
+                fn new_mp() -> Self::MpTy {
+                    new_mpfloat::<Self::FTy>()
+                }
+
+                fn run(this: &mut Self::MpTy, input: Self::RustArgs) -> Self::RustRet {
+                    this.assign(input.0);
+                    *this <<= input.1;
+                    prep_retval::<Self::FTy>(this, Ordering::Equal)
+                }
+            }
         }
     };
 }

crates/libm-test/src/precision.rs

@@ -585,8 +585,12 @@ fn int_float_common<F1: Float, F2: Float>(
     DEFAULT
 }
 
+#[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 {}
 
 impl MaybeOverride<(f32, f32, f32)> for SpecialCase {
     fn check_float<F: Float>(

crates/libm-test/tests/compare_built_musl.rs

@@ -89,8 +89,12 @@ libm_macros::for_each_function! {
         fdimf16,
         floorf128,
         floorf16,
+        ldexpf128,
+        ldexpf16,
         rintf128,
         rintf16,
+        scalbnf128,
+        scalbnf16,
         sqrtf128,
         sqrtf16,
         truncf128,

crates/util/src/main.rs

@@ -96,8 +96,12 @@ fn do_eval(basis: &str, op: &str, inputs: &[&str]) {
             | fdimf16
             | floorf128
             | floorf16
+            | ldexpf128
+            | ldexpf16
             | rintf128
             | rintf16
+            | scalbnf128
+            | scalbnf16
             | sqrtf128
             | sqrtf16
             | truncf128

etc/function-definitions.json

@@ -506,6 +506,18 @@
         ],
         "type": "f32"
     },
+    "ldexpf128": {
+        "sources": [
+            "src/math/ldexpf128.rs"
+        ],
+        "type": "f128"
+    },
+    "ldexpf16": {
+        "sources": [
+            "src/math/ldexpf16.rs"
+        ],
+        "type": "f16"
+    },
     "lgamma": {
         "sources": [
             "src/libm_helper.rs",
@@ -698,16 +710,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

@@ -73,6 +73,8 @@ jn
 jnf
 ldexp
 ldexpf
+ldexpf128
+ldexpf16
 lgamma
 lgamma_r
 lgammaf
@@ -103,6 +105,8 @@ round
 roundf
 scalbn
 scalbnf
+scalbnf128
+scalbnf16
 sin
 sincos
 sincosf

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,104 @@
+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>,
+{
+    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_bits((F::EXP_BIAS + sig_total_bits).cast() << F::SIG_BITS);
+
+    // std::println!("{exp_max} {exp_min} {n}");
+    // std::dbg!(x, exp_max, exp_min, n);
+
+    if n > exp_max {
+        x *= f_exp_max;
+        n -= exp_max;
+        // std::dbg!(11, x, n);
+        if n > exp_max {
+            x *= f_exp_max;
+            n -= exp_max;
+            // std::dbg!(12, x, n);
+            if n > exp_max {
+                n = exp_max;
+                // std::dbg!(13, x, n);
+            }
+        }
+    } 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;
+        // std::dbg!(21, x, n);
+        if n < exp_min {
+            x *= mul;
+            n += add;
+            // std::dbg!(22, x, n);
+            if n < exp_min {
+                n = exp_min;
+                // std::dbg!(23, x, n);
+            }
+        }
+    }
+
+    x * F::from_bits(F::Int::cast_from(F::EXP_BIAS as i32 + n) << F::SIG_BITS)
+}
+
+// DELETE
+
+extern crate std;
+
+#[test]
+fn testme() {
+    assert_eq!(scalbn::<f16>(f16::from_bits(0x6ecb), -1336428830), f16::ZERO);
+}
+
+#[test]
+fn testme2() {
+    // assert_eq!(scalbn(-f64::INFINITY, -2147033648), f64::ZERO);
+}

src/math/ldexpf128.rs

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

src/math/ldexpf16.rs

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