`<{f16,f32,f64,f128} as Rem>::rem` documented definition is misleading w.r.t. intermediate rounding

[traviscross]

Our documentation for impl Rem for {f16, f32, f64, f128} says:

The remainder has the same sign as the dividend and is computed as: x - (x / y).trunc() * y.

But that's not true. E.g.:

fn main() {
    let (x, y) = (11f64, 1.1f64);
    assert_eq!(x - (x / y).trunc() * y, x % y);
    //~^ PANIC
    // assertion `left == right` failed
    // left: 0.0
    // right: 1.0999999999999992
}

This mismatch creates a hazard when trying to correctly encode algorithms that rely on this semantic.

This tripped me up, e.g., when authoring:

• Fix {f16,f32,f64,f128}::div_euclid #133755

cc #133485 #107904

cc @cuviper @Noratrieb @tczajka @Neutron3529 @BartMassey

[traviscross]

Notably, and perhaps uncoincidentally, C's fmod has the same documented behavior...

The floating-point remainder of the division operation x / y calculated by this function is exactly the value x - n * y, where n is x / y with its fractional part truncated.

...and the same mismatch:

unsafe extern "C" {
    safe fn fmod(x: f64, y: f64) -> f64;
}

fn main() {
    let (x, y) = (11f64, 1.1f64);
    assert_eq!(x - (x / y).trunc() * y, fmod(x, y));
    //~^ PANIC
    // assertion `left == right` failed
    // left: 0.0
    // right: 1.0999999999999992
}

[quaternic]

Indeed, Rust's floating point % is actually equivalent to fmod. The documentation is only correct if you interpret the expression x - (x / y).trunc() * y as an exact formula. (The wording of "computed as" certainly doesn't help one towards that reading.)

[traviscross]

Yes, and the C docs state it even more clearly. They say, "it is exactly the value...".

[quaternic]

Yes, it is mathematically exact:

let x = 11_f64;  // == 11.0
let y = 1.1_f64; // == 1.100000000000000088817841970012523233890533447265625
                 // == 2476979795053773 / 2^51

// x / y == 11 / (2476979795053773 / 2^51)
//       == (11 * 2^51) / 2476979795053773
//       == 24769797950537728 / 2476979795053773
//       == 9 + 2476979795053771 / 2476979795053773 ( == 10 - 2 / 2476979795053773 )
// so the truncated quotient is 9
let iquot = 9;

// x - iquot * y == 11 - (9 * 2476979795053773 / 2^51)
//               == (11 * 2^51 - 9 * 2476979795053773) / 2^51
//               == 2476979795053771 / 2^51
//               == 1.099999999999999200639422269887290894985198974609375

assert_eq!(x % y, 1.099999999999999200639422269887290894985198974609375);

(sorry, the numbers get rather large with f64)

IMO, this is a bug in the documentation. Float % being equivalent to C's fmod is stated explicitly in https://rust-lang.github.io/rfcs/3514-float-semantics.html

[traviscross]

In my reading of the C23 (pre-)standard, section 7.12.10, either behavior would be acceptable:

The fmod functions return the value x − ny, for some integer n such that, if y is nonzero, the result has the same sign as x and magnitude less than the magnitude of y. If y is zero, whether a domain error occurs or the fmod functions return zero is implementation-defined.

That is, the standard is satisfied whether we choose 9 or 10 for n here.

[quaternic]

edit: I misinterpreted the C standard. See the next message.

Oh, I see. It's specified more loosely than the related remainder which is mandated by IEEE754.

So, fmod(11.0, 1.1) is allowed to return either of

• 11.0 - 9.0 * 1.1 == 1.0999999999999996
• 11.0 - 10.0 * 1.1 == 0.0

(in addition to the exact result of 1.0999999999999992_f64)

It's starting to look like a shaky foundation to be basing our operations on.

[tczajka]

That is, the standard is satisfied whether we choose 9 or 10 for n here.

So, fmod(11.0, 1.1) is allowed to return either of

• 11.0 - 9.0 * 1.1 == 1.0999999999999996
• 11.0 - 10.0 * 1.1 == 0.0

(in addition to the exact result of 1.0999999999999992_f64)

No, this is incorrect. The exact value is the only value that is allowed by the C standard (for an implementation that supports IEEE-754 with subnormals), as specified by F.10.7.1:

When subnormal results are supported, the returned value is exact and is independent of the current
rounding direction mode.

The C expression 11.0 - 10.0 * 1.1 is irrelevant. Section 7.12.10 says:

x − ny

That's a mathematical formula, not a C expression (as is clear by the lack of a * operator between n and y). For n=10 the value of the expression is negative which is not allowed.

[traviscross]

Agreeing and extending, IEEE 754-2008 says in 5.3.1:

When y ≠ 0, the remainder r = remainder(x, y) is defined for finite x and y regardless of the rounding-direction attribute by the mathematical relation r = x − y × n , where n is the integer nearest the exact number x/y ; whenever | n − x/y | = ½ , then n is even. Thus, the remainder is always exact. If r = 0, its sign shall be that of x. remainder(x, ∞) is x for finite x.

Then C23 F.10.7.1 says:

The double version of fmod behaves as though implemented by:

#include <math.h>
#include <fenv.h>
#pragma STDC FENV_ACCESS ON
double fmod(double x, double y)
{
    double result;
    result = remainder(fabs(x), (y = fabs(y)));
    if (signbit(result)) result += y;
    return copysign(result, x);
}

So, for 11 and 1.1, remainder will choose n = 10, resulting in us calculating:

remainder(11., 1.1) = exact!(11. - 1.1 * 10.) = -10f64.mul_add(1.1, -11.) = -0.0000000000000008881784197001252;
fmod(11., 1.1) = remainder(11., 1.1) + 1.1 = 1.0999999999999992;

[11happy]

Hello @traviscross , I am interested to work on this issue to fix div_euclid could you please assign this to me.
Thank you