Совсем недавно, помогая коллеге в решении вопроса о неповторяемости работы ряда тестов, я в очередной раз натолкнулся на задачу симуляции устройства, генерирующего последовательности случайных чисел. В этой статье я расскажу о том, какие сложности были обнаружены, какой подход к их разрешению был выбран, и как я проверял, что новое решение лучше предыдущего. Сразу отмечу, что вопрос создания и верификации случайных последовательностей очень тонкий, и почти никакое решение не может считаться в нём окончательным. Буду признателен за комментарии и исправления.
Вначале я кратко расскажу о предметной области аппаратных и псевдослучайных генераторов, об их характеристиках и требованиях к ним. Затем перейду к своему частному случаю — области программной симуляции вычислительных систем, и к конкретной задаче, которую нужно было решить.
Иногда самый надёжный способ получить случайное число — взять его из справочника. Источник изображения: www.flickr.com/photos/nothingpersonal/337684768/
Введение: ГСЧ, ГПСЧ, RNG, DRNG
Задача получения «случайной» последовательности чисел с помощью компьютера на удивление сложна. Интуитивно это можно обосновать тем фактом, что философия цифровых систем и архитектура современных ЭВМ направлены на достижение абсолютной точности, на изничтожение даже намёка на неопределённость в результате работы. Хранимые и передаваемые данные заверяются контрольными суммами, электрические цепи дублируются. Случайность изгнана, так как же её вернуть? И зачем?
На самом деле, многие практические задачи требуют для своего решения реализовывать алгоритмы, использующие достаточно длинные случайные последовательности. Имитационное моделирование, генетические алгоритмы, численное интегрирование и решение систем линейных уравнений (методы Монте Карло), компьютерные игры и современная криптография. В каждом случае имеется собственное понимание того, какие последовательности уже можно считать случайными, а какие ещё нет. Различаются и используемые генераторы случайных чисел (ГСЧ, англ. RNG — random number generator).
Казалось бы, природа предоставляет достаточно хаоса, куда тут человечеству за ним угнаться. Достаточно сделать оцифровывающее устройство для некоторого процесса (тепловые шумы, вспышки на солнце, кидание монетки), подключить его к компьютеру, и дело в шляпе.
Однако практическое применение «чистых» аппаратных ГСЧ ограничено следующими факторами. Во-первых, скорость генерации чисел может быть недостаточно высокой, а приложения могут потреблять их так быстро, что даже современные интерфейсы типа PCI-Express будут загружены под завязку. Во-вторых, статистические характеристики такой последовательности могут быть не самые лучшие — распределение не то или меняется со временем, ширина чисел меньше требуемой, тенденция выдавать нули чаще, чем единицы (англ. bias) и т.п. В-третьих, цена подобного качественного устройства может быть достаточно высокой, что будет ограничивать его область применения областями с самыми жёсткими требованиями к случайности. Наконец, энергопотребление и размеры ГСЧ могут быть неподходящими для мобильных устройств.
По этой причине часто используется комбинированный подход: первоначальное значение (seed) берётся из случайного источника, и вся последовательность генерируется из него с помощью неслучайной функции:
(a1, s1) = f(seed)
(a2, s2) = f(s1)
(a3, s3) = f(s2)
…
Здесь a1, a2,… — выходная псевдослучайная последовательность, а s1, s2,… — внутреннее состояние генератора. Результирующая конструкция называется генератор псевдослучайных чисел (ГПСЧ, англ. PRNG — pseudo-random number generator).
Качество последовательности, очевидно, зависит от выбора функции f(), а иногда также от значения seed (например, f(x) = a*x mod m выдаёт последовательность нулей при seed = 0). История знает множество случаев использования плохих функций. Дам лишь два примера.
- ai+1 = 65539 * ai mod 231, seed — нечётный. Это классическая RANDU, широко использовавшаяся в UNIX в 1960-1970 годах. Проваливает даже простейшие тесты (о них будет сказано далее).
- Дональд Кнут в [5] описывает свою попытку в 1959 году сделать супер-ГПСЧ с помощью алгоритма из 13 шагов, один запутаннее другого. На практике оказалось, что последовательность немедленно сошлась к одному числу, которое затем преобразовывалось само в себя.
Источники случайности в вычислительных платформах
Вернёмся к аппаратным ГСЧ. В архитектуре IBM PC они присутствуют уже некоторое время и предлагаются как в компонентах Intel, так и других вендоров. Несколько примеров ниже.
- В конце 20 века аппаратный ГСЧ был встроен в Intel 82802 Firmware Hub [3]. Работал он исправно, однако к нему были претензии по энергопотреблению, скорости генерации, а также по совместимости со всё время уменьшающимися размерами техпроцесса.
- По этой причине был спроектирован новый ГСЧ, находящийся прямо в ядре ЦПУ. Он был представлен в процессорах микроархитектуры Ivy Bridge в качестве инструкции RDRAND [4]. В дальнейшем он также начал использоваться для инструкции RDSEED. О различиях между этими двумя командами написано здесь. Пара статей на Хабре про RDRAND: habrahabr.ru/post/128666/ habrahabr.ru/post/145188/
- В современных платформах собственный ГСЧ находится в модуле TPM (англ. trusted platform module) и может быть использован как для нужд криптографии, так и для других целей. Подробнее о скорости работы и качестве выдаваемых чисел рассказывается в [8].
- Компания VIA также предлагает аппаратный ГСЧ в составе своих процессоров VIA C3 [6].
Общие требования к ГСЧ
Так как я работаю над программными моделями аппаратуры, возникает задача моделирования присутствующих в них ГСЧ. Модель, т.е. в конечном счёте реализующая генератор функция на Си, должна удовлетворять требованиям, выдвигаемым на физические ГСЧ. Если же это не так, то отличия должны быть поняты и объяснены особенностями процесса моделирования — нельзя оставлять такой вопрос на самотёк. В конце концов, всякая модель — это лишь приближение реальности, важно, чтобы оно было достаточно хорошо для практических нужд, а его ограничения были осознаны.
- Случайность. Самый расплывчатый критерий, и к нему мы ещё вернёмся. Но уже понятно, что стоит с умом подойти к выбору функции и проверить её тестами.
- Распределение. Будем считать, что требуемое распределение — равномерное U[0; 2N -1], где N=64. Современные генераторы выдают именно столько бит.
- Большой период. Всякая последовательность от ГПСЧ зацикливается из-за ограниченности числа различных внутренних состояний. Это не будет заметно, если её период будет очень велик. Наша цель — величина порядка 263 или выше.
- Скорость генерации. RDRAND может выдавать десятки миллионов чисел в секунду. Симуляция обычно работает значительно медленнее реального железа. Не будем гнаться за большой скоростью генерации, однако и слишком сложные функции по возможности выбирать не стоит — не хочется создавать bottleneck.
- Неугадываемость. Этот аспект важен для криптографических применений. Наша модель не используется для таких нужд. Кроме того, криптографически стойкие программные ГПСЧ обычно очень медленны.
Есть ещё одно требование, специфичное для симуляции.
- Повторяемость. Состояние симуляции в любой момент её работы может быть сохранено на диск в точку сохранения (англ. checkpoint), перенесено на другую машину и там восстановлено. Кроме того, необходимо поддержать возможность обращения течения времени, и ход всех процессов, в том числе выходящих из ГПСЧ чисел, должен повторяться. Это означает, что симуляция должна хранить значение внутреннего состояния генератора.
Проблемы с rand()
Самый первый и вообще никуда не годящийся генератор выглядел примерно так:
// init seed
int seed = ...;
srand(seed);
uint64 bad_rand() {
return rand();
}
Здесь нет даже намёка на повторяемость — хотя seed явно задаётся в начале работы, и последовательность будет всегда одна и та же, при откате симуляции до предыдущего состояния она не будет повторяться, так как состояние генератора в процессе работы нам неподконтрольно. Кроме того, у bad_rand() есть и другие проблемы, аналогичные тем, что описаны далее. Мы не будем больше к нему возвращаться в дальнейшем.
Затем код для генерации 64 битного случайного числа был переписан так:
uint64 unix_rand(uint64 *state) {
srand(*state);
uint32 r1 = rand();
srand(r1);
uint32 r2 = rand();
*state = r2;
uint64 result = ((uint64)r2 << 32) | r1;
return result;
}
Поскольку стандартный rand() из Libc возвращает int, который в худшем случае имеет ширину 32 бита, использовались два последовательных вызова и комбинация их в одно 64-битное число. Так как rand() не возвращает значение своего внутреннего состояния, приходилось задавать его каждый раз с помощью srand().
Этот подход всё ещё был плохим по ряду причин.
- Зависимость старшей половины бит от младшей. В качестве seed для старших 32 бит выступают младшие, т.е. части числа связаны функциональной зависимостью, которую можно обнаружить, т.е. показать неслучайность.
- RAND_MAX на моей системе был равен 231-1. Т.е. на самом деле rand() выдаёт только 31 бит! А в конечном числе меняются только 62 бита, 31-й и 63-й же всегда нули.
- Самый серьёзный недостаток — неопределённость функции rand(). Вообще нельзя сказать, будет ли RAND_MAX одинаковым на всех системах. Нельзя даже утверждать, что реализация rand() будет одинаковой на Windows и Linux, не будет зависеть от разрядности ОС, версии Libc, компилятора или чего-то ещё. Одна и та же симуляция, запущенная на разных хозяйских системах, будет выдавать разные результаты. Кошмар для отладки!
Новое решение
Наконец, мы привели код к следующему виду:
const uint64 a1 = ...;
const uint64 c1 = ...;
const uint64 a2 = ...;
const uint64 c2 = ...;
uint64 two_lcg_rnd(uint64 *state1, uint64 *state2) {
uint64 r1 = a1 * *state1 + c1;
uint64 r2 = a2 * *state2 + c2;
*state1 = r1;
*state2 = r2;
uint64 result = (r2 & 0xffffffff00000000ull) | (r1 >> 32);
return result;
}
- Использованы два линейных конгруэнтных генератора (англ. LCG) вида f(x) = a*x + c mod 264. Сильные и слабые стороны этого типа LCG хорошо изучены, а вычисления в целых числах по модулю степени двойки выполняются быстро.
- Константы a1, a2, c1, c2 выбраны так, чтобы периоды LCG были велики. Интересующиеся легко найдут их значения в Интернете.
- Для того, чтобы генераторы были независимы, a1 != a2, c1 != c2, а в качестве внутреннего состояния используются две независимые переменные. Но эти предосторожности не гарантирует независимости; требуется эмпирическая проверка.
- У каждого из возвращаемых LCG 64-битных чисел используются только старшие 32 бита. Младшие разряды LCG имеют значительно меньший период, поэтому они отброшены.
Использование собственной функции для ГПСЧ избавляет нас от самой трудноуловимой ошибки — недетерминизма симуляции. Но насколько случайным является её выдача?
Эмпирическая проверка
Принцип эмпирического анализа Г(П)СЧ — вычисление для входной последовательности статистики, определяющей истинность некоторой гипотезы о характере её поведения. Например, можно проверять равнораспределённость пар (троек, четвёрок...) соседних чисел, отсутствие периодов монотонного роста и спада и т.п. Для истинного ГСЧ ни одна из подобных гипотез не должна быть подтверждена с достаточной уверенностью. Чем больше различных статистик проходит конкретный ГПСЧ, тем выше его качество. Интересующимся сравнением качества ГСЧ из разных популярных программ я рекомендую прочитать [1].
Существует несколько пакетов для исследования ГСЧ, в том числе Dieharder [2] и NIST SP 800-22 [7]. Я использовал пакет TestU01 [1]. В нём предопределены три группы тестов, различающихся по числу проверяемых в них статистик: SmallCrush, Crush и BigCrush. Каждый из входящих в них подтестов выдаёт число p-value из полуинтервала [0,1). Для прохождения подтеста оно должно быть одновременно далеко как нуля, так и от единицы.
Я использовал средний по размеру набор Crush, который во всех моих экспериментах исполнялся от получаса до двух часов. Программы я запускал на такой системе:
ЦПУ: Intel® Core(TM) i5-4300U CPU @ 1.90GHz
uname -a: CYGWIN_NT-6.3 hostname 1.7.31(0.272/5/3) 2014-07-25 11:26 x86_64 Cygwin
Для второго генератора, основанного на rand(), результаты плачевны: 139 тестов из 144 закончились неудачно:
========= Summary results of Crush =========
Version: TestU01 1.2.3
Generator: unix_rand
Number of statistics: 144
Total CPU time: 00:46:32.40
The following tests gave p-values outside [0.001, 0.9990]:
(eps means a value < 1.0e-300):
(eps1 means a value < 1.0e-15):
Test p-value
----------------------------------------------
1 SerialOver, t = 2 eps
2 SerialOver, t = 4 eps
3 CollisionOver, t = 2 eps
4 CollisionOver, t = 2 eps
5 CollisionOver, t = 4 eps
6 CollisionOver, t = 4 eps
7 CollisionOver, t = 8 eps
8 CollisionOver, t = 8 eps
9 CollisionOver, t = 20 eps
10 CollisionOver, t = 20 eps
11 BirthdaySpacings, t = 2 eps
12 BirthdaySpacings, t = 3 eps
13 BirthdaySpacings, t = 4 eps
14 BirthdaySpacings, t = 7 eps
15 BirthdaySpacings, t = 7 eps
16 BirthdaySpacings, t = 8 eps
17 BirthdaySpacings, t = 8 eps
18 ClosePairs NP, t = 2 3.2e-157
18 ClosePairs mNP, t = 2 3.2e-157
18 ClosePairs mNP1, t = 2 eps
18 ClosePairs mNP2, t = 2 eps
18 ClosePairs NJumps, t = 2 eps
19 ClosePairs NP, t = 3 3.2e-157
19 ClosePairs mNP, t = 3 3.2e-157
19 ClosePairs mNP1, t = 3 eps
19 ClosePairs mNP2, t = 3 eps
19 ClosePairs NJumps, t = 3 eps
19 ClosePairs mNP2S, t = 3 eps
20 ClosePairs NP, t = 7 1.8e-79
20 ClosePairs mNP, t = 7 1.8e-79
20 ClosePairs mNP1, t = 7 eps
20 ClosePairs mNP2, t = 7 eps
20 ClosePairs NJumps, t = 7 eps
20 ClosePairs mNP2S, t = 7 eps
21 ClosePairsBitMatch, t = 2 6.9e-65
22 ClosePairsBitMatch, t = 4 7.5e-143
23 SimpPoker, d = 16 eps
24 SimpPoker, d = 16 eps
25 SimpPoker, d = 64 eps
26 SimpPoker, d = 64 eps
27 CouponCollector, d = 4 eps
28 CouponCollector, d = 4 eps
29 CouponCollector, d = 16 eps
30 CouponCollector, d = 16 eps
31 Gap, r = 0 eps
32 Gap, r = 27 eps
33 Gap, r = 0 eps
34 Gap, r = 22 eps
35 Run of U01, r = 0 eps
36 Run of U01, r = 15 eps
37 Permutation, r = 0 eps
38 Permutation, r = 15 eps
39 CollisionPermut, r = 0 eps
40 CollisionPermut, r = 15 eps
41 MaxOft, t = 5 eps
41 MaxOft AD, t = 5 3.2e-157
42 MaxOft, t = 10 eps
42 MaxOft AD, t = 10 1.8e-79
43 MaxOft, t = 20 eps
43 MaxOft AD, t = 20 1 - eps1
44 MaxOft, t = 30 eps
44 MaxOft AD, t = 30 1 - eps1
45 SampleProd, t = 10 1 - eps1
46 SampleProd, t = 30 1 - eps1
47 SampleMean eps
48 SampleCorr 1 - eps1
49 AppearanceSpacings, r = 0 1 - eps1
50 AppearanceSpacings, r = 20 1 - eps1
51 WeightDistrib, r = 0 eps
52 WeightDistrib, r = 8 eps
53 WeightDistrib, r = 16 eps
54 WeightDistrib, r = 24 eps
55 SumCollector eps
56 MatrixRank, 60 x 60 eps
57 MatrixRank, 60 x 60 eps
58 MatrixRank, 300 x 300 eps
60 MatrixRank, 1200 x 1200 eps
62 Savir2 eps
63 GCD, r = 0 eps
64 GCD, r = 10 eps
65 RandomWalk1 H (L = 90) eps
65 RandomWalk1 M (L = 90) eps
65 RandomWalk1 J (L = 90) eps
65 RandomWalk1 R (L = 90) eps
65 RandomWalk1 C (L = 90) eps
66 RandomWalk1 H (L = 90) eps
66 RandomWalk1 M (L = 90) eps
66 RandomWalk1 J (L = 90) eps
66 RandomWalk1 R (L = 90) eps
66 RandomWalk1 C (L = 90) eps
67 RandomWalk1 H (L = 1000) eps
67 RandomWalk1 M (L = 1000) eps
67 RandomWalk1 J (L = 1000) eps
67 RandomWalk1 R (L = 1000) eps
67 RandomWalk1 C (L = 1000) eps
68 RandomWalk1 H (L = 1000) eps
68 RandomWalk1 M (L = 1000) eps
68 RandomWalk1 J (L = 1000) eps
68 RandomWalk1 R (L = 1000) eps
68 RandomWalk1 C (L = 1000) eps
69 RandomWalk1 H (L = 10000) eps
69 RandomWalk1 M (L = 10000) eps
69 RandomWalk1 J (L = 10000) eps
69 RandomWalk1 R (L = 10000) eps
69 RandomWalk1 C (L = 10000) eps
70 RandomWalk1 H (L = 10000) eps
70 RandomWalk1 M (L = 10000) eps
70 RandomWalk1 J (L = 10000) eps
70 RandomWalk1 R (L = 10000) eps
70 RandomWalk1 C (L = 10000) eps
71 LinearComp, r = 0 1 - eps1
71 LinearComp, r = 0 1 - eps1
73 LempelZiv 1 - eps1
74 Fourier3, r = 0 eps
75 Fourier3, r = 20 eps
76 LongestHeadRun, r = 0 eps
76 LongestHeadRun, r = 0 1 - eps1
77 LongestHeadRun, r = 20 eps
77 LongestHeadRun, r = 20 1 - eps1
78 PeriodsInStrings, r = 0 eps
79 PeriodsInStrings, r = 15 eps
80 HammingWeight2, r = 0 eps
82 HammingCorr, L = 30 eps
83 HammingCorr, L = 300 eps
84 HammingCorr, L = 1200 eps
85 HammingIndep, L = 30 eps
86 HammingIndep, L = 30 eps
87 HammingIndep, L = 300 eps
88 HammingIndep, L = 300 eps
89 HammingIndep, L = 1200 eps
90 HammingIndep, L = 1200 eps
91 Run of bits, r = 0 eps
91 Run of bits, r = 0 eps
92 Run of bits, r = 20 eps
92 Run of bits, r = 20 1 - eps1
93 AutoCor, d = 1 1 - eps1
94 AutoCor, d = 1 eps
95 AutoCor, d = 30 1 - eps1
96 AutoCor, d = 10 1 - 5.2e-10
----------------------------------------------
All other tests were passed
Для нового генератора, основанного на двух независимых LCG, картина лучше:
========= Summary results of Crush =========
Version: TestU01 1.2.3
Generator: two_lcg_rnd
Number of statistics: 144
Total CPU time: 00:46:13.21
The following tests gave p-values outside [0.001, 0.9990]:
(eps means a value < 1.0e-300):
(eps1 means a value < 1.0e-15):
Test p-value
----------------------------------------------
1 SerialOver, t = 2 eps
2 SerialOver, t = 4 eps
8 CollisionOver, t = 8 1 - 8.0e-13
13 BirthdaySpacings, t = 4 1.2e-82
15 BirthdaySpacings, t = 7 eps
16 BirthdaySpacings, t = 8 eps
17 BirthdaySpacings, t = 8 eps
----------------------------------------------
All other tests were passed
Провалились лишь 7 тестов из 144. Неплохо, хотя для криптографических нужд я этот ГПСЧ не посоветую использовать. Отмечу, что отдельные «половинки» этого генератора также не прошли те же самые типы подтестов.
Из любопытства я прогнал TestU01 (самый длинный BigCrush) на последовательности, выдаваемой инструкцией RDRAND на том же ЦПУ:
========= Summary results of BigCrush =========
Version: TestU01 1.2.3
Generator: rdrand
Number of statistics: 160
Total CPU time: 15:58:06.92
The following tests gave p-values outside [0.001, 0.9990]:
(eps means a value < 1.0e-300):
(eps1 means a value < 1.0e-15):
Test p-value
----------------------------------------------
1 SerialOver, r = 0 eps
2 SerialOver, r = 22 eps
76 RandomWalk1 J (L=1000, r=0) 1.7e-4
----------------------------------------------
All other tests were passed
Поиск ответа на вопрос, почему даже на «истинно» случайном генераторе TestU01 выдаёт несколько неудач, оставим для будущих исследователей, хотя у меня есть гипотеза на этот счёт.
Заключение
Тема генерации псевдослучайных чисел и проверки их на случайность невероятно захватывающая. Всем заинтересовавшимся я рекомендую почитать, если вы ещё этого не сделали, соответствующую главу из второго тома Д. Кнута [5]. Вопрос использования ГСЧ для нужд симуляции, и не только вычислительных систем, разбирается в [9].
Литература
- Pierre L'Ecuyer and Richard Simard. 2007. TestU01: A C library for empirical testing of random number generators. ACM Trans. Math. Softw. 33, 4, Article 22 (August 2007). DOI=10.1145/1268776.1268777 simul.iro.umontreal.ca/testu01/tu01.html
- Robert G. Brown, Dirk Eddelbuettel, David Bauer. Dieharder: A Random Number Test Suite Version 3.31.1 www.phy.duke.edu/~rgb/General/dieharder.php
- Intel Corporation. Intel® 810 Chipset Design Guide, June 1999. download.intel.com/design/chipsets/designex/29065701.pdf
- Intel® Digital Random Number Generator (DRNG) Software Implementation Guide — software.intel.com/en-us/articles/intel-digital-random-number-generator-drng-software-implementation-guide
- Дональд Кнут. Искусство программирования. Том 2. Получисленные алгоритмы. Третье издание — 2011 — Вильямс. ISBN 978-5-8459-0081-4,
- Cryptography Research, Inc. Evaluation summary: VIA C3 Nehemiah random number generator — 2003 — www.cryptography.com/public/pdf/VIA_rngsum.pdf
- National Institute of Standards and Technology. A Statistical Test Suite for Random and Pseudorandom Number Generators for Cryptographic Applications — 2010 csrc.nist.gov/publications/nistpubs/800-22-rev1a/SP800-22rev1a.pdf
- Alin Suciu, Tudor Carean. Benchmarking the True Random Number Generator of TPM Chips. arxiv.org/ftp/arxiv/papers/1008/1008.2223.pdf
- Handbook of Simulation. Principles, Methodology, Advances, Applications, and Practice / под ред. J. Banks. — John Wiley & Sons, Inc., 1998. — ISBN 0-471-13403-1. — URL: books.google.com/books?id=dMZ1Zj3TBgAC
Автор: Atakua
