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

Люди ужасно плохо справляются с придумыванием случайных чисел. Я хотел научиться быстро генерировать «достаточно случайные» числа. Мне не нужно было что-то совершенное, просто способ придумывания случайных цифр за полминуты. Поискав онлайн, я нашёл старый пост в 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;
}

Получив последовательность, мы можем вывести её при помощи команд 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);

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})}
}

> 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)
}

> 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]

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

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

> say (1..30).grep(*.&orbits-for-mod == 1)
(3 6 11 15 18 23 27)

Поэтому «хорошие» варианты для 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
}

Затем мы проверяем период при помощи того же трюка с ==, что и раньше.

> 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, но его можно использовать разными способами; подробности см. в файле.