The drand48() and erand48() functions return values of type double. The
full 48 bits of r(n+1) are loaded into the mantissa of the returned
value, with the exponent set such that the values produced lie in the
interval [0.0, 1.0).
I changed the last line of the procedure like this to get the documented behavior:
The random function uses drand48:
I changed the last line of the procedure like this to get the documented behavior: