From 678480926f5ed77dcb26a64a0ae0b6ac8ed36f24 Mon Sep 17 00:00:00 2001
From: Travis Cross <tc@traviscross.com>
Date: Mon, 2 Dec 2024 13:22:48 +0000
Subject: [PATCH] Fix `{f16,f32,f64,f128}::rem_euclid`

The current implementation of `rem_euclid` for floating point numbers
violates the invariant, stated in the documentation, that:

```rust
a.rem_euclid(b) ~= a - b * a.div_euclid(b)
```

In a 2001 paper[^1], Daan Leijen (who notably later created the Koka
programming language) provides the correct formulation of this (and
of `div_euclid`) in "Algorithm E":

    q_E = q_T - I
    r_E = r_T + I * b
    where
        I = if r_T >= 0 then 0 else if b > 0 then 1 else -1
        q_T = trunc(a / b)
        r_T = a - b * q_T
        a is a dividend, a real number
        b is a divisor, a real number

In section 1.5 of the paper, he gives a proof of correctness.

To encode this in Rust, we might think to use `a % b` for
`r_T` (remainder of truncated division).  After all, we document[^2]
that `a % b` is computed as...

```rust
a - b * (a / b).trunc()
```

However, as it turns out, we do not currently compute `Rem` in this
way, as can be seen trivially with:

```rust
let (x, y) = (11f64, 1.1f64);
assert_eq!(x - (x / y).trunc() * y, x % y); //~ PANIC
```

For the moment, therefore, we've encoded `r_T` in the literal way.

As we know the maxim, from Knuth, to...

> Beware of bugs in the above code; I have only proved it correct, not
> tried it.

...we have additionally subjected our encoding of this formulation to
fuzzing.  It seems to hold up against the desired invariants.

This is of course a breaking change.  But the current implementation
is broken, and libs-api has signaled openness to fixing it.

[^1]: "Division and Modulus for Computer Scientists",
    Daan Leijen, 2001,
    <https://www.microsoft.com/en-us/research/wp-content/uploads/2016/02/divmodnote-letter.pdf>

[^2]: https://doc.rust-lang.org/std/ops/trait.Rem.html#impl-Rem-for-f64
---
 library/std/src/f128.rs | 13 +++++++++++--
 library/std/src/f16.rs  | 10 ++++++++--
 library/std/src/f32.rs  | 10 ++++++++--
 library/std/src/f64.rs  | 10 ++++++++--
 4 files changed, 35 insertions(+), 8 deletions(-)

diff --git a/library/std/src/f128.rs b/library/std/src/f128.rs
index e93e915159e40..9a07f294b05e6 100644
--- a/library/std/src/f128.rs
+++ b/library/std/src/f128.rs
@@ -309,8 +309,17 @@ impl f128 {
     #[unstable(feature = "f128", issue = "116909")]
     #[must_use = "method returns a new number and does not mutate the original value"]
     pub fn rem_euclid(self, rhs: f128) -> f128 {
-        let r = self % rhs;
-        if r < 0.0 { r + rhs.abs() } else { r }
+        if rhs.is_infinite() {
+            return self;
+        }
+        // FIXME(#133755): Though `self % rhs` is documented to be
+        // equivalent to this, it is not, and the distinction matters
+        // here.
+        let r = self - rhs * (self / rhs).trunc();
+        if r < 0.0 {
+            return if rhs > 0.0 { r + rhs } else { r - rhs };
+        }
+        r
     }
 
     /// Raises a number to an integer power.
diff --git a/library/std/src/f16.rs b/library/std/src/f16.rs
index 5b7fcaa28e064..48e0a229b34b8 100644
--- a/library/std/src/f16.rs
+++ b/library/std/src/f16.rs
@@ -309,8 +309,14 @@ impl f16 {
     #[unstable(feature = "f16", issue = "116909")]
     #[must_use = "method returns a new number and does not mutate the original value"]
     pub fn rem_euclid(self, rhs: f16) -> f16 {
-        let r = self % rhs;
-        if r < 0.0 { r + rhs.abs() } else { r }
+        // FIXME(#133755): Though `self % rhs` is documented to be
+        // equivalent to this, it is not, and the distinction matters
+        // here.
+        let r = self - rhs * (self / rhs).trunc();
+        if r < 0.0 {
+            return if rhs > 0.0 { r + rhs } else { r - rhs };
+        }
+        r
     }
 
     /// Raises a number to an integer power.
diff --git a/library/std/src/f32.rs b/library/std/src/f32.rs
index 7cb285bbff5f7..1654ae648b68d 100644
--- a/library/std/src/f32.rs
+++ b/library/std/src/f32.rs
@@ -285,8 +285,14 @@ impl f32 {
     #[inline]
     #[stable(feature = "euclidean_division", since = "1.38.0")]
     pub fn rem_euclid(self, rhs: f32) -> f32 {
-        let r = self % rhs;
-        if r < 0.0 { r + rhs.abs() } else { r }
+        // FIXME(#133755): Though `self % rhs` is documented to be
+        // equivalent to this, it is not, and the distinction matters
+        // here.
+        let r = self - rhs * (self / rhs).trunc();
+        if r < 0.0 {
+            return if rhs > 0.0 { r + rhs } else { r - rhs };
+        }
+        r
     }
 
     /// Raises a number to an integer power.
diff --git a/library/std/src/f64.rs b/library/std/src/f64.rs
index 47163c272de32..8487290c25a26 100644
--- a/library/std/src/f64.rs
+++ b/library/std/src/f64.rs
@@ -285,8 +285,14 @@ impl f64 {
     #[inline]
     #[stable(feature = "euclidean_division", since = "1.38.0")]
     pub fn rem_euclid(self, rhs: f64) -> f64 {
-        let r = self % rhs;
-        if r < 0.0 { r + rhs.abs() } else { r }
+        // FIXME(#133755): Though `self % rhs` is documented to be
+        // equivalent to this, it is not, and the distinction matters
+        // here.
+        let r = self - rhs * (self / rhs).trunc();
+        if r < 0.0 {
+            return if rhs > 0.0 { r + rhs } else { r - rhs };
+        }
+        r
     }
 
     /// Raises a number to an integer power.