c Random number package with simple test program.
c
c James L. Blue
c Leader, Mathematical Modeling Group
c Mathematics and Computational Sciences Division
c Information Technology Laboratory
c National Institute of Standards and Technology
c Gaithersburg, MD 20899
c blue@nist.gov
c
c
c This package provides a portable method for choosing N integers from
c the first M integers, with no duplications. There are four routines.
c
c 1. A simple test program
c 2. An initialization program, init. Init should be called with a different
c    argument for each call; the simplest method is to call it with argument
c    1 on the first call, 2 on the second, and so on.
c 3. A program for choosing N of the first M integers.
c 4. The underlying uniform random number generator.
c
c
c test program
c      Correct results are
c
c   uni   0.3564443
c   uni   0.3584030
c  choosing 3 out of 20 with seed 12345
c         12345           1           9          13
c  nofm requires n .le. m; have n and m          21          20
c
        program uni4
        integer irnd(100)
c
c
        call init(1)
        print *, ' uni ', uni(0)
        print *, ' uni ', uni(0)
c
        print *, 'choosing 3 out of 20 with seed 12345'
        call init(12345)
        call nofm(3, 20, irnd)
        print *, 12345, irnd(1), irnd(2), irnd(3)
c
        call nofm(21, 20, irnd)
c
        stop
        end
c
c
c init:
c   Initialize the random number generator, uni. As a seed,
c   uni requires an odd number in the range 1 to 2147483647, inclusive.
c   Allow the caller to use the positive integers in order as initializers.
c
c   Also, since very small seeds start off poorly, do a few random
c   numbers to get over the transient, then use the last for a new seed.
c
        subroutine init(iseed)
        data mbig / 2147483647 /
        i = iabs(iseed)
        i = min(i, mbig)
        if (mod(i, 2) .eq. 0) i = mbig - i
        x = uni(i)
        do 100 j = 1, 10
          x = uni(0)
100     continue
        i = x * mbig
        x = uni(i)
        return
        end
c
c nofm:
c   Select n items randomly from a total of m items.
c   Return an array of the ones chosen.
c   Method:
c     Algorithm S, p. 137 of
c     D.E. Knuth,
c     "The Art of Computer Programming,
c      Volume 2: Seminumerical Algorithms",
c     Addison-Wesley, 1969.
c
c  In the absence of roundoff, the algorithm always chooses n out of m.
c  With roundoff, it might come up one short. This happens rarely, but
c  if it does, just start over. (For example, in 10,000,000 different
c  choosings of 19 out of 20, there were 3 startovers.)
c
c  Sample values returned: If choosing 3 out of 20 after calling init(12345),
c  the items chosen are 1, 9, and 13.
c
      subroutine nofm(n, m, irnd)
        integer irnd(n)
c
      if (n .gt. m) then
          print *, 'nofm requires n .le. m; have n and m', n, m
          stop
      endif
      nn = m
      ihave = 0
 10   do 100 it = 0, nn-1
        x = uni(0)
        if ((nn - it) * x .lt. (n - ihave)) then
          ihave = ihave + 1
          irnd(ihave) = it + 1
          if (ihave .eq. n) goto 500
        endif
100   continue
500   continue
      if (ihave .ne. n) then
          print *, ' start over: got', ihave, ' wanted', n
          goto 10
      endif
      return
      end
c
c
c
      real function uni(jd)
c***begin prologue  uni
c***date written   810915
c***revision date  940629 (JLB)
c***authors
c   Blue, James L., NIST
c   Kahaner, David K., NIST
c   Marsaglia, George, Florida State University
c
c***purpose  This routine generates quasi uniform real random numbers in the
c            range 0 to 1, and can be used on any computer with which allows
c            integers at least as large as 2147483647 (in practice, any
c            computer with 32 or more bits).
c
c   use
c       first time....
c                   z = uni(jd)
c                     here jd is any  n o n - z e r o  integer.
c                     this causes initialization of the program
c                     and the first random number to be returned as z.
c       subsequent times...
c                   z = uni(0)
c                     causes the next random number to be returned as z.
c
c
c   machine dependencies...
c      mdig = 32; the machine must have at least 32 binary digits available
c              for representing integers, including the sign bit.
c      mbig = 2147483647; the machine must allow positive and negative
c              integers of this size.
c
c   remark
c        This program gives repeatable results on different computers, within
c        roundoff, since internal calculations are in integers.
c
c***reference  Marsaglia G., "a Current View of Random Number Generators, "
c              Proceedings Computer Science and Statistics: 16th
c              Symposium on the Interface, Elsevier, Amsterdam, 1985.
c***routines called: none
c***end prologue  uni
      integer ihist(17)
      integer i, j, j0, j1, jd, jseed, k0, k1, n, m, mbig, mdig
      logical notyet
c
      save i, j, ihist, mbig, notyet
c
c     data mdig, mbig, m / 16,      32767,   256 /
      data mdig, mbig, m / 32, 2147483647, 65536 /
      data notyet / .true. /
c
c***first executable statement  uni
      if (jd .eq. 0) go to 3
c
c  mbig is the largest positive integer, 2**(mdig-1) - 1
c  m = 2**(mdig/2)
c
c  We use the simple congruential generator x(n) = [9069*x(n-1)]mod(mm),
c  with mm = (2**(mdig-1)), to fill up the history array with integers in
c  the range 0 to mbig, inclusive.
c  Note that mm = (m**2)/2 = mbig + 1.
c  We use m to avoid overflow in the calculation, as follows:
c  Any integer j, 0<=j<=mbig, can be expressed as j1*m + j0, with j1=j/m
c  (discarding fractions) and j0=[j]mod m; then 0<=j0<m and 0<=j1<m.
c  To multiply two numbers modulus mm, we do
c  [j*k]mod(mm) = [j1*k1*m*m + (j0*k1+j1*k0)*m + j0*k0]mod(mm)
c      The first term is 0 modulus mm, since m*m = 2*mm.
c      We let j0*k0 = n1*m + n0, and combine n1*m with the 2nd term, giving
c              [(n1+j0*k1+j1*k0)*m + n0]mod(mm)
c      The final term [n0]mod(mm) is just n0 = [j0*k0]mod(m)
c      The remaining term is
c              [(n1+j0*k1+j1*k0)*m]mod(m**2/2)
c      Write this as
c              [a*m]mod(m**2/2)
c      And expand a = (a/(m/2))*(m/2) + [a]mod(m/2). Then the first term is
c              [(a/(m/2))*(m/2)*m]mod(m**2/2) = 0
c      and the second term is just
c              [a]mod(m/2)*m
c      since it is clearly already less than m**2/2.
c
      jseed = min0(iabs(jd), mbig)
      if ( mod(jseed, 2).eq.0 ) jseed=jseed-1
      k0 = mod(9069, m)
      k1 = 9069/m
      j0 = mod(jseed, m)
      j1 = jseed/m
      do 2 i = 1, 17
        n = j0*k0
        j1 = mod(n/m+j0*k1+j1*k0, m/2)
        j0 = mod(n, m)
        ihist(i) = j0+m*j1
    2 continue
      i = 5
      j = 17
      notyet = .false.
c
c  Begin main loop here.
c
    3 if (notyet) then
        print *, 'first call to uni must have nonzero argument'
        stop
      endif
      k = ihist(i) - ihist(j)
      if (k .lt. 0) k = k + mbig
      ihist(j) = k
      i = i-1
      if (i .eq. 0) i = 17
      j = j-1
      if (j .eq. 0) j = 17
      uni = float(k) / float(mbig)
      return
      end
c end of program