[Перевод] Генератор случайных чисел, который можно запустить в голове

Уровень сложности Сложный Время на прочтение 8 мин Количество просмотров 4.2K Блог компании RUVDS.com Программирование *Алгоритмы *Математика * Обзор Перевод Автор оригинала: Hillel Wayne

Люди ужасно плохо справляются с придумыванием случайных чисел. Я хотел научиться быстро генерировать «достаточно случайные» числа. Мне не нужно было что-то совершенное, просто способ придумывания случайных цифр за полминуты. Поискав онлайн, я нашёл старый пост в Usenet, написанный Джорджем Марсалья:

Выберите двухразрядное число, допустим, 23. Оно будет вашим «порождающим значением» (seed).

Создайте новое двухразрядное число: количество десяток плюс шесть, умноженное на количество единиц.

Пример последовательности: 23 –> (2 + 6 * 3) = 20 –> (2 + 6 * 0) = 02 –> 12 –> 13 –> 19 –> 55 –> 35 –> …

Его период будет порядком множителя (6) в группе остатков, простых относительно модуля, 10 (в данном случае 59).

«Случайными цифрами» будет количество единиц двухразрядных чисел, то есть 3,0,2,2,3,9,5,… то есть члены последовательности mod 10.

Больше всего Марсалья известен своим набором тестов diehard-генераторов случайных чисел (RNG), так что он в этом понимает (здесь и далее под RNG я имею в виду генератор псевдослучайных чисел (PRNG)). Мне стало любопытно, почему это работает и как он выбрал 6.

Мы будем писать на Raku, языке для гремлинов. На случай, если вы тоже гремлин, под спойлерами я буду объяснять все странные особенности.

▍ Введение

Последовательность периодична, то есть если мы итеративно будем применять её, то рано или поздно получим тот же элемент. Давайте начнём с функции («подпрограммы»), создающей всю последовательность:

my sub next-rng(Int $start, Int $unitmult = 6, —> List) { my @out; my $next = $start; repeat while $next !(elem) @out { @out.append($next); $next = sum($next.polymod(10) Z* $unitmult,1); }; return @out; }
Объяснение Raku — чрезвычайно странный язык, но я постараюсь объяснить максимально понятно.

• @ и $ — это сигилы для «позиционных» (напоминающих списки) и «скалярных» переменных. Если определить позиционную переменную, не присвоив ей ничего, то по умолчанию она становится пустым массивом.
• (elem) проверяет на членство, а ! можно использовать для отрицания любого инфиксного оператора. (elem) — это ASCII-версия, Raku может принимать и ∈.
• polymod разделяет число на остаток и частное, например, 346.polymod(10) = (6 34). Этот метод может принимать несколько параметров, так что можно выполнять операции наподобие num.polymod(60, 60, 24), чтобы получать часы-минуты-секунды.
• Z — это метаоператор «zip», применяющий инфиксный оператор поэлементно между двумя списками. (4, 6) Z* (6, 1) = (4*6, 6*1).

Получив последовательность, мы можем вывести её при помощи команд put или say, поведение которых немного различается. Я не буду вдаваться в подробности.

put next-rng(1);
01 06 36 39 57 47 46 40 04 24 26 38 51 11 07 42 16 37 45 34 27 44 28 50 05 30 03 18 49 58 53 23 20 02 12 13 19 55 35 33 21 08 48 52 17 43 22 14 25 32 15 31 09 54 29 56 41 10 01
Помните, что нужные нам случайные числа — это последняя цифра. То есть RNG работает так: 1 -> 6 -> 6 -> 9 -> 7 -> 7 -> …

▍ Исследуем свойства

Если RNG равномерный, то каждая цифра должна встречаться в последовательности одинаковое количество раз. Мы можем проверить это, преобразовав последние цифры в мультимножество, или «bag».

say bag(next-rng(1) <<%>> 10);
Объяснение

• <<op>> — это гипероператор, «отображающий» внутренний оператор по обоим спискам, также рекурсивно входя в список списков. Например, результатом ((1, 2), 3) <<+>> 10 будет ((11, 12), 13)! Гипероператоры имеют множество других странных свойств, которые делают их и полезными, и непонятными.
• Объекты класса Bag подсчитываю количество элементов в чём-то. bag((4, 5, 4)){4} = 2. Как ни странно, они могут содержать только скалярные переменные, но не массивы, списки или что-то подобное.

Bag(0(5) 1(6) 2(6) 3(6) 4(6) 5(6) 6(6) 7(6) 8(6) 9(5))
Похоже, распределение достаточно равномерное, хотя вероятность получить 0 или 9 чуть ниже.

Свою следующую идею я взял из тестов diehard. Цитата из Википедии:

Пересекающиеся перестановки (Overlapping Permutations)

Анализируются последовательности пяти последовательных случайных чисел. 120 возможных перестановок должны получаться со статистически эквивалентной вероятностью.

В датасете всего 54 пятичисловых последовательностей, так что вместо этого я применю тест к двухчисловым «переходам». Сделаю я это, выводя сетку 10 х 10, где (i, j)-тый индекс (от 0) — это количество переходов от i последней цифры к j. Например, последовательность содержит переход 28 -> 50 и никакие другие переходы вида X8 -> Y0, поэтому в ячейке (8, 0) должно быть значение 1.

sub successions-grid(@orbit) { my @pairs = (|@orbit , @orbit[0]).map(* % 10).rotor(2 => -1); for ^10 -> $x {put ($x X ^10).map({@pairs.grep($_).elems})} }
Объяснение

• | @f, $x конкатенирует $x напрямую в @f. Без | это был бы список из двух элементов.
• * в `* % 10` — это whatever, странный маленький оператор, выполняющий множество действий в разных контекстах, но обычно в этом случае он преобразует выражение в замыкание. Обычно. Это то же самое, что написать map({$_ % 10})1.
• rotor(2 => -1) получает два элемента, затем возвращается на один элемент назад, затем получает ещё два и так далее. Результатом [1, 2, 3, 4].rotor(2 => -1) будет [(1, 2), (2, 3), (3, 4)]. Также можно выполнить rotor(2), чтобы получить [(1, 2), (3, 4)], или rotor(1 => 1), чтобы получить [1, 3]. Rotor — очень крутая штука.
• ^10 — это всего лишь 0..9. Наконец-то что-то простое!
• X — это метаоператор векторного произведения. То есть если $x = 2, то $x X ^4 даст ((2 0), (2 1), (2 2), (2 3)). И да, этот оператор может стать намного страньше.
• grep(foo) возвращает список всех элементов, выполняя умное сравнение foo, а .elems — это количество элементов в списке. То есть @pairs.grep($_).elems — это количество элементов списка, соответствующих $_. Чтобы разобраться в этом, мне понадобилось слишком много времени

1. На самом деле, -> $x {$x % 10}, но это достаточно близко

> successions-grid(next-rng(1, 6)) 0 1 1 1 1 1 0 0 0 0 1 1 0 0 0 0 1 1 1 1 0 0 1 1 1 1 1 1 0 0 1 1 1 1 0 0 0 0 1 1 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 1 1 0 0 0 0 1 1 1 1 0 0 1 1 1 1 1 1 0 0 1 1 1 1 0 0 0 0 1 1 0 0 0 0 1 1 1 1 1 0
Из таблицы видно, что некоторые переходы невозможны. Если я сгенерировал 0, то сразу после него не смогу получить 6. Очевидно, это неидеальный RNG, но чего-то особо качественного я и не ждал.

▍ Почему 6?

Что, если я умножу последнюю цифру не на 6, а на 4?

> say next-rng(1, 4); 01 04 16 25 22 10
Ну, не знаю, вообще, мне нравится RNG, который никогда не выдаёт 3. Конкретные последовательности называются орбитами, а их длины — периодами. Давайте рассмотрим все возможные орбиты, которые можно получить при использовании 4 в качестве множителя:

sub orbits-for-mod(int $mult, $top = 20) { my &f = &next-rng.assuming(*, $mult); (1..$top).map(&f).unique(as => &set) }
Объяснение

• & — это сигил для «вызываемого объекта» или функции подпрограммы. Метод .assuming выполняет частичное применение функции, а передача * заставляет его частично применить второй параметр.1
• map возвращает последовательность списков, которую мы передаём unique. as => &set преобразует каждую последовательность из map в множество и сравнивает на уникальность их, а не исходные списки. Но в окончательном результате элементы используются до преобразования.

  Если это непонятно, то можно привести более простой пример: [-1, 1].unique(as => &abs) возвращает [-1], а [1, -1].unique(as => &abs) возвращает [1].

1. Дэниел Соквелл (codesections) любезно согласился прочитать первый черновик этого поста и сообщил мне об assume. Спасибо ему!

> say orbits-for-mod(4, 38).map(*.gist).join(«n»); [1 4 16 25 22 10] [2 8 32 11 5 20] [3 12 9 36 27 30] [6 24 18 33 15 21] [7 28 34 19 37 31] [13] [14 17 29 38 35 23] [26]
Объяснение Цитирую Дэниела Соквелла:

.map(*.gist).join(«n») здесь используется только для красивого вывода. cycles-for-mod возвращает Seq из Array; преобразование каждого Array при помощи .gist превращает его в строку, окружённую квадратными скобками, а .join(«n») помещает переход на новую строку после каждой из строк.

Если выбрать в качестве начального значения 13, то случайными числами будет последовательность 3, 3, 3, 3, 3, 3.

По очевидным причинам 4 никогда не должно быть множителем. На самом деле, чтобы множитель выдавал «хороший» RNG, он должен иметь ровно одну орбиту.

> say (1..30).grep(*.&orbits-for-mod == 1) (3 6 11 15 18 23 27)
Объяснение

• .& применяет функцию высшего порядка в качестве метода. grep(*.&f) — это эквивалент grep({f($_)}).
• &orbits-for-mod возвращает список. == приводит оба входных значения в числа, а приведение списка в число возвращает количество элементов. То есть мы проверяем, содержит ли возвращённый список только один элемент (иначе говоря, только одну орбиту). (Если вы не хотите выполнять сравнение без приведения, то используйте или ===, или eqv.)

Такой способ довольно медленный и проверяет только орбиты, начинающиеся с числа вплоть до 20. Поэтому он упустит орбиту 26 -> 26 для x=4. Обе эти проблемы я устраню позже.
Поэтому «хорошие» варианты для n — это 6, 11 и 18.

Стоит отметить, что если получится число из трёх разрядов, то с первыми двумя разрядами нужно обращаться как с одним числом. В случае n=11 число 162 приводит к значению 16 + 22, а не к 6 + 22 (или к 6 + 1 + 22).

▍ Почему это работает?

Для меня была непонятна вот эта часть объяснения:

Его период будет порядком множителя (6) в группе остатков, простых относительно модуля, 10 (в данном случае 59).

Поговорив с друзьями и почитав Википедию, я начал что-то понимать. Мы вычисляем в голове RNG умножения с переносом (multiply with carry, MWC) с константами a=x и c=10. Этот выбор имеет удобное свойство: если MWC(x) = y, то 10y mod (10n-1) = x!

MWC: 01 -> 06 -> 36 -> … -> 41 -> 10 -> 01 10y%59: 01 -> 10 -> 41 -> … -> 36 -> 06 -> 01
И это очень здорово! Мне проще математически рассуждать об 10y mod 59, чем об «умножении последней цифры на шесть и прибавлении первой цифры». Например, понятно, почему RNG генерирует 0 и 9 немного реже, чем другие цифры: какой бы множитель мы ни выбрали, генерируемая последовательность будет идти от 1 до 10n-2, «исключая» 10n-1 (которое заканчивается на 9) и 10n (которое заканчивается на 0).

«Умножение и деление с остатком» также называют RNG Лемера.

▍ Ищем более качественные RNG

А как работают другие числа? Мы знаем, что хороший множитель создаёт только одну орбиту, а выше я показал код для вычисления этого. К сожалению, этот алгоритм в худшем случае имеет сложность O(n²). [Если множитель n имеет одну орбиту, тогда мы выполним next-rng примерно для 10n-2 чисел, и функция выполнит итерации 10n-2 раз (потому что ей нужно обойти каждое число в орбите). Если бы я сделал так, чтобы код пропускал числа, которые я уже видел в орбите, то время исполнения сократилось бы до O(n).] Если воспринимать алгоритм MWC, как «Лемера наоборот», то можно получить более качественный способ: если n — хороший множитель, то период орбиты, начинающейся с 1, должен быть 10n-2.

Метод Лемера также даёт нам более быстрый способ вычисления орбиты:

sub oneorbit(x) { 10, * *10% (10*x — 1) … 1 }
Объяснение

• Указав x вместо $x в качестве параметра, мы можем использовать в теле x вместо $x.
• … — это оператор последовательности. Он может выполнять множество разных действий, но для нас важно то, что если написать 10, &f … 1, он начнёт с 10, а затем будет применять &f, пока наконец не сгенерирует 1.
• In * *10%[etc]: первая * — это Whatever, а вторая * — обычное умножение. Затем это преобразуется в функцию -> $a {$a * 10 % (10*x — 1)}.

На самом деле этот код создаёт орбиту в обратном порядке, но нас интересует только период, так что это не важно.
Затем мы проверяем период при помощи того же трюка с ==, что и раньше.

> say (1..100).grep({oneorbit($_) == 10*$_-2}); (2 3 6 11 15 18 23 27 38 39 42 50 51 62 66 71)
Теперь мне понятно, почему Марсалья выбрал 6: большинство программистов знает таблицу умножения на 6 и понимает, что она никогда не вернётся к трёхразрядному числу, так что этап сложения очень прост. Орбита состоит всего из 58 чисел, и мы не получим последовательностей из одинаковых чисел, но если вам нужно быстро взять несколько случайных цифр, то это вполне приемлемо.

Если вам нужно больше случайности, то у меня есть пара кандидатов. 50 имеет период 498 и невероятно легко вычисляется. Если последняя цифра чётная, то не нужно выполнять никаких переносов: 238 -> 423!

Тем не менее, последовательность для 50 не выглядит такой же случайной, как другие последовательности. В определённый момент она генерирует 9 чётных чисел, за которыми идут 8 нечётных. Не используйте её для имитации бросков монеты.

Ещё одно интересное число — это 18. Оно имеет приличный период (178) и обладает всеми возможными переходами цифр:

> successions-grid(next-rng(1, 18)) 1 2 2 2 2 2 2 2 1 1 2 2 2 2 2 2 1 1 2 2 2 2 2 2 1 1 2 2 2 2 2 2 1 1 2 2 2 2 2 2 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 1 2 2 2 2 2 2 1 1 2 2 2 2 2 2 1 1 2 2 2 2 2 2 1 1 2 2 2 2 2 2 1 1 2 2 2 2 2 2 2 1
Его недостаток в том, что придётся учить таблицу умножения на 18. Это не очень сложно: я освоил её примерно за десять минут практики. Но я всё равно пока неидеально справляюсь с выполнением всего этапа MWC, но могу стабильно получать следующую случайную цифру примерно каждые пять секунд. Меня это вполне устраивает.

Использованный мной для исследований код на Raku можно посмотреть здесь. Он настроен как CLI, но его можно использовать разными способами; подробности см. в файле.

Скидки, итоги розыгрышей и новости о спутнике RUVDS — в нашем Telegram-канале 🚀

Теги:

• ruvds_переводы
• генератор случайных чисел
• гспч
• random number generator
• raku
• генерация случайных чисел

Хабы:

• Блог компании RUVDS.com
• Программирование
• Алгоритмы
• Математика

Источник