Конгруэнтность

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

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

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

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

ЛКГ используют рекуррентную формулу вида \[s_{i+1} = 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_k, s_{k+1}, s_{k+2}\}\):

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

\(c\neq 0\)

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

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

Теорема

ЛКГ, задаваемый рекуррентным соотношением \[s_{i+1} = 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^{e+1},\] то \[x^p = 1 \mod p^{e+1},\] \[x^p \neq 1 \mod p^{e+2}.\]

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

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

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

Лемма 2

Пусть число \(m\) раскладывается на простые множители \[m=p_1^{e_1} \ldots p_t^{e_t}.\] Тогда длинна периода \(\lambda\) линейной конгруэнтной последовательности, определяемой параметрами \((s_0, a, c, m)\), является наименьшим общим кратным длин \(\lambda_j\) периодов линейных конгруэнтных последовательностей \((s_0 \mod p_j^{e_j}, a \mod p_j^{e_j}, c \mod p_j^{e_j}, p_j^{e_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_{n-1} + 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+\lambda}\) для всех достаточно больших \(n\), то \(y_n = y_{n+\lambda}\) и \(z_n = z_{n+\lambda}\), следовательно, \(\lambda\) кратно \(\lambda_1\) и \(\lambda_2\), и следовательно \(\lambda \ge \lambda'\). С другой стороны, \(y_n = y_{n+\lambda_1}\) и \(z_n = z_{n+\lambda_2}\) для достаточно больших \(n\) и следовательно \(x_n = x_{n+\lambda'}\), и тогда \(\lambda \le \lambda'\).

\(\not\Delta\)

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

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

Лемма 3

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

\(\Delta\):

Положим \(\lambda=p^e\) и при этом \(a\neq 1 \mod p\). Тогда \[\frac{a^\lambda-1}{a-1} = 0 \mod p^e,\] iff \[a^\lambda-1 = 0 \mod p^e,\] и \[a^\lambda = a^{p^e} = 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^{p^g} = 1 \mod p^{f+g},\] \[a^{p^g} \neq 1 \mod p^{f+g+1},\] и тогда \[\frac{a^{p^g}-1}{a-1} = 0 \mod p^{g}, \tag{*}\] \[\frac{a^{p^g}-1}{a-1} \neq 0 \mod p^{g+1},\tag{**}\] в частности, это верно для \(g=e\). Для последовательности \((0, a, 1, p^e)\), справедливо \[s_n = \frac{a^n-1}{a-1} \mod p^e.\] Но это значит, что период этой последовательности есть \(\lambda\). То есть \(s_n = 0\) iff \(n\) кратно \(\lambda\). Но поскольку \[\frac{a^{p^e}-1}{a-1} = 0 \mod p^{e},\] то \(p^e\) кратно \(\lambda\). Поскольку \(p\) – простое, это возможно только если \(\lambda=p^g\) для некоторого \(g\). Но \(\lambda\) – наименьшее число, для которого верно \[\frac{a^\lambda-1}{a-1} = 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^{e-f}.\] По теореме Эйлера, \(a^{\varphi(m)} = 1 \mod m\), поэтому \(\lambda\) является делителем \(\varphi(p^{e-f})\), \[\varphi(p^{e-f}) = p^{e-f-1}(p-1).\] Eсли \(x_0\) и \(p^e\) взаимно простые, \(p^f=1\), \(f=0\) и \[\varphi(p^{e}) = p^{e-1}(p-1).\]

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

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

Теорема²

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

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

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

---

1. Schrage algorithm↩︎

2. Gauss, C.F. Disquisitiones Arithmeticæ (1801), §90–92↩︎