Статья предоставлена (c) Nikitine Valeri F. 2000, web: algorithm.narod.ru
Самая простая последовательность, которую можно предложить для
реализации генератора равномерного распределения:
I(j+1)=a*I(j)(mod m)
при соответствующем выборе констант. Константы были предложены
Park и Miller:
a=75=16807, m=231-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);
}
Использование маски связано с тем, что специфика алгоритма не
позволяет устанавливать счетчик в нуль. Но, как показывает опыт,
большинство пользователей счетчиков делают именно так. Маска гарантирует,
что установленный счетчик не будет нулем. Если вы очень уверены в
том, что человек не допустит подобной ошибки после вашего предупреждения,
то можете убрать из программы все инструкции, связанные с маской.
|