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.

Since we document[^2] that `a % b` is computed as...

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

...we can use `a % b` to compute `r_T`.

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

Author: traviscross

library/std/src/f128.rs

@@ -310,7 +310,10 @@ impl f128 {
     #[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 r < 0.0 {
+            return if rhs > 0.0 { r + rhs } else { r - rhs };
+        }
+        r
     }
 
     /// Raises a number to an integer power.

library/std/src/f16.rs

@@ -310,7 +310,10 @@ impl f16 {
     #[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 }
+        if r < 0.0 {
+            return if rhs > 0.0 { r + rhs } else { r - rhs };
+        }
+        r
     }
 
     /// Raises a number to an integer power.

library/std/src/f32.rs

@@ -286,7 +286,10 @@ impl f32 {
     #[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 }
+        if r < 0.0 {
+            return if rhs > 0.0 { r + rhs } else { r - rhs };
+        }
+        r
     }
 
     /// Raises a number to an integer power.

library/std/src/f64.rs

@@ -286,7 +286,10 @@ impl f64 {
     #[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 }
+        if r < 0.0 {
+            return if rhs > 0.0 { r + rhs } else { r - rhs };
+        }
+        r
     }
 
     /// Raises a number to an integer power.