Intel® Fortran Compiler
Build applications that can scale for the future with optimized code designed for Intel® Xeon® and compatible processors.
Announcements
Welcome to the Intel Community. If you get an answer you like, please mark it as an Accepted Solution to help others. Thank you!
26730 Discussions

## Random-Number generator question about single precision

Beginner
187 Views

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.

5 Replies
Valued Contributor III
187 Views

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.

Black Belt
187 Views

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

Valued Contributor II
187 Views

You should be using MKL random number generators --

Beginner
187 Views

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.

Black Belt
187 Views

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 .