Полезное для программистов:

Фриланс
Новости
Статьи
   
Рубрики:


Минимальный генератор Парка-Миллера с перетасовкой и без

Поиск:
Самая простая последовательность, которую можно предложить для реализации генератора равномерного распределения:
I(j+1)=a*I(j)(mod m)

при соответствующем выборе констант. Константы были предложены Park и Miller:
a=(7^5)=16807, m=(2^31)-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, где он начинает повторяться.







Просмотров: 2702

 

 

Новые статьи:


Популярные:
  1. Как сделать цикличным проигрывание MIDI-файла?
  2. Создание AVI файла из рисунков
  3. Как устройство "отключить в данной конфигурации"?
  4. Kто в данный момент присоединен через Сеть?
  5. Как узнать количество доступной памяти?
  6. Как реализовать в RichEdit разноцветный текст?
  7. Как скрыть свое приложение от ProcessViewer
  8. Как программно нажать/скрыть/показ кнопку "Start"?
  9. Модуль работы с ресурсами в PE файлах
10. Функции вызова диалоговых окон выбора
11. Проверка граматики средствами Word'а из Delphi.
12. Модуль для упрощенного вызова сообщений
13. Функции для записи и чтение своих данных в, ЕХЕ- файле
14. Рекурсивный просмотр директорий
15. Network Traffic Monitor
16. Разные модули
17. Универсальная функция для обращения к любым экспортируем функциям DLL
18. Библиотека от VladS
19. Протектор для UPX'а
20. Еще об ICQ, сообщения по контакт листу?
21. Использование открытых интерфейсов
22. Теория и практика использования RTTI
23. Работа с TApplication
24. Примеры использования Drag and Drop для различных визуальных компонентов
25. Что такое порт? Правила для работы с портами
26. Симфония на клавиатуре
27. Загрузка DLL
28. Исправление автоинкремента
29. Взаимодействие с чужими окнами
30. Проверить дубляжи в столбце


 

 

 
 
На главную