# Bug in floating-point range?

**URL:** <https://discourse.julialang.org/t/bug-in-floating-point-range/45706>\
**Category:** General Usage\
**Created:** [August 28, 2020, 8:53pm UTC](https://discourse.julialang.org/t/bug-in-floating-point-range/45706 "2020-08-28T20:53:19Z")\
**Posts on this page:** 13\
**Page:** 1

<div class="post-metadata">

**Author:** ![Mikhail\_Kagalenko](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mikhail_kagalenko/32/13257_2.png) [@Mikhail\_Kagalenko](https://discourse.julialang.org/u/Mikhail_Kagalenko)\
**Post date:** [August 28, 2020, 8:53pm UTC](https://discourse.julialang.org/t/bug-in-floating-point-range/45706/1 "2020-08-28T20:53:19Z")

</div>

```julia
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)

```

---

<div class="post-metadata">

**Author:** ![dpsanders](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpsanders/32/3573_2.png) [@dpsanders](https://discourse.julialang.org/u/dpsanders)\
**Post date:** [August 28, 2020, 9:12pm UTC](https://discourse.julialang.org/t/bug-in-floating-point-range/45706/2 "2020-08-28T21:12:03Z")

</div>

What’s the bug?

---

<div class="post-metadata">

**Author:** ![Henrique\_Becker](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/henrique_becker/32/15443_2.png) [@Henrique\_Becker](https://discourse.julialang.org/u/Henrique_Becker)\
**Post date:** [August 28, 2020, 9:13pm UTC](https://discourse.julialang.org/t/bug-in-floating-point-range/45706/3 "2020-08-28T21:13:41Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![Mikhail\_Kagalenko](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mikhail_kagalenko/32/13257_2.png) [@Mikhail\_Kagalenko](https://discourse.julialang.org/u/Mikhail_Kagalenko)\
**Post date:** [August 28, 2020, 9:13pm UTC](https://discourse.julialang.org/t/bug-in-floating-point-range/45706/4 "2020-08-28T21:13:50Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![Mikhail\_Kagalenko](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mikhail_kagalenko/32/13257_2.png) [@Mikhail\_Kagalenko](https://discourse.julialang.org/u/Mikhail_Kagalenko)\
**Post date:** [August 28, 2020, 9:16pm UTC](https://discourse.julialang.org/t/bug-in-floating-point-range/45706/5 "2020-08-28T21:16:30Z")

</div>

Maybe it will be clearer in this way:

```julia
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

```

---

<div class="post-metadata">

**Author:** ![dpsanders](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpsanders/32/3573_2.png) [@dpsanders](https://discourse.julialang.org/u/dpsanders)\
**Post date:** [August 28, 2020, 9:17pm UTC](https://discourse.julialang.org/t/bug-in-floating-point-range/45706/6 "2020-08-28T21:17:34Z")

</div>

There does indeed seem to be some fishiness here:

```julia
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

```

---

<div class="post-metadata">

**Author:** ![dpsanders](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpsanders/32/3573_2.png) [@dpsanders](https://discourse.julialang.org/u/dpsanders)\
**Post date:** [August 28, 2020, 9:19pm UTC](https://discourse.julialang.org/t/bug-in-floating-point-range/45706/7 "2020-08-28T21:19:07Z")

</div>

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

```julia
]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

```

---

<div class="post-metadata">

**Author:** ![Henrique\_Becker](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/henrique_becker/32/15443_2.png) [@Henrique\_Becker](https://discourse.julialang.org/u/Henrique_Becker)\
**Post date:** [August 28, 2020, 9:32pm UTC](https://discourse.julialang.org/t/bug-in-floating-point-range/45706/8 "2020-08-28T21:32:37Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![Mikhail\_Kagalenko](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mikhail_kagalenko/32/13257_2.png) [@Mikhail\_Kagalenko](https://discourse.julialang.org/u/Mikhail_Kagalenko)\
**Post date:** [August 28, 2020, 9:33pm UTC](https://discourse.julialang.org/t/bug-in-floating-point-range/45706/9 "2020-08-28T21:33:11Z")

</div>

Looks like it

---

<div class="post-metadata">

**Author:** ![ffevotte](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ffevotte/32/6587_2.png) [@ffevotte](https://discourse.julialang.org/u/ffevotte)\
**Post date:** [August 28, 2020, 9:36pm UTC](https://discourse.julialang.org/t/bug-in-floating-point-range/45706/10 "2020-08-28T21:36:55Z")

</div>

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
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 `BigFloat`s, but not `Complex{BigFloat}`s:

```julia
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 `BigFloat`s. But I haven’t looked more into it yet…

---

<div class="post-metadata">

**Author:** ![ffevotte](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ffevotte/32/6587_2.png) [@ffevotte](https://discourse.julialang.org/u/ffevotte)\
**Post date:** [August 28, 2020, 10:02pm UTC](https://discourse.julialang.org/t/bug-in-floating-point-range/45706/11 "2020-08-28T22:02:04Z")

</div>

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](https://github.com/JuliaLang/julia/pull/37281))_. But here is a simple fix that seems do the trick:

```julia
@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](https://github.com/JuliaLang/julia/blob/master/base/range.jl#L686_L690)

```julia
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

```

---

<div class="post-metadata">

**Author:** ![mbauman](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mbauman/32/31082_2.png) [@mbauman](https://discourse.julialang.org/u/mbauman)\
**Post date:** [August 28, 2020, 10:29pm UTC](https://discourse.julialang.org/t/bug-in-floating-point-range/45706/12 "2020-08-28T22:29:30Z")

</div>

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 `LinRange`s are particularly susceptible to this. They don’t hit the theoretically “best” intermediate values as robustly as the step ranges:

```julia
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
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 `BigFloat`s, how it compares when you use a `step` keyword instead of a `length`)

---

<div class="post-metadata">

**Author:** ![ffevotte](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ffevotte/32/6587_2.png) [@ffevotte](https://discourse.julialang.org/u/ffevotte)\
**Post date:** [August 29, 2020, 12:58pm UTC](https://discourse.julialang.org/t/bug-in-floating-point-range/45706/13 "2020-08-29T12:58:42Z")

</div>

> [@mbauman](#):
>
> the big issue that @ffevotte identifies with complex ranges. This is clearly a bug and should be reported if it hasn’t already.

> <https://github.com/JuliaLang/julia/issues/37276>
>
> 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.

> <https://github.com/JuliaLang/julia/pull/37281>
>
> 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?