Following expression gives inconsistent results for different values of n.
r = 0.0:pi/n:pi
I expect it to be length n+1 and include pi. For instance, for n=24 and n=26 the behavior is as I expected, but for n=25 it does not include pi.
I guess rounding error in pi/n might be the reason. I am using range(0.0,stop=pi,length=n+1) to avoid this but I wanted to bring this up in case this is an unknown issue.
Unfortunately, this is a basic issue with floating-point arithmetic itself:
julia> 25*(pi/25)
3.1415926535897936
julia> float(pi)
3.141592653589793
There’s no principle way to know that you intended the floating-point value pi/25 to be 1/25 of pi. The correct way to express this range is with range as you are doing. This will soon be expressible as range(0, pi, length=26) which is slightly shorter and nicer.
Thanks! Actually I noticed this while translating some of my code from Matlab. Apparently Matlab has a way of dealing with this: it works as expected for any n in Matlab. Do you plan to depreciate the colon syntax and/or range(0.0,stop=pi,step=pi/n)? If not, it might be nice to fix this.
Matlab just assumes that you meant to hit the endpoint exactly if you’re within a few (5 IIRC) ulps, which is… of highly questionable correctness. It’s possible this can be fixed in a more disciplined way, but it’s fairly tricky.
See also the discussion in ranges with ambiguous lengths · Issue #20532 · JuliaLang/julia · GitHub and Issues · JuliaLang/julia · GitHub Matlab handles this by rounding the endpoint to a few ulps.
For example, 0:(pi/25)*(1 + 3e-16):pi in Matlab also has length 26. The tradeoff with this is that the n-th point is no longer equal to initial + (n-1)*step. If you round too much, you can also increase the variation in the step size (i.e. diff(diff(somerange))), but Matlab’s endpoint rounding is small enough that this variation doesn’t seem worse than what we get anyway from roundoff errors in the cases I’ve checked.
range(0, stop=pi, length=n+1) is the “Julian way” to go then. Thanks guys.