/*  TAUSWORTHE random number generator */
/*--- EDIT # 0111	10 May 1982   11:03:08	DB1:[21,111]RAND.C;5  */
/*--- PREVIOUS EDIT	6 May 1982   13:47:52	DB1:[21,111]RAND.C;4  */
/*
	.title random
;+
; index	Random number generator
;
; Usage
;	long int random();
;
;	double frandom();
;
;	VOID randinit(seed);	Don't give a seed other than zero!!
;		You don't need to do this to start, only if you wish
;		to re-start the generator.
;
; in
;	none
;
; out
;	returns the random number
;
; Description
;	Random returns a 28-bit random number.
;	Frandom returns a 28-bit random number in the range 0.0--1.0
;	Randinit initializes the generator. The seed should be
;	zero, to use the default seed. You should NOT select your
;	own seed, as the default one was very carefully chosen.
;
; Bugs
;	You do not need to call randinit initially, only if you
;	want to re-start the generator after pulling out some 
;	random numbers. The generator always goes thru the same
;	sequence after randinit with the same seed.
;
;	If you want to start at a random spot, you should get and
;	ignore the first "n" numbers. You might select "n" to be
;	the current time of day in msec, or something.
;
; status
;	none
;
; Author
;	 Published by Robert C. Tausworthe
;	 Coded in C for DECUS (and myself) by Ray Van Tassle
;
; internal
;	This is a Tausworthe random number generator. Numbers, formed as
;	non-overlapping, adjacent 28-bit words taken from the
;	bit stream produced by the formula:
;	"a(m+532) = a(m+37) + a(m) (modulo 2)"
;	will not repeat within the projected age of the
;	solar system, will show no ensemble correlation, will exhibit
;	uniform distribution of adjacent numbers up to 19 dimensions,
;	and will not deviate from random runs-up and
;	runs-down behavior.
;
;	The generator based on this polynomial has period 1.4E160,
;	has virtually no (average) correlation between any numbers
;	separated by distances less then 5.0E158, and has adjacencies
;	up to 19 dimensions uniformly distributed. Runs-up and runs-down
;	statistics up to length 16 are impeccable. The period and
;	maximum correlation distance are, in fact, so great that the
;	generator would have to produce numbers at a nano-second
;	rate for more than 1.0E142 years before nonrandom
;	distribution or correlation effects would be noticeable as nonrandom.
;	(And I, for one, can't stay awake that long).
;
;	Almost 4.0E14 numbers would have to be examined to detect
;	deviations in runs-up and runs-down statistics as nonrandom.
;	The original author claims that all tests run on it
;	(on the Caltech PDP-10 and the JPL Univac 1108) validate
;	the properties claimed by the theory. In the theory the
;	only factor left to chance is the specification of the initial
;	seed. The stated uniformity, zero-correlation, and runs
;	statistics are all based on the single assumption that the
;	seed be chosen randomly.
;	(Consider the recursiveness of randomly picking a random
;	seed for a random number generator!)
;	Of course, the default value was not randomly chosen, but
;	was chosen specifically to look random except for the first word
;	and, certainly to the extent of the test run, this appears
;	to have worked beautifully. The number was chosen as
;	41,475,557 for the L=28 case from theoretical results
;	published in Ahrens & Dieter.
;
;	It was also demonstrated that the generator is also capable
;	(as is every known random number generator) of producing
;	numbers with 3-sigma variations from randomness over a few
;	thousand samples when the wrong seed is supplied. Moral: don't
;	pick your own seed, use the default one.
;
;

; Updates
;	 Date		Vers 	Who	Description
;	 8 Apr 1982	0001	RVT	Initial coding
;
;-
*/
#define TRYOUT		/* create a testing MAIN program */
#undef TRYOUT

#ifndef decus		/* Whitesmith's otherwise DECUS */
#include <std.h>
#else
#include <stdio.h>
#endif


/* by Ray Van Tassle 3 Feb 1982
*/

#define TWO28 (0X10000000l)	/* 2**28    */
#define MSK28 (0X0fffffffl)	/* 2**28 - 1	*/
#define MULTIPLIER 41475557l	/* the seed */

int rn_index = -1;	/* index for next random number */
long int w[19] = 0;	/* array of 19 28-bit random numbers */
long dsl = 0;

/***********************************/
/* start things off */
randinit(starter)
int starter;

{
   register int i;
   long int start;

   start = starter;
   if (starter == 0)
      start = MULTIPLIER;

/* The seed is set up with a linear congruential multiplier */
   w[0] = start;
   for (i = 1; i <= 18; i++) {
      w[i] = (w[i-1] * MULTIPLIER) & MSK28;
   }
   rn_index = 0;
}


/*************************************/
/* Returns a long int random number.
 * Range 0 to 2**28-1.
 *
 * The routine that really gets the next random number.
 * They are generated 19 at a time.
 * If we have used up all 19, generate the next set.
*/
long int random()
{
   register int i;
   int i1, i2;

   if (rn_index < 0)	/* User didn't init, so I will */
      randinit(0);

   if (rn_index > 18) {
/* generate the next 532 bits sequence */
      for (i = 0; i <= 18; i++) {
	 i1 = i+1; i2 = i+2;
	 if (i1 > 18)
	    i1 = i1-19;
	 if (i2 > 18)
	    i2 = i2-19;
	 dsl = (w[i1] << 9) | (w[i2] >> (28-9));
	 dsl = dsl & MSK28;
	 w[i] = w[i] ^ dsl;
      }

      rn_index = 0;
   }
   return (w[rn_index++]);
}


/***************************/
/* Get a double floating random number, 
 *     range 0.0-1.0
*/
double frandom()
{
   return ((double)random() / TWO28);
}


#ifdef TRYOUT
/********************/
/****** Test driver program ******/
/* The 1st 40 numbers generated are below, listed left to right.
 (The 1st 19 are the generated seed, the others are from the full
algorithm proper)
Long int:
41475557.	184889049.	214711581.	51164145.	212599701.
214881609.	257740621.	51570913.	241995333.	138728633.
166592125.	187600593.	99102197.	221617449.	114948781.
159159745.	238129829.	88645273.	77067229.	136670844.
58660024.	95149704.	185309288.	17207422.	89149227.
177015936.	168715241.	130415992.	62402524.	82861248.
192243063.	248676142.	247514630.	267549035.	173662056.
252485687.	179250589.	239799218.	101457609.	79824345.

float:
0.1545085	0.6887654	0.7998630	0.1906013	0.7919956
0.8004964	0.9601586	0.1921166	0.9015029	0.5168044
0.6206040	0.6988667	0.3691845	0.8255893	0.4282176
0.5929163	0.8871027	0.3302294	0.2870978	0.5091386
0.2185256	0.3544603	0.6903309	0.0641026	0.3321068
0.6594357	0.6285133	0.4858374	0.2324675	0.3086822
0.7161612	0.9263908	0.9220638	0.9966978	0.6469416
0.9405825	0.6677605	0.8933217	0.3779590	0.2973689
*/

main(argc, argv)
int argc;
char *argv[];
{
   int cflag;		/* number of random numbers to type */
   int fflag;		/* FALSE means long integer, TRUE means double */
   char c;
   char *cp;
   register int i;

   fflag = cflag = 0;
   printf ("\n");
   for (i=1; i<argc; ++i) {
      cp = argv[i];
      if (*cp == '-') {
	 argv[i] = NULL;
	 ++cp;
	 while (c = *cp++)
	    switch (tolower(c)) {

	 case 'f':
	    ++fflag;
	    break;

	 case 'c':
#ifndef decus
	    btoi(cp, 5, &cflag, 10);
#else
	    cflag = atoi(cp);
#endif
	    break;
	 }
      }
   }

   if (cflag == 0)
      cflag = 40;
   for (i = 1; i <= cflag; i++) {
      if (fflag)
	 printf ("%11.7f\n", frandom());
      else
	 printf ("%11ld.\n", random());
   }
}
#endif