[BUG] Rounding error that break the ensurance of f32::div_euclid

Considering the following code:

let a = 1.10000001430511474609375_f32; // equivlent to let a = 1.1f32; but 1.1 is not the exact value.
assert_eq!(a, 1.1);
let b = 11f32;
let c = b.div_euclid(a); // expect 9, since 1.10000001430511474609375 * 10 > 11
println!("{c}"); // I'm sorry, but c is 10 here.

Although div/rem with float types is strongly discouraged, the doc said that,

This computes the integer n such that self = n * rhs + self.rem_euclid(rhs)

Here, b.div_euclid(a) * a + b.rem_euclid(a) is 12.099999, not 11. (rem_euclid yields the correct result.)

I wonder whether it should be fixed and is it worth creating an issue.

And what the possible solution it could be? Modify the document? Using a very ugly implementation so that b.div_euclid(a) + b.rem_euclid(a) equals to b? Or just mark the div_euclid/rem_euclid for float types as unsafe since they will meet some corner cases very often?


Related API:

Rust, failed

let a=11f64;
let b=1.1f64;
println!("{} {}", a/b, a%b)
// 10 1.0999999999999992
let a=11f64;
let b=1.1f64;
println!("{} {}", a.div_euclid(b), a.rem_euclid(b))
// 10 1.0999999999999992

python

>>> divmod(11,1.1)
(9.0, 1.0999999999999992)
8 Likes

Yep it's definitely a bug. The documentation of div_euclid states unambiguously:

Precision

The result of this operation is guaranteed to be the rounded infinite-precision result.

So a.div_euclid(b) should be 9 not 10.

Nitpick: the actual value of 1.1f32 is 1.10000002384185791015625 not 1.10000001430511474609375, but the answer is still supposed to be 9 not 10.

The implementation rounds to the nearest representable value before truncating, so no wonder it's not equivalent to only truncating as claimed by the docs.

6 Likes

div_euclid can be done correctly by using an integer division on appropriately scaled mantissas (requires dividing u128 by u64 for f64), with a special case when the divisor doesn't fit in u128 that requires some care but it can be done in a constant number of operations.

It could also be done in floating point by setting a non-default rounding mode. There is also a special case required for large answers (sometimes those have to be rounded up).

I think rem_euclid is correctly implemented.

Perhaps not.

In some rare cases, with roundings, even self < self.div_euclid(small_number) * small_number sometimes holds. According to document, it should be self.div_euclid(small_number) * small_number < self.

#![feature(f16, f128, float_next_up_down)]
fn main() {
    assert_eq!(1f64, 1f64.next_down().div_euclid(1e-16) * 1e-16);
    // assert are passed, since `self.div_euclid(small)` uses incorrect roundings.
}

The roundings cannot be avoided, and thus the result might not be what we want.

I strongly doubt that, whether those div/rem should be unsafe.

I think we could copy Python implemention: cpython/Objects/floatobject.c at 3.13 · python/cpython · GitHub.

PR's up: breaking change: Copy Python implementation for `float::div_euclid` by tesuji · Pull Request #133485 · rust-lang/rust · GitHub

May need an ACP though?

My interpretation is that self >= n * rhs and self == n * rhs + self.rem_euclid(rhs) should hold in real numbers, i.e. in infinite precision. The * and + in these formulas stand for exact multiplication and addition, not rounded multiplication and addition.

It's certainly possible (and not overly difficult) to implement it that way, as I described above.

This also doesn't look correctly rounded, the code doesn't even claim that in the comments, it's just doing some rough approximation of div.

2 Likes

This is worth discussed. Since we said "euclid", one may expect self >= n * rhs always holds no matter whether the rounding is performed.

No matter what properties it should hold, we could conclude that, self == n * rhs + self.rem_euclid(rhs) might not always hold due to various conditions. Thus it might be better to introduce an unsafe here, to alert users, that function may not do what exactly the users want to do.

The ACP mainly for safetyness. I strongly suspect div_euclid and rem_euclid on floating types are unsafe.

Plus, your conclusion, divmod(-11, -1.1) == 9 is incorrect, it should be 10. Unless you modify implementations for integers (signed only, unsigned unaffected), you cannot just modify divmod(-11, -1.1) == 9.

They aren't. Safety in Rust only refers to whether the function has some usage requirement (usually preconditions) that must be followed to avoid the unsound possibility for UB. A logically incorrect numeric result is not unsafe.

6 Likes

Maybe we should modify the defination?

In this case we may provide a divrem(_euclid) function:

fn divrem(self, rhs:Self) -> (Self, Self) {
    let rem = self % rhs;
    // for _euclid version, adding the following line:
    // if rem<0 {rem += rhs.abs()}
    (((self - rem) / rhs).round(), rem)
}

For this defination, we could ensure self.div_euclid(rhs) == (self - self.rem_euclid(rhs)) / rhs

I have no idea whether a .round() should be called for (self - rem) / rhs, it seems unnecessary but I cannot confirm.

For integers, self.div_euclid(rhs) == (self - self.rem_euclid(rhs)) / rhs is self.div_euclid(rhs) * rhs == self - self.rem_euclid(rhs), that is just what div_euclid is defined for integers.

According to Python, I think it is reasonable for the result of div_euclid(-11,-1.1) == 9:

In [2]: divmod(-11.0,-1.1)
Out[2]: (9.0, -1.0999999999999992)

In [3]: 9*(-1.1) + -1.0999999999999992
Out[3]: -11.0

But, if you talked about Rust rem_euclid then you're right, it shall be 10.0.

assert_eq!(f64::rem_euclid(-11.0, -1.1) + 10.0*(-1.1), -11.0);

I'm not sure which behavior is more reasonable.


Anyway I'm not an expert in this field. And I cannot prove that the Python implementation is correct for all floating-point numbers. I marked my PR above for discussion-only.

Maybe we could have both, divrem(self, div: Self) -> (Self, Self), which uses python's style, and divrem_euclid(self, div: Self) -> (Self, Self), which uses Rust's style.

In this case, self.div_euclid(div) should be self.divrem(div).0, and in a rare case we want a faster path, self.divrem(div).0 could be called.

As you can see, we might need additional ACP for those additional apis (divrem/divrem_euclid)


What's more, you could try divmod(-109,-11) in python compares to (-109).div_euclid(-11) in Rust.

Since it should be consist that (-109f32).div_euclid(-11f32) yields the same result as (-109i32).div_euclid(-11i32), you could figure that, the python way might not correct.

What's more, rust's defination provide a way for rounding: if you want to calculate a / b and rounding the result up, where both a and b are positive, you could just calculate (-a).div_euclid(-b). For python, no such path exists.

Using div_euclid to implement div_ceil when div_euclid is (logically) implemented in terms of div_trunc seems especially silly and roundabout, especially because it only works for one quadrant.

The python Euclid approach would be a.div_euclid(b) + if a.rem_euclid(b) > 0 { 1 } else { 0 }, if I'm not mistaken.

1 Like

It seems that we could rely on % operations:

fn check(a:(i64,i32), b:(i64,i32)) {
    assert_eq!(((a.0 as f64 * 2f64.powi(a.1)) / 2f64.powi(a.1)) as i64, a.0);
    assert_eq!(((b.0 as f64 * 2f64.powi(b.1)) / 2f64.powi(b.1)) as i64, b.0);
    println! {"{a0} * 2^{a1mb1} % {b0} == {}", ((a.0 as f64 * 2f64.powi(a.1)) % (b.0 as f64 * 2f64.powi(b.1))) * 2f64.powi(-b.1), a0 = a.0, b0=b.0, a1mb1=a.1-b.1}
}
fn main() {
    let a = (581941187211157i64, 970);
    let b = (751376156658013i64, -1000);
    check(a,b); // got 581941187211157 * 2^1970 % 751376156658013 == 210428852289837
}

the result is verified in Pari/GP. Thus it might be only div_euclid that contains BUG.

Floating point numbers are usually used to represent numbers where we neither know nor need to know what they are exactly. Insuring precise cut-offs with floating-points is hard or impossible. We may fix it in this case, but it will crop up elsewhere. We should just drop the guarantee in the docs. If sharply predictable cut-offs are needed, some other types of number should be used, e.g. integers, fixed-points, fractions.

2 Likes

I disagree because the guarantee is very feasible to uphold and I'm working on making it happen.

IEEE-754 actually recommends such a guarantee for all basic functions such as f32::sin, f32::ln, etc, but that would be a much more difficult thing to get and costly for performance.

5 Likes

Existing issue for this: f64::div_euclid and f64::rem_euclid yield inconsistent results · Issue #107904 · rust-lang/rust · GitHub

I don't see f32 mentioned there, so it would be useful to add examples of that as well.

5 Likes

Actually, all of f32, f64 and f128 failed to pass 11.0.div_euclid(1.1) * 1.1 + 11.0.rem_euclid(1.1) == 11.0 (approximately)

f16 passes test with 1.1 but failed with 1.2.

Most platforms support rounding modes other than nearest even. The most obvious fix would be to simply perform the division in a round down rounding mode.

The real problem is that Rust explicitly does not support non-default rounding modes, and applying the fix would trigger instant UB.

2 Likes

This does not quite work: the resulting integer needs to be rounded to the nearest representable number (not down).

Sure, if you want to do it and you are confident it is going to be fixed once and for all, that sounds great.