/* rand_com[b].c
subtract two random numbers modulo moduli1-1 and shuffle
see
L'Ecuyer, Comm. of the ACM 1988 vol. 31
Numerical Recipes in C, 2nd edition, pp. 278-86
shuffling -- Knuth, vol. II
Copyright (c) 1994, 1995 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 init_rand_comb(long *seed1, long *seed2) ;
long get_init_rand(int) ;
long rand_comb(void);
long genr_rand_diff(void);
long genr_rand(long a, long x, long modulus, long q, long r) ;
/* first generator */
#define MOD1 2147483563L /* modulus */
#define MULT1 40014L /* multiplier */
/* modulus=multiplier*quotient+remainder */
#define Q1 53668L /* quotient =[modulus/multiplier] */
#define R1 12211L /* remainder */
/* second generator */
#define MOD2 2147483399L
#define MULT2 40692L
#define Q2 52774L
#define R2 3791L
#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 804307721L
#define GENR1(init_rand) genr_rand(MULT1, init_rand, MOD1, Q1, R1)
#define GENR2(init_rand) genr_rand(MULT2, init_rand, MOD2, Q2, R2)
#define IMPOSSIBLE_RAND (-1)
#define STARTUP_RANDS 16 /* throw away initial draws */
#define DIM_RAND 150 /* size of array shuffled */
static long rand1, rand2 ;
static long rand_num = IMPOSSIBLE_RAND ;
static long rands[DIM_RAND];
/* initialize generators with seeds
fill rands[] with initial values */
void init_rand_comb(long *seed1, long *seed2)
{
extern long rand1, rand2 ;
extern long rand_num;
extern long rands[] ;
int i ;
if (*seed1 <= 0 || *seed1 > MAX_VALUE1)
*seed1 = get_init_rand(MAX_VALUE1);
if (*seed2 <= 0 || *seed2 > MAX_VALUE2)
*seed2 = get_init_rand(MAX_VALUE2);
/* seed the routine */
rand1 = *seed1;
rand2 = *seed2;
genr_rand_diff() ;
for (i = 1; i < STARTUP_RANDS; i++) /* throw some away */
genr_rand_diff() ;
/* fill the array of randum number values */
for (i = 0; i < DIM_RAND; i++)
rands[i] = genr_rand_diff() ;
/* initialize rand_num for shuffling */
rand_num = rands[DIM_RAND-1] ;
}
/* get a long initial seed for generator
assumes that rand returns a short integer */
long get_init_rand(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 difference between random numbers
assumes 0 < rand1 < MOD1
0 < rand2 < MOD2
output 1 <= rand_num <= MOD_COMB */
long genr_rand_diff(void)
{
extern long rand1, rand2;
long rand_new ;
rand1 = GENR1(rand1) ;
rand2 = GENR2(rand2) ;
rand_new = rand1 - rand2 ;
if (rand_new <= 0)
rand_new += MOD_COMB ;
assert(rand_new >= 1 && rand_new <= MOD_COMB) ;
return rand_new ;
}
/* 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 < init_rand < modulus
a * a <= modulus
[a*x/a*q]-[a*x/modulus] <= 1
(for only one addition of modulus below) */
long genr_rand(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 from rands[] and replace it*/
long rand_comb(void)
{
extern long rand_num ;
extern long rands[] ;
int i ;
/* if not initialized, do it now */
if (rand_num == IMPOSSIBLE_RAND) {
rand_num = 1 ;
init_rand_comb(&rand_num, &rand_num) ;
}
/* rand_num from previous call determines next rand_num from
rands[] */
i = (int) (((double)DIM_RAND*rand_num)/(double)(MAX_VALUE)) ;
rand_num = rands[i] ;
/* replace random number used */
rands[i] = genr_rand_diff();
return rand_num ;
}
#if defined(TESTING)
/* Test the generator */
#include <stdio.h>
void main(void)
{
long seed1=1, seed2=1 ;
int i ;
init_rand_comb(&seed1, &seed2);
printf("Seeds for random number generator are %ld %ld\n",
seed1, seed2) ;
i = STARTUP_RANDS + DIM_RAND ;
do {
rand_comb();
i++ ;
} while (i < 9999) ;
printf("On draw 10000, random number should be %ld\n", EXP_VAL) ;
printf("On draw %d, random number is %ld\n", i+1, rand_comb()) ;
}
#endif TESTING
/* End of File */