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

Add fmodf16 #469

Merged
merged 2 commits into from
Jan 24, 2025
Merged

Add fmodf16 #469

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
566e8e1
Add a generic version of `fmod`
tgross35 Jan 24, 2025
856e314
Add `fmodf16` using the generic implementation
tgross35 Jan 24, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion crates/libm-macros/src/shared.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -47,7 +47,7 @@ const ALL_OPERATIONS_NESTED: &[(FloatTy, Signature, Option<Signature>, &[&str])]
FloatTy::F16,
Signature { args: &[Ty::F16, Ty::F16], returns: &[Ty::F16] },
None,
&["copysignf16", "fdimf16", "fmaxf16", "fminf16"],
&["copysignf16", "fdimf16", "fmaxf16", "fminf16", "fmodf16"],
),
(
// `(f32, f32) -> f32`
Expand Down
1 change: 1 addition & 0 deletions crates/libm-test/benches/icount.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -111,6 +111,7 @@ main!(
icount_bench_fmin_group,
icount_bench_fminf_group,
icount_bench_fmod_group,
icount_bench_fmodf16_group,
icount_bench_fmodf_group,
icount_bench_frexp_group,
icount_bench_frexpf_group,
Expand Down
1 change: 1 addition & 0 deletions crates/libm-test/benches/random.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -131,6 +131,7 @@ libm_macros::for_each_function! {
| fmaxf16
| fminf128
| fminf16
| fmodf16
| rintf128
| rintf16
| roundf128
Expand Down
17 changes: 17 additions & 0 deletions crates/libm-test/src/mpfloat.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -152,6 +152,7 @@ libm_macros::for_each_function! {
floorf16,
fmod,
fmodf,
fmodf16,
frexp,
frexpf,
ilogb,
Expand Down Expand Up @@ -525,6 +526,22 @@ impl MpOp for crate::op::lgammaf_r::Routine {
}
}

// No fmodf128 yet
impl MpOp for crate::op::fmodf16::Routine {
type MpTy = (MpFloat, MpFloat);

fn new_mp() -> Self::MpTy {
(new_mpfloat::<Self::FTy>(), new_mpfloat::<Self::FTy>())
}

fn run(this: &mut Self::MpTy, input: Self::RustArgs) -> Self::RustRet {
this.0.assign(input.0);
this.1.assign(input.1);
let ord = this.0.rem_assign_round(&this.1, Nearest);
prep_retval::<Self::RustRet>(&mut this.0, ord)
}
}

/* stub implementations so we don't need to special case them */

impl MpOp for crate::op::nextafter::Routine {
Expand Down
1 change: 1 addition & 0 deletions crates/libm-test/tests/compare_built_musl.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -93,6 +93,7 @@ libm_macros::for_each_function! {
fmaxf16,
fminf128,
fminf16,
fmodf16,
rintf128,
rintf16,
roundf128,
Expand Down
1 change: 1 addition & 0 deletions crates/util/src/main.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -100,6 +100,7 @@ fn do_eval(basis: &str, op: &str, inputs: &[&str]) {
| fmaxf16
| fminf128
| fminf16
| fmodf16
| rintf128
| rintf16
| roundf128
Expand Down
13 changes: 11 additions & 2 deletions etc/function-definitions.json
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -437,16 +437,25 @@
"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"
},
"fmodf16": {
"sources": [
"src/math/fmodf16.rs",
"src/math/generic/fmod.rs"
],
"type": "f16"
},
"frexp": {
"sources": [
"src/libm_helper.rs",
Expand Down
1 change: 1 addition & 0 deletions etc/function-list.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -63,6 +63,7 @@ fminf128
fminf16
fmod
fmodf
fmodf16
frexp
frexpf
hypot
Expand Down
77 changes: 2 additions & 75 deletions src/math/fmod.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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)
}
87 changes: 2 additions & 85 deletions src/math/fmodf.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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)
}
5 changes: 5 additions & 0 deletions src/math/fmodf16.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +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 fmodf16(x: f16, y: f16) -> f16 {
super::generic::fmod(x, y)
}
84 changes: 84 additions & 0 deletions src/math/generic/fmod.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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;
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think this was mistranslated from the previous loop. It fails to account for the sign bit. I believe it should be

Suggested change
let i = ix << F::EXP_BITS;
let i = ix << (F::EXP_BITS + 1);

I discovered this because I translated this implementation to Wasm by hand for use in the Scala.js-to-Wasm backend:
scala-js/scala-js#5147
My test cases with subnormals failed with the formula as is, but do pass with the added + 1.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks for finding this, I think you are correct. I'll make sure the failures get tested with the fix. I'll double check with all the cases from scala-js/scala-js#5147, but do you happen to know off the top of your head which ones failed?

Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Among the tests I wrote, the following fail without the + 1s.

Tests are test(expectedResult, x, y):

For Double (binary64):

test(2.0696e-320, 2.1, 3.123e-320)
test(1.772535e-318, 2.1, 2.253547e-318)

For Float (binary32):

test(8.085e-42f, 2.1f, 8.858e-42f)
test(6.1636e-40f, 2.1f, 6.39164e-40f)
test(4.77036e-40f, 5.5f, 6.39164e-40f)
test(-5.64734e-40f, -151.189f, 6.39164e-40f)

Basically they're cases where x is a normal but y is a subnormal. If they're both subnormals, the two errors "cancel out". And if x is a subnormal but y is a normal, then x < y and so the earlier branch exits early before getting here.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks for the detailed info, opened #535 to look at this.

ex -= i.leading_zeros() as i32;
ix <<= -ex + 1;
} else {
ix &= F::Int::MAX >> F::EXP_BITS;
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

There is something similar here. I believe the intent was also to have the + 1 :

Suggested change
ix &= F::Int::MAX >> F::EXP_BITS;
ix &= F::Int::MAX >> (F::EXP_BITS + 1);

However in this branch, it doesn't make any difference. The only bit for which it can change the result is anyway or'ed back to 1 on the very next line. But it seems more clearly correct with the + 1, as it then creates a mask that exactly matches the mantissa bits, rather than "the mantissa bits plus arbitrarily the low-order bit of the exponent".

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)
}
2 changes: 2 additions & 0 deletions src/math/generic/mod.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,7 @@ mod fdim;
mod floor;
mod fmax;
mod fmin;
mod fmod;
mod rint;
mod round;
mod scalbn;
Expand All @@ -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;
Expand Down
Loading