Palli
November 23, 2024, 5:38pm
8
Right (though that last one a bit of a tough one), also he reported more, e.g.:
opened 12:38AM - 21 Apr 23 UTC
bug
maths
float16
correctness bug ⚠
Julia has a [`div`](https://docs.julialang.org/en/v1/base/math/#Base.div) functi… on that is used to implement [floor division](https://docs.julialang.org/en/v1/base/math/#Base.fld) (`fld`) and [ceil division](https://docs.julialang.org/en/v1/base/math/#Base.cld) (`cld`). I think I’ve found a bug that causes all of these division functions to return incorrect results.
Floor division is documented to return the largest integer less than or equal to `x / y`. It should never return a number that is *greater* than `x / y`.
But it does:
```julia
julia> 514 / Float16(0.75)
Float16(685.5)
julia> div(514, Float16(0.75)) # should round down, but rounds up instead
Float16(686.0)
julia> fld(514, Float16(0.75)) # likewise
Float16(686.0)
julia> fld(514, Float16(0.75)) ≤ 514 / Float16(0.75)
false
```
Similarly, ceil division should never return a number that is *smaller* than regular division, but it does:
```julia
julia> 515 / Float16(0.75)
Float16(686.5)
julia> cld(515, Float16(0.75)) # should round up, but rounds down instead
Float16(686.0)
julia> cld(515, Float16(0.75)) ≥ 515 / Float16(0.75)
false
```
This behavior is not limited to 16-bit floats. Here’s a case where `fld` produces an incorrect result for `Float32` inputs:
```julia
julia> 4_194_307 / Float32(0.75) # = 5592409.5
5.5924095f6
julia> fld(4_194_307, Float32(0.75)) # = 5592410, incorrectly rounded up
5.59241f6
julia> fld(4_194_307, Float32(0.75)) ≤ 4_194_307 / Float32(0.75)
false
```
And here’s the same for `cld`:
```julia
julia> 4_194_308 / Float32(0.75) # = 5592410.5
julia> cld(4_194_308, Float32(0.75)) # = 5592410, incorrectly rounded down
5.59241f6
julia> cld(4_194_308, Float32(0.75)) ≥ 4_194_308 / Float32(0.75)
false
```
The equivalent operations in Python produce the correct results:
```python
# For 16-bit floats:
>>> np.float16(514) / np.float16(0.75) # regular division
685.5
>>> np.float16(514) // np.float16(0.75) # floor division
685.0
# For 32-bit floats:
>>> np.float32(4_194_307) / np.float32(0.75) # regular division
5592409.5
>>> np.float32(4_194_307) // np.float32(0.75) # floor division
5592409.0
```
Examples of this incorrect behavior are not hard to find – for most floats, you can find a divisor that will make either `fld` or `cld` return the wrong answer.
Here are some examples for `Float16` where either `fld` or `cld` is incorrect:
- `cld(1, Float16(0.000999)) < 1 / Float16(0.000999)`
- `cld(2, Float16(0.001999)) < 2 / Float16(0.001999)`
- `cld(3, Float16(0.002934)) < 3 / Float16(0.002934)`
- `cld(4, Float16(0.003998)) < 4 / Float16(0.003998)`
- `fld(5, Float16(0.004925)) > 5 / Float16(0.004925)`
And here are some for `Float32`:
- `fld(5, Float32(6.556511e-7)) > 5 / Float32(6.556511e-7)`
- `fld(10, Float32(1.3113022e-6)) > 10 / Float32(1.3113022e-6)`
- `fld(11, Float32(1.4305115e-6)) > 11 / Float32(1.4305115e-6)`
- `cld(16, Float32(2.8014183e-6)) < 16 / Float32(2.8014183e-6)`
- `cld(17, Float32(2.2053719e-6)) < 17 / Float32(2.2053719e-6)`
For simplicity I’ve presented examples where the first argument is an integer; this bug also occurs for non-integral inputs.
A divisor producing the wrong result can be found for over 51% of all possible 16-bit floats. I have not evaluated how widespread this is for `Float32`, but the results above suggest that it is similarly easy to find failures there too.
I’ve tracked the invocations down to [this definition](https://github.com/JuliaLang/julia/blob/72aec423c2ab9f80c249d63fdd68b35833cfd7ed/base/div.jl#L370) of `div`, which has existed at least as far back as [Julia 1.6](https://github.com/JuliaLang/julia/blob/f9720dc2ebd6cd9e3086365f281e62506444ef37/base/div.jl#L279):
```julia
div(x::T, y::T, r::RoundingMode) where {T<:AbstractFloat} =
convert(T, round((x - rem(x, y, r)) / y))
```
Most recent update in August.
The advantage of mapping ÷ to fld is that then % needs to be mod, which is not so uncommon to use with a negative first argument. Both rem and mod give you a useful cyclic behavior for positive numbers. If you go past zero mod keeps being cyclic, whereas rem makes a sign shift and goes over to a cycle with negative numbers. I’m sure there’s some context where rem is useful with negative first argument, but I’m yet to find it. In my experience mod is always the operator you want and if you carel…
I’m not generally enthusiastic about Python but this is something they got right, and Julia (and many other languages) sadly didn’t.
Julia does (or well tries) to do same as C and C++. It’s not clear to me Python’s way is better, either definition works, assuming people know which is used.
Also a (minor? i.e. not correctness, only stackoverflow) bug was fixed in OrderedCollections.jl, and the PR fix, approved in Feb 2021, got forgotten for 4 years, and I reported it some months before that with Bug in OrderedDict (unlike in LittleDict, and Dict) · Issue #65 · JuliaCollections/OrderedCollections.jl · GitHub
1 Like