Bug in floating-point range?

General Usage
1 
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)
2 Likes
2

What’s the bug?

3

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

4

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

1 Like
5

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
1 Like
6

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
1 Like
7

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
1 Like
8

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

9

Looks like it

10

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…

8 Likes
11

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
3 Likes
12

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)

8 Likes
13
3 Likes