Add frexpf16, frexpf128, ilogbf16, and ilogbf128

crates/libm-macros/src/shared.rs

@@ -106,6 +106,13 @@ const ALL_OPERATIONS_NESTED: &[(FloatTy, Signature, Option<Signature>, &[&str])]
         None,
         &["fma"],
     ),
+    (
+        // `(f16) -> i32`
+        FloatTy::F16,
+        Signature { args: &[Ty::F16], returns: &[Ty::I32] },
+        None,
+        &["ilogbf16"],
+    ),
     (
         // `(f32) -> i32`
         FloatTy::F32,
@@ -120,6 +127,13 @@ const ALL_OPERATIONS_NESTED: &[(FloatTy, Signature, Option<Signature>, &[&str])]
         None,
         &["ilogb"],
     ),
+    (
+        // `(f128) -> i32`
+        FloatTy::F128,
+        Signature { args: &[Ty::F128], returns: &[Ty::I32] },
+        None,
+        &["ilogbf128"],
+    ),
     (
         // `(i32, f32) -> f32`
         FloatTy::F32,
@@ -162,6 +176,13 @@ const ALL_OPERATIONS_NESTED: &[(FloatTy, Signature, Option<Signature>, &[&str])]
         Some(Signature { args: &[Ty::F64, Ty::MutF64], returns: &[Ty::F64] }),
         &["modf"],
     ),
+    (
+        // `(f16, &mut c_int) -> f16` as `(f16) -> (f16, i32)`
+        FloatTy::F16,
+        Signature { args: &[Ty::F16], returns: &[Ty::F16, Ty::I32] },
+        Some(Signature { args: &[Ty::F16, Ty::MutCInt], returns: &[Ty::F16] }),
+        &["frexpf16"],
+    ),
     (
         // `(f32, &mut c_int) -> f32` as `(f32) -> (f32, i32)`
         FloatTy::F32,
@@ -176,6 +197,13 @@ const ALL_OPERATIONS_NESTED: &[(FloatTy, Signature, Option<Signature>, &[&str])]
         Some(Signature { args: &[Ty::F64, Ty::MutCInt], returns: &[Ty::F64] }),
         &["frexp", "lgamma_r"],
     ),
+    (
+        // `(f128, &mut c_int) -> f128` as `(f128) -> (f128, i32)`
+        FloatTy::F128,
+        Signature { args: &[Ty::F128], returns: &[Ty::F128, Ty::I32] },
+        Some(Signature { args: &[Ty::F128, Ty::MutCInt], returns: &[Ty::F128] }),
+        &["frexpf128"],
+    ),
     (
         // `(f32, f32, &mut c_int) -> f32` as `(f32, f32) -> (f32, i32)`
         FloatTy::F32,

crates/libm-test/benches/random.rs

@@ -133,6 +133,10 @@ libm_macros::for_each_function! {
         | fminf16
         | fmodf128
         | fmodf16
+        | frexpf128
+        | frexpf16
+        | ilogbf128
+        | ilogbf16
         | rintf128
         | rintf16
         | roundf128

crates/libm-test/src/mpfloat.rs

@@ -153,8 +153,12 @@ libm_macros::for_each_function! {
         fmodf16,
         frexp,
         frexpf,
+        frexpf128,
+        frexpf16,
         ilogb,
         ilogbf,
+        ilogbf128,
+        ilogbf16,
         jn,
         jnf,
         ldexp,
@@ -299,43 +303,6 @@ macro_rules! impl_op_for_ty {
                 }
             }
 
-            impl MpOp for crate::op::[<frexp $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);
-                    let exp = this.frexp_mut();
-                    (prep_retval::<Self::FTy>(this, Ordering::Equal), exp)
-                }
-            }
-
-            impl MpOp for crate::op::[<ilogb $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);
-
-                    // `get_exp` follows `frexp` for `0.5 <= |m| < 1.0`. Adjust the exponent by
-                    // one to scale the significand to `1.0 <= |m| < 2.0`.
-                    this.get_exp().map(|v| v - 1).unwrap_or_else(|| {
-                        if this.is_infinite() {
-                            i32::MAX
-                        } else {
-                            // Zero or NaN
-                            i32::MIN
-                        }
-                    })
-                }
-            }
-
             impl MpOp for crate::op::[<jn $suffix>]::Routine {
                 type MpTy = MpFloat;
 
@@ -466,6 +433,48 @@ macro_rules! impl_op_for_ty_all {
                     prep_retval::<Self::RustRet>(&mut this.0, ord)
                 }
             }
+
+            impl MpOp for crate::op::[<frexp $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 {
+                    // Implementation taken from `rug::Float::to_f32_exp`.
+                    this.assign(input.0);
+                    let exp = this.get_exp().unwrap_or(0);
+                    if exp != 0 {
+                        *this >>= exp;
+                    }
+
+                    (prep_retval::<Self::FTy>(this, Ordering::Equal), exp)
+                }
+            }
+
+            impl MpOp for crate::op::[<ilogb $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);
+
+                    // `get_exp` follows `frexp` for `0.5 <= |m| < 1.0`. Adjust the exponent by
+                    // one to scale the significand to `1.0 <= |m| < 2.0`.
+                    this.get_exp().map(|v| v - 1).unwrap_or_else(|| {
+                        if this.is_infinite() {
+                            i32::MAX
+                        } else {
+                            // Zero or NaN
+                            i32::MIN
+                        }
+                    })
+                }
+            }
         }
     };
 }

crates/libm-test/src/test_traits.rs

@@ -384,3 +384,15 @@ impl_tuples!(
     (f32, f32);
     (f64, f64);
 );
+
+#[cfg(f16_enabled)]
+impl_tuples!(
+    (f16, i32);
+    (f16, f16);
+);
+
+#[cfg(f128_enabled)]
+impl_tuples!(
+    (f128, i32);
+    (f128, f128);
+);

crates/libm-test/tests/compare_built_musl.rs

@@ -95,6 +95,10 @@ libm_macros::for_each_function! {
         fminf16,
         fmodf128,
         fmodf16,
+        frexpf128,
+        frexpf16,
+        ilogbf128,
+        ilogbf16,
         rintf128,
         rintf16,
         roundf128,

crates/util/src/main.rs

@@ -102,6 +102,10 @@ fn do_eval(basis: &str, op: &str, inputs: &[&str]) {
             | fminf16
             | fmodf128
             | fmodf16
+            | frexpf128
+            | frexpf16
+            | ilogbf128
+            | ilogbf16
             | rintf128
             | rintf16
             | roundf128

etc/function-definitions.json

@@ -466,16 +466,32 @@
     "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"
     },
+    "frexpf128": {
+        "sources": [
+            "src/math/frexpf128.rs",
+            "src/math/generic/frexp.rs"
+        ],
+        "type": "f128"
+    },
+    "frexpf16": {
+        "sources": [
+            "src/math/frexpf16.rs",
+            "src/math/generic/frexp.rs"
+        ],
+        "type": "f16"
+    },
     "hypot": {
         "sources": [
             "src/libm_helper.rs",
@@ -492,16 +508,32 @@
     "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"
     },
+    "ilogbf128": {
+        "sources": [
+            "src/math/generic/ilogb.rs",
+            "src/math/ilogbf128.rs"
+        ],
+        "type": "f128"
+    },
+    "ilogbf16": {
+        "sources": [
+            "src/math/generic/ilogb.rs",
+            "src/math/ilogbf16.rs"
+        ],
+        "type": "f16"
+    },
     "j0": {
         "sources": [
             "src/libm_helper.rs",

etc/function-list.txt

@@ -67,10 +67,14 @@ fmodf128
 fmodf16
 frexp
 frexpf
+frexpf128
+frexpf16
 hypot
 hypotf
 ilogb
 ilogbf
+ilogbf128
+ilogbf16
 j0
 j0f
 j1

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/frexpf128.rs

@@ -0,0 +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 frexpf128(x: f128) -> (f128, i32) {
+    super::generic::frexp(x)
+}

src/math/frexpf16.rs

@@ -0,0 +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 frexpf16(x: f16) -> (f16, i32) {
+    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_SAT {
+        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_SAT 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

@@ -6,6 +6,8 @@ mod floor;
 mod fmax;
 mod fmin;
 mod fmod;
+mod frexp;
+mod ilogb;
 mod rint;
 mod round;
 mod scalbn;
@@ -20,6 +22,8 @@ pub use floor::floor;
 pub use fmax::fmax;
 pub use fmin::fmin;
 pub use fmod::fmod;
+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/ilogbf128.rs

@@ -0,0 +1,5 @@
+/// Extract the binary exponent of `x`.
+#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
+pub fn ilogbf128(x: f128) -> i32 {
+    super::generic::ilogb(x)
+}

src/math/ilogbf16.rs

@@ -0,0 +1,5 @@
+/// Extract the binary exponent of `x`.
+#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
+pub fn ilogbf16(x: f16) -> i32 {
+    super::generic::ilogb(x)
+}

src/math/mod.rs

@@ -349,6 +349,8 @@ cfg_if! {
         mod fmaxf16;
         mod fminf16;
         mod fmodf16;
+        mod frexpf16;
+        mod ilogbf16;
         mod rintf16;
         mod roundf16;
         mod sqrtf16;
@@ -362,6 +364,8 @@ cfg_if! {
         pub use self::fmaxf16::fmaxf16;
         pub use self::fminf16::fminf16;
         pub use self::fmodf16::fmodf16;
+        pub use self::frexpf16::frexpf16;
+        pub use self::ilogbf16::ilogbf16;
         pub use self::rintf16::rintf16;
         pub use self::roundf16::roundf16;
         pub use self::sqrtf16::sqrtf16;
@@ -379,6 +383,8 @@ cfg_if! {
         mod fmaxf128;
         mod fminf128;
         mod fmodf128;
+        mod frexpf128;
+        mod ilogbf128;
         mod rintf128;
         mod roundf128;
         mod sqrtf128;
@@ -392,6 +398,8 @@ cfg_if! {
         pub use self::fmaxf128::fmaxf128;
         pub use self::fminf128::fminf128;
         pub use self::fmodf128::fmodf128;
+        pub use self::frexpf128::frexpf128;
+        pub use self::ilogbf128::ilogbf128;
         pub use self::rintf128::rintf128;
         pub use self::roundf128::roundf128;
         pub use self::sqrtf128::sqrtf128;