Меню

Линейный конгруэнтный генератор формула



Генерация случайных чисел

Дата создания: 2009-05-03 15:27:38
Последний раз редактировалось: 2012-02-08 06:53:25

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

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

Генерация псевдослучайных чисел

Обычно используется генерация псевдослучайных чисел. Т.е. числа не совсем случайные. Мы уже рассматривали функцию rand(). Если Вы напишите программу использующую данную функцию, то при каждом запуске программы, rand() будет генерировать одну и ту же последовательность чисел. Последовательность сгенерированная rand() определяется начальным числом (seed). Сначала задаётся начальное число, затем, по определённой формуле вичисляются все остальные числа последовательности. Зная начальное число и формулу по которой рассчитываются числа, можно вычислить следующее число.

Такие алгоритмы называются детерминистическими (предопределёнными). Т.е. в них последовательность чисел создаётся на основе начальных данных. При этом последовательность чисел — псевдослучайная: т.е. зная начальное число, можно узнать и все следующие. Но пользователь навряд ли их узнает, числа будут казаться случайными.

Одним из предопределённых (детерминистических) алгоритмов создания случайных чисел является линейно-конгруэнтный.

Установка начального значения (для функции rand(), которую мы использовали до сих пор) выглядит примерно так:

Функция srnd устанавливает начальное значение i.

Линейный конгруэнтный (linear congruential) метод генерации случайных чисел

Существует много методов генерации случайных чисел. Линейно-конгруэнтный лишь один из них. Метод довольно старый — 1950х годов. Разработал его Деррик Лемер.

Для реализации алгоритма необходиом задать четыре параметра:

Диапазон значений m, при этом m > 0.
Множитель a (0 = 2, b >= 1.

Вичислим шестой элемент с помощью этой формулы (мы ведь уже знаем пятый):

X5+1 = (a*1*X5 + (a*1-1)*c/b) % m = (2*1*5 + (2*1-1)*3/1) % 10 = 13 % 10 = 3

Всё росто, правда?

У последовательности созданной с помощью линейного конгруэнтного метода и определённой целыми параметрами m, a, c и X, период равен числу m когда выполняются следующие условия:
1.Наибольший общий делитель c и m равен 1.
2.b — кратно любому простому числу, являющемуся делителем m.
3.Если m кратно 4, тогда b также кратно 4.

Выбор m

Период не может быть больше числа m. Следовательно m должно быть довольно большим.

Примеры (x — целое число):

Очень часто для m выбирают одно из простых чисел Мерсенна [3] . Часто число 2 31 — 1 используется, когда вычисления ведутся с 32-ух битными данными.

Множитель a

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

Инкремент c

Данный параметр можно выбирать довольно произвольно. Очень часто его задают в виде нуля, но при этом уменьшается длина периода и X != 0.

Начальное значение (seed)

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

Вот так вот! Генерация последовательных чисел — это самая настоящая мистификация!

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

Вообще говоря, с помощью арифметических методов нельзя построить по настоящему случайную последовательность чисел. Как невесело шутил наш друг — Джон Фон Нейман:

«Anyone who uses arithmetic methods to produce random numbers is in a state of sin.» (C) Джон Фон Нейман.

Грустно всё это. Вот так живёшь-живёшь. Находясь в состоянии неведения, думаешь: «А вот круто бы было создать генератор случайных чисел на основе линейного когнруэнтного метода!». А оказывается что по настоящему случайную последовательность таким образом создать нельзя. Крах надежд! Депрессия и вот ты уже на пороге первой стадии алкоголизма. Этот груз невозможности реализовать свои желания будет давить на тебя всю оставшуюся жизнь. Что-то я отвлёкся. Продолжим.

Реализация генерации случайных чисел

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

Для последовательности рассмотренной выше, генератор будет выглядеть так:

X объявляется как глобальная переменная:

Это самый примитивный вариант генератора. В реальности такое чудовище использовать нельзя! Но мы будем это делать!

В данной реализации мы никак не отслеживаем возможность переполнения переменных. Вместо типа int используйте тип _int64 — это целый восьмибайтный тип.

Простые числа

Простые числа — числа которые можно разделить без остатка только на единицу и на само число. Например: 11, 3, 7.

Читайте также:  Судовые дизель генераторы конструкция

Код Грея

[1] Код Грея представляет собой последовательность чисел, где следующее число отличается от предыдущей одним битом. Т.е. чтобы угадать последовательность, нужно смотреть не на десятичные числа, а на их двоичные эквиваленты.

Вот возрастающая последовательность чисел от 0 до 7 в двоичном коде:

Используя код Грея, получим такую последовательность:

Простые числа Мерсенна

[3] На данный момент известно 44 числа Мерсенна. Из названия понятно, что данные числа — простые, т.е. делятся только на единицу и на самих себя. Кроме того данные числа должны подчиняться следующей формуле:

Где n также простое число. Для первых девяти простых чисел Мерсенна n следующие: <2, 3, 5, 7, 13, 17, 19, 31, 61>.

Поразрядные логические операции

[2] Мы уже рассматривали логические операции && (И) и || (ИЛИ). Существуют также поразрядные логические операции & (поразрядное И) и | (поразрядное ИЛИ). Работате операция следующим образом:

То есть в данной операции сравниваются соответствующие позиции двух чисел и если они равны, то в результируюем числе в данной позиции пишется 1, в противном случае — ноль.

Степень числа

Другие генераторы случайных чисел

Кроме линейного конгруэнтного генератора существует множество других. Например Вихрь Мерсенна, изобретённый двумя японскими учёными (непомню как из зовут) в 1997. У него очень большой (очень большой это мягко сказано) период. Кстати, вихрь Мерсенна использует линейный конгруэнтный генератор для установления начального значения (seed).

Генераторы случайных чисел применяются при шифровании. Но здесь используются специальные генераторы — криптографически защищённые генераторы псевдослучайных чисел. Например блочный шифр (block cipher), потоковый шифр (stream cipher), который был разработан на основе шифра Вернама.

Упражнения:

Реализуйте генераторы случайных чисел на основе линейного когруэнтного метода. В качестве параметров воспользуйтесь следующими:

  • m = 2 32 ; a = 1664525; c = 1013904223. Данные параметры использовались в примере реализации генератора в одной старой книжке по Фортрану.
  • m = 2 32 ; a= 214013; c = 2531011; Данные параметры используется в реализации метода в Visual C++. При этом берутся старшие биты 30..16. Берутся именно верхние биты, т.к. в линейном конгруэнтном методе нижние биты гораздо менее случайны. Так как мы пока не умеем брать конкретные биты, можете использовать всё число.
  • В качестве m выберите одно из чисел Мерсенна. Начните с малых (2 3 -1,2 5 -1). Попробуйте подставить 2 31 -1 и 2 61 -1.

Данный материал поможет вам понять как работают генераторы случайных чисел. К этому материалу будет продолжение: мы напишем нормальный генератор и рассмотрим различные приёмы работы с генераторами псевдослучайных чисел.

Источник

Конгруэнтные генераторы. Линейные конгруэнтные генераторы.

Все конгруэнтные генераторы используют для получения нового состояния рекуррентную формулу вида \[s_ = f(s_i) \mod m,\] где \(m\) – модуль, являющийся положительной целой степенью простого числа. Выходная последовательность обычно совпадает с последовательностью состояний.

Семейство конгруэнтных генераторов определяется общим видом функции \(f: S\to S\) , например, линейные конгруэнтные генераторы используют линейную функцию \(f\) . Конкретный генератор в семействе определяется конкретными параметрами функции \(f\) .

Линейные конгруэнтные генераторы (ЛКГ)

ЛКГ используют рекуррентную формулу вида \[s_ = a\cdot s_i + c \mod m,\] где \(m\) – модуль, \(s_i\) – состояния, \(a\) – “множитель”, \(c\) – “приращение”. Так же необходимо задать \(s_0\) – начальное состояние.

Конкретный генератор из семейства задаётся значениями \(a\) , \(c\) , \(m\) .

Если \(c = 0\) , такие генераторы называются мультипликативными, или генераторами Лемера.

Длина периода

Удачный выбор параметров ЛКГ даёт известную и достаточно большую длину периода. Однако, на практике, неудачный выбор параметров сделать гораздо проще, чем удачный. Так, например, при \(a=1, c=1\) , ЛКГ превращается в обычный счётчик по модулю \(m\) , который явно не похож на случайную последовательность.

Исторически, неудачный выбор параметров приводил к крайне сомнительным результатам. Одним из примеров является генератор RANDU, разработанный в IBM и широко использовавшийся в 1970-х годах (в частности, компиляторами FORTRAN). В нём использовались параметры \(a=65539=2^<16>+3, c=0, m=2^<31>\) , и чётным начальным значением. Хотя этот метод даёт равномерно распределённые значения, они крайне очевидно не являются независимыми. Действительно, рассмотрим последовательность 3 значений \(\, s_\>\) :

\[s_ = (2^<16>+3)s_ = (2^<16>+3)^2 s_\] \[s_ = (2^<16>+3)s_ = (2^<32>+6 \cdot 2^<16>+9) s_\] Т.к. \(2^ <32>\mod 2^31 = 0\) , \[s_ = (6 \cdot 2^<16>+9) s_\] Учитывая \(s_ = (2^<16>+3) s_\) , можем записать \[s_ = [6 (2^ <16>+ 3) — 9] s_\] \[s_ = 6 s_ — 9 s_\]

То есть каждое третье значение однозначно определяется двумя предыдущими. По сути, это уравнение описывает плоскость в 3-мерном пространстве, и если рассматривать каждые три последовательных числа как координаты точки, то все эти точки будут лежать на такой плоскости (технически, из-за модульной арифметики, это семейство 15 плоскостей)

Читайте также:  Генератор рато 6000 д

Нормированное распределение точек в RANDU

Существует три варианта выбора параметров. Рассмотрим каждый из них отдельно.

\(m\) простое, \(c=0\)

Этот случай соответствует генератору Лемера. Если \(a\) является примитивным элементом мультипликативной группы целых по модулю \(m\) , то период генератора очевидно равен \(m-1\) . Начальное состояние \(s_0\neq 0\) .

Один из недостатков простого модуля в том, что при использовании двоичной арифметики при вычислениях, в отличие от случая \(m=2^n\) , необходим явный шаг вычисления модуля, и двойной объем памяти при умножении.

Второй недостаток заключается в том, что если необходимо получить равномерно распределённую последовательность бит, числа в диапазоне \([1,m)\) не всегда тривиально можно преобразовать.

В качестве \(m\) часто выбираются числа Мерсенна, например \(2^<31>-1\) или \(2^<61>-1\) .

Для вычисления в ограниченном числе разрядов можно применить алгоритм Шража 1 . Представим \(m\) в виде \[m = qa+r,\] \(q = \left\lfloor \frac \right\rfloor\) , \(r = m \mod a\) .

Затем рассчитаем \[as \mod m = a(s \mod q) − r \left\lfloor \frac \right\rfloor\]

Поскольку \(s\mod q , первый член строго меньше чем \(am/a = m\) . Если выбрать \(a\) таким образом, чтобы \(r \le q\) (т.е. \(r/q \le 1\) ), то второй член так же меньше \(m\) : \[r \left\lfloor \frac \right\rfloor \le r\frac ≤ x

Тогда оба этих члена можно вычислить в том же числе разрядов, в которых представимо \(m-1\) , а их разность лежит в пределах \([1−m, m−1]\) , который легко привести к диапазону \([0, m−1]\) , просто прибавляя \(m\) , если первый член меньше второго.

\(m=2^n\) , \(c=0\)

Если \(m=2^n\) , то при использовании двоичной арифметики в \(n\) разрядах, вычисление резко упрощается. Действительно, поскольку беззнаковые вычисления в \(n\) разрядах по сути являются арифметикой \(\mod n\) , для вычисления достаточно одного умножения.

Есть, однако, обратная сторона.

Возьмём ЛКГ \[s_ = as_i \mod 2^n\] и рассмотрим младшие \(l\) бит состояния \(s_i\) : \[t_i = s_i \mod 2^l.\] Тогда \[t_ = (as_i \mod 2^n) \mod 2^l\] \[t_ = at_i \mod 2^l,\] и тогда младшие \(l\) бит тоже образуют линейную конгруэнтную последовательность с меньшим модулем. А поскольку период меньше модуля, то в лучшем случае, младший бит никогда не меняется, второй по старшинству бит меняется каждый раз, и так далее. Таким образом младшие биты оказываются значительно менее “случайными”, чем старшие.

\(c\neq 0\)

При корректном выборе параметров и \(c\neq 0\) , длина периода может достигать \(m\) .

Теорема Халла-Добелла

ЛКГ, задаваемый рекуррентным соотношением \[s_ = a s_i + c \mod m,\] \(c \neq 0\) достигает периода \(m\) iff:

  1. \(m\) и \(c\) взаимно простые
  2. \(a-1\) делится на все простые множители \(m\)
  3. \(a-1\) делится на 4, если \(m\) делится на 4.

Лемма 1

Пусть \(p\) – простое число, а \(e\) – положительное целое число, такое, что \(p^e > 2\) . Если \[x = 1 \mod p^e,\] \[x\neq 1 \mod p^,\] то \[x^p = 1 \mod p^,\] \[x^p \neq 1 \mod p^.\]

\(\Delta\) : По условию теоремы, \(x = 1 + qp^e\) для некого целого \(q\) . Тогда \[x^p = (1+qp^e)^p = 1 + C_p^1 qp^e + \ldots + C_p^q^p^+q^pp^,\] где \(C_n^k = \frac\) – биномиальный коэффициент. \[x^p = 1 + qp^(1 + C_p^2 q p^ + \ldots + C_p^q^p^+q^ p^).\] Выражения вида \[C_p^k q^ p^\] должны делиться на \(p^\) . Действительно, для \(1 \(C_p^k\) делится на \(p\) .

Последний член \(q^ p^\) так же делится на \(p\) , поскольку \(e(p-1) > 1\) при \(p^e>2\) (что легко проверить перебором, при \(e=1\) , \(p>2\) , при \(e=2\) , \(p\ge 2\) ).

Таким образом, все члены в скобках кроме \(1\) делятся на \(p\) и тогда \[x^p = 1+qp^ \mod p^.\] \(\not \Delta\)

Пусть число \(m\) раскладывается на простые множители \[m=p_1^ \ldots p_t^.\] Тогда длинна периода \(\lambda\) линейной конгруэнтной последовательности, определяемой параметрами \((s_0, a, c, m)\) , является наименьшим общим кратным длин \(\lambda_j\) периодов линейных конгруэнтных последовательностей \((s_0 \mod p_j^, a \mod p_j^, c \mod p_j^, p_j^)\) , \(1 \le j \le t\) .

\(\Delta\) : Используя индукцию по \(t\) , достаточно доказать лемму для \((s_0, a, c, m_1m_2)\) , где \(m_1\) и \(m_2\) – взаимно простые. Тогда \(\lambda\) будет НОК \(\lambda_1\) последовательности \((s_0 \mod m_1, a \mod m_1, c \mod m_1, m_1)\) и \(\lambda_2\) последовательности \((s_0 \mod m_2, a \mod m_2, c \mod m_2, m_2)\) . Отметим, что если элементы этих последовательностей обозначены соответственно \(x_i\) , \(y_i\) , \(z_i\) , то справедливы равенства \[y_n = x_n \mod m_1,\] \[z_n = x_n \mod m_2.\] Действительно, \[y_0 = x_0 \mod m_1,\] \[y_1 = (a y_0 + c) \mod m_1 = (a x_0 + c) \mod m_1,\] по индукции, \[y_n = (a x_ + c) \mod m_1 = x_n \mod m_1.\]

Аналогично для \(z_n\) .

Тогда \(\forall k\) , \[x_n = x_k \mod m_1m_2 \iff\] \[y_n = y_k \mod m_1,\] \[z_n = z_k \mod m_2.\]

Читайте также:  Усилитель заряда для генератора

Пусть теперь \(\lambda’\) – НОК( \(\lambda_1\) , \(\lambda_2\) ). Докажем, что \(\lambda’ = \lambda\) . Так как \(x_n = x_\) для всех достаточно больших \(n\) , то \(y_n = y_\) и \(z_n = z_\) , следовательно, \(\lambda\) кратно \(\lambda_1\) и \(\lambda_2\) , и следовательно \(\lambda \ge \lambda’\) . С другой стороны, \(y_n = y_\) и \(z_n = z_\) для достаточно больших \(n\) и следовательно \(x_n = x_\) , и тогда \(\lambda \le \lambda’\) .

Из леммы 2 следует, что теорему достаточно доказать для случая, когда \(m\) – степень простого числа, поскольку \[m = p_1^ \ldots p_t^ = \lambda = НОК(\lambda_1,\ldots,\lambda_t)\le\lambda_1\ldots\lambda_t\le p_1^ \ldots p_t^\] iff \(\lambda_j = p_j^\) для \(1\le j \le t\) .

Поэтому положим \(m=p^e\) , где \(p\) – простое, а \(e\) – натуральное. Для \(a=1\) утверждение теоремы достаточно очевидно, поэтому рассмотрим случай \(a>1\) . Период имеет длину \(m\) iff каждое целое \(0\le x встречается в последовательности ровно один раз. Действительно, никакое значение не может встретиться более одного раза и значений \(\mod m\) ровно \(m\) . Тогда без ограничения общности можем положить \(s_0 = 0\) . Рассмотрим общую формулу для \(s_n\) : \[\begin s_ & = as_+c \mod m \\& = a^2s_+ac+c \mod m \\& = a^3s_+a^2c+ac+c \mod m \\& = \ldots \\& = a^ks_+c\sum_^a^k \mod m \\& = a^ks_+c\frac \mod m. \end\] Тогда при \(s_0=0\) имеем \[s_n = c\frac \mod m.\] Если \(c\) и \(m\) не являются взаимно простыми, то \(s_n \neq 1\) , то есть условие (1) теоремы является необходимым. Период по определению есть наименьшее положительное \(n\) , для которого \(s_n = s_0 = 0\) . Тогда теорема сводится к следующей лемме:

Пусть \(1 , где \(p\) – простое число. Если \(\lambda\) есть наименьшее натуральное число, для которого \[\frac = 0 \mod p^e,\] то \(\lambda = p^e\) iff \(a=1 \mod p\) при \(p>2\) и \(a=1\mod 4\) при \(p=2\) .

Положим \(\lambda=p^e\) и при этом \(a\neq 1 \mod p\) . Тогда \[\frac = 0 \mod p^e,\] iff \[a^\lambda-1 = 0 \mod p^e,\] и \[a^\lambda = a^ = 1 \mod p^e,\] но по теореме Ферма \(a^n = a \mod n\) , что приводит к противоречию. Если \(p=2\) и \(a\neq 1 \mod 4\) , то \(a=3\mod 4 \neq 1 \mod 2^e\) для \(e>1\) и \[a^ <2^e>= 1 \mod 2^e.\]

Следовательно, чтобы \(\lambda=m\) , условия (2,3) теоремы необходимы. То есть, необходимо, чтобы \(a=1+qp^f\) , где \(p^f>2\) и \(q\) не кратно \(p\) .

Остаётся доказать, что это условие достаточно. Применяя лемму 1, находим, что \(\forall g\ge 0\) : \[a^ = 1 \mod p^,\] \[a^ \neq 1 \mod p^,\] и тогда \[\frac-1> = 0 \mod p^, \tag<*>\] \[\frac-1> \neq 0 \mod p^,\tag<**>\] в частности, это верно для \(g=e\) . Для последовательности \((0, a, 1, p^e)\) , справедливо \[s_n = \frac \mod p^e.\] Но это значит, что период этой последовательности есть \(\lambda\) . То есть \(s_n = 0\) iff \(n\) кратно \(\lambda\) . Но поскольку \[\frac-1> = 0 \mod p^,\] то \(p^e\) кратно \(\lambda\) . Поскольку \(p\) – простое, это возможно только если \(\lambda=p^g\) для некоторого \(g\) . Но \(\lambda\) – наименьшее число, для которого верно \[\frac = 0 \mod p^e,\] и тогда индуктивные рассуждения для \((*)\) и \((**)\) показывают, что \(\lambda = p^e\) . \(\not \Delta\)

Следствия теоремы

Из доказательства теоремы Халла-Добелла прямо следуют так же выражения максимального периода для случая \(c=0\) .

Действительно, если \(c=0\) , то \(x_n = a^n x_0 \mod m\) . По лемме 2, период любой ЛКП зависит от длин последовательностей при \(m=p^e\) . Если \(x_n = a^n x_0 \mod p^e,\) и \(a\) кратно \(p\) , период будет иметь длину не более \(e\) . Поэтому рассмотрим случай, когда \(a\) и \(p\) – взаимно простые. Тогда период \(\lambda \in \mathbb N\) – это наименьшее число, такое, что \[x_0 = a^\lambda x_0 \mod p^e.\] Если теперь \(p^f = НОД(x_0, p^e)\) , то это требование эквивалентно \[a^\lambda = 1 \mod p^.\] По теореме Эйлера, \(a^ <\varphi(m)>= 1 \mod m\) , поэтому \(\lambda\) является делителем \(\varphi(p^)\) , \[\varphi(p^) = p^(p-1).\] Eсли \(x_0\) и \(p^e\) взаимно простые, \(p^f=1\) , \(f=0\) и \[\varphi(p^) = p^(p-1).\]

Это позволяет определить \(\lambda(m)\) : \[\lambda(m) = \left\lbrace\begin 1, && m=2 \\ 2, && m=4 \\ 2^, && m=2^e, e\ge3 \\ p^(p-1), && m=p^e, p>2 \\ p^(p-1), && m=p^e, p>2 \\ НОК(\lambda(p_1^),\ldots,\lambda(p_n^)), && m = p_1^\ldots p_n^ \end\right.\tag<*>\]

Вышеперечисленные рассуждения можно кратко сформулировать следующим образом.

Максимальным периодом, возможным при \(c=0\) является \(\lambda\) , определённое в \((\text<*>)\) . Этот период достигается, если:

  1. \(x_0\) и \(m\) – взаимно простые
  2. \(a\) является первообразным элементом мультипликативной группы по модулю \(m\) .

В случае \(m=2^e\) , имеем \[\lambda = p^ = m/4,\] а требования на множитель сводятся к \(a = 3 \mod 8\) или \(a=5\mod 8.\)

Источник