A few billion digits may not be enough to reproduce the problem if you
are using rand(..., "double") which is the default. I found it happens
more frequently when using rand(..., "single").
That said, I am not sure it is a bug. It may be expected
behavior for the Mersenne Twister to potentially return values at
either one or the other extreme of the range (0.0, 1.0). I can't find
the web page today, but I do remember reading a Matlab documentation
page that said for certain values of a and b, the _expression_ x = rand * a
+ b may return either a or b. I realize this is Octave and not Matlab,
but it is possible the implementation is similar enough that the same
bounds can occur.