Компактные генераторы равномерного распределения случайных чисел
Страница 4. Алгоритм Л


 

Алгоритм Л'Экюера, комбинирующий две последовательности

Если требуется число вызовов, превышающее по порядку 108, то для этого случая Л'Экюер рекомендует комбинировать вывод двух последовательностей с близкими, но отличающимися константами. В его исследованиях неплохой результат был получен для значений:
m1=2147483563, a1=40014, q1=53668, r1=12211;
m2=2147483399, a2=40692, q2=52774, r2=3791.
Причем для современных компьютеров период генерируемой последовательности становится недостижимым (длина оценивается по порядку как 1018).

/* Алгоритм Л'Экюера, комбинирующий две последовательности. */

/* предыдущие константы и функции уже были объявлены выше */
#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;
/* инициализация случайной последовательности (первая установка коэффициентов */
if(dummy<=0 || !iy) {
if(dummy<0) dummy=-dummy;
else if(dummy==0) dummy=1;
dummy2=dummy;
for(j=NTAB+NWUP-1;j>=0;j--) {
k=dummy/IQ1;
if((dummy=IA1*(dummy-k*IQ1)-IR1*k)<0) dummy+=IM1;
if(j<NTAB) iv[j]=dummy;
}
iy=iv[0];
}
/* генерируем две последовательности. */
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;
iy=iv[j=iy/NDIV]-dummy2;iv[j]=dummy;
if(iy<1) iy+=IMM1;
if((temp=AM*iy)>RNMX) return(RNMX);
else return(temp);
}
 
« Предыдущая статья   Следующая статья »