RAND_GEN
Generates a rank-1 array of random numbers. The output array entries are positive and less than 1 in value.
Required Argument
X — Rank-1 array containing the random numbers. (Output)
Optional Arguments
IRND = IRND (Output)
Rank-1 integer array. These integers are the internal results of the Generalized Feedback Shift Register (GFSR) algorithm. The values are scaled to yield the floating-point array X. The output array entries are between 1 and 231‑ 1 in value.
ISTATE_IN = ISTATE_IN (Input)
Rank-1 integer array of size 3p + 2, where p = 521, that defines the ensuing state of the GFSR generator. It is used to reset the internal tables to a previously defined state. It is the result of a previous use of the “ISTATE_OUT=” optional argument.
ISTATE_OUT = ISTATE_OUT (Output)
Rank-1 integer array of size 3p + 2 that describes the current state of the GFSR generator. It is normally used to later reset the internal tables to the state defined following a return from the GFSR generator. It is the result of a use of the generator without a user initialization, or it is the result of a previous use of the optional argument “ISTATE_IN=” followed by updates to the internal tables from newly generated values. Example 2 illustrates use of ISTATE_IN and ISTATE_OUT for setting and then resetting RAND_GEN so that the sequence of integers, irnd, is repeatable.
IOPT = IOPT(:) (Input[/Output])
Derived type array with the same precision as the array x; used for passing optional data to RAND_GEN. The options are as follows:
Packaged Options for RAND_GEN |
|---|
Option Prefix = ? | Option Name | Option Value |
s_, d_ | Rand_gen_generator_seed | 1 |
s_, d_ | Rand_gen_LCM_modulus | 2 |
s_, d_ | Rand_gen_use_Fushimi_start | 3 |
IOPT(IO) = ?_options(?_rand_gen_generator_seed, ?_dummy)
Sets the initial values for the GFSR. The present value of the seed, obtained by default from the real-time clock as described below, swaps places with
iopt(IO + 1)%idummy. If the seed is set before any current usage of RAND_GEN, the exchanged value will be zero.
IOPT(IO) = ?_options(?_rand_gen_LCM_modulus, ?_dummy)
IOPT(IO+1) = ?_options(modulus, ?_dummy)
Sets the initial values for the GFSR. The present value of the LCM, with default value k = 16807, swaps places with iopt(IO+1)%idummy.
IOPT(IO) = ?_options(?_rand_gen_use_Fushimi_start, ?_dummy)
Starts the GFSR sequence as suggested by Fushimi (1990). The default starting sequence is with the LCM recurrence described below.
FORTRAN 90 Interface
Generic: CALL RAND_GEN (X [, …])
Specific: The specific interface names are S_RAND_GEN and D_RAND_GEN.
Description
This GFSR algorithm is based on the recurrence
where a ⊕ b is the exclusive OR operation on two integers a and b. This operation is performed until SIZE(x) numbers have been generated. The subscripts in the recurrence formula are computed modulo 3p. These numbers are converted to floating point by effectively multiplying the positive integer quantity
by a scale factor slightly smaller than 1./(huge(1)). The values p = 521 and q = 32 yield a sequence with a period approximately
The default initial values for the sequence of integers {xt} are created by a congruential generator starting with an odd integer seed
obtained by the Fortran 90 real-time clock routine:
CALL SYSTEM_CLOCK(COUNT=count,CLOCK_RATE=CLRATE)
An error condition is noted if the value of CLRATE=0. This indicates that the processor does not have a functioning real-time clock. In this exceptional case a starting seed must be provided by the user with the optional argument “iopt=” and option number ?_rand_generator_seed. The value v is the current clock for this day, in milliseconds. This value is obtained using the date routine:
CALL DATE_AND_TIME(VALUES=values)
and converting values(5:8) to milliseconds.
The LCM generator initializes the sequence {xt} using the following recurrence:
The default value of k = 16807. Using the optional argument “iopt=” and the packaged option number ?_rand_gen_LCM_modulus, k can be given an alternate value. The option number ?_rand_gen_generator_seed can be used to set the initial value of m instead of using the asynchronous value given by the system clock. This is illustrated in Example 2. If the default choice of m results in an unsatisfactory starting sequence or it is necessary to duplicate the sequence, then it is recommended that users set the initial seed value to one of their own choosing. Resetting the seed complicates the usage of the routine.
This software is based on Fushimi (1990), who gives a more elaborate starting sequence for the {xt} . The starting sequence suggested by Fushimi can be used with the option number ?_rand_gen_use_Fushimi_start. Fushimi’s starting process is more expensive than the default method, and it is equivalent to starting in another place of the sequence with period 2p.
Fatal and Terminal Error Messages
See the messages.gls file for error messages for RAND_GEN. These error messages are numbered 521‑528; 541‑548.
Examples
Example 1: Running Mean and Variance
An array of random numbers is obtained. The sample mean and variance are computed. These
values are compared with the same quantities computed using a stable method for the running means and variances, sequentially moving through the data. Details about the running mean and variance are found in Henrici (1982, pp. 21‑23).
use rand_gen_int
implicit none
! This is Example 1 for RAND_GEN.
integer i
integer, parameter :: n=1000
real(kind(1e0)), parameter :: one=1e0, zero=0e0
real(kind(1e0)) x(n), mean_1(0:n), mean_2(0:n), s_1(0:n), s_2(0:n)
! Obtain random numbers.
call rand_gen(x)
! Calculate each partial mean.
do i=1,n
mean_1(i) = sum(x(1:i))/i
end do
! Calculate each partial variance.
do i=1,n
s_1(i)=sum((x(1:i)-mean_1(i))**2)/i
end do
mean_2(0)=zero
mean_2(1)=x(1)
s_2(0:1)=zero
! Alternately calculate each running mean and variance,
! handling the random numbers once.
do i=2,n
mean_2(i)=((i-1)*mean_2(i-1)+x(i))/i
s_2(i) = (i-1)*s_2(i-1)/i+(mean_2(i)-x(i))**2/(i-1)
end do
! Check that the two sets of means and variances agree.
if (maxval(abs(mean_1(1:)-mean_2(1:))/mean_1(1:)) <= &
sqrt(epsilon(one))) then
if (maxval(abs(s_1(2:)-s_2(2:))/s_1(2:)) <= &
sqrt(epsilon(one))) then
write (*,*) 'Example 1 for RAND_GEN is correct.'
end if
end if
end
Output
Example 1 for RAND_GEN is correct.
Example 2: Seeding, Using, and Restoring the Generator
use rand_gen_int
implicit none
! This is Example 2 for RAND_GEN.
integer i
integer, parameter :: n=34, p=521
real(kind(1e0)), parameter :: one=1.0e0, zero=0.0e0
integer irndi(n), i_out(3*p+2), hidden_message(n)
type(s_options) :: iopti(2)=s_options(0,zero)
character*34 message, returned_message
! This is the message to be hidden.
message = 'SAVE YOURSELF. WE ARE DISCOVERED!'
! Start the generator with a known seed.
iopti(1) = s_options(s_rand_gen_generator_seed,zero)
iopti(2) = s_options(123,zero)
call rand_gen(x, iopt=iopti)
! Save the state of the generator.
call rand_gen(x, istate_out=i_out)
! Get random integers.
call rand_gen(y, irnd=irndi)
! Hide text using collating sequence subtracted from integers.
do i=1, n
hidden_message(i) = irndi(i) - ichar(message(i:i))
end do
! Reset generator to previous state and generate the previous
! random integers.
call rand_gen(x, irnd=irndi, istate_in=i_out)
! Subtract hidden text from integers and convert to character.
do i=1, n
returned_message(i:i) = char(irndi(i) - hidden_message(i))
end do
! Check the results.
if (returned_message == message) then
write (*,*) 'Example 2 for RAND_GEN is correct.'
end if
end
Output
Example 2 for RAND_GEN is correct.
Example 3: Generating Strategy with a Histogram
We generate random integers but with the frequency as in a histogram with nbins slots. The generator is initially used a large number of times to demonstrate that it is making choices with the same shape as the histogram. This is not required to generate samples. The program next generates a summary set of integers according to the histogram. These are not repeatable and are representative of the histogram in the sense of looking at 20 integers during generation of a large number of samples.
use rand_gen_int
use show_int
implicit none
! This is Example 3 for RAND_GEN.
integer i, i_bin, i_map, i_left, i_right
integer, parameter :: n_work=1000
integer, parameter :: n_bins=10
integer, parameter :: scale=1000
integer, parameter :: total_counts=100
integer, parameter :: n_samples=total_counts*scale
integer, dimension(n_bins) :: histogram= &
(/4, 6, 8, 14, 20, 17, 12, 9, 7, 3 /)
integer, dimension(n_work) :: working=0
integer, dimension(n_bins) :: distribution=0
integer break_points(0:n_bins)
real(kind(1e0)) rn(n_samples)
real(kind(1e0)), parameter :: tolerance=0.005
integer, parameter :: n_samples_20=20
integer rand_num_20(n_samples_20)
real(kind(1e0)) rn_20(n_samples_20)
! Compute the normalized cumulative distribution.
break_points(0)=0
do i=1,n_bins
break_points(i)=break_points(i-1)+histogram(i)
end do
break_points=break_points*n_work/total_counts
! Obtain uniform random numbers.
call rand_gen(rn)
! Set up the secondary mapping array.
do i_bin=1,n_bins
i_left=break_points(i_bin-1)+1
i_right=break_points(i_bin)
do i=i_left, i_right
working(i)=i_bin
end do
end do
! Map the random numbers into the 'distribution' array.
! This is made approximately proportional to the histogram.
do i=1,n_samples
i_map=nint(rn(i)*(n_work-1)+1)
distribution(working(i_map))= &
distribution(working(i_map))+1
end do
! Check the agreement between the distribution of the
! generated random numbers and the original histogram.
write (*, '(A)', advance='no') 'Original: '
write (*, '(10I6)') histogram*scale
write (*, '(A)', advance='no') 'Generated:'
write (*, '(10I6)') distribution
if (maxval(abs(histogram(1:)*scale-distribution(1:))) &
<= tolerance*n_samples) then
write(*, '(A/)') 'Example 3 for RAND_GEN is correct.'
end if
! Generate 20 integers in 1, 10 according to the distribution
! induced by the histogram.
call rand_gen(rn_20)
! Map from the uniform distribution to the induced distribution.
do i=1,n_samples_20
i_map=nint(rn_20(i)*(n_work-1)+1)
rand_num_20(i)=working(i_map)
end do
call show(rand_num_20,&
'Twenty integers generated according to the histogram:')
end
Output
Example 3 for RAND_GEN is correct.
