Случайные числа и детерминистичная симуляция
Совсем недавно, помогая коллеге в решении вопроса о неповторяемости работы ряда тестов, я в очередной раз натолкнулся на задачу симуляции устройства, генерирующего последовательности случайных чисел. В этой статье я расскажу о том, какие сложности были обнаружены, какой подход к их разрешению был выбран, и как я проверял, что новое решение лучше предыдущего. Сразу отмечу, что вопрос создания и верификации случайных последовательностей очень тонкий, и почти никакое решение не может считаться в нём окончательным. Буду признателен за комментарии и исправления.
Вначале я кратко расскажу о предметной области аппаратных и псевдослучайных генераторов, об их характеристиках и требованиях к ним. Затем перейду к своему частному случаю — области программной симуляции вычислительных систем, и к конкретной задаче, которую нужно было решить.
Иногда самый надёжный способ получить случайное число — взять его из справочника. Источник изображения: 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™ i5–4300U CPU @ 1.90GHzuname -a: CYGWIN_NT-6.3 hostname 1.7.31(0.272/5/3) 2014–07–25 11:26×86_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×60 eps 57 MatrixRank, 60×60 eps 58 MatrixRank, 300×300 eps 60 MatrixRank, 1200×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
