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