From e9d7e6c37ce5729af98577c96b80f56a8f6b5f34 Mon Sep 17 00:00:00 2001 From: Trevor Gross Date: Fri, 3 Jan 2025 04:34:21 +0000 Subject: [PATCH] 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 | 14 +++ crates/libm-test/benches/random.rs | 4 + crates/libm-test/src/mpfloat.rs | 63 +++++------ crates/libm-test/src/precision.rs | 4 + crates/libm-test/tests/compare_built_musl.rs | 4 + crates/util/src/main.rs | 4 + etc/function-definitions.json | 28 +++++ etc/function-list.txt | 4 + src/math/generic/mod.rs | 2 + src/math/generic/scalbn.rs | 104 +++++++++++++++++++ src/math/ldexpf128.rs | 4 + src/math/ldexpf16.rs | 4 + src/math/mod.rs | 8 ++ src/math/scalbn.rs | 33 +----- src/math/scalbnf.rs | 29 +----- src/math/scalbnf128.rs | 4 + src/math/scalbnf16.rs | 4 + 17 files changed, 230 insertions(+), 87 deletions(-) create mode 100644 src/math/generic/scalbn.rs create mode 100644 src/math/ldexpf128.rs create mode 100644 src/math/ldexpf16.rs create mode 100644 src/math/scalbnf128.rs create mode 100644 src/math/scalbnf16.rs diff --git a/crates/libm-macros/src/shared.rs b/crates/libm-macros/src/shared.rs index 80bd3e907..66872bf3d 100644 --- a/crates/libm-macros/src/shared.rs +++ b/crates/libm-macros/src/shared.rs @@ -134,6 +134,13 @@ const ALL_OPERATIONS_NESTED: &[(FloatTy, Signature, Option, &[&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, &[&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, diff --git a/crates/libm-test/benches/random.rs b/crates/libm-test/benches/random.rs index 4d050e817..a31f4f2ad 100644 --- a/crates/libm-test/benches/random.rs +++ b/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 diff --git a/crates/libm-test/src/mpfloat.rs b/crates/libm-test/src/mpfloat.rs index a404f227b..127064a4e 100644 --- a/crates/libm-test/src/mpfloat.rs +++ b/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::[]::Routine { - type MpTy = ]::Routine as MpOp>::MpTy; - - fn new_mp() -> Self::MpTy { - ]::Routine as MpOp>::new_mp() - } - - fn run(this: &mut Self::MpTy, input: Self::RustArgs) -> Self::RustRet { - ]::Routine as MpOp>::run(this, input) - } - } - - impl MpOp for crate::op::[]::Routine { - type MpTy = MpFloat; - - fn new_mp() -> Self::MpTy { - new_mpfloat::() - } - - fn run(this: &mut Self::MpTy, input: Self::RustArgs) -> Self::RustRet { - this.assign(input.0); - *this <<= input.1; - prep_retval::(this, Ordering::Equal) - } - } - impl MpOp for crate::op::[]::Routine { type MpTy = (MpFloat, MpFloat); @@ -475,6 +452,34 @@ macro_rules! impl_op_for_ty_all { prep_retval::(&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::[]::Routine { + type MpTy = ]::Routine as MpOp>::MpTy; + + fn new_mp() -> Self::MpTy { + ]::Routine as MpOp>::new_mp() + } + + fn run(this: &mut Self::MpTy, input: Self::RustArgs) -> Self::RustRet { + ]::Routine as MpOp>::run(this, input) + } + } + + impl MpOp for crate::op::[]::Routine { + type MpTy = MpFloat; + + fn new_mp() -> Self::MpTy { + new_mpfloat::() + } + + fn run(this: &mut Self::MpTy, input: Self::RustArgs) -> Self::RustRet { + this.assign(input.0); + *this <<= input.1; + prep_retval::(this, Ordering::Equal) + } + } } }; } diff --git a/crates/libm-test/src/precision.rs b/crates/libm-test/src/precision.rs index 800425f12..ba7a33a04 100644 --- a/crates/libm-test/src/precision.rs +++ b/crates/libm-test/src/precision.rs @@ -585,8 +585,12 @@ fn int_float_common( 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( diff --git a/crates/libm-test/tests/compare_built_musl.rs b/crates/libm-test/tests/compare_built_musl.rs index f009816c9..26a446561 100644 --- a/crates/libm-test/tests/compare_built_musl.rs +++ b/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, diff --git a/crates/util/src/main.rs b/crates/util/src/main.rs index 889823d2e..3e91905a8 100644 --- a/crates/util/src/main.rs +++ b/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 diff --git a/etc/function-definitions.json b/etc/function-definitions.json index d3810b940..a6ef79b82 100644 --- a/etc/function-definitions.json +++ b/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", diff --git a/etc/function-list.txt b/etc/function-list.txt index 41bb4e06b..0aef383c8 100644 --- a/etc/function-list.txt +++ b/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 diff --git a/src/math/generic/mod.rs b/src/math/generic/mod.rs index d3df650e1..c7741cb46 100644 --- a/src/math/generic/mod.rs +++ b/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; diff --git a/src/math/generic/scalbn.rs b/src/math/generic/scalbn.rs new file mode 100644 index 000000000..39a749273 --- /dev/null +++ b/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(mut x: F, mut n: i32) -> F +where + u32: CastInto, + F::Int: CastFrom, + F::Int: CastFrom, +{ + 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::::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::from_bits(0x6ecb), -1336428830), f16::ZERO); +} + +#[test] +fn testme2() { + // assert_eq!(scalbn(-f64::INFINITY, -2147033648), f64::ZERO); +} diff --git a/src/math/ldexpf128.rs b/src/math/ldexpf128.rs new file mode 100644 index 000000000..b35277d15 --- /dev/null +++ b/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) +} diff --git a/src/math/ldexpf16.rs b/src/math/ldexpf16.rs new file mode 100644 index 000000000..8de6cffd6 --- /dev/null +++ b/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) +} diff --git a/src/math/mod.rs b/src/math/mod.rs index 53d06974c..7fe1606ac 100644 --- a/src/math/mod.rs +++ b/src/math/mod.rs @@ -346,7 +346,9 @@ cfg_if! { mod fabsf16; mod fdimf16; mod floorf16; + mod ldexpf16; mod rintf16; + mod scalbnf16; mod sqrtf16; mod truncf16; @@ -355,7 +357,9 @@ cfg_if! { pub use self::fabsf16::fabsf16; pub use self::fdimf16::fdimf16; pub use self::floorf16::floorf16; + pub use self::ldexpf16::ldexpf16; pub use self::rintf16::rintf16; + pub use self::scalbnf16::scalbnf16; pub use self::sqrtf16::sqrtf16; pub use self::truncf16::truncf16; } @@ -368,7 +372,9 @@ cfg_if! { mod fabsf128; mod fdimf128; mod floorf128; + mod ldexpf128; mod rintf128; + mod scalbnf128; mod sqrtf128; mod truncf128; @@ -377,7 +383,9 @@ cfg_if! { pub use self::fabsf128::fabsf128; pub use self::fdimf128::fdimf128; pub use self::floorf128::floorf128; + pub use self::ldexpf128::ldexpf128; pub use self::rintf128::rintf128; + pub use self::scalbnf128::scalbnf128; pub use self::sqrtf128::sqrtf128; pub use self::truncf128::truncf128; } diff --git a/src/math/scalbn.rs b/src/math/scalbn.rs index 00c455a10..f809dad51 100644 --- a/src/math/scalbn.rs +++ b/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) } diff --git a/src/math/scalbnf.rs b/src/math/scalbnf.rs index 73f4bb57a..57e7ba76f 100644 --- a/src/math/scalbnf.rs +++ b/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) } diff --git a/src/math/scalbnf128.rs b/src/math/scalbnf128.rs new file mode 100644 index 000000000..c1d2b4855 --- /dev/null +++ b/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) +} diff --git a/src/math/scalbnf16.rs b/src/math/scalbnf16.rs new file mode 100644 index 000000000..2209e1a17 --- /dev/null +++ b/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) +}