/*	RandComb.c
	combination multiplicative random number generator
		subtract two random numbers modulo moduli1-1
	see
		L'Ecuyer, Comm. of the ACM 1988 vol. 31
		Numerical Recipes in C, 2nd edition, pp.  278-86
		Dwyer, C/C++ Users Journal, June 1995
		Dwyer and Williams, 2001
	Copyright (c) 1994, 1995, 2001 by Gerald P. Dwyer, Jr.				*/

#include	<time.h>
#include	<stdlib.h>
#include	<float.h>
#include	<assert.h>

#define	TESTING

#define	TRUE (-1)
#define	FALSE 0

void InitRandComb(long *seed1, long *seed2) ;
long GetInitRand(int) ;
long RandComb(void) ;
long GenrRand(long a, long x, long modulus, long q, long r) ;

						/* first generator -- smaller of two moduluses */
#define	MULT1		 16493L				/* multiplier */
#define	MOD1		 2147483563L		/* modulus */
						/* modulus=multiplier*quotient+remainder */
#define	Q1			 130205L				/* quotient = [modulus/multiplier] */
#define	R1		 	 12498L				/* remainder */

						/* second generator */
#define	MULT2		 43049L				/* multiplier */
#define	MOD2		 2147483629L		/* modulus */
						/* modulus=multiplier*quotient+remainder */
#define	Q2			 49884L				/* quotient = [modulus/multiplier] */
#define	R2			 27313L				/* remainder */


#define	MOD_COMB		(MOD1-1)

#define	MIN_VALUE1	      1
#define	MAX_VALUE1  (MOD1-1)
#define	MIN_VALUE2	      1
#define	MAX_VALUE2  (MOD2-1)
#define	MAX_VALUE		( (MOD1<MOD2) ? MAX_VALUE1 : MAX_VALUE2)
#define	EXP_VAL		634305488L


static long rand1   = (-1) ;
static long rand2   = (-1) ;
static long RandNum = (-1) ;


/* initialize generators with seeds */
void InitRandComb(long *seed1, long *seed2)
{
	extern long rand1, rand2 ;
	extern long RandNum ;

	if (*seed1 <= 0 || *seed1 > MAX_VALUE1)
		*seed1 = GetInitRand(MAX_VALUE1) ;
	if (*seed2 <= 0 || *seed2 > MAX_VALUE2)
		*seed2 = GetInitRand(MAX_VALUE2) ;

						/* seed the routine */
	rand1 = *seed1 ;
	rand2 = *seed2 ;

	RandNum = rand1 - rand2 ;
	if (RandNum <= 0)
		RandNum += MOD_COMB ;
}


/* get a long initial seed for generator
	assumes that rand returns a short integer */
long GetInitRand(int max_value)
{
	long seed ;

	srand((unsigned int)time(NULL)) ;    /* initialize system generator */
	do {
		seed  = ((long)rand())*rand() ;
		seed += ((long)rand())*rand() ;
	} while (seed > max_value) ;

	assert(seed > 0) ;

	return seed ;
}


/* generate the next value in sequence from generator
	uses approximate factoring
	residue = (a * x) mod modulus
			= a*x - [(a*x)/modulus]*modulus
	where
		[(a*x)/modulus] = integer less than or equal to (a*x)/modulus
	approximate factoring avoids overflow associated with a*x and
		uses equivalence of above with
	residue = a * (x - q * k) - r * k + (k-k1) * modulus
	where
		modulus = a * q + r
		q  = [modulus/a]
		k  = [x/q]		(= [ax/aq])
		k1 = [a*x/m]
	assumes
		a, m > 0
		0 < InitRand < modulus
		a * a <= modulus
		[a*x/a*q]-[a*x/modulus] <= 1
			(for only one addition of modulus below) */
long	GenrRand(long a, long x, long modulus, long q, long r)
{
	long k, residue ;

	k = x / q ;
	residue = a * (x - q * k) - r * k ;
	if (residue < 0)
		residue += modulus ;

	assert(residue >= 1 && residue <= modulus-1) ;

	return residue ;
}


/* get a random number */
long RandComb(void)
{
	extern long RandNum, rand1, rand2 ;

	rand1 = GenrRand(MULT1, rand1, MOD1, Q1, R1) ;
	rand2 = GenrRand(MULT2, rand2, MOD2, Q2, R2) ;
	RandNum = rand1 - rand2 ;
	if (RandNum <= 0)
		RandNum += MOD_COMB ;

	assert(RandNum >= 1 && RandNum <= MOD_COMB) ;

	return RandNum ;
}


#if	defined(TESTING)
/* Test the generator */
#include     <stdio.h>
void main(void)
{
	extern long rand1, rand2, RandNum ;

	long seed1=1, seed2=1 ;
	int i ;

	InitRandComb(&seed1, &seed2) ;
	printf("Seeds for random number generator are %ld   %ld\n",
			seed1, seed2) ;
	i = 1 ;
	do {
		RandComb() ;
		i++ ;
	} while (i < 10000) ;

	printf("On draw 10000, random number should be %ld\n", EXP_VAL) ;
	printf("On draw %5d, random number is        %ld\n", i, RandComb()) ;

	seed1 = 1 ;
	seed2 = 1 ;
	InitRandComb(&seed1, &seed2) ;
	printf("rand1 %12ld rand2 %12ld RandNum %12ld\n", rand1, rand2, RandNum) ;
	for (i = 0; i < 5; i++) {
		RandComb() ;
		printf("rand1 %12ld rand2 %12ld RandNum %12ld\n", rand1, rand2, RandNum) ;
	}
}
#endif TESTING