mwcrans - multiply with carry pseudo-random number generators
cc [ flag ... ] file ... -lsunmath -lm [ library ... #include <sunmath.h> int i_mwcran_(void);
unsigned int u_mwcran_(void);
long i_lmwcran_(void);
unsigned long u_lmwcran_(void);
long long i_llmwcran_(void);
unsigned long long u_llmwcran_(void);
float r_mwcran_(void);
double d_mwcran_(void);
void i_mwcrans_(int *x, const int *n, const int *l,
const int *u);
void u_mwcrans_(unsigned *x, const int *n,
const unsigned *l, const unsigned *u);
void i_lmwcrans_(long *x, const int *n, const long *l,
const long *u);
void u_lmwcrans_(unsigned long *x, const int *n,
constg unsigned long *l, const unsigned long *u);
void i_llmwcrans_(long long *x, const int *n,
const long long *l, const long long *u);
void u_llmwcrans_(unsigned long long *x, const int *n,
const unsigned long long *, const unsigned long long *u);
void r_mwcrans_(float *x, const int *n, const float *l,
const float *u);
void d_mwcrans_(double *x, const int *n, const double *l,
const double *u);
void i_init_mwcrans_(void);
void smwcran_(const int *seed);
void i_set_mwcrans_(const int *p);
void i_get_mwcrans_(int *p);
These functions generate sequences of pseudo-random numbers using the multiply-with-carry algorithm developed by Marsaglia. Given a multiplier M and the initial value of a seed X and carry C (all 32 bit integers), the algorithm generates a new seed and carry as follows:
1. Z = X*M + C (here Z is a 64 bit integer) 2. new seed X = lower 32 bits of Z 3. new carry C = upper 32 bits of Z
The period of the sequence of seeds and carries is M*2**31 - 1, provided that both (M*2**32 - 1) and (M*2**31 - 1) are prime.
The functions listed above internally use two 32 bit multiply-with-carry generators denoted by mwcran0 and mwcran1. The multiplier used in mwcran0 is 526533 (0x808C5), and the multiplier used in mwcran1 is 557325 (0x8810D). The periods of both of these generators are approximately 2**50. The following descriptions refer to these generators.
u_mwcran_() invokes mwcran0 and returns a pseudo-random 32 bit unsigned integer between 0 and 2**32 - 1.
i_mwcran_() returns a pseudo-random 32 bit signed integer between 0 and 2**31 - 1 by calling u_mwcran_() and masking off the most significant bit of the result.
u_llmwcran_() invokes both mwcran0 and mwcran1 and concatenates the two 32 bit values together to return a pseudo-random 64 bit unsigned integer between 0 and 2**64 - 1. The period of this random number generator is approximately 2**100.
i_llmwcran_() returns a pseudo-random 64 bit signed integer between 0 and 2**63 - 1 by calling u_llmwcran_() and masking off the most significant bit of the result.
The functions i_lmwcran_() and u_lmwcran_() return pseudo-random long integers. The ranges of these numbers depend on the width of the long integer type: if a long int is 32 bits wide (i.e., the ILP32 data model), these functions are the same as i_mwcran_() and u_mwcran_(), respectively; if a long int is 64 bits wide (the LP64 data model), these functions are the same as i_llmwcran_() and u_llmwcran_(), respectively.
r_mwcran_() returns a pseudo-random single precision floating point number in the range [0,1). It produces this number by calling mwcran0 repeatedly to generate a sequence of random bits, which is then interpreted as a binary fraction 0.xxxxxxxxxxxxxxx... and truncated to the float format. The resulting sequence of floating point numbers is uniformly distributed in the range [0,1).
d_mwcran_() returns a pseudo-random double precision floating point number in the range [0,1). It produces this number by calling mwcran0 and mwcran1 repeatedly to generate a sequence of random bits, which is then interpreted as a binary fraction 0.xxxxxxxxxxxxxxx... and truncated to the double format. The resulting sequence of floating point numbers is uniformly distributed in the range [0,1).
i_mwcrans_(x, n, l, u) invokes mwcran0 repeatedly to generate *n pseudo-random 32 bit integers x[0], x[1], ... x[*n - 1] uniformly distributed between *l and *u. If [*l,*u] = [0,0x7fffffff], then i_mwcrans_() yields the same *n random numbers as would be generated by calling i_mwcran_() *n times.
u_mwcrans_(x, n, l, u) invokes mwcran0 repeatedly to generate *n pseudo-random 32 bit unsigned integers x[0], x[1], ... x[*n - 1] uniformly distributed between *l and *u. If [*l,*u] = [0,0xffffffff], then u_mwcrans_() yields the same *n random numbers as would be generated by calling u_mwcran_() *n times.
i_llmwcrans_(x, n, l, u) invokes mwcran0 and mwcran1 repeatedly to generate *n pseudo-random 64 bit integers x[0], x[1], ... x[*n - 1] uniformly distributed between *l and *u. If [*l,*u] = [0,0x7fffffffffffffff], then i_llmwcrans_() yields the same *n random numbers as would be generated by calling i_llmwcran_() *n times.
u_llmwcrans_(x, n, l, u) invokes mwcran0 and mwcran1 repeatedly to generate *n pseudo-random 64 bit unsigned integers x[0], x[1], ... x[*n - 1] uniformly distributed between *l and *u. If [*l,*u] = [0,0xffffffffffffffff], then u_llmwcrans_() yields the same *n random numbers as would be generated by calling u_llmwcran_() *n times.
The functions i_lmwcrans_() and u_lmwcrans_() generate arrays of pseudo-random long integers. The ranges of the numbers they generate depend on the width of the long integer type: if a long int is 32 bits wide (i.e., the ILP32 data model), these functions are the same as i_mwcrans_() and u_mwcrans_(), respectively; if a long int is 64 bits wide (the LP64 data model), these functions are the same as i_llmwcrans_() and u_llmwcrans_(), respectively.
r_mwcrans_(x, n, l, u) invokes mwcran0 repeatedly to generate *n pseudo-random single precision floating point numbers x[0], x[1], ... x[*n - 1] uniformly distributed between *l and *u (up to a truncation error). If *l is zero and *u is the largest single precision floating point number less than 1, then r_mwcrans_() yields the same *n random numbers as would be generated by calling r_mwcran_() *n times.
d_mwcrans_(x, n, l, u) invokes mwcran0 and mwcran1 repeatedly to generate *n pseudo-random double precision floating point numbers x[0], x[1], ... x[*n - 1] uniformly distributed between *l and *u (up to a truncation error). If *l is zero and *u is the largest double precision floating point number less than 1, then d_mwcrans_() yields the same *n random numbers as would be generated by calling d_mwcran_() *n times.
i_init_mwcrans_() initializes the seeds and carries for mwcran0 and mwcran1 to default values.
smwcran_(m) initializes the seeds and carries for mwcran0 and mwcran1 with a scrambled value of *m according to the following formulas.
For mwcran0:
X0 = MWCRAN_SEED0 + (*m)*0x110005
C0 = MWCRAN_CARRY0 + (*m)*0x110005
For mwcran1:
X1 = MWCRAN_SEED1 + (*m)*0x100021
C1 = MWCRAN_CARRY1 + (*m)*0x100021
Here MWCRAN_SEED0, MWCRAN_CARRY0, MWCRAN_SEED1, and MWCRAN_CARRY1 are the default initial values established by i_init_mwcrans_() for the seeds and carries of mwcran0 and mwcran1, respectively. In particular, calling smwcran_(m) with *m equal to zero is equivalent to calling i_init_mwcrans_(). (This function provides a simple way to obtain different sequences of random numbers. For more precise control of the seeds and carries, use i_set_mwcrans_().)
i_get_mwcrans_(p) and i_set_mwcrans_(p) respectively get and set the state table of mwcran0 and mwcran1. The table is an array of four integers p[0]-p[3] containing the seeds and carries for both generators: p[0] and p[1] contain the seed and carry for mwcran0, and p[2] and p[3] contain the seed and carry for mwcran1.
The following example verifies that u_mwcran_() and u_mwcrans_() are consistent and tests the uniformity of the distribution of the numbers they generate:
/* sample program to check u_mwcran*() */
#include <sunmath.h>
int main() {
unsigned x[1000],y[1000],i,lb,ub;
int n,hex1[16],hex2[16],hex3[16],seed,j;
double t1,t2,t3,v;
seed = 40;
/* generate 1000 random unsigned ints with seed initialized to 40 */
smwcran_(&seed);
for (i=0;i<1000;i++) x[i] = u_mwcran_();
/* generate the same 1000 random unsigned ints by u_mwcrans_() */
n = 1000; lb = 0; ub = 0xffffffff;
smwcran_(&seed); /* reset seed */
u_mwcrans_(y,&n,&lb,&ub);
/* verify that x and y are indeed equal */
for (n=0,i=0;i<1000;i++) n |= (x[i]-y[i]);
if(n==0) printf("Array x is equal to array y.\n\n");
else printf("***Error: array x is not equal to array y.\n\n");
/* Check uniformity by counting the appearance of hexadecimal digits
of 1000 random numbers. The expected value is 500. The counting
is repeated three times with different seed values. */
/* initialize the counting array hex1[],hex2[],hex3[] */
n = 1000; for (i=0;i<16;i++) hex1[i]=hex2[i]=hex3[i]=0;
/* count the hex digits of 1000 random numbers with seed=1 */
seed = 1; smwcran_(&seed);
u_mwcrans_(x,&n,&lb,&ub);
for (i=0;i<n;i++) {
j = x[i]; hex1[j&0xf] += 1; hex1[(j>>4)&0xf] += 1;
hex1[(j>>8)&0xf] += 1; hex1[(j>>12)&0xf] += 1;
hex1[(j>>16)&0xf] += 1; hex1[(j>>20)&0xf] += 1;
hex1[(j>>24)&0xf] += 1; hex1[(j>>28)&0xf] += 1;
}
/* count the hex digits of 1000 random number with seed=2 */
seed = 2; smwcran_(&seed);
u_mwcrans_(x,&n,&lb,&ub);
for (i=0;i<n;i++) {
j = x[i]; hex2[j&0xf] += 1; hex2[(j>>4)&0xf] += 1;
hex2[(j>>8)&0xf] += 1; hex2[(j>>12)&0xf] += 1;
hex2[(j>>16)&0xf] += 1; hex2[(j>>20)&0xf] += 1;
hex2[(j>>24)&0xf] += 1; hex2[(j>>28)&0xf] += 1;
}
/* count the hex digits of 1000 random number with seed=3 */
seed = 3; smwcran_(&seed);
u_mwcrans_(x,&n,&lb,&ub);
for (i=0;i<n;i++) {
j = x[i]; hex3[j&0xf] += 1; hex3[(j>>4)&0xf] += 1;
hex3[(j>>8)&0xf] += 1; hex3[(j>>12)&0xf] += 1;
hex3[(j>>16)&0xf] += 1; hex3[(j>>20)&0xf] += 1;
hex3[(j>>24)&0xf] += 1; hex3[(j>>28)&0xf] += 1;
}
/* Compute the Chi-square of each test */
t1 = t2 = t3 = 0.0;
for (i=0;i<16;i++) {
v = hex1[i]-500; t1 += v*v/500.0;
v = hex2[i]-500; t2 += v*v/500.0;
v = hex3[i]-500; t3 += v*v/500.0;
}
/* print out result */
printf("Expected value of each hex digit's appearance is 500.\n");
printf("Observed result with seed=1,2, and 3:\n");
printf(" Hex digit Number of appearances with\n");
printf(" seed=1 seed=2 seed=3\n");
for (i=0;i<16;i++) {
printf(" %01X: %4d %4d %4d\n",
i,hex1[i],hex2[i],hex3[i]);
}
printf(" ----------------------------\n");
printf("Chi-square value%7.2g %8.2g %8.2g\n",t1,t2,t3);
printf("Note: A reasonable range of the Chi-square value is\n");
printf(" within [7.26, 25.00] which corresponds to the 5\n");
printf(" percent and 95 percent points of the Chi-square\n");
printf(" distribution with degree of freedom equal to 15.\n");
return 0;
}
The output from the preceding program reads:
Array x is equal to array y.
Expected value of each hex digit's appearance is 500.
Observed result with seed=1,2, and 3:
Hex digit Number of appearances with
seed=1 seed=2 seed=3
0: 514 493 521
1: 529 507 480
2: 477 495 493
3: 495 541 517
4: 518 504 486
5: 496 464 484
6: 467 488 484
7: 511 487 540
8: 517 499 525
9: 500 489 490
A: 492 506 511
B: 475 504 482
C: 499 504 504
D: 485 514 493
E: 520 531 515
F: 505 474 475
----------------------------
Chi-square value 9.5 11 11
Note: A reasonable range of the Chi-square value is
within [7.26, 25.00] which corresponds to the 5
percent and 95 percent points of the Chi-square
distribution with degree of freedom equal to 15.
See attributes(7) for descriptions of the following attributes:
|
addrans(3SUNMATH), lcrans(3SUNMATH), shufrans(3SUNMATH), attributes(7)
Because these functions share the internal generators mwcran0 and mwcran1, they are not independent of one another. For example, calling i_init_mwcran_() followed by d_mwcran_() will produce a different result than calling i_init_mwcran_() followed by i_mwcran_() and then d_mwcran_(). (The generators are independent of those in other threads, however.)
The random numbers generated by mwcran0 and mwcran1 pass Marsaglia's Diehard test.