Fix incorrect results from floating-point division (fld, cld, div)

[lvlte]

This addresses #49450 (more details in this comment).

[adienes]

unfortunately does not solve the whole issue

julia> x = 1.125899906841625e15
1.125899906841625e15

julia> y = 0.1
0.1

julia> div(x, y, RoundDown)
1.125899906841625e16

julia> div(big(x), big(y), r)
1.1258999068416249e16

[lvlte]

@adienes @Seelengrab Float64 cannot represent the value 1.1258999068416249e16.

julia> x, y = 1.125899906841625e15, 0.1
(1.125899906841625e15, 0.1)

julia> prevfloat(x/y), x/y, nextfloat(x/y)
(1.1258999068416248e16, 1.125899906841625e16, 1.1258999068416252e16)

As already mentioned in the issue :

If eps(x/y) > 1, x/y rounds to an unsafe integer which can't be floored/ceiled if it needs to because x/y ± 1 is not representable (except for 2^53 - 1).

[adienes]

Float64 can represent 1.1258999068416248e16 though, which would be the correct result under RoundDown

[lvlte]

Float64 can represent 1.1258999068416248e16 though, which would be the correct result under RoundDown

Could you at least explain why do you think it should be the correct result ?

fld is defined as the "Largest integer less than or equal to x / y", and 1.1258999068416248e16 is not the largest integer less than or equal to x / y because the remainder is not significant :

julia> 1.125899906841625e15 / 0.1
1.125899906841625e16

julia> (1.125899906841625e15 - 1.125899906841625e15 % 0.1) / 0.1
1.125899906841625e16

julia> 1.125899906841625e15 / 0.1 - (1.125899906841625e15 % 0.1) / 0.1
1.125899906841625e16

Getting 1.1258999068416248e16 would implies the remainder is significantly greater than it really is.

Also, 1.1258999068416248e16 is off by one just as 1.125899906841625e16, but the latter is closer to the exact result of x divided by y, and rounding in that direction is consistent with the behavior of floating-point rounding when there is a tie, which is to round to the nearest integer.

[adienes]

so your concern is that since the condition is impossible (since the largest integer not exceeding x/y is not in fact representable in the return type), we should be allowed to violate the rounding conditions?

I don't mean that question as sophistry; it is a fair argument and one that has appeared a few times to attempt to justify behavior of accidental integer overflow e.g. in #34325, also discussed in #50486 . but I think due to the high performance needs of functions like div and fld we won't reach the same conclusions as in those issues.

there are basically three paths forward as I see it:

• when the conditions in the docstring cannot be satisfied, throw a DomainError. this is the approach taken by, for example, gcd and lcm
• promote to a type with higher precision that can represent the correct value
• amend the docstring to make a promise it can actually keep, and then amend the implementation to match that new specification

unfortunately the first two paths would be disastrous for performance and type-stability and are non-starters, so by process of elimination that leaves only the third. I would suggest an amendment to the docstring --- and, I believe, this is what most discussion participants in #49450 were implicitly assuming already --- to something like

return the largest integer representable in T less than or equal to x / y

and under this definition, a type-stable and fast /(::Float64, ::Float64) is at least possible, although it is not the current implementation

[lvlte]

You are missing the point : the floating-point division x/y returns an integer and its remainder, even normalized, is not significant enough compared to the result, so for consistency the fld and cld for the same type should give the same result. Using more precision should remain at the discretion of the user when using these functions. I don't think we should tweak every operation result considering the same operation with types allowing more precision. This would be like trying to fix all cases where x - x % y equals x in Float64 with x ± eps(x) because the results are different with BigFloat, a completely arbitrary choice which is not requested by the user and that could break consistency between operations.

I'm not comfortable with the idea that abs(cld(x, y) - fld(x, y)) equals 2 (in this example) as it would imply for all cases where eps(x/y) > 1 that the absolute difference should be equal to eps(x/y), which would make the results dependent on the magnitude of x/y, instead of having a consistent value that is either 0 or 1.

return the largest integer representable in T less than or equal to x / y.

1.125899906841625e16 still satisfies the premise.

Also, the current implementation gives the same results for those case where eps(x/y) > 1. This PR solves the issue for all other cases where cld, fld, div, give an inexact result and where the exact result is representable.

You didn't explain why 1.1258999068416248e16 should be preferred, and what is violating the rounding conditions.

[adienes]

here's what I'm looking at:

julia> x, y
(1.125899906841625e15, 0.1)

julia> big(x) / big(y)
1.125899906841624937500000000055459109526978519729750209567102918209446384572968e16

so to my eyes the "correct" integer would be 11258999068416249 . but this is not representable in Float64. the reason I suggest to round downwards to 11258999068416248 rather than upwards to 11258999068416250 is that is usually what I would assume flooring division (aka RoundDown) does

as a side note:

You are missing the point ... You didn't explain

I would invite you to use phrasing that more clearly reaffirms that we are collaborating together to solve the bug, rather than that we are arguing 🙂

[lvlte]

so your concern is that [since ...] we should be allowed to violate the rounding conditions?

doesn't sound like "we are collaborating together to solve the bug" but rather like "we are arguing" about the legitimacy of a PR that doesn't improve anything or even asks to break things because my concern is. Again, the implementation I propose fixes cld, fld, div, giving an inexact result where the exact result is representable, and for those edge cases where eps(x/y) > 1, it doesn't imply being allowed to violate anything but just matches the current implementation, which i think is correct and I explained why.

[adienes]

well I can appreciate the fervent defense and I do apologize for having given off the wrong impression during review.

this might be a fun one to discuss on triage, so I'll tag for review on the next meeting

[Seelengrab]

FWIW, the other partial fix in #49637 for Float16 & Float32 as far as I can tell makes the same decision that @adienes mentions, in particular see #49637 (comment)

So it would be good to be consistent and have Float64 behave the same.

[lvlte]

The comment you mention doesn't involve x,y such that eps(x/y) > 1, but yes it reminds us that having x/y rounded to an integer doesn't mean the mathematical result is an integer. To better illustrate my point, i would like to insist on the consistency between the quotient and remainder, and between the operations that involve them :

Example 1 eps(x/y) ≤ 1 (no disagreement, this PR just fixes inexact results) :

julia> x, y = 1.0, 0.1
(1.0, 0.1)

julia> a, r = divrem(x, y, RoundDown)    # accurate
(9.0, 0.09999999999999995)

julia> a + r/y == x/y, a + r/y, x/y      # accurate, consistent
(true, 10.0, 10.0)

Example 2 eps(x/y) > 1 (current behavior, which this PR currently doesn't affect) :

julia> x, y = 1.125899906841625e15, 0.1
(1.125899906841625e15, 0.1)

julia> a, r = divrem(x, y, RoundDown)    # inaccurate (intrinsic roundoff)
(1.125899906841625e16, 0.03750000000005546)

julia> a + r/y == x/y, a + r/y, x/y      # inaccurate (intrinsic roundoff), consistent
(true, 1.125899906841625e16, 1.125899906841625e16)

Now from example 2, let say because the remainder is not zero we force a to be rounded further down to the previous representable integer :

julia> a = 1.1258999068416248e16         # inaccurate (arbitrary roundoff)
1.1258999068416248e16

julia> a + r/y == x/y, a + r/y, x/y      # inaccurate (arbitrary roundoff), inconsistent
(false, 1.1258999068416248e16, 1.125899906841625e16)

[lvlte]

Regarding the behavior of div with Float16/Float32, I don't think the partial fix actually makes the decision of forcing rounding further down/up in those situations where eps(x/y) > 1 and where the remainder is not significant. I mean, currently div() might behave like that with the partial fix in certain cases, but not in others (which is problematic, I'm just discovering this), eg. if we compute fld as Float32(round(Float64(x) / Float64(y), RoundDown)) step by step, we get :

julia> x, y = 2.096153f6, 0.1f0
(2.096153f6, 0.1f0)

julia> round(Float64(x) / Float64(y), RoundDown)
2.0961529e7

julia> Float32(2.0961529e7)    # rounded down
2.0961528f7

julia> x, y = nextfloat(x), y + 2*eps(y)
(2.0961531f6, 0.10000002f0)

julia> round(Float64(x) / Float64(y), RoundDown)
2.0961527e7

julia> Float32(2.0961527e7)    # rounded up
2.0961528f7

The rounding for non-representable numbers is done by Float32() here and it can go both ways if not set explicitly.

If the decision was made to round down/up further for non representable numbers, we should have a partial fix that makes the "final" rounding explicit :

div(x::Float32, y::Float32, r::RoundingMode) = Float32(round(Float64(x) / Float64(y), r), r)
div(x::Float16, y::Float16, r::RoundingMode) = Float16(round(Float32(x) / Float32(y), r), r)

Not only is this not the case, but the comment that comes with that fix tells us cases where eps(x/y) > 1 are not addressed :

# Vincent Lefèvre: "The Euclidean Division Implemented with a Floating-Point Division and a Floor"
# https://inria.hal.science/inria-00070403
# Theorem 1 implies that the following are exact if eps(x/y) <= 1
div(x::Float32, y::Float32, r::RoundingMode) = Float32(round(Float64(x) / Float64(y), r))

Also, #49450 was closed as completed after the Float16/Float32 fix, and reopened shortly after only to address inexact results with Float64 and mixed types involving Float64, not to address those cases where eps(x/y) > 1 with Float16/Float32. So I'm wondering (even more) why the concern of solving unsolvable cases suddenly become a necessity in the meantime, without anybody mentioning it in the issue, or even just asking (given the comment above) "what about those cases where eps(x/y) > 1" ?

Now, if we were to leave those edge cases aside for Float64 as they are for Float16/Float32, the new fix proposed here should work seamlessly for all of them (I mean without the partial fix, but I haven't tested). In the other case, if eps(x/y) > 1 ... can still be done once in the generic div() definition.

[adienes]

as a note for triage: it does not seem like there is a whole ton of precedent to draw from. but I guess there is one in Common Lisp:

(format t "~a~%" (ffloor 1.125899906841625d15 0.1d0)) ;; 1.125899906841625d16

so they currently do the same thing Julia does here (and round upwards < eps, rather than downwards > eps)