NAME
mwcrans - multiply with carry pseudo-random number genera-
tors
SYNOPSIS
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, const
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 *l, 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);
DESCRIPTION
These functions generate sequences of pseudo-random numbers
using the multiply-with-carry algorithm developed by Mar-
saglia. Given a multiplier M and the initial value of a
seed X and carry C (all 32 bit integers), the algorithm gen-
erates 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 concaten-
ates 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_(), respec-
tively.
r_mwcran_() returns a pseudo-random single precision float-
ing 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 frac-
tion 0.xxxxxxxxxxxxxxx... and truncated to the float format.
The resulting sequence of floating point numbers is uni-
formly distributed in the range [0,1).
d_mwcran_() returns a pseudo-random double precision float-
ing 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 gen-
erate *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 gen-
erate *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 repeat-
edly 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 repeat-
edly 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 gen-
erate *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 repeat-
edly to generate *n pseudo-random double precision floating
point numbers x[0], x[1], ... x[*n - 1] uniformly distri-
buted 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 pre-
cise 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.
EXAMPLE
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.
ATTRIBUTES
See attributes(5) for descriptions of the following attri-
butes:
______________________________________________________________
| ATTRIBUTE TYPE | ATTRIBUTE VALUE |
|____________________|________________________________________|
| Availability | SPROsunms, SPROlang (32 bit libraries)|
| | SP64sunms, SP64lang (64 bit libraries)|
| Interface Stability| Evolving |
| MT-Level | MT-Safe |
|____________________|________________________________________|
SEE ALSO
addrans(3M), lcrans(3M), shufrans(3M), attributes(5)
NOTES
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.