Генератор случайных чисел Лемера

Генератор случайных чисел Лемера (названный в честь D. H. Lehmer), иногда также называемый генератором случайных чисел Парка-Миллера (после Стивена К. Парк и Кит У. Miller) — тип линейного конгруэнтного генератора (LCG), который работает в мультипликативной группе целых чисел по модулю n. Общая формула является

$X_{k+1}=a\cdot X_{k}{\bmod {m}},$

где модуль m является простым числом или степенью простого числа, множитель a является элементом с высоким мультипликативным порядком по модулю m (например, примитивный корень по модулю n), а семя X0 является coprime к m.

Другими названиями являются мультипликативный линейный конгруэнтный генератор (MLCG) и мультипликативный конгруэнтный генератор (MCG).

Параметры общего использования править

В 1988 году Парк и Миллер предложили ГСЧ Лемера с определенными параметрами m = 2³¹ − 1 = 2 147 483 647 (простое число Мерсенна M31) и a = 7⁵ = 16 807 (примитивный корень по модулю M31), теперь известный как MINSTD. Хотя позже MINSTD подвергся критике со стороны Марсальи и Салливана (1993), он все еще используется сегодня (в частности, в CarbonLib и minstd_rand0) C++11. Парк, Миллер и Стокмайер ответили на критику (1993), сказав:

Учитывая динамичный характер области, неспециалистам трудно принимать решения о том, какой генератор использовать. «Дайте мне что-то, что я могу понять, реализовать и портировать… это не должно быть современным, просто убедитесь, что это достаточно хорошо и эффективно». Наша статья и связанный с ней минимальный стандартный генератор были попыткой ответить на этот запрос. Пять лет спустя мы не видим необходимости изменять наш ответ, кроме как предложить использование множителя a = 48271 вместо 16807.

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

Sinclair ZX81 и его преемники используют ГСЧ Лемера с параметрами m = 2¹⁶ + 1 = 65 537 (простое число Ферма F4) и a = 75 (примитивный корень по модулю F4). Генератор случайных чисел CRAY RANF представляет собой ГСЧ Лемера с модулем мощности двух m = 2⁴⁸ и a = 44 485 709 377 909. Научная библиотека GNU включает в себя несколько генераторов случайных чисел формы Лемера, включая MINSTD, RANF и печально известный генератор случайных чисел IBM RANDU.

Выбор модуля править

Чаще всего модуль выбирается в качестве простого числа, что делает выбор семени coprime тривиальным (сойдет любой 0 < X0 < m). Это дает выход наилучшего качества, но вносит некоторую сложность реализации, и диапазон вывода вряд ли будет соответствовать желаемому применению; преобразование в желаемый диапазон требует дополнительного умножения.

Использование модуля m, который является степенью двойки, создает особенно удобную компьютерную реализацию, но обходится дорого: период не более m/4, а нижние биты имеют периоды короче этого. Это связано с тем, что самые низкие k бит образуют генератор по модулю-2k сами по себе; биты более высокого порядка никогда не влияют на биты более низкого порядка. Значения Xi всегда нечётные (бит 0 никогда не меняется), биты 2 и 1 чередуются (нижние 3 бита повторяются с периодом 2), нижние 4 бита повторяются с периодом 4 и так далее. Поэтому приложение, использующее эти случайные числа, должно использовать наиболее значимые биты; сокращение до меньшего диапазона с использованием операции по модулю с чётным модулем приведет к катастрофическим результатам.

Для достижения этого периода множитель должен удовлетворять ≡ ± 3 (mod 8), а семя X0 должно быть нечётным.

Использование составного модуля возможно, но генератор должен быть засеян значением, равным m, иначе период будет значительно сокращен. Например, модуль F5 = 232 + 1 может показаться привлекательным, так как выходы могут быть легко сопоставлены с 32-битным словом 0 ≤ Xi — 1 < 232. Однако семя X0 = 6700417 (которое делит 232 + 1) или любое кратное приведет к выходу с периодом всего 640.

Более популярной реализацией для больших периодов является комбинированный линейный конгруэнтный генератор; объединение (например, путем суммирования их выходов) нескольких генераторов эквивалентно выходу одного генератора, модуль которого является произведением модулей компонентных генераторов. Хотя периоды будут иметь общий делитель 2, модули могут быть выбраны так, что являются единственным общим делителем и результирующим периодом (m1 — 1) (m2 — 1)···(mk — 1)/2k-1.: 744 Одним из примеров этого является генератор Вихмана-Хилла.

Отношение к LCG править

Хотя Lehmer RNG можно рассматривать как частный случай линейного конгруэнтного генератора с c = 0, это частный случай, который подразумевает определенные ограничения и свойства. В ЧАСТНОСТИ, для Lehmer RNG начальное семя X0 должно быть взаимно простым с модулем m, который не требуется для LCG в целом. Выбор модуля m и множителя a также более ограничен для Lehmer RNG. В отличие от LCG, максимальный период Lehmer RNG равен m — 1, и он такой, когда m является простым и a является примитивным корнем по модулю m.

С другой стороны, (дискретные логарифмы (для основания a или любого примитивного корня по модулю m)) в представляют собой линейную конгруэнтную последовательность по модулю .

Осуществление править

Простой модуль требует вычисления произведения двойной ширины и явного шага уменьшения. Если используется модуль чуть меньше степени 2 (популярны простые числа Мерсенна 231 — 1 и 261 — 1, а также 232 — 5 и 264 — 59), сокращение по модулю m = 2e — d может быть реализовано дешевле, чем общее деление двойной ширины с использованием тождества 2e ≡ d (mod m).

Произведение на две части, умножает верхнюю часть на d и складывает их вместе: (ax mod 2e) + d ⌊ax/2e⌋. Затем m можно вычитать до тех пор, пока результат не окажется в диапазоне. Число ограничено ad/m, которое легко может быть ограниченной единице, если d мало и выбрано a < m/d. (Это условие также гарантирует, что d [2e] и a < m/d). (Это условие также гарантирует, что d ⌊ax/2e⌋ является произведением одинарной ширины; если оно не выполняется, необходимо вычислить произведение двойной ширины).

Когда модуль является простым числом Мерсенна (d = 1), процедура особенно проста. Умножение на d не только тривиально, но и условное вычитание может быть заменено безусловным сдвигом и сложением. Чтобы увидеть это, обратите внимание, что алгоритм гарантирует, что x × 0 (mod m), что означает, что x = 0 и x = m оба невозможны. Это позволяет избежать необходимости рассматривать эквивалентные e-bit представления состояния; только значения, где высокие биты являются ненулевыми, уменьшают потребность.

Низкие биты продукта ax не могут представлять значение больше m, а высокие биты никогда не будут иметь значение больше a — 1 ≤ m — 2. Таким образом, на первом этапе сокращения получается значение не более m + a — 1 ≤ 2m — 2 = 2e+1 — 4. Это (e + 1)-битное число, которое может быть больше m (то есть может иметь битовое множество e), но высокая половина не более 1, и если это так, то низкие биты e будут строго меньше m. Таким образом, независимо от того, равен ли высокий бит 1 или 0, второй шаг сокращения (сложение половин) никогда не переполнит e битов, и сумма будет желаемым значением.

Если d > 1, условного вычитания также можно избежать, но процедура более сложная. Фундаментальная проблема модуля, подобного 232-5, заключается в обеспечении того, чтобы мы производили только одно представление для таких значений, как 1 ≡ 232-4. Решение состоит в том, чтобы временно добавить d, чтобы диапазон возможных значений был от d до 2e — 1, и уменьшить значения, большие, чем биты e, таким образом, чтобы никогда не генерировать представления меньше d. Наконец, вычитание временного смещения даёт желаемое значение.

Начнем с предположения, что у нас есть частично уменьшенное значение y, ограниченное так, что 0 ≤ y < 2m = 2e+1 — 2d. В этом случае один шаг вычитания со смещением даст 0 ≤ y′ = ((y + d) mod 2e) + d ⌊(y + d)/2e⌋ — d < m. Чтобы убедиться в этом, рассмотрим два случая:

0 ≤ y < m = 2^e − d

В этом случае y + d < 2e и y′ = y < m, как и требуется.
m ≤ y < 2m

В этом случае 2e ≤ y + d < 2e+1 является (e + 1)-битовым числом, и ⌊(y + d)/2e⌋ = 1. Таким образом, y′ = (y + d) — 2e + d — d = y — 2e + d = y — m < m, как и требуется. Поскольку старшая часть умножения равна d, сумма равна по крайней мере d, и вычитание смещения никогда не приводит к переполнению.

(Для конкретного случая генератора Лемера нулевое состояние или его изображение y = m никогда не будет иметь место, поэтому смещение d — 1 будет работать точно так же, если это более удобно. Это уменьшает смещение до 0 в простом случае Мерсенна, когда d = 1.)

Сокращение более крупного топора продукта до менее чем 2m = 2e + 1 — 2d может быть выполнено одним или несколькими шагами сокращения без смещения.

Если ad ≤ m, то достаточно одного дополнительного шага сокращения. Поскольку x < m, ax < am ≤ (a — 1)2e, и один шаг сокращения преобразует это не более чем в 2e — 1 + (a — 1)d = m + ad — 1. Это в пределах 2 м, если ad — 1 < m, что является первоначальным предположением.

Если ad > m, то на первом этапе сокращения можно получить сумму, превышающую 2m = 2e+1 — 2d, что слишком много для конечного этапа сокращения. (Это также требует умножения на d, чтобы получить произведение больше, чем биты e, как упоминалось выше.) Однако до тех пор, пока d2 < 2e, первое сокращение будет производить значение в диапазоне, необходимом для применения в предыдущем случае двух этапов сокращения.

Метод Шраге править

Если произведение двойной ширины недоступно, метод Шраге, также называемый методом аппроксимированной факторизации, может быть использован для вычисления ax mod m, но это обходится дорого:

Модуль должен быть представлен в виде целого числа со знаком; арифметические операции должны допускать диапазон ±m. Выбор множителя a ограничен. Мы требуем, чтобы m mod a ≤ γm/a, обычно достигается путем выбора a ≤ √m. Требуется одно деление (с остатком) на итерацию.

Хотя этот метод популярен для портативных реализаций на языках высокого уровня, в которых отсутствуют операции двойной ширины: 744 на современных компьютерах деление на константу обычно реализуется с использованием умножения двойной ширины, поэтому этого метода следует избегать, если проблема в эффективности. Даже в высокоуровневых языках, если множитель a ограничен √m, то ax с двойной шириной можно вычислить, используя два умножения с одной шириной, и уменьшить, используя методы, описанные выше.

Для использования метода Шраге, первый множитель m = qa + r, то есть предварительно вычислить вспомогательные константы r = m mod a и q = qm/a q = (m-r)/a. Затем каждую итерацию вычисляем ax ≡ a(x mod q) — r Χx/qΧ (mod m).

Равенство существует потому, что

{\begin{aligned}aq&=a\cdot (m-r)/a\\&=(a/a)\cdot (m-r)\\&=(m-r)\\&\equiv -r{\pmod {m}}\\\end{aligned}}

мы разберем x = (x mod q) + q∠x/q, мы получим:

{\begin{aligned}ax&=a(x{\bmod {q}})+aq\lfloor x/q\rfloor \\&=a(x{\bmod {q}})+(m-r)\lfloor x/q\rfloor \\&\equiv a(x{\bmod {q}})-r\lfloor x/q\rfloor {\pmod {m}}\\\end{aligned}}

Причина, по которой он не переполняется, заключается в том, что оба термина меньше m. Поскольку x mod q < q ≤ m/a, первый член строго меньше am/a = m и может быть вычислен с помощью произведения одной ширины.

Если a выбрано так, что r ≤ q (и, следовательно, r/q ≤ 1), то второй член также меньше m: r ≤ rx/q = x(r/q) ≤ x(1) = x < m. Таким образом, разница лежит в диапазоне [1-m, m-1] и может быть уменьшена до [0, m-1] с помощью одного условного сложения.

Этот метод может быть расширен, чтобы позволить отрицательное r (-q ≤ r < 0), изменив конечное сокращение на условное вычитание.

Метод также может быть расширен, чтобы позволить большее a, применяя его рекурсивно: 102 Из двух вычитанных членов для получения конечного результата только второй (r × x / q) рискует переполнить. Но это само по себе модульное умножение на константу времени компиляции r, и может быть реализовано той же техникой. Поскольку каждый шаг в среднем удваивает размер множителя (0 ≤ r < a, среднее значение (a-1)/2), это, по-видимому, потребует одного шага на бит и будет впечатляюще неэффективным. Тем не менее, каждый шаг также делит х на постоянно растущее частное q = qm/a, и быстро достигается точка, где аргумент равен 0 и рекурсия может быть завершена.

Пример кода C99 править

Используя C-код, Park-Miller RNG можно записать следующим образом:

uint32_t lcg_parkmiller(uint32_t *state)
{
   return *state = (uint64_t)*state * 48271 % 0x7fffffff;
}

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

Чтобы избежать 64-битного деления, сделайте сокращение вручную:

uint32_t lcg_parkmiller(uint32_t *state)
{
   uint64_t продукт = (uint64_t)*состояние * 48271;
   uint32_t x = (продукт и 0x7fffffff) + (продукт " 31);
   x = (x & 0x7fffffff) + (x " 31);
   return *state = x;
}

Чтобы использовать только 32-битную арифметику, используйте метод Шраге:

uint32_t lcg_parkmiller(uint32_t *state)
{
   // Предварительно вычисленные параметры для метода Шраге
   const uint32_t M = 0x7fffffff;
   const uint32_t A = 48271;
   const uint32_t Q = M / A; // 44488
   const uint32_t R = M % A; // 3399
   uint32_t div = *state / Q; // max: M / Q = A = 48 271
   uint32_t rem = *state % Q; // max: Q — 1 = 44 487
   int32_t s = rem * A; // max: 44,487 * 48,271 = 2,147,431,977 = 0x7fff3629
   int32_t t = div * R; // max: 48,271 * 3,399 = 164,073,129
   int32_t result = s — t;

   if (result < 0)
      result += M;
   return *state = result;
}

или использовать два 16×16-битных умножения

uint32_t lcg_parkmiller(uint32_t *state)
{
   const uint32_t A = 48271;
   uint32_t low = (*state & 0x7fff) * A; // max: 32,767 * 48,271 = 1,581,695,857 = 0x5e46c371
   uint32_t high = (*state > 15) * A; // max: 65 535 * 48 271 = 3 163 439 985 = 0xbc8e4371
   uint32_t x = low + ((high & 0xffff) « 15) + (high» 16); // max: 0x5e46c371 + 0x7fff8000 + 0xbc8e = 0xde46ffff
   x = (x & 0x7fffffff) + (x > 31);
   return *state = x;
}

Другой популярный генератор Лемера использует простой модуль 232-5:

uint32_t lcg_rand(uint32_t *state)
{
   return *state = (uint64_t)*state * 279470273u % 0xfffffffb;
}

Это также может быть записано без 64-битного деления

uint32_t lcg_rand(uint32_t *state)
{
   uint64_t product = (uint64_t)*state * 279470273u;
   uint32_t x;

   // Не требуется, потому что 5 * 279470273 = 0x5349e3c5 помещается в 32 бита.
   // product = (product & 0xffffffff) + 5 * (product " 32);
   // Множитель, превышающий 0x33333333 = 858,993,459, нуждается в этом.
   // Результат умножения укладывается в 32 бита, но сумма может быть 33 бита.
   product = (product & 0xffffffff) + 5 * (uint32_t) (product " 32);
   product += 4;

   // Эта сумма гарантированно равна 32 битам.
   x = (uint32_t)product + 5 * (uint32_t)(product > 32);
   return *state = x — 4;
}

Многие другие генераторы Lehmer обладают хорошими свойствами. Следующий генератор Лемера по модулю 2128 требует 128-битной поддержки от компилятора и использует множитель, вычисленный L’Ecuyer. [1] Он имеет период 2126:

static unsigned __int128 state;

// Состояние должно быть засеяно нечётным значением.
void seed(unsigned __int128 seed)
{
	state = seed << 1 | 1;
}

uint64_t next(void)
{
    // GCC не может писать 128-битные литералы, поэтому мы используем выражение
	const unsigned __int128 mult =
		(unsigned __int128)0x12e15e35b500f16e << 64 |
		0x2e714eb2b37916a5;
	state *= mult;
	return state >> 64;
}

Генератор вычисляет нечётное 128-битное значение и возвращает верхние 64 бита. Этот генератор проходит тест BigCrush от TestU01, но не проходит тест TMFn от PractRand. Этот тест был разработан, чтобы точно уловить дефект этого типа генератора: поскольку модуль равен мощности 2, период наименьшего бита на выходе составляет всего 262, а не 2126. Линейные конгруэнтные генераторы с модулем мощности 2 имеют аналогичное поведение.

Следующая основная процедура улучшает скорость вышеупомянутого кода для целочисленных рабочих нагрузок (если объявление константы разрешено оптимизировать из цикла вычислений компилятором):

uint64_t next(void)
{
	uint64_t result = state >> 64;
    // GCC не может писать 128-битные литералы, поэтому мы используем выражение
	const unsigned __int128 mult =
		(unsigned __int128)0x12e15e35b500f16e << 64 |
		0x2e714eb2b37916a5;
	state *= mult;
	return result;
}

Однако, поскольку умножение отложено, оно не подходит для хеширования, так как первый вызов просто возвращает верхние 64 бита начального состояния.