/*

Random number package with simple test program.
James L. Blue
Leader, Mathematical Modeling Group
Mathematics and Computational Sciences Division
Information Technology Laboratory
National Institute of Standards and Technology
Gaithersburg, MD 20899
blue@nist.gov

This package provides a portable method for choosing N integers from
the first M integers, with no duplications. There are four routines.

1. A simple test program
2. An initialization program, init. Init should be called with a different
   argument for each call; the simplest method is to call it with argument
   1 on the first call, 2 on the second, and so on.
3. A program for choosing N of the first M integers.
4. The underlying uniform random number generator.

*/

#include <stdio.h>

#define mbig  2147483647

extern void init(int iseed);
extern void nofm(int n, int m, int *irnd);
extern float uni(int jd);

/* test program
     Correct results are

uni    0.3564443
uni    0.3584030
choosing 3 out of 20 with seed 12345
       12345            1            9           13
nofm requires n <= m; have n 21 and m 20
*/

main()
{
    int irnd[100];

    init(1);
    printf("uni %12.7f\n", uni(0));
    printf("uni %12.7f\n", uni(0));

    printf("choosing 3 out of 20 with seed 12345\n");
    init(12345);
    nofm(3, 20, irnd);
    printf("%12d %12d %12d %12d\n", 12345, irnd[0], irnd[1], irnd[2]);

    nofm(21, 20, irnd);

    exit(0);
}

/* init:
   Initialize the random number generator, uni. As a seed,
   uni requires an odd number in the range 1 to 2147483647, inclusive.
   Allow the caller to use the positive integers in order as initializers.

   Also, since very small seeds start off poorly, do a few random
   numbers to get over the transient, then use the last for a new seed.
*/
void
init(int iseed)
{
    int i, j;
    float x;

    i = iseed;
    if (i < 0) i = -i;
    if (i > mbig) i = mbig;
    if ((i % 2) == 0)
        i = mbig - i;
    x = uni(i);
    for (j = 0; j < 10; j++)
        x = uni(0);
    i = x * mbig;
    x = uni(i);
}

/*
 nofm:
   Select n items randomly from a total of m items.
   Return an array of the ones chosen.
   Method:
     Algorithm S, p. 137 of
     D.E. Knuth,
     "The Art of Computer Programming,
      Volume 2: Seminumerical Algorithms",
     Addison-Wesley, 1969.

  In the absence of roundoff, the algorithm always chooses n out of m.
  With roundoff, it might come up one short. This happens rarely, but
  if it does, just start over. (For example, in 10,000,000 different
  choosings of 19 out of 20, there were 3 startovers.)

  Sample values returned: If choosing 3 out of 20 after calling init(12345),
  the items chosen are 1, 9, and 13.
*/

void
nofm(int n, int m, int *irnd)
{
    int ihave, it, nn;
    float x;

    if (n > m)
    {   printf("nofm requires n <= m; have n %d and m %d\n", n, m);
        exit(4);
    }
    nn = m;
    ihave = 0;
    do
    {
        for (it = 0; it < nn; it++)
        {   x = uni(0);
            if ((nn - it) * x < (n - ihave))
            {   irnd[ihave++] = it + 1;
                if (ihave == n) break;
            }
        }
        if (ihave != n)
            printf(" start over: got %d, wanted %d\n", ihave, n);
    }
    while (ihave != n);
}


/* uni
***date written   810915
***revision date  940629 (JLB)
***authors
   Blue, James L., Applied and Computational Mathematics Division, NIST
   Kahaner, David K., Applied and Computational Mathematics Division, NIST
   Marsaglia, George, Florida State University

***purpose  This routine generates quasi uniform real random numbers in the
        range 0 to 1, and can be used on any computer which allows
        integers at least as large as 2147483647 (in practice, any
        computer with 32 or more bits).

   use
       first time....
           z = uni(jd)
             here jd is any  n o n - z e r o  integer.
             this causes initialization of the program
             and the first random number to be returned as z.
       subsequent times...
           z = uni(0)
             causes the next random number to be returned as z.


   machine dependencies...
      mdig = 32; the machine must have at least 32 binary digits available
          for representing integers, including the sign bit.
      mbig = 2147483647; the machine must allow positive and negative
          integers of this size.

   remark
    This program gives repeatable results on different computers, within
    roundoff, since internal calculations are in integers.

***reference  Marsaglia G., "a Current View of Random Number Generators, "
          Proceedings Computer Science and Statistics: 16th
          Symposium on the Interface, Elsevier, Amsterdam, 1985.
***routines called: none
***end prologue  uni
*/

float
uni(int jd)
{
    int k;
    static int i, j, ihist[17], notyet = 1;

    if (jd != 0)
    {   int j0, j1, jseed, k0, k1, n, m;
    /*
      Fill up the history array, taking care not to overflow.

      mbig is the largest positive integer, 2**(mdig-1) - 1
      m = 2**(mdig/2)

      We use the simple congruential generator x(n) = [9069*x(n-1]mod(mm),
      with mm = (2**(mdig-1)), to fill up the history array with integers in
      the range 0 to mbig, inclusive.
      Note that mm = (m**2)/2 = mbig + 1.
      We use m to avoid overflow in the calculation, as follows:
      Any integer j, 0<=j<=mbig, can be expressed as j1*m + j0, with j1=j/m
      (discarding fractions) and j0=[j]mod m; then 0<=j0<m and 0<=j1<m.
      To multiply two numbers modulus mm, we do
      [j*k]mod(mm) = [j1*k1*m*m + (j0*k1+j1*k0)*m + j0*k0]mod(mm)
          The first term is 0 modulus mm, since m*m = 2*mm.
          We let j0*k0 = n1*m + n0, and combine n1*m with the 2nd term, giving
              [n1+j0*k1+j1*k0)*m + n0]mod(mm)
          The final term [n0]mod(mm) is just n0 = [j0*k0]mod(m)
          The remaining term is
              [n1+j0*k1+j1*k0)*m]mod(m**2/2)
          Write this as
              [a*m]mod(m**2/2)
          And expand a = (a/(m/2))*(m/2) + [a]mod(m/2). Then the first term is
              [a/(m/2))*(m/2)*m]mod(m**2/2) = 0
          and the second term is just
              [a]mod(m/2)*m
          since it is clearly already less than m**2/2.
    */
        jseed = jd;
        if (jseed < 0) jseed = -jseed;
        if (jseed > mbig) jseed = mbig;
        if ((jseed % 2) == 0) jseed--;
        m = 65536;
        k0 = 9069 % m;
        k1 = 9069/m;
        j0 = jseed % m;
        j1 = jseed/m;
        for (i = 0; i < 17; i++)
        {    n = j0*k0;
             j1 = (n/m+j0*k1+j1*k0) % (m/2);
             j0 = n % m;
             ihist[i] = j0+m*j1;
        }
        i =  4;
        j = 16;
        notyet = 0;
    }

    /*  Begin main loop here. */

    if (notyet)
    {   printf("first call to uni must have nonzero argument");
        exit(2);
    }
    k = ihist[i] - ihist[j];
    if (k < 0) k += mbig;
    ihist[j] = k;
    if (--i < 0) i = 16;
    if (--j < 0) j = 16;
    return (float)(k) / (float)(mbig);
}
/* end of program */