The more I’ve thought about this, the more I’ve come to believe that the current Julia (and C, C++, Fortran, Python, Perl and Ruby) implementations are wrong. Here’s why:
The rand
function does not actually draw from the infinite set of points on the (open or closed) interval. It draws from the finite set 0:eps(T):1-eps(T)
. (Which has the convenient property that the elements can be represented in type T
.)
We might want to think about rand
as drawing from the infinite set (the open or closed interval) and then projecting (rounding) onto a finite set of floating-point numbers. But for that metaphor to work with the above-mentioned set, the interval must be (-eps(T)/2, 1-eps(T)/2)
. (The open or the closed interval. Mathematically, there’s no difference. The probability of hitting an endpoint in a finite number of draws is zero.)
Setting the least significant bit to 1, as suggested by @rfourquet would be equiivalent to drawing from eps(T)/2:eps(T):1-eps(T)/2
, which in turn would be equivalent to drawing from the (open or closed) continuous interval (0,1) and then correctly rounding to the nearest member of that set. This would not cost any extra time, and is the correct fix, i.m.o.