Случайные числа и детерминистичная симуляция

в 7:21, , рубрики: PRNG, rand, randu, simulation, srand, Блог компании Intel, ГПСЧ, ГСЧ, криптография, Программирование, случайные числа

Совсем недавно, помогая коллеге в решении вопроса о неповторяемости работы ряда тестов, я в очередной раз натолкнулся на задачу симуляции устройства, генерирующего последовательности случайных чисел. В этой статье я расскажу о том, какие сложности были обнаружены, какой подход к их разрешению был выбран, и как я проверял, что новое решение лучше предыдущего. Сразу отмечу, что вопрос создания и верификации случайных последовательностей очень тонкий, и почти никакое решение не может считаться в нём окончательным. Буду признателен за комментарии и исправления.

Вначале я кратко расскажу о предметной области аппаратных и псевдослучайных генераторов, об их характеристиках и требованиях к ним. Затем перейду к своему частному случаю — области программной симуляции вычислительных систем, и к конкретной задаче, которую нужно было решить.

Случайные числа и детерминистичная симуляция
Иногда самый надёжный способ получить случайное число — взять его из справочника. Источник изображения: 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). История знает множество случаев использования плохих функций. Дам лишь два примера.

  1. ai+1 = 65539 * ai mod 231, seed — нечётный. Это классическая RANDU, широко использовавшаяся в UNIX в 1960-1970 годах. Проваливает даже простейшие тесты (о них будет сказано далее).
  2. Дональд Кнут в [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().
Этот подход всё ещё был плохим по ряду причин.

  1. Зависимость старшей половины бит от младшей. В качестве seed для старших 32 бит выступают младшие, т.е. части числа связаны функциональной зависимостью, которую можно обнаружить, т.е. показать неслучайность.
  2. RAND_MAX на моей системе был равен 231-1. Т.е. на самом деле rand() выдаёт только 31 бит! А в конечном числе меняются только 62 бита, 31-й и 63-й же всегда нули.
  3. Самый серьёзный недостаток — неопределённость функции 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].

Литература

  1. 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
  2. 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
  3. Intel Corporation. Intel® 810 Chipset Design Guide, June 1999. download.intel.com/design/chipsets/designex/29065701.pdf
  4. Intel® Digital Random Number Generator (DRNG) Software Implementation Guide — software.intel.com/en-us/articles/intel-digital-random-number-generator-drng-software-implementation-guide
  5. Дональд Кнут. Искусство программирования. Том 2. Получисленные алгоритмы. Третье издание — 2011 — Вильямс. ISBN 978-5-8459-0081-4,
  6. Cryptography Research, Inc. Evaluation summary: VIA C3 Nehemiah random number generator — 2003 — www.cryptography.com/public/pdf/VIA_rngsum.pdf
  7. 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
  8. Alin Suciu, Tudor Carean. Benchmarking the True Random Number Generator of TPM Chips. arxiv.org/ftp/arxiv/papers/1008/1008.2223.pdf
  9. 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

Источник

* - обязательные к заполнению поля


https://ajax.googleapis.com/ajax/libs/jquery/3.4.1/jquery.min.js