function rand (r)
c
c this pseudo-random number generator is portable amoung a wide
c variety of computers. rand(r) undoubtedly is not as good as many
c readily available installation dependent versions, and so this
c routine is not recommended for widespread usage. its redeeming
c feature is that the exact same random numbers (to within final round-
c off error) can be generated from machine to machine. thus, programs
c that make use of random numbers can be easily transported to and
c checked in a new environment.
c the random numbers are generated by the linear congruential
c method described, e.g., by knuth in seminumerical methods (p.9),
c addison-wesley, 1969. given the i-th number of a pseudo-random
c sequence, the i+1 -st number is generated from
c x(i+1) = (a*x(i) + c) mod m,
c where here m = 2**22 = 4194304, c = 1731 and several suitable values
c of the multiplier a are discussed below. both the multiplier a and
c random number x are represented in double precision as two 11-bit
c words. the constants are chosen so that the period is the maximum
c possible, 4194304.
c in order that the same numbers be generated from machine to
c machine, it is necessary that 23-bit integers be reducible modulo
c 2**11 exactly, that 23-bit integers be added exactly, and that 11-bit
c integers be multiplied exactly. furthermore, if the restart option
c is used (where r is between 0 and 1), then the product r*2**22 =
c r*4194304 must be correct to the nearest integer.
c the first four random numbers should be .0004127026,
c .6750836372, .1614754200, and .9086198807. the tenth random number
c is .5527787209, and the hundredth is .3600893021 . the thousandth
c number should be .2176990509 .
c in order to generate several effectively independent sequences
c with the same generator, it is necessary to know the random number
c for several widely spaced calls. the i-th random number times 2**22,
c where i=k*p/8 and p is the period of the sequence (p = 2**22), is
c still of the form l*p/8. in particular we find the i-th random
c number multiplied by 2**22 is given by
c i = 0 1*p/8 2*p/8 3*p/8 4*p/8 5*p/8 6*p/8 7*p/8 8*p/8
c rand= 0 5*p/8 2*p/8 7*p/8 4*p/8 1*p/8 6*p/8 3*p/8 0
c thus the 4*p/8 = 2097152 random number is 2097152/2**22.
c several multipliers have been subjected to the spectral test
c (see knuth, p. 82). four suitable multipliers roughly in order of
c goodness according to the spectral test are
c 3146757 = 1536*2048 + 1029 = 2**21 + 2**20 + 2**10 + 5
c 2098181 = 1024*2048 + 1029 = 2**21 + 2**10 + 5
c 3146245 = 1536*2048 + 517 = 2**21 + 2**20 + 2**9 + 5
c 2776669 = 1355*2048 + 1629 = 5**9 + 7**7 + 1
c
c in the table below log10(nu(i)) gives roughly the number of
c random decimal digits in the random numbers considered i at a time.
c c is the primary measure of goodness. in both cases bigger is better.
c
c log10 nu(i) c(i)
c a i=2 i=3 i=4 i=5 i=2 i=3 i=4 i=5
c
c 3146757 3.3 2.0 1.6 1.3 3.1 1.3 4.6 2.6
c 2098181 3.3 2.0 1.6 1.2 3.2 1.3 4.6 1.7
c 3146245 3.3 2.2 1.5 1.1 3.2 4.2 1.1 0.4
c 2776669 3.3 2.1 1.6 1.3 2.5 2.0 1.9 2.6
c best
c possible 3.3 2.3 1.7 1.4 3.6 5.9 9.7 14.9
c
c input argument --
c r if r=0., the next random number of the sequence is generated.
c if r.lt.0., the last generated number will be returned for
c possible use in a restart procedure.
c if r.gt.0., the sequence of random numbers will start with the
c seed r mod 1. this seed is also returned as the value of
c rand provided the arithmetic is done exactly.
c
c output value --
c rand a pseudo-random number between 0. and 1.
c
c ia1 and ia0 are the hi and lo parts of a. ia1ma0 = ia1 - ia0.
data ia1, ia0, ia1ma0 /1536, 1029, 507/
data ic /1731/
data ix1, ix0 /0, 0/
c
if (r.lt.0.) go to 10
if (r.gt.0.) go to 20
c
c a*x = 2**22*ia1*ix1 + 2**11*(ia1*ix1 + (ia1-ia0)*(ix0-ix1)
c + ia0*ix0) + ia0*ix0
c
iy0 = ia0*ix0
iy1 = ia1*ix1 + ia1ma0*(ix0-ix1) + iy0
iy0 = iy0 + ic
ix0 = mod (iy0, 2048)
iy1 = iy1 + (iy0-ix0)/2048
ix1 = mod (iy1, 2048)
c
10 rand = ix1*2048 + ix0
rand = rand / 4194304.
return
c
20 ix1 = amod(r,1.)*4194304. + 0.5
ix0 = mod (ix1, 2048)
ix1 = (ix1-ix0)/2048
go to 10
c
end