molpit
Login:
Password:
remember
DIA00136 О генераторе случайных чисел

В одной из бесед PIT поднял серьезный вопрос о том, насколько качественный генератор случайных чисел я использую для программ. Как правило я пользовался стандартным генератором, который описан в начале второй главы книги Н. Вирта «Программирование на Обероне: шаг дальше Паскаля и Модулы» [1]. Что показательно, в книге это самая первая рассмотренная программа. Я решил разобраться в этом алгоритме подробнее и поэтому приведу здесь вольный перевод фрагментов книги Вирта.

<...> Концепции случайности и алгоритма (рецепта вычислений) ­— конечно несовместимы. Мы можем лишь надеяться получить алгоритм способный производить очень большие последовательности чисел таким путем, что их паттерн будет неразличим. <...>

<...> Метод, предложенный Д.Г. Лемером в 1949 году — мультипликативный линейный конгруэнтный метод — выдержал испытания временем.

Математические основы этого метода являются красивым примером элегантности и простоты. Все, что требуется — это разумный выбор двух целочисленных параметров: коэффициента m и множителя a для простой рекуррентной формулы:

 z_{n+1} = a \, z_n \; \mathbf{mod} \; m

На Компонентном Паскале это выражение запишется как:

z := (a * z) MOD m;

Последовательность должна быть начата с неким значением z_1, называемым seed. Оказывается, что выбор m и a критичен. Если m = 2^{31} - 1 = 2147483647 и a = 7^5 = 16807, то все числа между первым и m-тым возникнут лишь однажды в последовательности, определяемой уравнением выше, и составят более двух миллиардов случайных чисел, которые проходят сильные статистические тесты.

Несложно увидеть, что произведение a * z может выйти за границы целого, которое может быть представлено 32-мя битами. Существует один умный прием умножения, предложенный Л. Шраге, который решает эту проблему. Мы введем две новые константы q = \lfloor m / a \rfloor = 127 773 и r = m \, \mathbf{mod} \, a = 2 836. Тогда вычисления будут выглядеть так:

gamma := a * (z MOD q) - r * (z DIV q);
IF gamma > 0 THEN
    z := gamma
ELSE
    z := gamma + m
END;

На практике большие целые требуются редко. Поэтому мы нормализуем:

 r = z_n / m

Таким образом, получили действительное число r_n, заключенное в интервале 0 < r_n \leq 1.

В итоге мы оформляем процедуру для вычисления в виде:

PROCEDURE Uniform* (): REAL;
    CONST a = 16807; m = 2147483647; q = m DIV a; r = m MOD a;
    VAR gamma: INTEGER;
BEGIN
    gamma := a * (z MOD q) - r * (z DIV q);
    IF gamma > 0 THEN
        z := gamma
    ELSE
        z := gamma + m
    END;
    RETURN z / m (* value of the function *)
END Uniform;

<...>

На этом достаточно вольный перевод фрагментов завершим. В Блэкбоксе код этого генератора находится в файле: Obx/Mod/Random.odc.

z должен быть задан заранее, тогда процедура Uniform может быть выполнена многократно, создавая новые случайные числа. Желательно, чтобы она задавалась также каждый раз разная, поэтому для затравки я использую таймер компьютерного времени. Это гарантирует, что дважды запущенные приложения всегда будут иметь разные наборы случайных чисел.

Этот генератор пока меня не подводил, то есть решения выглядели всегда физично. Хотя данный алгоритм признан кроптографически неустойчивым, для ряда физических задач это не важно. Тут целесообразно привести цитату из Калиткина:

«Строго говоря, закономерность псевдослучайных чисел должна быть незаметна по отношению к требуемому частному применению. Поэтому в несложных или удачно сформулированных задачах можно использовать последовательности не очень хорошего качества, но при этом необходимы специальные проверки» [2].

В случае, когда потребуется более изощренные генераторы, уже реализована целая коллекция в виде подсистемы Rng, которая содержит 28 разных генераторов случайных чисел, а также специальные инструменты для оценки их качества [3].

Источники:

1. Reiser M., Wirth N. Programming in Oberon: steps beyond Pascal and Modula. New York, N.Y. : Reading, Mass: ACM Press ; Addison-Wesley Pub. Co, 1992. 320 p.

2. Калиткин Н.Н. Численные методы. Москва: Наука, 1978. 150 с.

3. Коллекция компонентов Гельмута Цинна

Ivan Denisov 19 Aug 2014 21:22
© International Open Laboratory for Advanced Science and Technology — MOLPIT, 2009–2021