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

Генератор случайных чисел Лемера (названный в честь Д. Х. Лемера ), иногда также называемый генератором случайных чисел Парка – Миллера (в честь Стивена К. Парка и Кейта В. Миллера), представляет собой тип линейного конгруэнтного генератора (LCG), который работает в мультипликативной группе целых чисел по модулю. н . Общая формула

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

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

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

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

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

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

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

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

Выбор модуля

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

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

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

Использование составного модуля возможно, но генератор должен быть заполнен значением, взаимно простым с m , иначе период будет значительно уменьшен. Например, модуль F ₅ = 2 ³² + 1 может показаться привлекательным, поскольку выходные данные можно легко преобразовать в 32-битное слово 0 ≤ X _i − 1 < 2 ³² . Однако начальное число X ₀ = 6700417 (которое делит 2 ³² + 1) или любое кратное число приведет к выходу с периодом всего 640.

Более популярная реализация для больших периодов — комбинированный линейный конгруэнтный генератор ; объединение (например, путем суммирования их выходов) нескольких генераторов эквивалентно выходу одного генератора, модуль которого является произведением модулей составляющих генераторов. и чей период является наименьшим общим кратным составляющих периодов. Хотя периоды имеют общий делитель 2, модули можно выбрать так, чтобы это был единственный общий делитель, и результирующий период будет равен ( m ₁ - 1)( m ₂ - 1)···( m _k - 1)/ 2 ^{к -1} . Одним из примеров этого является генератор Вихмана – Хилла .

Отношение к LCG

Хотя ГСЧ Лемера можно рассматривать как частный случай линейного конгруэнтного генератора с c = 0 , это особый случай, который подразумевает определенные ограничения и свойства. В частности, для ГСЧ Лемера начальное начальное число X ₀ должно быть взаимно простым с модулем m , что в целом не требуется для LCG. Выбор модуля m и множителя a также является более ограничительным для ГСЧ Лемера. В отличие от LCG, максимальный период ГСЧ Лемера равен m − 1, и он таков, когда m простое число, а a является примитивным корнем по модулю m .

С другой стороны, дискретные логарифмы (по основанию a или любому примитивному корню по модулю m ) X _k в $\mathbb {Z} _{m}$ представляют собой линейную конгруэнтную последовательность по модулю Эйлера. $\varphi (m)$.

Выполнение

Простой модуль требует вычисления произведения двойной ширины и явного шага приведения. Если используется модуль чуть меньше степени 2 ( простые числа Мерсенна 2 ³¹ − 1 и 2 ⁶¹ − 1 популярен, как и 2 ³² − 5 и 2 ⁶⁴ − 59), приведение по модулю m = 2 ^и − d можно реализовать дешевле, чем обычное деление двойной ширины, используя тождество 2 ^и ≡ d (мод м ) .

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

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

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

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

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

0 ≤ у < м = 2 ^и - д

В этом случае y + d < 2 ^и и y ′ = y < m , по желанию.

м ≤ у < 2 м

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

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

большего продукта Уменьшение оси до менее 2 м = 2 ^{и +1} − 2 d может быть выполнено одним или несколькими шагами уменьшения без смещения.

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

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

Метод Шраге

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

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

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

Чтобы использовать метод Шраге, сначала разложите m = qa + r , т.е. предварительно вычислите вспомогательные константы r = m mod a и q = ⌊ m / a ⌋ = ( 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   ⌊ x / q ⌋ ≤ rx / q = x ( r / q ) ≤ x (1 ) знак равно Икс < м . Таким образом, разница лежит в диапазоне [1− m , m −1] и может быть уменьшена до [0, m −1] одним условным добавлением.

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

Этот метод также можно расширить, чтобы обеспечить большее значение a, применяя его рекурсивно. Из двух членов, вычтенных для получения окончательного результата, только второй ( r   ⌊ x / q ⌋ ) рискует переполниться. Но это само по себе является модульным умножением на константу времени компиляции r и может быть реализовано тем же методом. Поскольку каждый шаг в среднем уменьшает вдвое размер множителя (0 ≤ r < a , среднее значение ( a −1)/2), это может потребовать один шаг на бит и быть крайне неэффективным. Однако каждый шаг также делит x на постоянно возрастающее частное q = ⌊ m / a ⌋ , и быстро достигается точка, в которой аргумент равен 0, и рекурсия может быть прекращена.

Пример кода C99

Используя код C , ГСЧ Парка-Миллера можно записать следующим образом:

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 product = (uint64_t)*state * 48271;
	uint32_t x = (product & 0x7fffffff) + (product >> 31);

	x = (x & 0x7fffffff) + (x >> 31);
	return *state = x;
}

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

uint32_t lcg_parkmiller(uint32_t *state)
{
	// Precomputed parameters for Schrage's method
	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;
}

Другой популярный генератор Лемера использует простой модуль 2. ³² −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;

	// Not required because 5 * 279470273 = 0x5349e3c5 fits in 32 bits.
	// product = (product & 0xffffffff) + 5 * (product >> 32);
	// A multiplier larger than 0x33333333 = 858,993,459 would need it.

	// The multiply result fits into 32 bits, but the sum might be 33 bits.
	product = (product & 0xffffffff) + 5 * (uint32_t)(product >> 32);

	product += 4;
	// This sum is guaranteed to be 32 bits.
	x = (uint32_t)product + 5 * (uint32_t)(product >> 32);
	return *state = x - 4;
}

Многие другие генераторы Лемера обладают хорошими свойствами. Следующий модуль-2 ¹²⁸ Генератор Лемера требует 128-битной поддержки со стороны компилятора и использует множитель, вычисленный Л'Экуйером. Имеет период 2 ¹²⁶ :

static unsigned __int128 state;

/* The state must be seeded with an odd value. */
void seed(unsigned __int128 seed)
{
	state = seed << 1 | 1;
}

uint64_t next(void)
{
	// GCC cannot write 128-bit literals, so we use an expression
	const unsigned __int128 mult =
		(unsigned __int128)0x12e15e35b500f16e << 64 |
		0x2e714eb2b37916a5;
	state *= mult;
	return state >> 64;
}

Генератор вычисляет нечетное 128-битное значение и возвращает его старшие 64 бита.

Этот генератор проходит BigCrush от TestU01 , но не проходит тест TMFn от PractRand . Этот тест был разработан, чтобы точно выявить дефект генератора этого типа: поскольку модуль представляет собой степень 2, период младшего бита на выходе равен только 2. ⁶² , а не 2 ¹²⁶ . Линейные конгруэнтные генераторы с модулем степени 2 ведут себя аналогично.

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

uint64_t next(void)
{
	uint64_t result = state >> 64;
	// GCC cannot write 128-bit literals, so we use an expression
	const unsigned __int128 mult =
		(unsigned __int128)0x12e15e35b500f16e << 64 |
		0x2e714eb2b37916a5;
	state *= mult;
	return result;
}

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

Ссылки

• Лемер, Д.Х. (1949). «Математические методы в больших вычислительных устройствах». Материалы второго симпозиума по крупномасштабной цифровой вычислительной технике . стр. 141–146 . МР   0044899 . (журнальная версия: Анналы вычислительной лаборатории Гарвардского университета , Том 26 (1951)).
• Стив Парк, Генераторы случайных чисел

Внешние ссылки

• Простые числа, меньшие степени двойки, могут быть полезны для выбора модулей. Часть Prime Pages .