/*
 * -------------------------------------------------------------------------- 
 * This is a Turbo Pascal unit for random number generation.  The use of this 
 * unit is recommended as a replacement for the Turbo Pascal function Random  
 * and procedure Randomize - particularly in simulation applications where    
 * the statistical 'goodness' of the random number generator is important.    
 *                                                                            
 * The generator upon which this unit is based is a so-called 'Lehmer random  
 * number generator' which returns a pseudo-random real number uniformly      
 * distributed between 0.0 and 1.0.  The period is (m - 1) where              
 * m = 2,147,483,647 and the smallest and largest possible values are (1 / m) 
 * and 1 - (1 / m) respectively.  For more details see -                      
 *                                                                            
 *        "Random Number Generators: Good Ones Are Hard To Find"              
 *                   Steve Park & Keith Miller                                
 *           Communications of the ACM, October, 1988                         
 *                                                                            
 * This unit actually supplies 256 different "streams" of random numbers.     
 * The procedure Select_Stream can be used to switch between streams.         
 *                                                                            
 * Note that a math co-processor is assumed.  To use this unit without a math 
 * co-processor, change all instances of 'double' to 'real' and change the    
 * compiler directive from $N+ to $N-.                                        
 *                                                                            
 * Note that as of 1-08-91 the multiplier used in this unit has been changed  
 * from the previous "minimal standard" 16807 to a new value of 48271.  To    
 * use this unit in its old (16807) form change the constants 'a', 'check'    
 * and 'a256' as indicated in the implementation section comments.            
 *                                                                            
 * -------------------------------------------------------------------------- 
 * Author          : Steve Park                                               
 * Language        : Turbo Pascal, 5.0                                        
 * Latest Revision : 2-19-92                                                  
 * Reference       : Lecture Notes on Simulation, by Steve Park               
 * -------------------------------------------------------------------------- 
 * C Translation By : Tracey A. Beauchat
 * -------------------------------------------------------------------------- 
 */


#include <stdlib.h>

#include "rngs.h"

/* Private constants */

#define    m      2147483647 /* DO NOT CHANGE THIS VALUE */
#define    max    (m - 1)
#define    a      48271      /*change to 16807 to use the "minimal standard" */
#define    check  399268537  /*change to 1043618065 if a = 16807             */
#define    a256   22925      /*jump multiplier, change to 36563 if a = 16807 */

#define    MaxStreamIndex 255

/* Private global variables */

static int s;                               /* stream index   */
static long int seed[MaxStreamIndex+1];     /* array of seeds */

/* Callable functions */

void initialize(void)
{
   Select_Stream(0);
   Randomize(123456789);
} /* initialize */

double Random(void)
{
   static long int q = m / a;
   static long int r = m % a;

   long int t;

   t  = a * (seed[s] % q) - r * (seed[s] / q);
   seed[s] = (t > 0) ? t : t + m;
   t = seed[s] - 1;
   return((double)t / (double)m);
}  /* Random */


/* This is a private procedure used to */
/* "plant" a sequence of 255 seeds     */
/* seperated one from the next by      */
/* 8,367,782 steps.                    */
static void Plant_Seeds(void)
{                 
   static long int quotient  = m / a256;
   static long int remainder = m % a256;
   
   int      j;
   long int t;

    for (j = 1; j <= MaxStreamIndex; j++) {
       t = a256 * (seed[j - 1] % quotient)
                  - remainder * (seed[j - 1] / quotient);
       seed[j] = (t > 0) ? t : t + m;
    }
} /* Plant_Seeds */


void

Select_Stream(int stream)
{
   if ((s < 0) || (s > MaxStreamIndex)) {
     fprintf(stderr, "Stream index is out of range (0..%d)\n", MaxStreamIndex);
      exit(0);
   }
   s = stream;
} /* Select_Stream */


void

Randomize(long int x)
{
   int  ok;
   char buf[25];

   if (x == m)                     /* error trap on this value      */
      x = 0;

   if (x > 0)                      /* do nothing                    */
      x = x;
   else if (x == 0) {              /* request user input            */
      do {
          printf("\nenter a positive integer seed (9 digits or less) >> ");
          gets(buf);
	  x = atoi(buf);
          ok = ((0 < x) && (x < m));
          if (!ok) 
            printf("%cinput out of range ... try again\n", 0x07);
      } while (!ok);

      printf("\n");
   }   

/* Remove this for now: */

/*   else { */                         /* seed from the clock           */
/*      do {
        x = Current_Time();
      } while ((x < 0) || (x > m));
   }
*/

   seed[0] = x;                       /* this is the seed for stream 0 */
   Plant_Seeds();                     /* - now, plant the other seeds  */
} /* Randomize */


void

Get_Seed(long int *x)
{
   *x = seed[s];
} /* Get_Seed */


void

Put_Seed(long int x)
{
   if ((x > 0) && (x < m)) seed[s] = x;
} /* Put_Seed */

/* Current_Time removed.....*/

/*Test_Random removed...*/