1 | /* SPDX-License-Identifier: MIT OR Apache-2.0 */ |
2 | use super::super::{CastFrom, Float, Int, MinInt}; |
3 | |
4 | #[inline ] |
5 | pub 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) |
50 | fn 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: ::Int = 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. |
61 | fn 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(default:x); |
66 | } |
67 | x |
68 | } |
69 | |