Bug in floating-point range?

Mikhail_Kagalenko · August 28, 2020, 8:53pm

julia> rng=range(big(1)/100, stop=big(10)/100, length=10)
10-element LinRange{BigFloat}:
 0.010,0.020,0.030,0.040,0.050,0.060,0.070,0.080,0.090,0.10

julia> Float64(maximum(abs.(angle.(exp.(im*rng[1:9]))-rng[1:9])))
1.0795210693868056e-78

julia> Float64(maximum(abs.(angle.(exp.(im*rng))-rng)))
4.440892098500626e-18


julia> versioninfo()
Julia Version 1.5.1
Commit 697e782ab8 (2020-08-25 20:08 UTC)
Platform Info:
  OS: Linux (i686-pc-linux-gnu)
  CPU: Intel(R) Pentium(R) 4 CPU 3.00GHz
  WORD_SIZE: 32
  LIBM: libopenlibm
  LLVM: libLLVM-9.0.1 (ORCJIT, prescott)

dpsanders · August 28, 2020, 9:12pm

What’s the bug?

Henrique_Becker · August 28, 2020, 9:13pm

From 1 to 9 you have 9 elements. Not 10. The two lines do not do the same, if you were assuming that.

Mikhail_Kagalenko · August 28, 2020, 9:13pm

Different accuracy of the result depending on the length of the range. Large errors for 10 elements are not just in the last element.

Mikhail_Kagalenko · August 28, 2020, 9:16pm

Maybe it will be clearer in this way:

ulia> Float64.(abs.(angle.(exp.(im*rng[1:9]))-rng[1:9]))
9-element Array{Float64,1}:
 0.0
 0.0
 0.0
 0.0
 0.0
 8.636168555094445e-78
 0.0
 8.636168555094445e-78
 0.0

julia> Float64.(abs.(angle.(exp.(im*rng))-rng))
10-element Array{Float64,1}:
 0.0
 1.1102230246251566e-18
 1.1102230246251566e-18
 1.1102230246251566e-18
 2.220446049250313e-18
 2.220446049250313e-18
 3.3306690738754695e-18
 1.1102230246251566e-18
 4.440892098500626e-18
 0.0

dpsanders · August 28, 2020, 9:17pm

There does indeed seem to be some fishiness here:

julia> collect(im * rng[1:9])
9-element Array{Complex{BigFloat},1}:
 0.0 + 0.009999999999999999999999999999999999999999999999999999999999999999999999999999995im
 0.0 + 0.01999999999999999999999999999999999999999999999999999999999999999999999999999999im
 0.0 + 0.02999999999999999999999999999999999999999999999999999999999999999999999999999985im
 0.0 + 0.04000000000000000000000000000000000000000000000000000000000000000000000000000052im
 0.0 + 0.05000000000000000000000000000000000000000000000000000000000000000000000000000011im
 0.0 + 0.0599999999999999999999999999999999999999999999999999999999999999999999999999997im
 0.0 + 0.07000000000000000000000000000000000000000000000000000000000000000000000000000091im
 0.0 + 0.07999999999999999999999999999999999999999999999999999999999999999999999999999996im
 0.0 + 0.09000000000000000000000000000000000000000000000000000000000000000000000000000009im

julia> collect(im * rng[1:10])
10-element Array{Complex{BigFloat},1}:
 0.0 + 0.009999999999999999999999999999999999999999999999999999999999999999999999999999995im
 0.0 + 0.01999999999999999888977697537484345957636833190917968749999999999999999999999987im
 0.0 + 0.0299999999999999988897769753748434595763683319091796875im
 0.0 + 0.03999999999999999888977697537484345957636833190917968750000000000000000000000013im
 0.0 + 0.04999999999999999777955395074968691915273666381835937499999999999999999999999987im
 0.0 + 0.06000000000000000222044604925031308084726333618164062500000000000000000000000047im
 0.0 + 0.06999999999999999666933092612453037872910499572753906250000000000000000000000082im
 0.0 + 0.08000000000000000111022302462515654042363166809082031250000000000000000000000035im
 0.0 + 0.08999999999999999555910790149937383830547332763671874999999999999999999999999961im
 0.0 + 0.1000000000000000000000000000000000000000000000000000000000000000000000000000002im

dpsanders · August 28, 2020, 9:19pm

This seems to reduce to how the iteration through the range happens:

]julia> collect(rng[1:10])[1:9] .- collect(rng[1:9])
9-element Array{BigFloat,1}:
 0.0
 0.0
 2.698802673467013945433234957125124865973750113886337932819907334427684938488258e-79
 0.0
 0.0
 0.0
 1.079521069386805578173293982850049946389500045554535173127962933771073975395303e-78
 0.0
 0.0

Henrique_Becker · August 28, 2020, 9:32pm

So, LinRange compute the elements in a different way depending how how it is subscribed?

Mikhail_Kagalenko · August 28, 2020, 9:33pm

Looks like it

ffevotte · August 28, 2020, 9:36pm

Yep, there seems to be an issue with complex ranges not being as accurate as real ones.

We can extract from your example an even more minimal case:

julia> rng=range(big(1)/100, stop=big(10)/100, length=10)
10-element LinRange{BigFloat}:
 0.010,0.020,0.030,0.040,0.050,0.060,0.070,0.080,0.090,0.10

julia> im*rng
10-element LinRange{Complex{BigFloat}}:
0.0+0.010im,0.0+0.020im,0.0+0.030im,0.0+0.040im,0.0+0.050im,0.0+0.060im,0.0+0.070im,0.0+0.080im,0.0+0.090im,0.0+0.10im

# Both ways of computing the end points yield the same results:
julia> im*rng[1] == (im*rng)[1] && im*rng[end] == (im*rng)[end]
true

# But indexing any other element than the end points yields different results
julia> im*rng[9]
0.0 + 0.09000000000000000000000000000000000000000000000000000000000000000000000000000009im

julia> (im*rng)[9]
0.0 + 0.08999999999999999555910790149937383830547332763671874999999999999999999999999961im

In other words, im*rng is exactly what we expect it to be, but getindex does not work as accurately for both ranges.

A bit of inspection using the debugger seems to indicate that the culprit is a Base.lerpi function that is called internally, and is specialized for BigFloats, but not Complex{BigFloat}s:

julia> methods(Base.lerpi)
# 3 methods for generic function "lerpi":
[1] lerpi(j::Integer, d::Integer, a::Rational, b::Rational) in Base at rational.jl:467
[2] lerpi(j::Integer, d::Integer, a::BigFloat, b::BigFloat) in Base.MPFR at mpfr.jl:1033
[3] lerpi(j::Integer, d::Integer, a::T, b::T) where T in Base at range.jl:687

julia> Base.lerpi(8, 9, rng[1], rng[end])
0.09000000000000000000000000000000000000000000000000000000000000000000000000000009

julia> Base.lerpi(8, 9, im*rng[1], im*rng[end])
0.0 + 0.08999999999999999555910790149937383830547332763671874999999999999999999999999961im

So I guess this could be fixed by defining a specialized method for Base.lerpi working on Complex{BigFloat}s, probably in the same way as the one working on BigFloats. But I haven’t looked more into it yet…

ffevotte · August 28, 2020, 10:02pm

This is less elaborate than what gets done in the BigFloat case. And it would probably need more careful inspection in order to check that it does not beak anything nor slows down unrelated cases (EDIT: it does break other things; please see PR#37281). But here is a simple fix that seems do the trick:

@eval Base begin
    function lerpi(j::Integer, d::Integer, a::T, b::T) where T
        @_inline_meta
        # t is currently a Float64 independently of T.
        # If T is smaller than Float64, we want it to stay that way.
        # If T is larger than Float64, let's make sure j/d is computed with extra precision.
        t = T(j)/Float64(d)
        T((1-t)*a + t*b)
    end
end

For comparison, here is the current version of this function:
https://github.com/JuliaLang/julia/blob/master/base/range.jl#L686_L690

julia> rng=range(big(1)/100, stop=big(10)/100, length=10)
10-element LinRange{BigFloat}:
 0.010,0.020,0.030,0.040,0.050,0.060,0.070,0.080,0.090,0.10

julia> abs.(angle.(exp.(im*rng[1:9]))-rng[1:9])  .|> Float64
9-element Array{Float64,1}:
 0.0
 0.0
 2.698802673467014e-79
 5.397605346934028e-79
 0.0
 5.397605346934028e-79
 1.0795210693868056e-78
 0.0
 1.0795210693868056e-78

julia> abs.(angle.(exp.(im*rng))-rng)  .|> Float64
10-element Array{Float64,1}:
 0.0
 0.0
 5.397605346934028e-79
 0.0
 5.397605346934028e-79
 0.0
 0.0
 0.0
 1.0795210693868056e-78
 0.0

mbauman · August 28, 2020, 10:29pm

There are two things going on here — one is the big issue that @ffevotte identifies with complex ranges. This is clearly a bug and should be reported if it hasn’t already.

The other is significantly smaller (in magnitude) — subsetting ranges may indeed shift values by an ULP or so, and LinRanges are particularly susceptible to this. They don’t hit the theoretically “best” intermediate values as robustly as the step ranges:

julia> collect(LinRange(0, .8, 9))
9-element Vector{Float64}:
 0.0
 0.1
 0.2
 0.30000000000000004
 0.4
 0.5
 0.6000000000000001
 0.7000000000000001
 0.8

When you recompute a subset, we grab that imperfect intermediate value and use it as the new exact endpoint. This can end up changing the computation of the other values.

julia> collect(LinRange(0, .8, 9)[1:8])
8-element Vector{Float64}:
 0.0
 0.1
 0.2
 0.3
 0.4
 0.5000000000000001
 0.6
 0.7000000000000001

(Note how this compares with the default Float64 range collect(range(0, .8, length=9)) or if you’re using BigFloats, how it compares when you use a step keyword instead of a length)

ffevotte · August 29, 2020, 12:58pm

github.com/JuliaLang/julia

Inaccuracies for `LinRange`

opened 12:55PM - 29 Aug 20 UTC

ffevotte

ranges

This is a follow-up to discourse thread ["Bug in floating-point range?"](https:/…/discourse.julialang.org/t/bug-in-floating-point-range/45706?u=ffevotte). The most minimal example I could come up with to evidence the issue is the following one: ```julia function test_range(T) @show r = range(T(0), stop=T(1), length=11) @show val = r[2] ref = 1//10 err = (Rational{BigInt}(val) - ref) / ref |> abs |> Float64 @show err end ``` Indexing into a `LinRange{BigFloat}` gives very accurate results, with errors in the order of ulp(BigFloat): ```julia julia> test_range(BigFloat); r = range(T(0), stop = T(1), length = 11) = range(0.0, stop=1.0, length=11) val = r[2] = 0.1000000000000000000000000000000000000000000000000000000000000000000000000000002 err = 2.1590421387736112e-78 ``` However, the same calculation performed on a `LinRange{Complex{BigFloat}}` is much less accurate, with errors in the same order of magnitude as if 64-bit floats had been used: ```julia julia> test_range(Complex{BigFloat}); r = range(T(0), stop = T(1), length = 11) = range(0.0 + 0.0im, stop=1.0 + 0.0im, length=11) val = r[2] = 0.1000000000000000055511151231257827021181583404541015625 + 0.0im err = 5.551115123125783e-17 julia> test_range(Complex{Float64}); r = range(T(0), stop = T(1), length = 11) = range(0.0 + 0.0im, stop=1.0 + 0.0im, length=11) val = r[2] = 0.1 + 0.0im err = 5.551115123125783e-17 ``` <br/> As mentioned in the discourse thread above, I think I have identified `Base.lerpi` as the origin of the issue: ```julia julia> Base.lerpi(1, 10, big"0.", big"1.") 0.1000000000000000000000000000000000000000000000000000000000000000000000000000002 julia> Base.lerpi(1, 10, Complex(big"0."), Complex(big"1.")) |> real 0.1000000000000000055511151231257827021181583404541015625 ``` and I think I have a fix for this. I'll try and propose a PR as soon as I can.

github.com/JuliaLang/julia

More accurate Base.lerpi for higher-precision numbers

JuliaLang:master ← ffevotte:issue37276

opened 09:07PM - 29 Aug 20 UTC

ffevotte

+10 -1

This should fix #37276. Examples of things that are more accurate with this P…R include: ``` julia> LinRange(Complex{BigFloat}(0), Complex{BigFloat}(1), 11)[2] 0.1000000000000000000000000000000000000000000000000000000000000000000000000000002 + 0.0im julia> LinRange([0//1], [1//1], 11)[2] 1-element Vector{Rational{Int64}}: 1//10 ``` in comparison to the current master: ``` julia> LinRange(Complex{BigFloat}(0), Complex{BigFloat}(1), 11)[2] 0.1000000000000000055511151231257827021181583404541015625 + 0.0im julia> LinRange([0//1], [1//1], 11)[2] 1-element Array{Rational{Int64},1}: 3602879701896397//36028797018963968 ``` I think that the proposed implementation of a generic `Base.lerpi` makes it less useful to have a specific method for rationals (such as defined in [rational.jl:467](https://github.com/JuliaLang/julia/blob/master/base/rational.jl#L467)). One potential downside to the way things are done in this PR, is that there could potentially be a loss of accuracy w.r.t how things are currently done, when working on ranges of small (i.e. less than 64 bits) floating-point numbers. I'm not sure whether that was intended, but in the current implementation, the `j/d` quotient usually results in `t` being a `Float64`, meaning that computations involving `a` and `b` will be promoted to `Float64` too if the end points are smaller FP numbers. With the current implementation, all calculations use the same precision as the range end points. It would not be too difficult to restore the current master behavior for smaller FP number types, but I'm not sure whether it is worth making the implementation more complicated. In particular, I have no idea how many users expect and rely on the fact that points in a `LinRange{Float16}` are computed in double precision before being rounded to half precision. Would somebody have some insight on this?