Skip to content
This repository was archived by the owner on Apr 28, 2025. It is now read-only.

Commit 566e8e1

Browse files
tgross35tgross35
committed
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`.
1 parent 42a8dfb commit 566e8e1

File tree

6 files changed

+96
-162
lines changed

6 files changed

+96
-162
lines changed

‎etc/function-definitions.json

Lines changed: 4 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -437,13 +437,15 @@
437437
"fmod": {
438438
"sources": [
439439
"src/libm_helper.rs",
440-
"src/math/fmod.rs"
440+
"src/math/fmod.rs",
441+
"src/math/generic/fmod.rs"
441442
],
442443
"type": "f64"
443444
},
444445
"fmodf": {
445446
"sources": [
446-
"src/math/fmodf.rs"
447+
"src/math/fmodf.rs",
448+
"src/math/generic/fmod.rs"
447449
],
448450
"type": "f32"
449451
},

‎src/math/fmod.rs

Lines changed: 2 additions & 75 deletions
Original file line numberDiff line numberDiff line change
@@ -1,78 +1,5 @@
1+
/// Calculate the remainder of `x / y`, the precise result of `x - trunc(x / y) * y`.
12
#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
23
pub fn fmod(x: f64, y: f64) -> f64 {
3-
let mut uxi = x.to_bits();
4-
let mut uyi = y.to_bits();
5-
let mut ex = ((uxi >> 52) & 0x7ff) as i64;
6-
let mut ey = ((uyi >> 52) & 0x7ff) as i64;
7-
let sx = uxi >> 63;
8-
let mut i;
9-
10-
if uyi << 1 == 0 || y.is_nan() || ex == 0x7ff {
11-
return (x * y) / (x * y);
12-
}
13-
if uxi << 1 <= uyi << 1 {
14-
if uxi << 1 == uyi << 1 {
15-
return 0.0 * x;
16-
}
17-
return x;
18-
}
19-
20-
/* normalize x and y */
21-
if ex == 0 {
22-
i = uxi << 12;
23-
while i >> 63 == 0 {
24-
ex -= 1;
25-
i <<= 1;
26-
}
27-
uxi <<= -ex + 1;
28-
} else {
29-
uxi &= u64::MAX >> 12;
30-
uxi |= 1 << 52;
31-
}
32-
if ey == 0 {
33-
i = uyi << 12;
34-
while i >> 63 == 0 {
35-
ey -= 1;
36-
i <<= 1;
37-
}
38-
uyi <<= -ey + 1;
39-
} else {
40-
uyi &= u64::MAX >> 12;
41-
uyi |= 1 << 52;
42-
}
43-
44-
/* x mod y */
45-
while ex > ey {
46-
i = uxi.wrapping_sub(uyi);
47-
if i >> 63 == 0 {
48-
if i == 0 {
49-
return 0.0 * x;
50-
}
51-
uxi = i;
52-
}
53-
uxi <<= 1;
54-
ex -= 1;
55-
}
56-
i = uxi.wrapping_sub(uyi);
57-
if i >> 63 == 0 {
58-
if i == 0 {
59-
return 0.0 * x;
60-
}
61-
uxi = i;
62-
}
63-
while uxi >> 52 == 0 {
64-
uxi <<= 1;
65-
ex -= 1;
66-
}
67-
68-
/* scale result */
69-
if ex > 0 {
70-
uxi -= 1 << 52;
71-
uxi |= (ex as u64) << 52;
72-
} else {
73-
uxi >>= -ex + 1;
74-
}
75-
uxi |= sx << 63;
76-
77-
f64::from_bits(uxi)
4+
super::generic::fmod(x, y)
785
}

‎src/math/fmodf.rs

Lines changed: 2 additions & 85 deletions
Original file line numberDiff line numberDiff line change
@@ -1,88 +1,5 @@
1-
use core::f32;
2-
1+
/// Calculate the remainder of `x / y`, the precise result of `x - trunc(x / y) * y`.
32
#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
43
pub fn fmodf(x: f32, y: f32) -> f32 {
5-
let mut uxi = x.to_bits();
6-
let mut uyi = y.to_bits();
7-
let mut ex = ((uxi >> 23) & 0xff) as i32;
8-
let mut ey = ((uyi >> 23) & 0xff) as i32;
9-
let sx = uxi & 0x80000000;
10-
let mut i;
11-
12-
if uyi << 1 == 0 || y.is_nan() || ex == 0xff {
13-
return (x * y) / (x * y);
14-
}
15-
16-
if uxi << 1 <= uyi << 1 {
17-
if uxi << 1 == uyi << 1 {
18-
return 0.0 * x;
19-
}
20-
21-
return x;
22-
}
23-
24-
/* normalize x and y */
25-
if ex == 0 {
26-
i = uxi << 9;
27-
while i >> 31 == 0 {
28-
ex -= 1;
29-
i <<= 1;
30-
}
31-
32-
uxi <<= -ex + 1;
33-
} else {
34-
uxi &= u32::MAX >> 9;
35-
uxi |= 1 << 23;
36-
}
37-
38-
if ey == 0 {
39-
i = uyi << 9;
40-
while i >> 31 == 0 {
41-
ey -= 1;
42-
i <<= 1;
43-
}
44-
45-
uyi <<= -ey + 1;
46-
} else {
47-
uyi &= u32::MAX >> 9;
48-
uyi |= 1 << 23;
49-
}
50-
51-
/* x mod y */
52-
while ex > ey {
53-
i = uxi.wrapping_sub(uyi);
54-
if i >> 31 == 0 {
55-
if i == 0 {
56-
return 0.0 * x;
57-
}
58-
uxi = i;
59-
}
60-
uxi <<= 1;
61-
62-
ex -= 1;
63-
}
64-
65-
i = uxi.wrapping_sub(uyi);
66-
if i >> 31 == 0 {
67-
if i == 0 {
68-
return 0.0 * x;
69-
}
70-
uxi = i;
71-
}
72-
73-
while uxi >> 23 == 0 {
74-
uxi <<= 1;
75-
ex -= 1;
76-
}
77-
78-
/* scale result up */
79-
if ex > 0 {
80-
uxi -= 1 << 23;
81-
uxi |= (ex as u32) << 23;
82-
} else {
83-
uxi >>= -ex + 1;
84-
}
85-
uxi |= sx;
86-
87-
f32::from_bits(uxi)
4+
super::generic::fmod(x, y)
885
}

‎src/math/generic/fmod.rs

Lines changed: 84 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,84 @@
1+
/* SPDX-License-Identifier: MIT */
2+
/* origin: musl src/math/fmod.c. Ported to generic Rust algorithm in 2025, TG. */
3+
4+
use super::super::{CastFrom, Float, Int, MinInt};
5+
6+
#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
7+
pub fn fmod<F: Float>(x: F, y: F) -> F {
8+
let zero = F::Int::ZERO;
9+
let one = F::Int::ONE;
10+
let mut ix = x.to_bits();
11+
let mut iy = y.to_bits();
12+
let mut ex = x.exp().signed();
13+
let mut ey = y.exp().signed();
14+
let sx = ix & F::SIGN_MASK;
15+
16+
if iy << 1 == zero || y.is_nan() || ex == F::EXP_MAX as i32 {
17+
return (x * y) / (x * y);
18+
}
19+
20+
if ix << 1 <= iy << 1 {
21+
if ix << 1 == iy << 1 {
22+
return F::ZERO * x;
23+
}
24+
return x;
25+
}
26+
27+
/* normalize x and y */
28+
if ex == 0 {
29+
let i = ix << F::EXP_BITS;
30+
ex -= i.leading_zeros() as i32;
31+
ix <<= -ex + 1;
32+
} else {
33+
ix &= F::Int::MAX >> F::EXP_BITS;
34+
ix |= one << F::SIG_BITS;
35+
}
36+
37+
if ey == 0 {
38+
let i = iy << F::EXP_BITS;
39+
ey -= i.leading_zeros() as i32;
40+
iy <<= -ey + 1;
41+
} else {
42+
iy &= F::Int::MAX >> F::EXP_BITS;
43+
iy |= one << F::SIG_BITS;
44+
}
45+
46+
/* x mod y */
47+
while ex > ey {
48+
let i = ix.wrapping_sub(iy);
49+
if i >> (F::BITS - 1) == zero {
50+
if i == zero {
51+
return F::ZERO * x;
52+
}
53+
ix = i;
54+
}
55+
56+
ix <<= 1;
57+
ex -= 1;
58+
}
59+
60+
let i = ix.wrapping_sub(iy);
61+
if i >> (F::BITS - 1) == zero {
62+
if i == zero {
63+
return F::ZERO * x;
64+
}
65+
66+
ix = i;
67+
}
68+
69+
let shift = ix.leading_zeros().saturating_sub(F::EXP_BITS);
70+
ix <<= shift;
71+
ex -= shift as i32;
72+
73+
/* scale result */
74+
if ex > 0 {
75+
ix -= one << F::SIG_BITS;
76+
ix |= F::Int::cast_from(ex) << F::SIG_BITS;
77+
} else {
78+
ix >>= -ex + 1;
79+
}
80+
81+
ix |= sx;
82+
83+
F::from_bits(ix)
84+
}

‎src/math/generic/mod.rs

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -5,6 +5,7 @@ mod fdim;
55
mod floor;
66
mod fmax;
77
mod fmin;
8+
mod fmod;
89
mod rint;
910
mod round;
1011
mod scalbn;
@@ -18,6 +19,7 @@ pub use fdim::fdim;
1819
pub use floor::floor;
1920
pub use fmax::fmax;
2021
pub use fmin::fmin;
22+
pub use fmod::fmod;
2123
pub use rint::rint;
2224
pub use round::round;
2325
pub use scalbn::scalbn;

‎src/math/support/int_traits.rs

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -45,7 +45,9 @@ pub trait Int:
4545
+ ops::BitOrAssign
4646
+ ops::BitXorAssign
4747
+ ops::ShlAssign<i32>
48+
+ ops::ShlAssign<u32>
4849
+ ops::ShrAssign<u32>
50+
+ ops::ShrAssign<i32>
4951
+ ops::Add<Output = Self>
5052
+ ops::Sub<Output = Self>
5153
+ ops::Mul<Output = Self>

0 commit comments

Comments
 (0)