diff --git a/library/compiler-builtins/libm/etc/function-definitions.json b/library/compiler-builtins/libm/etc/function-definitions.json index b6653295c10c..866e9a439309 100644 --- a/library/compiler-builtins/libm/etc/function-definitions.json +++ b/library/compiler-builtins/libm/etc/function-definitions.json @@ -437,13 +437,15 @@ "fmod": { "sources": [ "src/libm_helper.rs", - "src/math/fmod.rs" + "src/math/fmod.rs", + "src/math/generic/fmod.rs" ], "type": "f64" }, "fmodf": { "sources": [ - "src/math/fmodf.rs" + "src/math/fmodf.rs", + "src/math/generic/fmod.rs" ], "type": "f32" }, diff --git a/library/compiler-builtins/libm/src/math/fmod.rs b/library/compiler-builtins/libm/src/math/fmod.rs index b68e6b0ea3c8..d9786b53d719 100644 --- a/library/compiler-builtins/libm/src/math/fmod.rs +++ b/library/compiler-builtins/libm/src/math/fmod.rs @@ -1,78 +1,5 @@ +/// Calculate the remainder of `x / y`, the precise result of `x - trunc(x / y) * y`. #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)] pub fn fmod(x: f64, y: f64) -> f64 { - let mut uxi = x.to_bits(); - let mut uyi = y.to_bits(); - let mut ex = ((uxi >> 52) & 0x7ff) as i64; - let mut ey = ((uyi >> 52) & 0x7ff) as i64; - let sx = uxi >> 63; - let mut i; - - if uyi << 1 == 0 || y.is_nan() || ex == 0x7ff { - return (x * y) / (x * y); - } - if uxi << 1 <= uyi << 1 { - if uxi << 1 == uyi << 1 { - return 0.0 * x; - } - return x; - } - - /* normalize x and y */ - if ex == 0 { - i = uxi << 12; - while i >> 63 == 0 { - ex -= 1; - i <<= 1; - } - uxi <<= -ex + 1; - } else { - uxi &= u64::MAX >> 12; - uxi |= 1 << 52; - } - if ey == 0 { - i = uyi << 12; - while i >> 63 == 0 { - ey -= 1; - i <<= 1; - } - uyi <<= -ey + 1; - } else { - uyi &= u64::MAX >> 12; - uyi |= 1 << 52; - } - - /* x mod y */ - while ex > ey { - i = uxi.wrapping_sub(uyi); - if i >> 63 == 0 { - if i == 0 { - return 0.0 * x; - } - uxi = i; - } - uxi <<= 1; - ex -= 1; - } - i = uxi.wrapping_sub(uyi); - if i >> 63 == 0 { - if i == 0 { - return 0.0 * x; - } - uxi = i; - } - while uxi >> 52 == 0 { - uxi <<= 1; - ex -= 1; - } - - /* scale result */ - if ex > 0 { - uxi -= 1 << 52; - uxi |= (ex as u64) << 52; - } else { - uxi >>= -ex + 1; - } - uxi |= sx << 63; - - f64::from_bits(uxi) + super::generic::fmod(x, y) } diff --git a/library/compiler-builtins/libm/src/math/fmodf.rs b/library/compiler-builtins/libm/src/math/fmodf.rs index 4de18195770d..4e95696e20d6 100644 --- a/library/compiler-builtins/libm/src/math/fmodf.rs +++ b/library/compiler-builtins/libm/src/math/fmodf.rs @@ -1,88 +1,5 @@ -use core::f32; - +/// Calculate the remainder of `x / y`, the precise result of `x - trunc(x / y) * y`. #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)] pub fn fmodf(x: f32, y: f32) -> f32 { - let mut uxi = x.to_bits(); - let mut uyi = y.to_bits(); - let mut ex = ((uxi >> 23) & 0xff) as i32; - let mut ey = ((uyi >> 23) & 0xff) as i32; - let sx = uxi & 0x80000000; - let mut i; - - if uyi << 1 == 0 || y.is_nan() || ex == 0xff { - return (x * y) / (x * y); - } - - if uxi << 1 <= uyi << 1 { - if uxi << 1 == uyi << 1 { - return 0.0 * x; - } - - return x; - } - - /* normalize x and y */ - if ex == 0 { - i = uxi << 9; - while i >> 31 == 0 { - ex -= 1; - i <<= 1; - } - - uxi <<= -ex + 1; - } else { - uxi &= u32::MAX >> 9; - uxi |= 1 << 23; - } - - if ey == 0 { - i = uyi << 9; - while i >> 31 == 0 { - ey -= 1; - i <<= 1; - } - - uyi <<= -ey + 1; - } else { - uyi &= u32::MAX >> 9; - uyi |= 1 << 23; - } - - /* x mod y */ - while ex > ey { - i = uxi.wrapping_sub(uyi); - if i >> 31 == 0 { - if i == 0 { - return 0.0 * x; - } - uxi = i; - } - uxi <<= 1; - - ex -= 1; - } - - i = uxi.wrapping_sub(uyi); - if i >> 31 == 0 { - if i == 0 { - return 0.0 * x; - } - uxi = i; - } - - while uxi >> 23 == 0 { - uxi <<= 1; - ex -= 1; - } - - /* scale result up */ - if ex > 0 { - uxi -= 1 << 23; - uxi |= (ex as u32) << 23; - } else { - uxi >>= -ex + 1; - } - uxi |= sx; - - f32::from_bits(uxi) + super::generic::fmod(x, y) } diff --git a/library/compiler-builtins/libm/src/math/generic/fmod.rs b/library/compiler-builtins/libm/src/math/generic/fmod.rs new file mode 100644 index 000000000000..93da6c51e7e1 --- /dev/null +++ b/library/compiler-builtins/libm/src/math/generic/fmod.rs @@ -0,0 +1,84 @@ +/* SPDX-License-Identifier: MIT */ +/* origin: musl src/math/fmod.c. Ported to generic Rust algorithm in 2025, TG. */ + +use super::super::{CastFrom, Float, Int, MinInt}; + +#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)] +pub fn fmod(x: F, y: F) -> F { + let zero = F::Int::ZERO; + let one = F::Int::ONE; + let mut ix = x.to_bits(); + let mut iy = y.to_bits(); + let mut ex = x.exp().signed(); + let mut ey = y.exp().signed(); + let sx = ix & F::SIGN_MASK; + + if iy << 1 == zero || y.is_nan() || ex == F::EXP_MAX as i32 { + return (x * y) / (x * y); + } + + if ix << 1 <= iy << 1 { + if ix << 1 == iy << 1 { + return F::ZERO * x; + } + return x; + } + + /* normalize x and y */ + if ex == 0 { + let i = ix << F::EXP_BITS; + ex -= i.leading_zeros() as i32; + ix <<= -ex + 1; + } else { + ix &= F::Int::MAX >> F::EXP_BITS; + ix |= one << F::SIG_BITS; + } + + if ey == 0 { + let i = iy << F::EXP_BITS; + ey -= i.leading_zeros() as i32; + iy <<= -ey + 1; + } else { + iy &= F::Int::MAX >> F::EXP_BITS; + iy |= one << F::SIG_BITS; + } + + /* x mod y */ + while ex > ey { + let i = ix.wrapping_sub(iy); + if i >> (F::BITS - 1) == zero { + if i == zero { + return F::ZERO * x; + } + ix = i; + } + + ix <<= 1; + ex -= 1; + } + + let i = ix.wrapping_sub(iy); + if i >> (F::BITS - 1) == zero { + if i == zero { + return F::ZERO * x; + } + + ix = i; + } + + let shift = ix.leading_zeros().saturating_sub(F::EXP_BITS); + ix <<= shift; + ex -= shift as i32; + + /* scale result */ + if ex > 0 { + ix -= one << F::SIG_BITS; + ix |= F::Int::cast_from(ex) << F::SIG_BITS; + } else { + ix >>= -ex + 1; + } + + ix |= sx; + + F::from_bits(ix) +} diff --git a/library/compiler-builtins/libm/src/math/generic/mod.rs b/library/compiler-builtins/libm/src/math/generic/mod.rs index 819781a219aa..68686b0b2550 100644 --- a/library/compiler-builtins/libm/src/math/generic/mod.rs +++ b/library/compiler-builtins/libm/src/math/generic/mod.rs @@ -5,6 +5,7 @@ mod fdim; mod floor; mod fmax; mod fmin; +mod fmod; mod rint; mod round; mod scalbn; @@ -18,6 +19,7 @@ pub use fdim::fdim; pub use floor::floor; pub use fmax::fmax; pub use fmin::fmin; +pub use fmod::fmod; pub use rint::rint; pub use round::round; pub use scalbn::scalbn; diff --git a/library/compiler-builtins/libm/src/math/support/int_traits.rs b/library/compiler-builtins/libm/src/math/support/int_traits.rs index cf19762e8169..b403c658cb69 100644 --- a/library/compiler-builtins/libm/src/math/support/int_traits.rs +++ b/library/compiler-builtins/libm/src/math/support/int_traits.rs @@ -45,7 +45,9 @@ pub trait Int: + ops::BitOrAssign + ops::BitXorAssign + ops::ShlAssign + + ops::ShlAssign + ops::ShrAssign + + ops::ShrAssign + ops::Add + ops::Sub + ops::Mul