Fix {f16,f32,f64,f128}::div_euclid
#133755
Conversation
rustbot
added
S-waiting-on-review
This comment has been minimized.
This comment has been minimized.
traviscross
force-pushed
the
TC/fix-rem_euclid
branch
2 times, most recently
from
1aea373 to
21349e5
Compare
This comment has been minimized.
This comment has been minimized.
traviscross
force-pushed
the
TC/fix-rem_euclid
branch
from
21349e5 to
6784809
Compare
| if r < 0.0 { | ||
| return if rhs > 0.0 { r + rhs } else { r - rhs }; | ||
| } | ||
| r |
There was a problem hiding this comment.
nit: is there a reason this has a return instead of
| if r < 0.0 { | |
| return if rhs > 0.0 { r + rhs } else { r - rhs }; | |
| } | |
| r | |
| if r < 0.0 { | |
| if rhs > 0.0 { r + rhs } else { r - rhs } | |
| } else { | |
| r | |
| } |
There was a problem hiding this comment.
For consistency (both stylistic and as it makes the mathematically analogy more clear), I followed the style of the existing div_euclid implementation:
pub fn div_euclid(self, rhs: f64) -> f64 {
let q = (self / rhs).trunc();
if self % rhs < 0.0 {
return if rhs > 0.0 { q - 1.0 } else { q + 1.0 };
}
q
}
This comment has been minimized.
This comment has been minimized.
traviscross
force-pushed
the
TC/fix-rem_euclid
branch
from
6784809 to
1921207
Compare
This comment has been minimized.
This comment has been minimized.
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
```
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/1.83.0/std/ops/trait.Rem.html#impl-Rem-for-f64
traviscross
force-pushed
the
TC/fix-rem_euclid
branch
from
1921207 to
0b3b949
Compare
|
In talking with @ehuss, he pointed out RFC 2169 where In looking at that, interestingly, the definition given for fn mod_euc(self, rhs: f64) -> f64 {
let r = self % rhs;
if r < 0.0 {
return if rhs > 0.0 { r + rhs } else { r - rhs }
}
r
}...matches exactly the original implementation in this PR (which, due to #133758, is now being adjusted slightly). I wonder what the history was in moving from the definition in the RFC, which seems more obvious and is in line with the 2001 Leijen paper, to the one now in the standard library. |
|
Looking at the git history, the current implementation goes all the way back to the original PR: Looking through the discussion there, and also in the RFC thread and the tracking issue... ...I don't immediately see any discussion about the mathematics to justify the change from the RFC specification to the current implementation, which is: pub fn rem_euclid(self, rhs: f64) -> f64 {
let r = self % rhs;
if r < 0.0 { r + rhs.abs() } else { r }
} |
|
Because,
So the changed expression would only differ when |
|
This PR should be rejected. |
|
It's a draft, so there's nothing to reject at this point. And yes, having analyzed it, I agree that |
traviscross
changed the title
{f16,f32,f64,f128}::rem_euclid{f16,f32,f64,f128}::div_euclid
|
(The state of this draft PR lags the current state of my analysis. I have a hypothesis for a fix for |
|
Closing in favor of: (...so as to fix the branch name.) |
DRAFT
We're still analyzing this.
cc #133485
See also:
<{f16,f32,f64,f128} as Rem>::remare not remainder of truncated division, as documented #133758"Division and Modulus for Computer Scientists",
Daan Leijen, 2001,
https://www.microsoft.com/en-us/research/wp-content/uploads/2016/02/divmodnote-letter.pdf
r? ghost