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.