В генераторе Парка-Миллера мы находили последовательность так:
I(j+1)=a*I(j)(mod m)
где a и m - особым образом выбранные константы.
Прямое приложение этого метода возможно на языках ассемблера, но языки высокого уровня могут при этом зафиксировать переполнение. Для обхода этого Scharge предложил метод частичной факторизации модуля. Модуль разлагается в выражение:
m=a*q+r
Если r и 0, то при этом величины a*(z mod q) и r*[z/q] всегда лежат в интервале 0,...,m-1. Для умножения (a*z)(mod m) при этом используется алгоритм:
t = a(z mod q)-r[z/q]
если t<0, то t += m.
(a*z)(mod m)=t.
Если требуется число вызовов, превышающее по порядку 108, то для этого случая L'Ecuyer рекомендует комбинировать вывод двух последовательностей с близкими, но отличающимися константами. В его исследованиях хороший результат был получен для значений:
m1=2147483563, a1=40014, q1=53668, r1=12211;
m2=2147483399, a2=40692, q2=52774, r2=3791.
При этом для современных компьютеров период генерируемой последовательности становится недостижимым (длина оценивается по порядку как 1018).
/* L'Ecuyer algorithm for uniform random generator
with practically endless period. Combining 2 sequences. */
/* previous constants, static variables and functions are valid as if
the programs on this page are determined in one module */
#define IM1 2147483563
#define IM2 2147483399
#undef AM
#define AM (1./IM1)
#define IMM1 (IM1-1)
#define IA1 40014
#define IA2 40692
#define IQ1 53668
#define IQ2 52774
#define IR1 12211
#define IR2 3791
#undef NDIV
#define NDIV (1+IMM1/NTAB)
float unirand2(void) {
int j;
long k;
static long dummy2=123456789;
static long iy=0;
static long iv[NTAB];
float temp;
/* initialize the random sequence (first set of coefficients, the
routine close to that in the function above */
if(dummy<=0 || !iy) {
/* avoid negative or zero seed */
if(dummy<0) dummy=-dummy;
else if(dummy==0) dummy=1;
dummy2=dummy;
/* after NWUP warmups, initialize shuffle table */
for(j=NTAB+NWUP-1;j>=0;j--) {
k=dummy/IQ1;
if((dummy=IA1*(dummy-k*IQ1)-IR1*k)<0) dummy+=IM1;
if(j }
/* first specimen from the table */
iy=iv[0];
}
/* regular work: generate 2 sequences. */
k=dummy/IQ1;
if((dummy=IA1*(dummy-k*IQ1)-IR1*k)<0) dummy+=IM1;
k=dummy2/IQ2;
if((dummy2=IA2*(dummy2-k*IQ2)-IR2*k)<0) dummy2+=IM2;
/* shuffle output combining 2 sequences */
iy=iv[j=iy/NDIV]-dummy2;iv[j]=dummy;
if(iy<1) iy+=IMM1;
/* return the result, as in the previous function */
if((temp=AM*iy)>RNMX) return(RNMX);
else return(temp);
}
8 8 8
| |
*
Обратная связь