> random() may return 1.0 exactly

That shouldn't be possible.  Although the code does assume C doubles have at least 53 bits of mantissa precision (in which case it does arithmetic that's exact in at least 53 bits - cannot round up to 1.0; but _could_ round up if the platform C double has less than 53 bits of precision).

> py> x = 0.9999999999999999
> py> for i in range(1, 1000000):
> ...     if int(i*x) == i:
> ...             print i
> ...             break
> ...
> 2049

Very surprising!  Which platform & Python is that?  The loop runs to completion on my box:

Python 2.7.9 (default, Dec 10 2014, 12:28:03) [MSC v.1500 64 bit (AMD64)] on win32