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 2 3 1− 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 {x t} 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 {x t} 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 2 p.

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.

.p>.F90CH5.DOC!EX1_RAND_GEN;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

Output

Example 1 for RAND_GEN is correct.

Additional Examples

.p>.F90CH5.DOC!EX2_RAND_GEN;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)

real(kind(1e0)) x(n), y(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.

.p>.F90CH5.DOC!EX3_RAND_GEN;Example 3: Generating Strategy with a Histogram

We generate random integers but with the frequency as in a histogram with n bins 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

! 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.

.p>.F90CH5.DOC!EX4_RAND_GEN;Example 4: Generating with a Cosine Distribution

We generate random numbers based on the continuous distribution function

Using the cumulative

we generate the samples by obtaining uniform samples u, 0 < u < 1 and solve the equation

These are evaluated in vector form, that is all entries at one time, using Newton's method:

An iteration counter forces the loop to terminate, but this is not often required although it is an important detail.

use rand_gen_int

use show_int

use Numerical_Libraries

IMPLICIT NONE

! This is Example 4 for RAND_GEN.

integer i, i_map, k

integer, parameter :: n_bins=36

integer, parameter :: offset=18

integer, parameter :: n_samples=10000

integer, parameter :: n_samples_30=30

integer, parameter :: COUNT=15

real(kind(1e0)) probabilities(n_bins)

real(kind(1e0)), dimension(n_bins) :: counts=0.0

real(kind(1e0)), dimension(n_samples) :: rn, x, f, fprime, dx

real(kind(1e0)), dimension(n_samples_30) :: rn_30, &

x_30, f_30, fprime_30, dx_30

real(kind(1e0)), parameter :: one=1e0, zero=0e0, half=0.5e0

real(kind(1e0)), parameter :: tolerance=0.01

real(kind(1e0)) two_pi, omega

! Initialize values of 'two_pi' and 'omega'.

two_pi=2.0*const((/'pi'/))

omega=two_pi/n_bins

! Compute the probabilities for each bin according to

! the probability density (cos(x)+1)/(2*pi), -pi<x<pi.

do i=1,n_bins

probabilities(i)=(sin(omega*(i-offset)) &

-sin(omega*(i-offset-1))+omega)/two_pi

end do

! Obtain uniform random numbers in (0,1).

call rand_gen(rn)

! Use Newton's method to solve the nonlinear equation:

! accumulated_distribution_function - random_number = 0.

x=zero; k=0

solve_equation: do

f=(sin(x)+x)/two_pi+half-rn

fprime=(one+cos(x))/two_pi

dx=f/fprime

x=x-dx; k=k+1

if (maxval(abs(dx)) <= sqrt(epsilon(one)) &

.or. k > COUNT) exit solve_equation

end do solve_equation

! Map the random numbers 'x' array into the 'counts' array.

do i=1,n_samples

i_map=int(x(i)/omega+offset)+1

counts(i_map)=counts(i_map)+one

end do

! Normalize the counts array.

counts=counts/n_samples

! Check that the generated random numbers are indeed

! based on the original distribution.

if (maxval(abs(counts(1:)-probabilities(1:))) &

<= tolerance) then

write (*,'(a/)') 'Example 4 for RAND_GEN is correct.'

end if

! Generate 30 random numbers in (-pi,pi) according to

! the probability density (cos(x)+1)/(2*pi), -pi<x<pi.

call rand_gen(rn_30)

x_30=0.0; k=0

solve_equation_30: do

f_30=(sin(x_30)+x_30)/two_pi+half-rn_30

fprime_30=(one+cos(x_30))/two_pi

dx_30=f_30/fprime_30

x_30=x_30-dx_30

if (maxval(abs(dx_30)) <= sqrt(epsilon(one))&

.or. k > COUNT) exit solve_equation_30

end do solve_equation_30

write(*,'(A)') 'Thirty random numbers generated ', &

'according to the probability density ',&

'pdf(x)=(cos(x)+1)/(2*pi), -pi<x<pi:'

call show(x_30)

end

Output

Example 4 for RAND_GEN is correct.

Visual Numerics, Inc.
Visual Numerics - Developers of IMSL and PV-WAVE
http://www.vni.com/
PHONE: 713.784.3131
FAX:713.781.9260