function random (t, n)
c
c this random number generator is portable amoung a wide variety of
c computers. it generates a random number between 0.0 and 1.0 accord-
c ing to the algorithm presented by bays and durham (toms, 2, 59,
c 1976). the motivation for using this scheme, which resembles the
c maclaren-marsaglia method, is to greatly increase the period of the
c random sequence. if the period of the basic generator (rand) is p,
c then the expected mean period of the sequence generated by random is
c given by new mean p = sqrt (pi*factorial(n)/(8*p)),
c where factorial(n) must be much greater than p in this asymptotic
c formula. generally, n should be 16 to maybe 32.
c
c input argument --
c n iabs(n) is the number of random numbers in an auxiliary table.
c note though that iabs(n)+1 is the number of items in array t.
c if n is positive and differs from its value in the previous
c invocation, then the table is initialized for the new value of
c n. if n is negative, iabs(n) is the number of items in an
c auxiliary table, but the tables are now assumed already to
c be initialized. this option enables the user to save the
c table t at the end of a long computer run and to restart with
c the same sequence. normally, random would be called at most
c once with negative n. subsequent invocations would have n
c positive and of the correct magnitude.
c
c input and output argument --
c t an array of iabs(n)+1 random numbers from a previous invocation
c of random. whenever n is positive and differs from the old
c n, the table is initialized. the first iabs(n) numbers are the
c table discussed in the reference, and the n+1 -st value is y.
c this array may be saved in order to restart a sequence.
c
c output value --
c random a random number between 0.0 and 1.0.
c
dimension t(1)
external rand
data nold, floatn / -1, -1.0 /
c
if (n.eq.nold) go to 20
c
nold = iabs(n)
floatn = nold
if (n.lt.0) dummy = rand (t(nold+1))
if (n.lt.0) go to 20
c
do 10 i=1,nold
t(i) = rand (0.)
10 continue
t(nold+1) = rand (0.)
c
20 j = t(nold+1)*floatn + 1.
t(nold+1) = t(j)
random = t(j)
t(j) = rand (0.)
c
return
end