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

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

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

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

/* предыдущие константы и функции уже были объявлены выше */
#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);
}