可以移植的随机数生成程序

ANSI C标准委员会的随机数生成算法：

I[j+1] = a * I[j] + c (mod m)  (c != 0)                  （1）

a = 1103515245,c = 12345,m = 2^32

优点：速度快、每次调用的运算量小

缺点：连续调用不能避免序列的相关性

 

Park和Miller的简单的乘法同余算法：

I[j+1] = a * I[j] (mod m)                                     （2）

选择的最低标准：（选择适当的选择a和m的值，可以与（1）中的c!=0生成程序一样好）

以下a和m的选择为Park和Miller提出：

a = 7^5 = 16807  m = 2^31 - 1 = 2147483647 （3）

理论证明是“完美的”可是高级语言直接实施（2）和（3）是不可能的。因为a和m-1的乘积超过了32位整数的最大值。虽然汇编中可以使用64位乘积寄存器实现，可是机器与机器之间的不可移植。于是Schrage提出了一个实现的诀窍：基于m的“近似因子分解”：

m = aq + r 即： q = [m/a],r = m mod a （方括号表示整数部分）    （4）

若r很小，特别是r < q且 0 < z < m - 1,则可以看出a(a mod q)和r[z/q]都处在0,...,m - 1范围内，且

az mod m = a(z mod q) - r[z / q]    如果>=0

                   a(z mod q) - r[a / q] + m 其他

对（4）式中的q = 127773和r = 2836

可以得出最低标准生成程序如下：

 

#include <stdio.h> #include <math.h> #include <windows.h> unsigned long next = 1; #define IA 16807 #define IM 2147483647 #define AM (1.0/IM) #define IQ 127773 #define IR 2836 #define MASK 123459876 float ran0( long *idum ) { long k; float ans; *idum ^= MASK; k = (*idum) / IQ; *idum = IA * (*idum - k * IQ) - IR * k; if( *idum < 0 ) { *idum += IM; } ans = (float)AM * (*idum); *idum ^= MASK; return ans; } int rad( void ) { next = next * 1103515245 + 12345; return (unsigned int)(next / 65536) % 32768; } void srad( unsigned int seed ) { next = seed; } int main( int argc, char **argv ) { long test; int i; unsigned long begin = ::GetTickCount(); for( i = 0; i < 1000; i++ ) { test = rad(); float k = ran0(&test); printf("%ld/n", test); } printf("Time = %dms", ::GetTickCount() - begin); return 0; }