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