Random-Number generator question about single precision

[pike__christopher]

I'm just curious about this:  I want a random uniform integer between and inclusive of 1 and 1**6.   To get this, I use a single precision real and the following code:

My code is 

integer, parameter ::  SinglePrecision = selected_real_kind(5) 

real (kind = SinglePrecision) ARandomReal

Call Random_Number(ARandomReal)

Print *, Ceiling(  (1.0-ARandomReal)*1**6)

Now, here is the question:  Suppose I asked for an integer between 1 and 1**11?   I'm sure the random number generator would fail.  A single precision real cannot be mapped one-to-one to such a large integer, but where would this go wrong?  I've done a quick test, and the numbers still look random, but I haven't dug deep into the tests for randomness.

Any theories on how Random_Number would go wrong.

[Arjen_Markus]

You probably mean 10**6 and 10**11.

The number of distinct random integers that you can produce this way is limited by the number of single-precision numbers that fit in the interval 0-1.0. If you want to go up to 10**11, then you will see more or less systematic gaps between the generated numbers, as the number of decimals is 5 to 6, so at most 10**5 distinct numbers in that interval. It is possible to be more precise than that, but this should give you a rough idea.

[jimdempseyatthecove]

From: https://software.intel.com/content/www/us/en/develop/documentation/fortran-compiler-developer-guide-and-reference/top/language-reference/a-to-z-reference/q-to-r/random-number.html

The RANDOM_NUMBER generator uses two separate congruential generators together to produce a period of approximately 10**18, and produces real pseudorandom results with a uniform distribution in [0, 1). It accepts two integer seeds, the first of which is reduced to the range [1, 2147483562]. The second seed is reduced to the range [1, 2147483398]. This means that the generator effectively uses two 31-bit seeds.

This seems to imply that internally, the result will be generated internally as a 62-bit integer result that is then partially used as the mantissa of the floating point REAL result. Typically a center portion of the 62-bit product (23 bits) is excised for the mantissa. It is unclear that DRAND or  DRANDM is called when a DOUBLEPRECISION scalar/array is used as an argument .OR. that the internally produced REAL is promoted to DOUBLEPRECISION.

While spreading the range, as Arjen reports is adequate for your purposes, when it is not, consider something like:

integer, parameter ::  SinglePrecision = selected_real_kind(5) 
integer, parameter ::  DoublePrecision = selected_real_kind(11) 
real (kind = SinglePrecision) ARandomReal
real (kind = DoublePrecision) ARandomDouble

Call Random_Number(ARandomReal)
Call Random_Number(ARandomDouble)
ARandomDouble = ARandomDouble * ARandomReal
Print *, Ceiling(  (1.0_DoublePrecision-ARandomDouble)*(10_DoublePrecision**11) )

*** The above ASSUMES Random_Number internally produces a SinglePrecision REAL. If it does not, then the above is not necessary.

Jim Dempsey

[JohnNichols]

You should be using MKL random number generators --

[pike__christopher]

Thanks for the answers.  These sort of confirm my guess, which is that the random integers I try to generate would develop gaps.   As for the MKL generator, I'm sort of making the assumption that the people who wrote fortran are smarter than me, so I use theirs. 

[mecej4]

The scholar George Marsaglia (1924-2011) wrote a famous paper on this property of some classes of pseudo-random numbers.

Marsaglia's Theorem: Random numbers fall mainly in the planes. See https://www.pnas.org/content/pnas/61/1/25.full.pdf .