Go to main content

man pages section 3: Extended Library Functions, Volume 4

Exit Print View

Updated: Thursday, June 13, 2019
 
 

mwcrans (3SUNMATH)

Name

mwcrans - multiply with carry pseudo-random number generators

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,
     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);

Description

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.

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(7) for descriptions of the following attributes:

ATTRIBUTE TYPE
ATTRIBUTE VALUE
Interface Stability
Committed
MT-Level
MT-Safe
Availability
system/library/math

See Also

addrans(3SUNMATH), lcrans(3SUNMATH), shufrans(3SUNMATH), attributes(7)

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.