Компактные генераторы равномерного распределения случайных чисел

Компактные генераторы равномерного распределения случайных чисел с хорошими статистическими свойствами, отличающиеся от быстрых и грязных, предназначенных для той же цели.

 

Минимальный генератор Парка-Миллера

Наиболее простая последовательность, которую можно предложить для реализации генератора равномерного распределения:
I(j+1)=a*I(j)(mod m)
при соответствующем выборе констант. Константы были предложены Парком и Миллером:
a=7⁵=16807, m=2³¹-1=2147483647
и протестированы в исследованиях Льюисом, Гудманом и Миллером (1969).

Реализация этого метода возможна на языках ассемблера, но языки высокого уровня могут зафиксировать переполнение. Для обхода этого Шарж предложил метод частичной факторизации модуля. Модуль разлагается в следующее выражение:
m=a*q+r
Если r<q и 0<z<m-1, то при этом величины 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.

В случае констант Парка-Миллера можно использовать q=12773 и r=2836.

/* Минимальный компактный генератор случайных чисел Парка и Миллера */

/* константы Льюиса-Гудмана-Миллера */
#define IA 16807
#define IM 2147483647
#define AM (1./IM)
/* константы Шаржа */
#define IQ 12773
#define IR 2836
/* Специальная маска (см. ниже) */
#define MASK 123456789

static long dummy;

/* начальное значение для всех генераторов */
void Seed(long dum) {dummy=dum;}

/* возвращает случайное число на промежутке от 0 до 1 */
float unirand0(void) {
 long k;
 float ans;
 dummy^=MASK;	/* избегаем dummy==0 */
 k=dummy/IQ;
 if((dummy=IA*(dummy-k*IQ)-IR*k)<0) dummy+=IM;
 ans=AM*dummy;
 dummy^=MASK;	/* восстанавливаем dummy */
 return(ans);
}

Использование маски вызвано тем, что специфика алгоритма не позволяет устанавливать счетчик в нуль. Но, как показывает опыт, большинство пользователей счетчиков делают именно так. Маска гарантирует, что установленный счетчик не будет нулем. Если вы очень уверены в том, что человек не допустит подобной ошибки после вашего предупреждения, то можете убрать из программы все инструкции, связанные с маской.

 

Улучшенная версия генератора Парка-Миллера

В дополнение к предыдущему алгоритму, производит перетасовку по методу Бейса и Дурхема, что позволяет уничтожить нежелательные корреляции между сериями сгенерированных случайных чисел.

/* Минимальный компактный генератор случайных чисел Бейса-Дурхема */

/* Предыдущие фукнции и константы были описаны выше. */
#define NTAB 32
#define NWUP 8
#define NDIV (1+(IM-1)/NTAB)
#define EPS 1.2e-7
#define RNMX (1.0-EPS)

float unirand1(void) {
 int j;
 long k;
 static long iy=0,iv[NTAB];
 float temp;
 /* инициализация */
 if(dummy<=0 || !iy) {
  /* следим, чтобы значение было положительным */
  if(dummy<0) dummy=-dummy;
  else if(dummy==0) dummy=1;
  for(j=NTAB+NWUP-1;j>=0;j--) {
   k=dummy/IQ;
   if((dummy=IA*(dummy-k*IQ)-IR*k)<0) dummy+=IM;
   if(j<NTAB) iv[j]=dummy;
  }
  /* первый элемент из таблицы */
  iy=iv[0];
 }
 /* генерируем новое число */
 k=dummy/IQ;
 if((dummy=IA*(dummy-k*IQ)-IR*k)<0) dummy+=IM;
 iy=iv[j=iy/NDIV];iv[j]=dummy;
 if((temp=AM*iy)>RNMX) return(RNMX);
 else return(temp);
}

Маска здесь не применяется, так как метод инициализации алгоритма отсекает ошибочные значения текущего счетчика.

Алгоритм Парка-Миллера с перетасовкой проходит все известные тесты, включая и те, где без перетасовки этот алгоритм дает сбой. Хорошее поведение наблюдается до значений числа вызовов счетчика порядка 10^8, где он начинает повторяться.

 

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

Если требуется число вызовов, превышающее по порядку 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);
}