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

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

Прямое приложение этого метода возможно на языках ассемблера, но языки высокого уровня могут при этом зафиксировать переполнение. Для обхода этого Scharge предложил метод частичной факторизации модуля. Модуль разлагается в выражение:
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.

/* Minimal portable random generator by Park and Miller */

/* Lewis-Goodman-Miller constants */
#define IA 16807
#define IM 2147483647
#define AM (1./IM)
/* Scharge constants */
#define IQ 12773
#define IR 2836
/* Special mask to be explained below */
#define MASK 123456789

static long dummy;

/* initial seed, for all the generators here */
void Seed(long dum) {dummy=dum;}

/* returns random uniformly distributed between 0 and 1 */
float unirand0(void) {
 long k;
 float ans;
 dummy^=MASK;	/* avoid dummy==0 */
 k=dummy/IQ;
 if((dummy=IA*(dummy-k*IQ)-IR*k)<0) dummy+=IM;
 ans=AM*dummy;
 dummy^=MASK;	/* restore unmasked dummy */
 return(ans);
}

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

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

/* Minimal random generator with Bays-Durham shuffle */

/* previous functions, variables and constants are described above. 
   Only additional follow */
#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;
 /* initialize */
 if(dummy<=0 || !iy) {
  /* avoid negative or zero seed */
  if(dummy<0) dummy=-dummy;
  else if(dummy==0) dummy=1;
  /* after NWUP warmups, initialize shuffle table */
  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;
  }
  /* first specimen from the table */
  iy=iv[0];
 }
 /* regular work: generate new number */
 k=dummy/IQ;
 if((dummy=IA*(dummy-k*IQ)-IR*k)<0) dummy+=IM;
 /* shuffle output */
 iy=iv[j=iy/NDIV];iv[j]=dummy;
 /* return */
 if((temp=AM*iy)>RNMX) return(RNMX);
 else return(temp);
}

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

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

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

/* 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<NTAB) iv[j]=dummy;
  }
  /* 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);
}