I guess this gets back to the definition of “uniformly distributed”? It does not mean that the set of all possible floating point numbers is equally sampled – that distribution would indeed have more small numbers (<0.5) than large numbers (>0.5). Rather, it means that the numbers in [0.25,0.5) are sampled half as often as those in [0.5,1.0), and those in [0.125,0.25) are sampled half as often as those in [0.25,0.5), and so on. Each binary exponent should be sampled half as often as the next larger exponent. In principle, that can continue all the way down to the exponent for tiny(1.0), and within each exponent, the full set of fractions could be equally sampled. Such a PRNG could then return the full set of normal floating point numbers from 0.0 up to nearest(1.0,-1.0) with uniform distribution. As I said previously, I do not know of any fortran compiler that implements such a PRNG for its random_number() intrinsic. The reason, as discussed in this thread, is that the simpler division approach, which return only 1/126 of those numbers, is “good enough” for most applications.