ranges terribly sensitive to floating point

[alanedelman]

This one is kind of funny

julia> v=.1:.1:.3
0.1:0.1:0.2

[JeffBezanson]

We are using the formula floor((stop-start)/step)+1 to determine the length of the range. Is there a better one we should use? (in base/range.jl)

[JeffBezanson]

Ah yes:

julia> .1 + 2*.1
0.30000000000000004

So our current formula is consistent with how we compute the elements. I suppose we'd have to change both. I don't suppose it would be ok to adjust the endpoint to 0.30000000000000004.

[kmsquire]

More specifically, the issue is here:

julia> a = (.3-.1)/.1
1.9999999999999998

Adding eps(x) works, but I'm not sure it's justified:

julia> eps(a)
2.220446049250313e-16

julia> a+eps(a)
2.0

[JeffBezanson]

In this particular case, it turns out you can get the exact values you want (0.1, 0.2, 0.3) by adjusting the step size to 0.09999999999999999. I don't know how this generalizes, but it seems promising.

[toivoh]

I'm not sure that you can reasonably expect to avoid this kind of off-by-one errors for floating point ranges. The definition is intrinsically sensitive to roundoff errors.

[timholy]

We are using the formula floor((stop-start)/step)+1 to determine the length of the range. Is there a better one we should use? (in base/range.jl)

How about floor((stop-start)/step + sqrt(eps(T)))+1?

[timholy]

...or maybe you want to avoid the sqrt for performance reasons, in which case a better choice might be

nf = (stop-start)/step + 1
nf += abs(nf)*eps(T)
n = ifloor(nf)

[alanedelman]

I'll bet most of you didn't know you could almost write a book on colon alone. :-)

There are two quesions for start:delta:stop

1. How to compute the length and the truestop
2. How to compute the numbers in between

Regarding 2 a few comments

Strategy a) [start+i_delta for i=0:n] <--- terrible strategy
Strategy b) Go halfway forwards from start and halfway backwards from stop
(middle may have an odd length -- and odd/even considerations)
Maybe MATLAB does this
Strategy c) Figure out how to compute the closest floating point number to
start + (i/n)_(stop-start) for all numerical ranges
so start and stop are in the sequence without any rounding
[I prefer c]

Regarding 1

You'll have to decide whether to preserve the user's stop without roundng,
or to modify it drastically == depending on whether stop "snaps" to the grid
well enough.
I'll probably regret this as not being good enough but as something quick and dirty
I'd compute

nf = (stop-start)/step + 1

and

n=round(nf)

if abs(n-nf) < eps(n) (* 3 or 4 or 5)

I would declare stop as official and not to be modified. Otherwise, I would
use start + (n-1)*delta as the official stop

[StefanKarpinski]

I think that dismissing [start+i*delta for i=0:n] as a terrible strategy is a bit undue. It has the enormous advantage of being simple, easy to understand, and fast.

[JeffBezanson]

Excellent, thank you.
In our representation the degrees of freedom are the step and the length, so I'm trying this logic:

    nf = (stop-start)/step + 1
    n = round(nf)
    if abs(n-nf) < eps(n)*3
        step = (stop-start)/(n-1)
    end

[stevengj]

In the case where you modify step, you also might want to change the rule to i < (n-1)-i ? start+i*step : end - ((n-1)-i)*step, according to Alan's strategy 2(b). The main purpose of this seems to be to ensure that the final point is exactly end. A slightly simpler rule for that purpose would be i == n-1 ? end : start + i*step.

Also, shouldn't the rule use n = floor(nf), so that you never step beyond end (plus roundoff)? e.g. length(1:1:1.8) should be 1, not 2.

[JeffBezanson]

Ok, I pushed a new implementation that is probably no worse than the old one. I think the best thing to do now is find more test cases, and if the behavior for any of them seems bad enough I can try further improvements.

[stevengj]

Should be something like:

if abs(n-nf) < eps(n)*3
    step = (stop-start)/(n-1)
else
    n = floor(nf)
end

I think. Your current patch gives

julia> [1:1:1.8]
2-element Float64 Array:
 1.0
 2.0