Community
cancel
Showing results for 
Search instead for 
Did you mean: 
Highlighted
31 Views

Random-Number generator question about single precision

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.

Chris
0 Kudos
5 Replies
Highlighted
Valued Contributor II
31 Views

You probably mean 10**6 and

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.

0 Kudos
Highlighted
31 Views

From: https://software.intel

From: https://software.intel.com/content/www/us/en/develop/documentation/fortran-compiler-developer-guide-...

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

0 Kudos
Highlighted
New Contributor I
31 Views

You should be using MKL

You should be using MKL random number generators --

0 Kudos
31 Views

Thanks for the answers.

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. 

Chris
0 Kudos
Highlighted
Black Belt
31 Views

The scholar George Marsaglia

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 .

0 Kudos