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);
这些函数使用 Marsaglia 开发的带进位相乘算法生成伪随机数序列。在提供了乘数 M 以及种子 X 和进位 C 的初始值(全部为 32 位整数)的情况下,该算法会按如下方式生成一个新种子和进位:
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
种子和进位的序列的周期为 M*2**31 - 1,前提是 (M*2**32 - 1) 和 (M*2**31 - 1) 都是素数。
上面列出的函数内部使用两个 32 位带进位相乘生成器,分别用 mwcran0 和 mwcran1 表示。mwcran0 中使用的乘数为 526533 (0x808C5),mwcran1 中使用的乘数为 557325 (0x8810D)。这两个生成器的周期都大约为 2**50。以下说明涉及这些生成器。
u_mwcran_() 调用 mwcran0 并返回介于 0 和 2**32 - 1 之间的伪随机 32 位无符号整数。
i_mwcran_() 通过调用 u_mwcran_() 并屏蔽结果的最高有效位,返回介于 0 和 2**31 - 1 之间的伪随机 32 位有符号整数。
u_llmwcran_() 同时调用 mwcran0 和 mwcran1 并将两个 32 位值连接在一起以返回介于 0 和 2**64 - 1 之间的伪随机 64 位无符号整数。此随机数生成器的周期大约为 2**100。
i_llmwcran_() 通过调用 u_llmwcran_() 并屏蔽结果的最高有效位,返回介于 0 和 2**63 - 1 之间的伪随机 64 位有符号整数。
函数 i_lmwcran_() 和 u_lmwcran_() 返回伪随机长整数。这些数的范围取决于长整数类型的宽度:如果长整数为 32 位宽(即 ILP32 数据模型),则这些函数分别与 i_mwcran_() 和 u_mwcran_() 相同;如果长整数为 64 位宽(LP64 数据模型),则这些函数分别与 i_llmwcran_() 和 u_llmwcran_() 相同。
r_mwcran_() 返回一个 [0,1) 范围中的伪随机单精度浮点数。此函数重复调用 mwcran0 生成一个随机位序列,然后将此序列解释为二进制小数 0.xxxxxxxxxxxxxxx... 并将其截断为浮点格式,从而生成此数。结果浮点数序列在 [0,1) 范围内均匀分布。
d_mwcran_() 返回一个 [0,1) 范围中的伪随机双精度浮点数。此函数重复调用 mwcran0 和 mwcran1 生成一个随机位序列,然后将此序列解释为二进制小数 0.xxxxxxxxxxxxxxx... 并将其截断为双精度格式,从而生成此数。结果浮点数序列在 [0,1) 范围内均匀分布。
i_mwcrans_(x, n, l, u) 重复调用 mwcran0 以生成 *n 个伪随机 32 位整数 x[0]、x[1]、...x[*n - 1] 在 *l 与 *u 之间均匀分布。如果 [*l,*u] = [0,0x7fffffff],则 i_mwcrans_() 获得的 *n 个随机数与通过调用 i_mwcran_() *n 次生成的随机数相同。
u_mwcrans_(x, n, l, u) 重复调用 mwcran0 以生成 *n 个伪随机 32 位无符号整数 x[0]、x[1]、...x[*n - 1] 在 *l 与 *u 之间均匀分布。如果 [*l,*u] = [0,0xffffffff],则 u_mwcrans_() 获得的 *n 个随机数与通过调用 u_mwcran_() *n 次生成的随机数相同。
i_llmwcrans_(x, n, l, u) 重复调用 mwcran0 和 mwcran1 以生成 *n 个伪随机 64 位整数 x[0]、x[1]、...x[*n - 1] 在 *l 与 *u 之间均匀分布。如果 [*l,*u] = [0,0x7fffffffffffffff],则 i_llmwcrans_() 获得的 *n 个随机数与通过调用 i_llmwcran_() *n 次生成的随机数相同。
u_llmwcrans_(x, n, l, u) 重复调用 mwcran0 和 mwcran1 以生成 *n 个伪随机 64 位无符号整数 x[0]、x[1]、...x[*n - 1] 在 *l 与 *u 之间均匀分布。如果 [*l,*u] = [0,0xffffffffffffffff],则 u_llmwcrans_() 获得的 *n 个随机数与通过调用 u_llmwcran_() *n 次生成的随机数相同。
函数 i_lmwcrans_() 和 u_lmwcrans_() 生成伪随机长整数的数组。它们生成的数的范围取决于长整数类型的宽度:如果长整数为 32 位宽(即 ILP32 数据模型),则这些函数分别与 i_mwcrans_() 和 u_mwcrans_() 相同;如果长整数为 64 位宽(LP64 数据模型),则这些函数分别与 i_llmwcrans_() 和 u_llmwcrans_() 相同。
r_mwcrans_(x, n, l, u) 重复调用 mwcran0 以生成 *n 个伪随机单精度浮点数 x[0]、x[1]、...x[*n - 1] 在 *l 与 *u(最高到发生截断错误)之间均匀分布。如果 *l 是零并且 *u 是小于 1 的最大单精度浮点数,则 r_mwcrans_() 获得的 *n 个随机数与通过调用 r_mwcran_() *n 次生成的随机数相同。
d_mwcrans_(x, n, l, u) 重复调用 mwcran0 和 mwcran1 以生成 *n 个伪随机双精度浮点数 x[0]、x[1]、...x[*n - 1] 在 *l 与 *u(最高到发生截断错误)之间均匀分布。如果 *l 是零并且 *u 是小于 1 的最大双精度浮点数,则 d_mwcrans_() 获得的 *n 个随机数与通过调用 d_mwcran_() *n 次生成的随机数相同。
i_init_mwcrans_() 将 mwcran0 和 mwcran1 的种子和进位初始化为缺省值。
smwcran_(m) 根据以下公式使用 *m 的杂乱值初始化 mwcran0 和 mwcran1 的种子和进位。
For mwcran0:
X0 = MWCRAN_SEED0 + (*m)*0x110005
C0 = MWCRAN_CARRY0 + (*m)*0x110005
For mwcran1:
X1 = MWCRAN_SEED1 + (*m)*0x100021
C1 = MWCRAN_CARRY1 + (*m)*0x100021
此处,MWCRAN_SEED0、MWCRAN_CARRY0、MWCRAN_SEED1 和 MWCRAN_CARRY1 是 i_init_mwcrans_() 分别为 mwcran0 和 mwcran1 的种子和进位建立的缺省初始值。需要特别指出的是,调用 smwcran_(m) 时 *m 等于零等效于调用 i_init_mwcrans_()。(此函数提供了一种获取不同随机数序列的简单方式。如需更加精确地控制种子和进位,可使用 i_set_mwcrans_()。)
i_get_mwcrans_(p) 和 i_set_mwcrans_(p) 分别获取和设置 mwcran0 和 mwcran1 的状态表。该表是四个整数 p[0]-p[3] 的数组,其中包含两个生成器的种子和进位:p[0] 和 p[1] 包含 mwcran0 的种子和进位,p[2] 和 p[3] 包含 mwcran1 的种子和进位。
以下示例验证 u_mwcran_() 和 u_mwcrans_() 是否一致,并测试它们生成的数的分布均匀性:
/* 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;
}
上一个程序的输出如下:
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(5):
|
addrans(3M)、lcrans(3M)、shufrans(3M)、attributes(5)
因为这些函数共享内部生成器 mwcran0 和 mwcran1,因此它们不是互相独立的。例如,调用 i_init_mwcran_() 之后再调用 d_mwcran_() 产生的结果不同于依次调用 i_init_mwcran_()、i_mwcran_() 和 d_mwcran_()。(但是,生成器独立于其他线程中的那些生成器。)
mwcran0 和 mwcran1 生成的随机数通过 Marsaglia 的 Diehard 测试。