breaking change: Copy Python implementation for `float::div_euclid`

[tesuji]

This is a breaking change trying to fix https://internals.rust-lang.org/t/bug-rounding-error-that-break-the-ensurance-of-f32-div-euclid/21917.

In summary, Rust float::div_euclid and Python divmod disagrees with each others:		Python	Rust	This PR
div_eulid(11.0, 1.1)	9.0	10.0	9.0
div_eulid(-11.0, 1.1)	-10.0	-9.0	-10.0
div_eulid(11.0, -1.1)	-10.0	-10.0	-10.0
div_eulid(-11.0, -1.1)	9.0	11.0	9.0
div_eulid(0.5, 1.1)	0.0	0.0	0.0
div_eulid(-0.5, 1.1)	-1.0	-1.0	-1.0
div_eulid(0.5, -1.1)	-1.0	0.0	-1.0
div_eulid(-0.5, -1.1)	0.0	1.0	0.0

Personally I think the Python's behavior is more correct. Because given real numbers a and b, q is euclidean division of a and b, and r is the eulidean remainder of them. q*b + r should be equal to a. Take the first example, 11.0 and 1.1. Currectly with Rust 10.0*1.1 + 1.0999998 > 11.0 + 1 > 11.0.

FIXME: The Python link uses "fmod " to get more exact remander. But this PR only uses raw % operator. Is this a problem in practice?

Reference:

• https://github.com/python/cpython/blob/3.13/Objects/floatobject.c#L662.
• https://github.com/rust-lang/rust/issues/107904

@rustbot label T-libs-impl

[rustbot]

r? @cuviper

rustbot has assigned @cuviper. They will have a look at your PR within the next two weeks and either review your PR or reassign to another reviewer.

Use r? to explicitly pick a reviewer

[Noratrieb]

https://doc.rust-lang.org/std/primitive.f32.html#method.div_euclid the documentation clearly states that this method is precise to a rounded infinite-precision result, so its current implementation is just wrong and the new implementation is current. libs-api can have a look, but this seems very correct to me.

[traviscross]

To inline the problem a bit further:

fn main() {
    let (lhs, rhs) = (11.0f64, 1.1f64);
    // From the docs of `div_euclid`:
    //
    // This computes the integer `n`...
    let n = lhs.div_euclid(rhs);
    let r = lhs.rem_euclid(rhs);
    dbg!(lhs, rhs, n, r);
    // lhs = 11.0
    // rhs = 1.1
    // n = 10.0
    // r = 1.0999999999999992
    //
    // ...such that `self = n * rhs + self.rem_euclid(rhs)`.
    assert_eq!(lhs, n * rhs + r);
    //~^ PANIC
    // assertion `left == right` failed
    //  left: 11.0
    // right: 12.1
}

Playground link

That is, the behavior contradicts what our documentation states (and what good mathematical sense would suggest) about the invariant that div_euclid and rem_euclid hold in relation to each other.

[cuviper]

https://doc.rust-lang.org/std/primitive.f32.html#method.div_euclid the documentation clearly states that this method is precise to a rounded infinite-precision result, so its current implementation is just wrong and the new implementation is current. libs-api can have a look, but this seems very correct to me.

The documentation is not necessarily wrong because rounding from an infinite result can include rounding up. It's definitely bad that div/rem aren't in sync though.

I'm worried about the license implications of blatantly copying Python too. If that's indeed a problem, we may need someone else to do a clean-room fix.

[traviscross]

Agreed. Over in...

• https://github.com/rust-lang/rust/issues/107904

...it was suggested that:

...it should be possible to achieve a consistent behavior by just calculating one euclidean function based on the other one. E.g.

pub fn rem_euclid(self, rhs: f64) -> f64 {
    self - self.div_euclid(rhs) * rhs
}

We could do something like that.

[cuviper]

On rounding:

• 1.1f32 is 1.10000002384185791015625
• Wolfram Alpha says 11 divided by that is 9.9999997832558418782008370876130822005470839295152332671071558192...
• This round up to 10f32! The closest smaller value would be 9.99999904632568359375.

[traviscross]

This paper has an extensive analysis of this problem:

• The Euclidean Division Implemented with a Floating-Point Division and a Floor, Vincent Lefèvre, 2005

There is a companion paper:

• The Euclidean Division Implemented with a Floating-Point Multiplication and a Floor, Vincent Lefèvre, 2005

[tczajka]

Personally I think the Python's behavior is more correct.

Are you going for "more correct" or actually correct in all cases? I believe that the proposed (Python) implementation only works correctly in some cases: if the integer result can be exactly represented (such as in the example of 9). But it doesn't necessarily round correctly when the result can't be exactly represented (very large integer results). More care (and a proof) is necessary to make sure of correct rounding in all cases.

[tczajka]

• 1.1f32 is 1.10000002384185791015625
• Wolfram Alpha says 11 divided by that is 9.9999997832558418782008370876130822005470839295152332671071558192...
• This round up to 10f32! The closest smaller value would be 9.99999904632568359375.

The documentation guarantees "the rounded infinite-precision result.". Therefore the correct behavior is unambiguously 9f32:

1. Compute floor(11 / 1.10000002384185791015625) = floor(9.99999978.....) = 9 in infinite precision.
2. Round the infinite-precision result, 9, to the nearest representable value, which is 9f32. No rounding necessary in this case because 9 can be exactly represented.

[Neutron3529]

Maybe we need an unsafe, since the function will not always yield results we want to have.

Plus, since it is documented, 0<=rem<divisor.abs(), div_eulid(-11.0, -1.1) should be 10 rather than 9, since with infinite precision, { -1.1f64 as f128 * 10f64 as f128 } < { -11.0f64 as f128 } = true, and -1.1f64 * 9 > -11.0. Otherwise, calling div_euclid for integers will yields inconsistant results.

[tczajka]

Python divmod doesn't implement Euclidean division, it implements flooring division which is different for negative numbers. (-0.5).div_euclid(-1.1) should be 1.0, not 0.0.

[traviscross]

Having analyzed this a bit, and having looked at the Rust implementation in core (but notably, not the code in this PR or that for any other implementation), and having read through the paper I cited above, my feeling is increasingly that any fix other than...

pub fn rem_euclid(self, rhs: f64) -> f64 {
    self - self.div_euclid(rhs) * rhs
}

...should be accompanied by a (preferably machine-checked) proof of correctness.