Round Off Error in Iterator Arithmetic

When I first started writing Julia (i.e., recently) I created vectors of intermediate results and operated upon them, that is complicated versions of:

temp0 = [k for k = 1:3]
3-element Vector{Int64}:
 1
 2
 3

temp1 = map(x -> 2*x + 1, temp0)
3-element Vector{Int64}:
 3
 5
 7

But then I realized that others were creating iterators, manipulating the iterators, but delaying the create ion of the vector as long as possible. It seems that this should be more efficient, especially when the vectors no longer fit in cache. And that the ability to do this might be an important design pattern in Julia. So I rewrote a big chunk of code to do complicated versions of:

temp0 = 1:3
1:3

temp1 = 2*temp0 .+ 1
3:2:7

collect(temp1)
3-element Vector{Int64}:
 3
 5
 7

It broke most of my code. After a bit of fiddling I discovered that the problem is roundoff that causes the iterator to terminate one iteration too soon. Here is an example.

n = 569
 569
start32 = Float32(153.97725105285645)
 153.97725f0
stop32 = Float32(246.612811088562)
 246.61281f0

range = stop32 - start32
 92.63556f0
step = range / (n-1)
 0.16309078f0
typeof(step)   # just to check that I didn't do something wrong
 Float32

iter = start32 : step : stop32
 153.97725f0:0.16309078f0:246.44972f0
length(collect(iter))
 568

Notice that the “length of the iterator” is one short. So, at the very least, this is a bug.

If you know nothing about Julia, you would expect this problem, to arise when the round off for the value pointed to by step yields a value enough higher than the exact step value so that n times the error in step is greater than eps(type) times stop. Perhaps I should write down a formal proof just to be sure, but I think that for random values of n, start, and stop this should occur somewhat often.

But if you use Julia you know that it doesn’t occur often. As near as I can tell it never happens if start, stop, and step are Float64's. I think this means that some code somewhere accounts for roundoff and fixes this problem in the common case.

The problem only occurs if you use some other float type. Of course, using smaller floats is really hot right now, especially for GPUs (I’d like 18-bit and 24-bit floats too).

Interestingly

typeof(iter)
StepRangeLen{Float32, Float64, Float64}

So the step size and the stop value are Float64? I can partly confirm this …

iter.step
0.163090780377388

though I can’t figure out how to access the rest of the fields in iter. If I manually specify the length of the iterator then the end point may be off by a bit but the mixed type problem disappears.

iterManual = StepRangeLen(start32,step,n)
153.97725f0:0.16309078f0:246.61282f0

typeof(iterManual)
StepRangeLen{Float32, Float32, Float32}

length(collect(iterManual))
569

Not sure if this is a sufficiently general workaround for my purposes.

After all that, I see that range computes the correct sequence using a rational representation of the grid points, which is reduced to a float as the last step, so that the step jitters from point to point but you always have the exact start, stop, and number of points.