NUM Числа и алгебра Алгоритм
Генераторы псевдослучайных чисел
Также: PRNG · pseudorandom number generator · ГПСЧ · генератор случайных чисел · линейный конгруэнтный генератор · LCG · RANDU · MINSTD · xorshift · PCG · splitmix64 · вихрь Мерсенна · Mersenne Twister · спектральный тест
Детерминированная программа, которая выдаёт числа, похожие на случайные. Плохой генератор выдаёт себя геометрией: тройки чисел RANDU лежат на 15 плоскостях, младшие биты линейного генератора повторяются с коротким периодом, а по двум числам java.util.Random восстанавливается всё будущее. Как устроены генераторы, как их проверяют и почему листинги этого каталога пользуются splitmix64.
В 1960-х годах в научной библиотеке IBM был генератор случайных чисел RANDU: умножить предыдущее число на 65539 и взять остаток от деления на 231. Им пользовались для моделирования много лет. Ниже — 6000 троек подряд идущих чисел RANDU, поставленных точками в куб. Похоже на облако. Нажмите «Повернуть к плоскостям».
Все тройки лежат на 15 параллельных плоскостях, а между плоскостями — пустота. Любой опыт, которому нужны три случайных числа сразу — точка в объёме, три шага блуждания, три параметра модели, — получает вместо случайных точек 15 листов. У splitmix64 с любого угла облако одинаковое.
Идея
Генератор псевдослучайных чисел — детерминированная программа: состояние и правило шага. Из одного и того же зерна он выдаёт одну и ту же последовательность, и рано или поздно она зациклится: состояний конечное число. Хороший генератор нельзя отличить от настоящей случайности никаким разумным тестом, и до его повторения не дойти. Плохой выдаёт себя структурой — плоскостями, периодичными битами, предсказуемостью.
Что такое зерно и почему в Python работает вихрь Мерсенна, рассказано в главе курса. Здесь — что у генераторов внутри, как их проверяют и как ломают.
Линейный конгруэнтный генератор
Самый старый и самый простой генератор предложил Деррик Лемер в 1949 году: xn+1 = (a·xn + c) mod m. Одно умножение и один остаток на число. Период не больше m, и он равен m ровно тогда, когда c взаимно просто с m, a − 1 делится на все простые делители m, а если m делится на 4, то и a − 1 делится на 4 — так в 1962 году доказали Томас Халл и Алан Добелл. Длинный период ещё ничего не гарантирует: RANDU и его плоскости это показывают.
Тройки подряд идущих чисел RANDU лежат на 15 плоскостях 9x − 6y + z = k·231, k = −5, …, 9.
Множитель 65539 = 216 + 3, поэтому 655392 = 232 + 6·216 + 9 = 232 + 6·65539 − 9. По модулю 231 слагаемое 232 пропадает, и a2 ≡ 6a − 9. Значит, xn+2 ≡ a2xn ≡ 6axn − 9xn ≡ 6xn+1 − 9xn, то есть 9xn − 6xn+1 + xn+2 делится на 231. Все числа лежат между 1 и 231 − 1, поэтому это выражение больше −6·231 и меньше 10·231, и частное k может быть только от −5 до 9 — пятнадцать значений.
Листинг проверяет это на 65 534 тройках подряд: все лежат на плоскостях, и заняты все 15. RANDU — крайний случай, но не исключение. Джордж Марсалья в 1968 году доказал, что у любого линейного генератора t-ки подряд идущих чисел лежат не больше чем на (t!·m)1/t гиперплоскостях. Для m = 231 и троек это около 2344 плоскостей у лучших множителей — и 15 у RANDU.
Спектральный тест
На плоскости то же видно на парах. Возьмём генератор с маленьким модулем 4096 и отметим все пары (xn, xn+1) за полный период. Двигайте множитель: при одних a точки лежат на трёх-четырёх прямых, при других прямых много, и они почти не видны.
Почему всегда есть прямые: если s1 + a·s2 делится на m, то s1xn + s2xn+1 по модулю m — одно и то же число для всех n, и все пары, поделённые на m, лежат на прямых s1u + s2v = целое плюс постоянный сдвиг. Соседние прямые такого семейства стоят на расстоянии 1/√(s12 + s22). Хуже всего семейство с самыми редкими прямыми, то есть кратчайший вектор (s1, s2). Его длину и измеряет спектральный тест, который Роберт Ковэю и Роберт Макферсон предложили в 1967 году. Длиннее (4/3)1/4·√m кратчайший вектор не бывает, поэтому отношение его длины к этой границе — оценка от 0 до 1.
| генератор | a | m | кратчайший вектор | оценка |
|---|---|---|---|---|
| RANDU | 65539 | 231 | (−32765, 32767) | 0,931 |
| MINSTD | 16807 | 231 − 1 | (−16807, 1) | 0,338 |
| MINSTD, 1993 | 48271 | 231 − 1 | (3399, 44488) | 0,896 |
| Numerical Recipes | 1664525 | 232 | (−54907, −43865) | 0,998 |
| java.util.Random | 25214903917 | 248 | (−4824621, 7847617) | 0,511 |
RANDU на плоскости выглядит хорошо: 0,93. Его беда — в объёме, где кратчайший вектор — (9, −6, 1) длиной меньше 11 при возможных примерно 1450. Поэтому тест проводят в нескольких измерениях сразу. MINSTD Стивена Парка и Кита Миллера с множителем 16807 на плоскости слаб: маленькое число x переходит в 16807·x, тоже относительно маленькое. В 1993 году авторы сами предложили множитель 48271, и оценка выросла с 0,34 до 0,90.
Кратчайший вектор находится быстро. Все решения s1 + a·s2 ≡ 0 образуют решётку, порождённую векторами (m, 0) и (−a, 1). Редукция Лагранжа — Гаусса — двумерный алгоритм Евклида: из длинного вектора вычитаем ближайшее кратное короткого и повторяем, пока длинный не перестанет становиться короче. Шагов O(log m), как у Евклида. Это вторая задача статьи. В трёх и больше измерениях то же делает редукция решётки, например алгоритм Ленстры — Ленстры — Ловаса.
Младшие биты
У генераторов с модулем 232 или 248 остаток почти ничего не стоит: это просто отбрасывание старших бит при переполнении. Цена — младшие биты. Ниже картинка из 16 384 чисел подряд: каждая точка — один бит очередного числа.
Если m = 2e, то бит номер k чисел линейного генератора повторяется с периодом не больше 2k+1.
Младшие k + 1 бит числа — это его остаток по модулю 2k+1. Поскольку 2k+1 делит m, остаток следующего числа зависит только от остатка предыдущего: xn+1 mod 2k+1 = (a·(xn mod 2k+1) + c) mod 2k+1. Это тот же генератор с модулем 2k+1, у него не больше 2k+1 состояний, и через столько шагов он обязан повториться.
Младший бит при нечётных a и c просто чередуется: 0, 1, 0, 1. Младший байт повторяется каждые 256 чисел — и в тесте листинга он «слишком ровный»: каждое значение ровно по 256 раз на 65 536 числах, чего у честных случайных чисел не бывает. Однажды это уже сбило замер на этом сайте: в статье об эмпирической сложности длинные числа, собранные из младших байтов такого генератора, вышли периодическими, и умножение Карацубы в Go пошло в 19 раз быстрее, чем на настоящих данных. Поэтому java.util.Random держит 48 бит состояния и отдаёт только старшие 32, а MINSTD берёт простой модуль 231 − 1, у которого этой беды нет.
Генераторы получше
xorshift. Джордж Марсалья в 2003 году предложил генераторы без умножения: состояние перемешивается сдвигами и исключающим ИЛИ — x ^= x << 13, x ^= x >> 17, x ^= x << 5. Это очень быстро, и простые тесты он проходит. Но каждый шаг линеен над битами: любой бит выхода — исключающее ИЛИ каких-то битов состояния. Алгоритм Берлекэмпа — Мэсси по такой последовательности битов находит линейную рекурренцию длины не больше 32 — у xorshift32 линейная сложность 32 вместо ожидаемых 500 на тысяче битов. Франсуа Паннетон и Пьер л’Экуйе в 2005 году показали слабости xorshift и на больших батареях тестов. Лечится это нелинейным шагом на выходе: умножить результат на константу или сложить два состояния.
PCG. Мелисса О’Нил в 2014 году взяла линейный генератор с модулем 264, у которого хороши старшие биты, и поставила на выход перестановку: из старших бит берётся величина сдвига или поворота для остальных. Решётка остаётся внутри, но наружу не видна. PCG — один из генераторов пакета math/rand/v2 в Go; глобальный генератор Go с версии 1.22 — ChaCha8.
splitmix64. Гай Стил, Даг Ли и Кристин Флуд в 2014 году описали генератор, у которого состояние — просто счётчик, растущий на нечётную константу γ, а всё качество — в функции перемешивания: два умножения и три сдвига с исключающим ИЛИ. k-е число равно mix(зерно + k·γ), его можно получить без всех предыдущих, и легко завести независимые потоки для параллельных задач. Им пользуются листинги этого каталога, где Python и Go должны выдать одни и те же случайные числа: пять строк дают одинаковый поток в Python, Go и в сценах на странице, а значит, листинг и сцена показывают один и тот же опыт.
Вихрь Мерсенна. Генератор Макото Мацумото и Такудзи Нисимуры 1998 года с периодом 219937 − 1 стоит в Python, R и многих других системах. Его состояние — 624 32-битных слова, около 2,5 КБ. Он линеен над битами, как xorshift, а каждое выданное число — обратимо перемешанное слово состояния, поэтому по 624 числам подряд состояние восстанавливается целиком, как объясняет глава курса.
Как проверяют генератор
Доказать, что генератор хороший, нельзя: тест может только поймать плохой. Листинг гоняет четыре простых теста по 65 536 числам каждого генератора. Первые три — хи-квадрат: младший байт должен принимать все 256 значений примерно поровну; пары младших четвёрок бит — все 256 сочетаний; тройки старших четырёх бит — все 4096 клеток куба. Отклонение пересчитано в z: у честного генератора |z| почти всегда меньше 3,3; большое положительное — слишком неровно, большое отрицательное — подозрительно ровно. Четвёртый тест — линейная сложность младшего бита на тысяче чисел, у случайных битов около 500.
| генератор | младший байт | пары младших бит | тройки старших бит | линейная сложность |
|---|---|---|---|---|
| RANDU | 276,8 | 825,4 | 60,0 | 1 |
| Numerical Recipes | −33,8 | 503,1 | 1,8 | 2 |
| MINSTD | −1,3 | 0,3 | 1,3 | 500 |
| java.util.Random | −0,9 | −1,5 | −1,6 | 501 |
| xorshift32 | 0,3 | −0,6 | 0,8 | 32 |
| splitmix64 | 0,6 | −0,4 | −0,8 | 499 |
У RANDU все числа нечётные — множитель и зерно нечётные, а модуль — степень двойки, — так что младший бит всегда 1, и тройки в кубе видят плоскости. Генератор из «Numerical Recipes» проваливает младшие биты: байт «слишком ровный», пары повторяются. xorshift32 проходит три теста и проваливает четвёртый. MINSTD, java.util.Random и splitmix64 проходят все четыре — что значит только то, что эти четыре теста их не поймали.
Настоящие батареи тестов больше. TestU01 Пьера л’Экуйе и Ришара Симара, вышедшая в 2007 году, в наборе BigCrush гоняет больше сотни тестов на сотнях миллиардов чисел. Но и батарея не последнее слово. В 1992 году Алан Ферренберг, Дэвид Ландау и Джоанна Вонг моделировали двумерную модель Изинга кластерным алгоритмом Вольфа с генератором R250, который считался хорошим, и получили систематически неверную энергию: алгоритм оказался чувствителен к тонким связям между битами, которых стандартные тесты не замечали. С тех пор проверенный приём — повторить важный расчёт с генератором другого устройства.
Предсказание
Генератор детерминирован: кто знает состояние, знает всё, что будет выдано. java.util.Random прячет 16 младших бит из 48, но по двум числам подряд их можно перебрать: 65 536 вариантов, одно умножение на проверку. Подходит, как правило, ровно один, и дальше каждое число известно заранее. Это первая задача статьи.
Ещё опаснее бедное зерно. В 2008 году выяснилось, что после правки в пакете OpenSSL дистрибутива Debian генератор ключей почти два года получал из всех источников случайности только номер процесса — около 32 тысяч вариантов. Все возможные ключи каждого типа можно было перебрать заранее. Поэтому для ключей, паролей и токенов нужны криптографические генераторы, которые берут случайность у операционной системы: secrets в Python и crypto/rand в Go — подробнее в главе курса о шифрах.
Код
Листинг заводит шесть генераторов и прогоняет на них четыре теста, проверяет плоскости RANDU, считает спектральный тест для пяти настоящих генераторов, восстанавливает состояние java.util.Random по двум числам и показывает, как splitmix64 сразу выдаёт стотысячное число.
Сколько это стоит
Каждое число — O(1): у линейного генератора одно умножение и остаток, у xorshift — три сдвига и три исключающих ИЛИ, у splitmix64 — сложение, два умножения и три сдвига. Память — от 4 байт у xorshift32 и 8 у splitmix64 до 2,5 КБ у вихря Мерсенна; период — от 231 − 2 у MINSTD до 219937 − 1. Спектральный тест на плоскости — O(log m) шагов редукции; в больших измерениях редукция решётки дороже, но для генератора её делают один раз.
Где работает
Моделирование. Метод Монте-Карло съедает миллиарды чисел, и здесь важно качество: плоскости RANDU прямо превращаются в неверные ответы.
Алгоритмы. Перемешивание, выборки из потока, случайный опорный элемент быстрой сортировки, случайные основания хешей.
Тесты. Генераторы тестов к задачам этого каталога создают входы из фиксированного зерна, поэтому каждый запуск даёт те же тесты, и ответы можно проверять.
Игры и генерация. Мир, уровень или раскладка строятся из одного числа-зерна: передал зерно — передал весь мир.
История
Джон фон Нейман в 1940-х годах предложил метод середины квадрата: возвести число в квадрат и взять средние цифры. Метод быстро зацикливается, и фон Нейман это знал. В 1949 году Деррик Лемер описал линейный конгруэнтный генератор, в 1962-м Халл и Добелл нашли условия полного периода. В 1960-х RANDU разошёлся с научной библиотекой IBM; в 1967 году Ковэю и Макферсон предложили спектральный тест, в 1968-м Марсалья доказал теорему о плоскостях, в 1969-м Дональд Кнут собрал теорию генераторов во втором томе «Искусства программирования». Дальше — MINSTD Парка и Миллера в 1988 году, вихрь Мерсенна в 1998-м, xorshift в 2003-м, батарея TestU01 в 2007-м, PCG и splitmix в 2014-м.
Задачи
Первая — предсказать java.util.Random по двум числам. Вторая — спектральный тест на плоскости для тысяч генераторов: редукция Лагранжа — Гаусса.
Следующее число Прямые генератора