libm/math/generic/
fmod.rs

1/* SPDX-License-Identifier: MIT OR Apache-2.0 */
2use crate::support::{CastFrom, Float, Int, MinInt};
3
4#[inline]
5pub fn fmod<F: Float>(x: F, y: F) -> F {
6    let _1 = F::Int::ONE;
7    let sx = x.to_bits() & F::SIGN_MASK;
8    let ux = x.to_bits() & !F::SIGN_MASK;
9    let uy = y.to_bits() & !F::SIGN_MASK;
10
11    // Cases that return NaN:
12    //   NaN % _
13    //   Inf % _
14    //     _ % NaN
15    //     _ % 0
16    let x_nan_or_inf = ux & F::EXP_MASK == F::EXP_MASK;
17    let y_nan_or_zero = uy.wrapping_sub(_1) & F::EXP_MASK == F::EXP_MASK;
18    if x_nan_or_inf | y_nan_or_zero {
19        return (x * y) / (x * y);
20    }
21
22    if ux < uy {
23        // |x| < |y|
24        return x;
25    }
26
27    let (num, ex) = into_sig_exp::<F>(ux);
28    let (div, ey) = into_sig_exp::<F>(uy);
29
30    // To compute `(num << ex) % (div << ey)`, first
31    // evaluate `rem = (num << (ex - ey)) % div` ...
32    let rem = reduction(num, ex - ey, div);
33    // ... so the result will be `rem << ey`
34
35    if rem.is_zero() {
36        // Return zero with the sign of `x`
37        return F::from_bits(sx);
38    };
39
40    // We would shift `rem` up by `ey`, but have to stop at `F::SIG_BITS`
41    let shift = ey.min(F::SIG_BITS - rem.ilog2());
42    // Anything past that is added to the exponent field
43    let bits = (rem << shift) + (F::Int::cast_from(ey - shift) << F::SIG_BITS);
44    F::from_bits(sx + bits)
45}
46
47/// Given the bits of a finite float, return a tuple of
48///  - the mantissa with the implicit bit (0 if subnormal, 1 otherwise)
49///  - the additional exponent past 1, (0 for subnormal, 0 or more otherwise)
50fn into_sig_exp<F: Float>(mut bits: F::Int) -> (F::Int, u32) {
51    bits &= !F::SIGN_MASK;
52    // Subtract 1 from the exponent, clamping at 0
53    let sat = bits.checked_sub(F::IMPLICIT_BIT).unwrap_or(F::Int::ZERO);
54    (
55        bits - (sat & F::EXP_MASK),
56        u32::cast_from(sat >> F::SIG_BITS),
57    )
58}
59
60/// Compute the remainder `(x * 2.pow(e)) % y` without overflow.
61fn reduction<I: Int>(mut x: I, e: u32, y: I) -> I {
62    x %= y;
63    for _ in 0..e {
64        x <<= 1;
65        x = x.checked_sub(y).unwrap_or(x);
66    }
67    x
68}
libm/math/generic/fmod.rs

libm/math/generic/
fmod.rs
