/* Listing 2
rand_por[t].c
a portable and reasonably fast random number
generator
multiplicative random number generator
see
L'Ecuyer - Comm. of the ACM, Oct. 1990, vol. 33.
Numerical Recipes in C, 2nd edition, pp. 278-86
L'Ecuyer and Cote, ACM Transactions on
Mathematical Software,
March 1991
Russian peasant algorithm -- Knuth, vol. II, pp.
442-43
Copyright (c) 1994, 1995 by Gerald P. Dwyer, Jr.
*/
#include
#include
#include
#include
#define TESTING
#define TRUE (-1)
#define FALSE 0
long init_rand_port(long seed) ;
long get_init_rand_port(void) ;
long genr_rand_port(long init_rand) ;
long rand_port(void) ;
double rand_rect_port(void) ;
long skip_ahead(long a, long init_rand, long modulu
s, long skip) ;
long mult_mod(long a, long x, long modulus) ;
#define MOD 2147483647L
/* modulus for generator */
#define MULT 41358L /* multiplier */
/* modulus = mult*quotient + remainder */
#define Q 51924L
/* int(modulus / multiplier) */
#define R 10855L /* remainder */
#define MAX_VALUE (MOD-1)
#define EXP_VAL 1285562981L
/* value for 10,000th draw */
#define IMPOSSIBLE_RAND (-1)
#define STARTUP_RANDS 16
/*throw away this number of initial random numbers */
static long rand_num = IMPOSSIBLE_RAND ;
/* initialize random number generator with seed */
long init_rand_port(long seed)
{
extern long rand_num ;
int i ;
if (seed < 1 || seed > MAX_VALUE)
/* if seed out of range */
seed = get_init_rand_port() ;
/* get seed */
rand_num = seed ;
for (i = 0; i < STARTUP_RANDS; i++)
/* and throw away */
rand_num = genr_rand_port(rand_num) ;
/* some initial ones */
return seed ;
}
/* get a long initial seed for generator
assumes that rand returns a short integer */
long get_init_rand_port(void)
{
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/modulus]
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_port(long init_rand)
{
long k, residue ;
k = init_rand / Q ;
residue = MULT * (init_rand - Q * k) - R * k ;
if (residue < 0)
residue += MOD ;
assert(residue >= 1 && residue <= MAX_VALUE) ;
return residue ;
}
/* get a random number */
long rand_port(void)
{
extern long rand_num ;
/* if not initialized, do it now */
if (rand_num == IMPOSSIBLE_RAND) {
rand_num = 1 ;
init_rand_port(rand_num) ;
}
rand_num = genr_rand_port(rand_num) ;
return rand_num ;
}
/* generates a value on (0,1) with mean of .5
range of values is [1/(MAX_VALUE+1), MAX_VALUE/
(MAX_VALUE+1)]
to get [0,1], use (double)(rand_port()-1)/(doubl
e)(MAX_VALUE-1) */
double rand_rect_port(void)
{
return (double)rand_port()/(double)
(MAX_VALUE+1) ;
}
/* skip ahead in recursion
residue = (a^skip * init) mod modulus
Use Russian peasant algorithm */
long skip_ahead(long a, long init_rand,
long modulus, long skip)
{
long residue = 1 ;
if (init_rand < 1 || init_rand > modulus-1 ||
skip < 0)
return -1 ;
while (skip > 0) {
if (skip % 2)
residue = mult_mod(a, residue, modulus) ;
a = mult_mod(a, a, modulus) ;
skip >>= 1 ;
}
residue =mult_mod(residue, init_rand, modulus) ;
assert(residue >= 1 && residue <= modulus-1) ;
return residue ;
}
/* calculate residue = (a * x) mod modulus for
arbitrary a and x
without overflow
assume 0 < a < modulus and 0 < x < modulus
use Russian peasant algorithm followed by
approximate factoring */
long mult_mod(long a, long x, long modulus)
{
long q, r, k, residue ;
residue = -modulus ;
/* to avoid overflow on addition */
while (a > SHRT_MAX) {
/* use Russian Peasant to reduce a */
if (a % 2) {
residue += x ;
if (residue > 0)
residue -= modulus ;
}
x += (x - modulus) ;
if (x < 0)
x += modulus ;
a >>= 1 ;
}
/* now apply approximate factoring to a
and compute (a * x) mod modulus */
q = modulus / a ;
r = modulus - a * q ;
k = x / q ;
x = a * (x - q * k) - r * k ;
while (x < 0)
x += modulus ;
/* add result to residue and take mod */
residue += x ;
if (residue < 0)
/* undo initial subtraction if necessary */
residue += modulus ;
assert(residue >= 1 && residue <= modulus-1) ;
return residue ;
}
#if defined(TESTING)
/* Test the generator */
#include
void main(void)
{
long seed ;
int i ;
seed = init_rand_port(1) ;
printf("Seed for random number generator is
%ld\n", seed) ;
i = STARTUP_RANDS ;
/* threw away STARTUP_RANDS */
do {
rand_port() ;
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_port()) ;
}
#endif /* TESTING */
리스트 3 Combination multiplicative
random number generator
/* rand_com[b].c
combination multiplicative random number
generator 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
#include
#include
#include
#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 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
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
