Different ways of declaring the same StepRangeLen are not equivalent

General Usage
1

Hey everyone,
This is my first post on here. I will do my best to follow the posted guidelines, but please let me know if I should change anything.

I noticed something slightly odd about the range function. Here is a minimal reproducible example.

npts = 21

y = 0.43693030000894467
r1 = range(start = -3y, stop = 3y, length = npts)
typeof(r1)
r2 = -3y:6y/(npts-1):3y
typeof(r2) 

r1 == r2

x = 2.0
r3 = range(start = -3x, stop = 3x, length = npts)
typeof(r3)
r4 = -3x:6x/(npts-1):3x
typeof(r4) 

r3 == r4

I have been using range with the syntax of r2 and r4 for a while now and encountered this issue just today. I noticed that for this specific value of y and this number of points, range behaves strangely. The (start:step:stop) syntax does not create the last point in the range, stop. The issue is easier to see if you compare:

collect(r1)
collect(r2)

Where collect(r2) creates a 20, not a 21-element Vector{Float64}.
I will probably simply switch how I am declaring StepRangeLen objects from now on, but does anybody have an explanation as to what is happening here?

Thank you for your help!

2

You’ve created a range where due to floating point precision in the step, it will not nicely end on 3y, thus skipping it.

Adding the floating point precision back to stop (even though eps(y)*npts is only 1.1657341758564144e-15 will demonstrate this:
r2 = -3y:6y/(npts-1):3y+eps(y)*npts.

If you want to have a fixed length, it’s best to use the range constructor.

edit: 7(eps(y)) is already enough here to make it to 21

1 Like
3

Ranges are tricky beasts, especially when it comes to floating point numbers. There are four key properties of a range (start, stop, number of points, and step size), but no matter how you cut it, one of those is fully dependent upon the other three.

So that’s fundamentally why there are multiple ways of specifying the range — and why you might get different results depending upon which way you go.

The start:step:stop formulation tries to fit as many steps as possible between start and stop. And the stop that you get back might be slightly less than the value you specified. But it’ll never be more. Depending upon the exact floating point values, it can sometimes be surprising that the stop you specified is a little short of where you expected the computed step to hit.

4 Likes
4

Thank you both, that makes a lot of sense!

5

Given the reasonable explanations here, this confused me for a good 5 minutes

julia> step(r1) === step(r2)
true

until I checked the internals:

julia> r1.step < r2.step
true

julia> r1.step
Base.TwicePrecision{Float64}(0.13107909000268325, 1.554312234475219e-16)

julia> r1.step.hi + r1.step.lo
0.13107909000268342

julia> r2.step
Base.TwicePrecision{Float64}(0.13107909000268342, 0.0)

and realized the relevant step method tossed out the extra precision in T(r.step):

julia> Float64(r1.step) === Float64(r2.step)
true

which is good for getting a familiar type to work with but bad for exact comparisons.

6

Just as a little tidbit, Julia does try pretty hard to do make start:step:stop hit the stop you specified with the step you meant, but only if “what you meant” could possibly round to the exact value you specified. For example, by the literally exact floating point values specified in the range 0.0:0.1:0.3 there should only be three values (0.0, 0.1, and 0.2). The exact value specified by 0.1 isn’t \frac{1}{10} but slightly more than that, and the exact value specified by 0.3 isn’t \frac{3}{10} but it’s slightly less… so the naive thing should miss 0.3. But there is a simple rational that rounds to 0.1 that can make it work, so Julia uses that instead! That’s what’s happening with the higher precision stuff @Benny is peeking at.

This whole “do what you might have meant” thing only works, however, if what you meant rounds exactly to the value you passed (and isn’t closer to any other floating point values). Doing more than one arithmetic computation in computing your step and/or endpoints yourself is a very easy way to get away from an exact round. Simplify your computations — specifically the “computation distance” between the values you use in the range — and things more often work as you expect:

julia> y = 0.43693030000894467
0.43693030000894467

julia> step = 6y/(npts-1)
0.13107909000268342

julia> r5 = -10step:step:10step
-1.3107909000268343:0.13107909000268342:1.3107909000268343

julia> length(r5)
21

julia> 10step == 3y
false

julia> println(10step); println(3y)
1.3107909000268343
1.310790900026834

Even better, of course, is to just specify the length directly with range if that’s really what you mean, but hopefully this helps you understand why a bit more.

2 Likes