//======================================================= file = unif.c =====
//=  Random number generator for uniform(0.0, 1.0)                          =
//===========================================================================
//=  Notes:                                                                 =
//=    1) This program is a "driver program" for the function rand_val().   =
//=       See the header of rand_val() for a full description.              =
//=    2) The function rand_val() is better than a z = rand() / RAND_MAX    =
//=       approach which is limited to 32767 different values.              =
//=-------------------------------------------------------------------------=
//= Example output:                                                         =
//=                                                                         =
//=   -------------------------------------------------- unif.c -----       =
//=   0.000008                                                              =
//=   0.131538                                                              =
//=   0.755605                                                              =
//=                                                                         =
//=   <SNIP SNIP>                                                           =
//=                                                                         =
//=   0.430814                                                              =
//=   0.691408                                                              =
//=   0.485973                                                              =
//=   ---------------------------------------------------------------       =
//=-------------------------------------------------------------------------=
//=  Build: gcc unif.c, bcc32 unif.c, cl unif.c                             =
//=-------------------------------------------------------------------------=
//=  Execute: unif                                                          =
//=-------------------------------------------------------------------------=
//=  Author: Kenneth J. Christensen                                         =
//=          University of South Florida                                    =
//=          WWW: http://www.csee.usf.edu/~christen                         =
//=          Email: christen@csee.usf.edu                                   =
//=-------------------------------------------------------------------------=
//=  History: KJC (01/21/99) - Genesis                                      =
//===========================================================================

//----- Include files -------------------------------------------------------
#include <stdio.h>          // Needed for printf()

//----- Function prototypes -------------------------------------------------
double rand_val(void);      // Function to generate uniform(0, 1)

//===========================================================================
//=  Main program                                                           =
//===========================================================================
void main(void)
{
  long i;      // Loop counter

  // Output a banner
  printf("-------------------------------------------------- unif.c -----\n");

  // Output 10000 random values between 0.0 and 1.0
  for (i=0; i<10000; i++)
    printf("%f \n", rand_val());

  // Output a tailer
  printf("---------------------------------------------------------------\n");
}

//=========================================================================
//= Multiplicative LCG for generating uniform(0.0, 1.0) random numbers    =
//=   - x_n = 7^5*x_(n-1)mod(2^31 - 1)                                    =
//=   - With x seeded to 1 the 10000th x value should be 1043618065       =
//=   - From R. Jain, "The Art of Computer Systems Performance Analysis," =
//=     John Wiley & Sons, 1991. (Page 443, Figure 26.2)                  =
//=========================================================================
double rand_val(void)
{
  const long  a =      16807;  // Multiplier
  const long  m = 2147483647;  // Modulus
  const long  q =     127773;  // m div a
  const long  r =       2836;  // m mod a
  static long x = 1;           // Random int value (seed is set to 1)
  long        x_div_q;         // x divided by q
  long        x_mod_q;         // x modulo q
  long        x_new;           // New x value

  // RNG using integer arithmetic
  x_div_q = x / q;
  x_mod_q = x % q;
  x_new = (a * x_mod_q) - (r * x_div_q);
  if (x_new > 0)
    x = x_new;
  else
    x = x_new + m;

  // Return a random value between 0.0 and 1.0
  return((double) x / m);
}