Generators of pseudorandom numbers

Датчики псевдослучайных чисел

Someone

Линейный конгруэнтный метод

Поскольку на нашем сайте можно найти датчик (генератор) псевдослучайных чисел RandAs, к нам иногда обращаются с вопросом о том, с помощью какого алгоритма можно генерировать такие числа. Этот вопрос подробно обсуждается в книге [1], поэтому здесь мы ограничимся только некоторыми простыми замечаниями.

Речь будет идти исключительно о генерации целых чисел отрезка [0;m-1], где m>1 - натуральное число.

Прежде всего, нужно понимать, что последовательность чисел, порождаемая любым алгоритмом, не может считаться случайной, если под случайной последовательностью понимать результат серии независимых испытаний, в которых все исходы равновероятны. Однако нам вполне достаточно того, что эта последовательность в интересующих нас приложениях выглядит как случайная. С другой стороны, в бесконечной случайной последовательности с вероятностью 1 должна бесконечно много раз встретиться любая наперёд заданная конечная последовательность целых чисел отрезка [0;m-1], например, последовательность, состоящая из миллиона нулей; никто не может "запретить" этой последовательности встретиться в самом начале наших "экспериментов". Так что обнаружение каких-либо закономерностей в просмотренном конечном отрезке бесконечной последовательности нельзя рассматривать как доказательство её неслучайности.

Разумеется, такие "отклонения от случайности" в наперёд заданном отрезке случайной последовательности маловероятны и встречаются редко, но они вполне возможны [Примечание 1].

Основная проблема, таким образом, состоит в том, чтобы найти алгоритм, порождающий последовательность целых чисел отрезка [0;m-1]

X0,X1,X2,…,Xk,…, (1)

использование которой в нашей программе давало бы такие же результаты, какие ожидаются при использовании истинно случайной последовательности. Это пожелание не является достаточно определённым. Обычно удовлетворяются тем, что некоторый набор статистических тестов даёт для нашей последовательности примерно такие же результаты, какие даёт теоретический расчёт для истинно случайной последовательности. Подробно этот вопрос рассматривается в уже упомянутой книге [1].

Обычно считается, что алгоритм датчика псевдослучайных чисел должен быть как можно более простым, чтобы выработка очередного псевдослучайного числа занимала мало времени. Это требование, однако, не выглядит чрезвычайно важным. Если исключить случай патологически сложного алгоритма, то время работы датчика вряд ли будет составлять сколько-нибудь существенную часть всего времени работы программы, так как практически в любой разумной программе выработка псевдослучайных чисел составляет малую часть всех вычислений. Исключением может быть, например, программа тестирования датчика псевдослучайных чисел, которая вырабатывает большое количество таких чисел, выполняя над ними лишь небольшое количество простых вычислений.

В качестве простого, но эффективного алгоритма выработки псевдослучайных последовательностей, в книге [1] рекомендуется и подробно изучается так называемый линейный конгруэнтный метод.

Этот метод состоит в следующем. Выбирается натуральное число m>1, которое далее называется модулем. Затем выбираются целые числа a (множитель) и c (приращение), удовлетворяющие условиям 0≤a<m и 0≤c<m. Начав с произвольного целого числа X0 (0≤X0<m), порождаем последовательность псевдослучайных чисел (1) по формуле [Примечание 2]

Следующее Xk, k≥0. (2)

Такая последовательность обязательно является периодической, так как число Xk может иметь не более m различных значений, причём, длина периода не может быть больше m. Многолетняя практика использования псевдослучайных последовательностей показывает, что крайне желательно, чтобы период псевдослучайной последовательности был настолько большим, чтобы реально используемая в программе часть последовательности составляла очень малую часть этого периода. А это означает, что модуль m должен быть очень большим числом.

Максимальный достижимый период последовательности (2) равен m. Такой период получается тогда и только тогда, когда параметры a, c и m удовлетворяют следующим условиям:

  1. число b=a-1 делится на все простые делители числа m;
  2. если m делится на 4, то и b делится на 4;
  3. число c взаимно просто с m.

Далее будем предполагать, что эти условия выполняются.

Заметим, что значение a=1 удовлетворяет условиям I и II, но брать такое a не следует: в этом случае

Xk=(X0+kc)mod m,

так что трудно рассчитывать на "случайность" такой последовательности. Поэтому далее мы будем предполагать, что a>1 и, соответственно, b>0. В этом случае, используя формулу суммы геометрической прогрессии, легко получить соотношение, обобщающее исходное соотношение (2):

Обобщение формулы (2). (3)

Условие b>0 вместе с условием максимальности периода сразу же исключают многие значения m. Например, если m=2^n1*p1*p2*…*pn2, где 2<p1<p2<…<pn2 - простые числа, а n1≤2, то единственное значение b<m, удовлетворяющее условиям I и II, - это b=0. Поэтому такие значения m рассматривать не будем.

Другое ограничение накладывается понятием мощности линейной конгруэнтной последовательности. Мощность определяется как наименьшее целое число s, для которого b^s делится на m. Из условия I следует, что такое s обязательно существует.

Предположим, что мощность равна s. Пользуясь для вычисления a^n=(1+b)^n формулой бинома Ньютона [Примечание 3], можно привести формулу (3) к виду

Преобразование формулы (3) (4)

или к виду

Преобразование формулы (3). (5)

Из этого выражения следует, что элементы линейной конгруэнтной последовательности выражаются через их номера многочленом степени s (по выбранному модулю m).

Практика показывает, что последовательности с небольшой мощностью "недостаточно случайны", поэтому желательно, чтобы мощность была "достаточно" большой; при этом большая мощность сама по себе ещё не гарантирует "хорошего" поведения псевдослучайной последовательности.

Практически это означает, что модуль m должен делиться на достаточно большую степень какого-нибудь простого числа p. Если m делится на p^n, где p - нечётное простое число, то мы можем рассчитывать на мощность s≥n, взяв b делящимся на p, но не делящимся на p^2; если m делится на 2^n, где n≥2, то можно рассчитывать на мощность s≥… [Примечание 4], взяв b делящимся на 4, но не делящимся на 8 (это означает, что множитель a при делении на 8 должен давать в остатке 5).

Поскольку компьютеры, как правило, используют двоичную систему счисления, чаще всего берут в качестве модуля степень числа 2: m=2^n. Это удобно ещё и потому, что при вычислении остатка при делении на 2^n можно обойтись без деления, просто отбросив старшие разряды числа. Как было только что сказано, в этом случае наибольшая мощность s=… получается, если множитель a при делении на 8 даёт в остатке 5.

Заметим также, что идея использовать числа последовательности (1) не подряд, а только каждое n-ое (n>1), вопреки интуитивному ожиданию, что в последовательности

X0,Xn,X2n,…,Xkn,… (6)

зависимость между соседними элементами меньше, чем в исходной последовательности (1), не оправдывается.

Действительно, из формулы (3) следует, что эта последовательность имеет такой же вид, что и исходная последовательность (1), только имеет параметры a^n mod m, (a^n-1) mod m и (a^n-1)c/b mod m. Ясно, что "прореженная" последовательность не может иметь больший период, чем последовательность (1). Из указанных формул, в частности, следует, что если m делится на p^t1, а b делится на p^t2, то a^n-1 делится на p^Min{t1,t2} и может делиться на большую степень числа p. Поэтому "прореженная" последовательность имеет мощность, не превосходящую мощности первоначальной последовательности. В действительности, как мы сейчас увидим, период (и мощность) последовательности (6) может быть меньше периода (и мощности) последовательности (1).

Обозначим d=НОД(n,m) [Примечание 5] и q=n/d. Проверим, что число q=(a^n-1)/b делится на d. Используя формулу бинома Ньютона для вычисления a^n=(1+b)^n, получим

q=(a^n-1)/b. (7)

Предположим, что d делится на p^s и не делится на p^(s+1), где p - простое число и s≥1. Тогда b делится на p, так как m делится на p, а b делится на все простые делители числа m.

Первое слагаемое суммы (7) делится на p^s.

Рассмотрим k-тое слагаемое k-тое слагаемое при k>1. Числитель этой дроби делится на p^t1, где t1, а знаменатель - на p^t2, где (учитывая, что p≥2, и что в первых двух суммах имеется только конечное число слагаемых, отличных от 0) t2 и, следовательно, t2≤k-1. Поэтому k-тый член суммы (7) делится на p^(t1-t2), где t1-t2. Поэтому все члены суммы (7) делятся на p^s. Таким образом, число q=(a^n-1)/b делится на все простые делители числа d в необходимых степенях; следовательно, оно делится и на d.

Пользуясь только что доказанной делимостью, проверим, что все члены последовательности (6) при делении на d дают один и тот же остаток. Для этого нужно доказать, что разность разность делится на d. Так как модуль m делится на d, и имея в виду формулу (3), достаточно проверить, что X(k+1)n-Xkn делится на d, а это сразу же следует из доказанного.

Отсюда следует также, что период последовательности (6) не превосходит m/d. Легко проверить, что этот период равен m/d. Действительно, по свойству наибольшего общего делителя существуют такие целые числа k и l, что выполняется равенство kn-lm=d. Тогда, как легко видеть, Xkn=Xd, X2kn=X2d и так далее. Поэтому последовательность (6) содержит все числа последовательности X0,Xd,X2d,…,Xkd,…, которых как раз m/d штук.

Обозначим r остаток, который дают при делении на d все члены последовательности (6). Считаем, естественно, что 0≤r<d. Тогда число Yk=(Xkn-r)/d является целым и 0≤Yk<m/d.

Элементы последовательности (6) выражаются через элементы последовательности

Y0,Y1,Y2,…,Yk,… (8)

формулой

Xkn=dYk+r. (9)

Последовательность (8) также является линейной конгруэнтной последовательностью:

Следующее Yk, k≥0, (10)

где mY=m/d, aY=a^n mod mY и сY=(a^n-1)(br+c)/bd mod mY.

Это ещё раз показывает, что "прореженная" последовательность не имеет никаких очевидных преимуществ по сравнению с исходной.

Однако такого рода "прореженные" последовательности при использовании датчиков псевдослучайных чисел возникают естественным образом. Нередко каждый цикл расчётов требует не одного псевдослучайного числа, а n>1. В таком случае элементы последовательности (1) используются группами:

X1,X2,…,Xn;….

Взяв первые элементы каждой группы, получим "прореженную" последовательность типа (6).

Спектральный тест

В книге [1] описано довольно большое число различных тестов для проверки датчиков псевдослучайных чисел. Наиболее важным там считается так называемый спектральный тест. Этот тест основан на дискретном преобразовании Фурье.

Предположим, что датчик псевдослучайных чисел используется следующим образом: каждый цикл расчётов использует N≥1 псевдослучайных чисел; в k-том цикле используются числа

Используемые числа, (11)

где n≥1 и 0≤t1<t2<…<tN [Примечание 6]. Например, в книге [1] спектральный тест даётся для случая n=1 (точнее, для случая взаимно простых n и m) и (t1,t2,…,tN)=(0,1,…,N-1); описанному выше разбиению на группы соответствуют значения n=N и (t1,t2,…,tN)=(1,2,…,N).

Между прочим, спектральный тест проверяет только "случайность" внутри групп (11), но не контролирует возможную "неслучайность" последовательности групп. Проверка последовательности (8) не выглядит достаточной для этого.

Применяя дискретное преобразование Фурье так, как это описано в книге [1], и имея в виду, что период последовательности (6) равен m/d, где, как и раньше, d=НОД(n,m), получим [Примечание 7]

f(s1,s2,…,sN).

Используя формулу (3) и обозначая s(x), приведём это выражение к виду

f(s1,s2,…,sN).

Подставив выражение (9), получим

f(s1,s2,…,sN).

Поскольку числа (8) - это "перемешанные" числа 0,1,2,…,m/d-1, последнее выражение можно записать также в виде

f(s1,s2,…,sN).

Последнее выражение представляет собой сумму геометрической прогрессии; используя соответствующую формулу [Примечание 8], получим (при e^(-2πis(a)d/m))

f(s1,s2,…,sN)=0,

если s(a)d/m - не целое число, так как в этом случае q^(m/d)=1, и

f(s1,s2,…,sN),

если s(a)d/m - целое число, так как в этом случае q=1.

Из этих вычислений следует, что f(s1,s2,…,sN)=0 во всех случаях, когда число s(a) не делится на m/d [Примечание 9]; если же s(a) делится на m/d, то Abs(f(s1,s2,…,sN))=1.

Таким образом, при использовании элементов линейной конгруэнтной последовательности (2) максимального периода группами (11) наименьшее волновое число равно

νN, (12)

где минимум берётся по всем наборам чисел (s1,s2,…,sN)≠(0,0,…,0), для которых число s(a) делится на m/d. Как и в книге [1], задача сводится к отысканию наименьшего значения квадратичной формы

Квадратичная форма

при целых значениях x1,x2,…,xN, хотя бы одно из которых не равно 0.

В книге [1] качество датчика рекомендуется оценивать по набору чисел [Примечание 10]

СN. (13)

Малые значения СN соответствуют "отсутствию" случайности, большие - "наличию".

Дальнейшие подробности и метод нахождения величин СN можно найти в книге [1], поэтому на этом мы закончим.

Результаты применения спектрального теста к датчику псевдослучайных чисел RandAs продемонстрированы в статье "Датчик псевдослучайных чисел RandAs".

Литература

[1] Д.Кнут. Искусство программирования для ЭВМ. Том2. Получисленные алгоритмы. "Мир", Москва, 1977.

Примечания

1) Необходимо учитывать ещё то обстоятельство, что человек очень легко находит "закономерности" в хаосе случайных точек, поэтому случайность предстаёт перед человеком под маской закономерности. Посмотрите, например, на звёздное небо. Другое дело, что "закономерность", обнаруженная на одном участке случайной последовательности, не работает на другом её участке. [Вернуться]
2) Запись "A=B mod C" означает, что число A равно остатку от деления числа B на число C. [Вернуться]
3) (a+b)^n, где Биномиальный коэффициент - биномиальные коэффициенты (0≤j≤n). Поскольку 0!=1, получаем Крайние коэффициенты. Считаем, что Нулевые коэффициенты при j<0 или j>n. При j>0 удобно пользоваться формулой Биномиальный коэффициент, из которой следует, что Биномиальный коэффициент является многочленом степени j от перемнной n. [Вернуться]
4) Квадратные скобки обозначают целую часть числа: [x] - наибольшее целое число, не превосходящее числа x. [Вернуться]
5) НОД(n,m) - наибольший общий делитель чисел n и m. [Вернуться]
6) Параметры n,t1,t2,…,tN рекомендуется выбирать так, чтобы группы, используемые в программе, не имели общих элементов. [Вернуться]
7) Exp(x)=e^x - показательная функция; i - мнимая единица. [Вернуться]
8) Сумма геометрической прогрессии, если q≠1, и Сумма геометрической прогрессии, если q=1. [Вернуться]
9) Так как числа a и m взаимно просты, то множитель a^t1 для делимости числа s(a) на m/d не играет роли и может быть опущен. [Вернуться]
10) Заметим, что всегда C1=2, поэтому можно считать, что N>1. [Вернуться]

Если у Вас возникнут какие-либо вопросы, Если вопросов будет много, мы, возможно, сделаем на сайте раздел вопросов и ответов.

Hosted by uCoz