Floating point modulo producing incorrect compuations

[Gip-Gip]

I tried this code:

fn ieee_mod(a: f32, b: f32) -> f32 {
  a - (a / b).round() * b
}

fn libc_mod(a: f32, b: f32) -> f32 {
  a % b
}


fn main() {
    let x = 1.5 * 3.14;
    let y = 0.5 * 3.14;
    
    assert_eq!(ieee_mod(x, y), libc_mod(x, y));
}

I expected the assert to pass.

Instead, this happened:

thread 'main' panicked at 'assertion failed: `(left == right)`
  left: `0.0`,
 right: `1.5699999`', bug.rs:14:5

Meta

This is most likely due to division being rounded, but mod not. This caused a bug in a sin table function I was writing.

This is standard libc behavior but it is still not correct and can cause bugs where none are expected.

rustc --version --verbose:

rustc 1.69.0 (84c898d65 2023-04-16)
binary: rustc
commit-hash: 84c898d65adf2f39a5a98507f1fe0ce10a2b8dbc
commit-date: 2023-04-16
host: x86_64-unknown-linux-gnu
release: 1.69.0
LLVM version: 15.0.7

Backtrace

stack backtrace:
   0: rust_begin_unwind
             at /rustc/84c898d65adf2f39a5a98507f1fe0ce10a2b8dbc/library/std/src/panicking.rs:579:5
   1: core::panicking::panic_fmt
             at /rustc/84c898d65adf2f39a5a98507f1fe0ce10a2b8dbc/library/core/src/panicking.rs:64:14
   2: core::panicking::assert_failed_inner
   3: core::panicking::assert_failed
   4: bug::main
   5: core::ops::function::FnOnce::call_once
note: Some details are omitted, run with `RUST_BACKTRACE=full` for a verbose backtrace.

[asquared31415]

Small note: the manual implementation should be trunc not round, but it's still producing the same result.

The mod impl says that these operations should be equivalent.

[konnorandrews]

After investigation, the answer given by Rust is congruent to 0.0 within the mod. However, the main concern to me is this behavior breaks the identity given for rem_euclid in the std docs quite badly.

The rem_euclid identity given is:

In particular, the return value r satisfies 0.0 <= r < rhs.abs() in most cases. However, due to a floating point round-off error it can result in r == rhs.abs(), violating the mathematical definition, if self is much smaller than rhs.abs() in magnitude and self < 0.0. This result is not an element of the function’s codomain, but it is the closest floating point number in the real numbers and thus fulfills the property self == self.div_euclid(rhs) * rhs + self.rem_euclid(rhs) approximately.

https://doc.rust-lang.org/std/primitive.f32.html#method.rem_euclid

Example playground: https://play.rust-lang.org/?version=stable&mode=debug&edition=2021&gist=838ccc4e48cd3ef6fa7c475ba1fae693

Output:

different answers in an absolute sense
[src/main.rs:6] test(x, y, ieee_remainder) = 0.0
[src/main.rs:7] test(x, y, mod_trunc) = 0.0
[src/main.rs:8] test(x, y, rust_std_mod) = 1.5707963
[src/main.rs:9] test(x, y, libm_fmodf) = 1.5707963
[src/main.rs:10] test(x, y, libm_remainderf) = -1.1920929e-7

the errors are within reason though not exact to the ieee definition
[src/main.rs:13] error(test(x, y, libm_remainderf), 0.0, y) = 1.1920929e-7
[src/main.rs:14] error(test(x, y, libm_fmodf), 0.0, y) = 1.1920929e-7

the identity given for rem_euclid fails for this case
[src/main.rs:17] test(x, y, rem_euclid_identity) = 6.283185
[src/main.rs:17] x = 4.712389

[konnorandrews]

Note also that the libc remainderf gives a much better result here, but Rust + LLVM doesn't use this implementation.

[konnorandrews]

Results on Windows 10 for those curious if there is a difference.

[src\main.rs:6] test(x, y, ieee_remainder) = 0.0
[src\main.rs:7] test(x, y, mod_trunc) = 0.0
[src\main.rs:8] test(x, y, rust_std_mod) = 1.5707963        
[src\main.rs:9] test(x, y, libm_fmodf) = 1.5707963
[src\main.rs:10] test(x, y, libm_remainderf) = -1.1920929e-7

the errors are within reason though not exact to the ieee definition
[src\main.rs:13] error(test(x, y, libm_remainderf), 0.0, y) = 1.1920929e-7
[src\main.rs:14] error(test(x, y, libm_fmodf), 0.0, y) = 1.1920929e-7

the identity given for rem_euclid fails for this case
[src\main.rs:17] test(x, y, rem_euclid_identity) = 6.283185
[src\main.rs:17] x = 4.712389

[quaternic]

Rem % for floats is equivalent to libc's fmod. Rust's documentation for it is indeed misleading. "Is computed as x - (x / y).trunc() * y" reads like that's an equivalent expression, when actually x % y is that expression as if computed without intermediate rounding.

Given that definition, I do believe % is giving the correct result. With x = 1.5 * 3.14_f32 and y = 0.5 * 3.14_f32:

4.710000038146972656250000000000 <- x
1.570000052452087402343750000000 <- y
1.569999933242797851562500000000 <- x - 2.0 * y
1.569999933242797851562500000000 <- x % y

There's an existing issue about rem_euclid being inconsistent with div_euclid: #107904

[cyqsimon]

I may have a related problem.

While working on https://github.com/imsnif/bandwhich/ I noticed that some computed test snapshots are different on Windows. Upon further investigation, it seems like it's due to the dreaded floating point math. Here's a sample:

println!("{:?}", 2.5f64.powi(39));

On Linux and MacOS this computes to 3308722450212111.0, whereas on Windows this computes to 3308722450212110.5. Similar problem exists for 2.5f64.powi(33) and 2.5f64.powi(34) and many others.

Tested on 1.70, 1.74, 1.77, and nightly; identical behivour between the versions.