f64::div_euclid and f64::rem_euclid yield inconsistent results

[Woyten]

I tried this code:

let numer = -0.7500000000000006f64;
let denom = 0.0833333333333334f64;
let div_euclid = -10.0;
let rem_euclid = 0.0;
assert_approx_eq!(numer.div_euclid(denom), div_euclid);
assert_approx_eq!(numer.rem_euclid(denom), rem_euclid);
assert_approx_eq!(div_euclid * denom + rem_euclid, numer - denom); // WRONG, should be numer

I expected to see this happen: div_euclid * denom + rem_euclid should yield numer. So either div_euclid should yield -9.0 or rem_euclid should yield 0.0833333333333334f64. Together they should satisfy div_euclid * denom + rem_euclid == numer.

Instead, this happened: div_euclid * denom + rem_euclid yields numer - denom. While both div_euclid and rem_euclid yield reasonable results in isolation, they are inconsistent which each other. The relation div_euclid * denom + rem_euclid == numer is not satisfied.

Meta

rustc --version --verbose:

rustc 1.68.0-nightly (4c83bd03a 2023-01-19)
binary: rustc
commit-hash: 4c83bd03a9d94af35c97a6b8b595d40e291af84a
commit-date: 2023-01-19
host: x86_64-unknown-linux-gnu
release: 1.68.0-nightly
LLVM version: 15.0.6

[edward-shen]

I can't seem to reproduce your issue:

fn main() {
    let numer = -0.7500000000000006f64;
    let denom = 0.0833333333333334f64;
    let div_euclid = -10.0;
    let rem_euclid = 0.0;
    assert!(numer.div_euclid(denom) - div_euclid < f64::EPSILON);
    assert!(numer.rem_euclid(denom) - rem_euclid < f64::EPSILON);
    assert!(div_euclid * denom + rem_euclid - numer < f64::EPSILON);
}

works fine for me: https://play.rust-lang.org/?version=nightly&mode=debug&edition=2021&gist=4b64e0e4da59fb03e821016c854f5955

[Woyten]

@edward-shen Thanks for your reply! I think you missed the .abs() term.

fn main() {
    let numer = -0.7500000000000006f64;
    let denom = 0.0833333333333334f64;
    let div_euclid = -10.0;
    let rem_euclid = 0.0;
    assert!((numer.div_euclid(denom) - div_euclid).abs() < f64::EPSILON);
    assert!((numer.rem_euclid(denom) - rem_euclid).abs() < f64::EPSILON);
    assert!((div_euclid * denom + rem_euclid - numer).abs() < f64::EPSILON);
}

Can you reproduce this example: https://play.rust-lang.org/?version=nightly&mode=debug&edition=2021&gist=585eb09aa751c09b8e79e8cec1ae15a7

[pzorin]

The problem seems to be that although numer == (-9.0) * denom and (numer / denom).trunc() == -9.0, in the real numbers we have |numer/denom|<9, because 2 trailing bits of 9.0*denom are truncated in f64, so the correction made in the implementation of div_euclid does not work in this case. On a related note, it would be good to clarify the promise div_euclid * denom + rem_euclid == numer made by div_euclid. It cannot hold for negative numer, because it is possible that rem_euclid.abs() > numer.abs(), so that the remainder will not have enough precision. (Update, since I am apparently not allowed to comment again): Consider as an example numer=-f64:MIN_POSITIVE and denom=1.0, then one might argue from the current description of f64:div_euclid that div_euclid should be -1, but then numer cannot be recovered precisely by addition of any remainder. This is one reason why the floating point remainder operation is usually defined to return a result with the same sign as the numerator. If a full precision remainder in a fixed fundamental domain modulo denom is desired, it can be chosen from the interval (-denom/2,denom/2]. This will still leave you with a rounding problem in the multiplication div_euclid * denom though, so in any case the promise should really be an inequality.

[Woyten]

@pzorin I think it is very important that the promise div_euclid * denom + rem_euclid == numer also holds for negative numbers. Otherwise, there would be no good reason to use {div,rem}_euclid in the first place. Note that there also is Div::div and Rem::rem but the reason to prefer {div,rem}_euclid is its strictly periodic behavior over positive and negative numbers. I wish this was the default case for Div::div and Rem::rem but that might be a topic for an RFC. :wink:

I haven't looked into the code yet but 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
}

[quaternic]

I've been looking into this for a bit now. Some observations:

For the division x % y there is an exact remainder r with |r| < |y| and the same sign as x, and Rust seems to compute precisely that by calling fmod, though its not clear from Rust's documentation if that's actually guaranteed.

The defining equation for x.{div/rem}_euclid(y) is x == n*y + r, where 0 <= r < |y|, and n is some integer

This can be exact, except for two edge cases:

• When |x| < |y|, and x < 0, the correct remainder would be x + |y|, which might not be exact, and in the worst case can round to |y|. In that case it might be most reasonable to return 0.0 instead of |y| to be able to make the promise that 0 <= r < |y|.
• The integer n can be large enough to not be exactly representable in a float. For the value returned by div_euclid, rounding error can be necessary whenever the integer quotient doesn't exactly fit in a float. I believe the most correct answer would be the exact integer that matches rem_euclid, rounded to the float type in some fashion.

Notably these two edge cases are disjoint. The remainder may need rounding only if |x| < |y| to begin with, so |n| <= 1.

If we do want to specify the result of div_euclid further, it might be reasonable to always round towards negative infinity. That would allow you to compute both an upper bound and a lower bound with

let r = x.rem_euclid(y);
let lb = x.div_euclid(y);
let ub = -x.div_euclid(-y);

This would satisfy x == n*y+r for some integer n with lb <= n <= ub

Additionally, it might be reasonable to never return -0.0 for the remainder, even though it compares >= 0.0, so the sign bit would always be positive.

I've written an implementation that does all of the above, but it's not currently as fast as it should be. The difficulty is that since fmod doesn't return information on the quotient it's problematic to recover one accurately afterwards, so I ended up reimplementing it. With the API having two separate functions, it's also not guaranteed you won't end up computing the same thing again for both rem_euclid and div_euclid.

[J-F-Liu]

I find a new pair of numbers yield inconsistent result:

fn main() {
    let s = 1.6470588235294124_f64;
    for a in [46.11764705882353_f64, 46.11764705882354_f64,46.11764705882355_f64] {
        println!("{:?} ÷ {:?} = {:?}⋯⋯{:?}", a, s, a.div_euclid(s), a.rem_euclid(s));
    }
}

The output is:

46.11764705882353 ÷ 1.6470588235294124 = 27.0⋯⋯1.6470588235293955
46.11764705882354 ÷ 1.6470588235294124 = 28.0⋯⋯1.6470588235294097
46.11764705882355 ÷ 1.6470588235294124 = 28.0⋯⋯4.440892098500626e-15

So 46.11764705882354.div_euclid(1.6470588235294124) should return 27.

[c00t]

This issue arises due to problems with the default rounding mode, which has been thoroughly analyzed in this paper: Vincent Lefèvre. The Euclidean Division Implemented with a Floating-Point Multiplication and a Floor. 2005. ⟨inria-00000159⟩

[tczajka]

The current documentation is clear about what should happen (and I think it is reasonable): exact values should be rounded according to usual floating point rounding rules.

rem_euclid does this correctly. div_euclid does not. So it's a bug in the implementation in div_euclid.

Consider as an example numer=-f64:MIN_POSITIVE and denom=1.0, then one might argue from the current description of f64:div_euclid that div_euclid should be -1, but then numer cannot be recovered precisely by addition of any remainder.

Yes, in this case div_euclid should be -1 and rem_euclid should be 1.0.

The equation numer == div_euclid * denom + rem_euclid holds in ideal real numbers, it doesn't have to hold after rounding (that would be impossible).

If a full precision remainder in a fixed fundamental domain modulo denom is desired, it can be chosen from the interval (-denom/2,denom/2]

That could be a useful operation, but that wouldn't be "Euclidean" division.

The difficulty is that since fmod doesn't return information on the quotient it's problematic to recover one accurately afterwards, so I ended up reimplementing it.

div_euclid can be faster than rem_euclid. Consider dividing 1e300 by 13. rem_euclid needs to do this precisely which requires essentially a bignum division loop (although it can be done in O(log exponent) rather than O(exponent) time), whereas div_euclid can be computed accurately in a single step without having to do the full division all the way down precisely.

[tczajka]

To implement this using integral division for f128 requires dividing u256 / u128 -> u128. Is it possible to utilize LLVM's i256 type for this?