Add a generic version of fmod

This can replace `fmod` and `fmodf`. As part of this change I was able
to replace some of the `while` loops with `leading_zeros`.

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"
     },

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

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

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<F: Float>(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)
+}

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;

library/compiler-builtins/libm/src/math/support/int_traits.rs

@@ -45,7 +45,9 @@ pub trait Int:
     + ops::BitOrAssign
     + ops::BitXorAssign
     + ops::ShlAssign<i32>
+    + ops::ShlAssign<u32>
     + ops::ShrAssign<u32>
+    + ops::ShrAssign<i32>
     + ops::Add<Output = Self>
     + ops::Sub<Output = Self>
     + ops::Mul<Output = Self>