First, call Random_seed(put = seed) requires seed be an array of 33 integers. I cannot understand how it can possibly use 33 integers. Isn’t this algorithm a markov chain, where the next output depends only on the previous output? If so, then only a single seed should be enough.
Second, the argument of Call Random_number is a 64 bit real which by convention is called RandomReal. However, in the past, I let RandomReal be a 32 bit real. The result was that the output would sometimes equal exactly zero. Whereas if RandomReal is 64 bits then the output is never equal to exactly zero. If Call Random_number generates exactly zero, then it breaks my code. From a theoretical point of view, I would argue that the output should never equal zero. So the moral to the story I guess is that RandomReal must be 64 bits. I didn’t see this in the documentation. Shouldn’t this be made clear in the documentation?
ps. My best guess about the first question is that the algorithm does not trust the programmer to create a “high quality” seed. So it lets the user enter 33 integers which individually might be low quality seeds (e.g. 5142), and it tries to combine them in some way to come up with a high quality seed.
Not quite, I’am afraid.
From the F2018 Standard Interpretation document:
16.9.156 RANDOM_NUMBER (HARVEST)
1 Description. Generate pseudorandom number(s).
2 Class. Subroutine.
3 Argument. HARVEST shall be of type real. It is an INTENT (OUT) argument. It may be a scalar or an array. It is assigned pseudorandom numbers from the uniform distribution in the interval 0 ≤ x < 1.
As you see, the argument may be a REAL entity of any kind and the result can be exact 0 while cannot be exact 1.
Good day msz59. So you are saying that Call Random_Number can equal zero even with a 64 bit output? I did a test of this. I repeatedly called Random_number with a 32 bit input, and I got zeros right away. But when I repeated the test with a 64 bit argument I never found a zero. I realize that zeros should occur less frequently with a 64 bit argument. I also realize that I only ran the test for a few hours. But still, I find it surprising that I have never encountered a zero with the 64 bit argument.
Also, from my very crude understanding of the Call Random_Number algorithm, it appears that if it did ever generate an exact zero, it would get stuck, and only generate zeros from that point forward. If you understand the algorithm better, maybe you can tell me that is wrong.
Here is a subroutine (and a test program that calls it) that sets the array of seeds used by random_number() given a single integer input, by using that integer as the initial state of an integer random number generator.
module set_random_seed_mod
implicit none
private
public :: set_random_seed, random_int_32_scalar
contains
pure elemental subroutine random_int_32_scalar(jsr,iran)
! generate random 32-bit integers
integer, intent(in out) :: jsr ! state of RNG
integer, intent(out) :: iran ! random integer
integer :: jz
jz = jsr
jsr = ieor(jsr, ishft(jsr, 13))
jsr = ieor(jsr, ishft(jsr, -17))
jsr = ieor(jsr, ishft(jsr, 5))
iran = jz + jsr
end subroutine random_int_32_scalar
!
subroutine set_random_seed(iseed)
! set the array of seeds used by the random_number() intrinsic
integer, intent(in) :: iseed
integer :: i, jseed, nseeds
integer, allocatable :: seeds(:)
jseed = iseed ! copy iseed to jseed since iseed is intent(in)
call random_seed(size=nseeds)
allocate (seeds(nseeds))
do i=1,nseeds
call random_int_32_scalar(jseed,seeds(i))
end do
call random_seed(put=seeds)
end subroutine set_random_seed
end module set_random_seed_mod
!
program xset_random_seed
use set_random_seed_mod, only: set_random_seed
use iso_fortran_env
implicit none
integer, parameter :: dp = kind(1.0d0), nrepetitions=2
integer :: iseed, irep
real(kind=dp) :: x(5)
print "('compiler_version: ',a)", compiler_version()
iseed = 0
do iseed=0,3
do irep=1,nrepetitions
call set_random_seed(iseed)
call random_number(x)
print "(i0,*(f7.4))", iseed, x
end do
end do
end program xset_random_seed
I’m sorry Beliavsky, but I don’t understand how this algorithm is useful. It isn’t too much bother to create an array of 33 integers for my seed. I just don’t understand what Call Random_Number wants with a dimension 32 seed array, when one seed number could theoretically do this job?
| “It is unlikely that 32 integers provided by the user by hand will have high entropy, which ideally they should.”
Agreed. But my hypothesis is that the code can create a higher entropty seed from 32 integers than it can from just a single user provided integer. Correct? So the makers of the algorithm must be using the 32 user provided integers to construct a better seed that any single integer by itself. Isn’t that the motivation for the 32 integers? I can’t think of any other good reason why they would require a seed array of size 32 (or 33, I forget).
ps. Any thoughts about my question of whether Random_Number with a 64 bit output ever equals zero?
Sixty-four bit FP number stores 53 bits of mantissa, 11 bits of exponent and 1 sign bit. The numbers in [0,1) range use all 53 bits of mantissa, 10 bits of exponent (they are all less than 2**0, so roughly half of all available exponents). Sign bit is always 0. So, there are more or less (say, with factor of 2 error) 2**63 different possible values, or 9.2E+18. You are probably unable to generate so many values in any thinkable time. A sample program
program rand
implicit none
integer, parameter :: size=1E7
real(kind(1d0)), allocatable :: tab(:)
integer :: i, loop=1000
allocate(tab(size))
do i=1,loop
call random_number(tab)
end do
end program rand
which actually computes 1e10 random values (1000x1e7) takes 35 secs to complete on 3.2GHZ i5-10505 Intel CPU. Scaled by 9.2e8 that gives over 1000 years. And that would only give you a chance (not a certainty) to get all possible values including 0. I think getting any particular value is as unlikely as getting zero.
Ok, I agree you are right. In practice I will never see a zero. Nonethless, I’m going to put a line in my code to say, if I get a zero then print a warning and stop.
If uniform random variate x is exactly zero, and this value breaks the code, two alternatives to stopping the code are setting x to tiny(x), which is 1.17549435E-38 and 2.2250738585072014E-308 for single and double precision, or requesting another variate.
Gfortran uses the xoshiro256** PRNG algorithm which is described here: Xorshift - Wikipedia. It has a very long cycle, which is why the internal state requires 8 integers (with default KIND int32). It has the remarkable property that it is possible to skip forward an arbitrary number of steps. This allows multiple PRNG streams to be generated in independent threads and/or on separate compute nodes that are guaranteed to not overlap with each other.
The fortran requirements are that the real numbers that are generated are in the range 0.0<=x<1.0, so a return value of exactly 0.0 is allowed. If your algorithm cannot use those numbers, then you need to test and generate a new value in its place. I think the way the numbers are generated, there are only 2**24 (=16 million) possible values returned for real32 arguments, so if you generate a few million numbers, you should expect to see 0.0 every so often.
edit: gfortran seed is 8 integers, not 33 integers.
Thanks Ron. It sounds like you know this well! So the algorithm is not a simple markov chain, as I first guessed. You re saying that it is like a markov chain where each new value depends on the previous 32 values? (I’m not sure what you would refer to this as in Markov lanaugage.)
I easily generate a million random numbers in my code, which goes on for days. I gave it a line that tells it to stop and print a warning if it produces a zero. I’ll be checking over coming days to see if it trips this wire.
The fortran standard does not specify the PRNG algorithm, so every compiler is free to choose its own. One of the above replies mentions NAG, and they use a different algorithm than gfortran does. The gfortran documentation (https://gcc.gnu.org/onlinedocs/gfortran/RANDOM_005fNUMBER.html) says that the period is 2^256-1 long. Since there are only 2^24 values generated (I think this is correct for real32, but I’m not certain), then each possible value will be generated about 2^232 times in that long period, and then it repeats. The wikipedia link above should be consulted for the details of the algorithm, but I think all the bits of all the 8 integers are modified when the internal state is changed.
edit: gfortran seed is 8 integers, not 33 integers.
those 32 integers (I guess your 33 is a typo) are probably counted in 8-bit values, as call random_seed(size = n)
returns n==8 (gfortran 13.2.0 x86_64-pc-linux-gnu)
do you know whether (and how) can one initialize such independent, non-overlapping PRNG streams? I guess it should be done with proper seeding, but how?
This was a mistake. You are right, gfortran seed is 8 integers, or 256 bits. I edited my original post to correct this mistake.
BTW, the NAG compiler seed is 630 integers. Each compiler is free to choose its own algorithm, which typically requires different numbers of bits of internal state storage.
In a parallel environment, this is done with the IMAGE_DISTINCT argument to the RANDOM_INIT() intrinsic as described here (RANDOM_INIT (The GNU Fortran Compiler)) for gfortran. This uses the jump algorithm mentioned in the wikipedia page to skip forward/backward. This is an inexpensive operation. In a single thread environment, I don’t think the programmer can jump easily. That might be a handy thing to be able to do.
I thought the wikipedia page had a good reference for the jump operation, but I’m not seeing it now. I did a google search, and I also don’t see a good reference for this in the first page. I know it involves computing a polynomial with modular arithmetic that depends on the number of jump steps, and then operating on the 256 bits of state, but I don’t have the original reference handy where I read this.
Well, the doc page only says that if that argument is .true., the seed is set to a processor-dependent value that is distinct from the seed set by a call to RANDOM_INIT in another image. Nothing about non-overlapping. After all, to get non-overlapping sequences the user would have to supply the length of the sequence.
BTW, the second part of the description of IMAGE_DISTINCT parameter:
If it is .false., the seed is set to a value that does depend which image called RANDOM_INIT.
seems to be wrong/incomplete. I’d rather guess a value that does not dependon which image called…*
This generator has a period of 2^{256} - 1, and when using multiple threads
up to 2^{128} threads can each generate 2^{128} random numbers before
any aliasing occurs.
It jumps ahead 2^128 steps for each image, so that gives 2^128 starting points within the long period, and each can do up to 2^128 steps without overlap (which they call aliasing).
I think this is right, the documentation seems to be incorrect.