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