/* newRandom.cp */

#include "newRandom.h"

/* 
	by 
	Robert M. Delaney
	Department of Physics
	Saint Louis University
	221 N. Grand Blvd.
	St. Louis, MO 63103
	(314) 658-2525
	
	revised 1/20/94 to eliminate assembly language.
	Remarkably the C program executes about twice as fast as my assembly language program!
	Using the C program on a 33MHz Radius Rocket took 3953 ticks to generate
	10,000,000 random doubles, or about 6.57 microseconds per random double.
	Using the original assembly language program took 7854 ticks, or about
	13.058 microseconds per random double.
	
	This pseudo-random number generator uses the method of Mitchell and Moore
	as described in "The Art of Computer Programming", Vol. 2, 2nd edition,
	p. 26, by Donald E. Knuth (Addison-Wesley). It is meant for scientific
	work, generating doubles (>= 0 and <1) rather than integers. The advantage
	of this method is its long period (greater than 3.6 x 10^16) and that
	no multiplication or division is required. Also, according to Knuth,
	"...it may well prove to be the very best source of random numbers for
	practical purposes."
	
	The loading of the Mitchell-Moore array, consisting of 55 8-byte integers,
	uses still another random number generator described on p. 349 of
	"Assembly Cookbook for the Apple II/IIe" by Don Lancaster (Howard W.
	Sams & Co.).
	
	The file should be compiled with ANSI (coprocessor option off) or
	ANSI-881 (coprocessor option on).
*/

#include <time.h>
#include <math.h>
#include "newRandom.h"


static unsigned long MMSeq[110]; /* will use as a sequence of 55 8-byte unsigned
							integers */
static long j,k;

static double scale1, scale2;


void InitArray(unsigned long seed)
{
	unsigned long *seedPtr = &seed;
	unsigned long saveSeed, saveSeed2;
	long	is, js;
	
	if(seed==0)
		seed = clock();
	
	for(js=0;js<220;js++)
	{
		// now shift seed left with 16 new random bits coming in from right
		// then the LSWord of seed will be a new random word to put in MM Array
		
		for(is=0;is<16;is++)
		{
			saveSeed = seed;
			saveSeed2 = seed;
			saveSeed>>=27;
			saveSeed&=1;  			// make bit 27 the LSBit and zero all other bits
			saveSeed2>>=30;
			saveSeed2&=1;	// make bit 30 the LSBit and zero all other bits
			saveSeed^=saveSeed2;	// EOR, the LSBit is the new random bit
			seed<<=1;
			seed|=saveSeed;	// seed has been shifted with one random bit coming in
		}
		
		((unsigned short*)MMSeq)[js] = ((unsigned short*)seedPtr)[1];
	}
	
	j = 46;
	k = 108; /* magic numbers */
	/* j and k are regarded as indices for a 55 element array of 8-byte integers */
	
	scale1 = pow(2,-32);
	scale2 = pow(2,-64);
	
}/* InitArray */


double randNum()
{
	unsigned long	saveLong;
	double			number;
	
	/*
		 add the two 8-byte unsigned integers pointed to by j and k
		 and put the result into container pointed to by k.
		 Use the k container to form the random double
	*/
	
	saveLong = MMSeq[k+1];
	MMSeq[k+1] = MMSeq[k+1] + MMSeq[j+1];
	MMSeq[k] = MMSeq[k] + MMSeq[j];
	if(MMSeq[k+1]<saveLong)
		MMSeq[k]++;	// add the carry
		
	number = scale1*MMSeq[k] + scale2*MMSeq[k+1];
	
	j-=2;
	k-=2;
	
	if (j == -2) j = 108;
	if (k == -2) k = 108;
	
	return(number);
	
}/* randNum */

