Генератор случайных чисел Лемера [1] (названный в честь Д. Х. Лемера ), иногда также называемый генератором случайных чисел Парка – Миллера (в честь Стивена К. Парка и Кейта В. Миллера), представляет собой тип линейного конгруэнтного генератора (LCG). который работает в мультипликативной группе целых чисел по модулю n . Общая формула
где модуль m является простым числом или степенью простого числа , множитель a является элементом высокого мультипликативного порядка по модулю m (например, примитивный корень по модулю n ), а семя X 0 взаимно просто с m .
Другие названия — мультипликативный линейный конгруэнтный генератор (MLCG) [2] и мультипликативный конгруэнтный генератор (MCG) .
В 1988 году Парк и Миллер [3] предложили ГСЧ Лемера с конкретными параметрами m = 2 31 − 1 = 2 147 483 647 ( простое число Мерсенна M 31 ) и a = 7 5 = 16 807 (примитивный корень по модулю M 31 ), теперь известный как МИНСТД . Хотя MINSTD позже подвергся критике со стороны Марсальи и Салливана (1993), [4] [5] он все еще используется сегодня (в частности, в CarbonLib и C++11 ) minstd_rand0
. Парк, Миллер и Стокмейер ответили на критику (1993), [6] сказав:
Учитывая динамичный характер данной области, неспециалистам сложно принять решение о том, какой генератор использовать. «Дайте мне что-нибудь, что я смогу понять, реализовать и портировать… это не обязательно должно быть самым современным, просто убедитесь, что это достаточно хорошо и эффективно». Наша статья и связанный с ней генератор минимальных стандартов были попыткой ответить на этот запрос. Пять лет спустя мы не видим необходимости менять наш ответ, кроме как предложить использовать множитель а = 48271 вместо 16807.
Эта пересмотренная константа используется в генераторе случайных чисел C ++11 .minstd_rand
Sinclair ZX81 и его преемники используют ГСЧ Лемера с параметрами m = 2 16 + 1 = 65 537 ( простое число Ферма F 4 ) и a = 75 (примитивный корень по модулю F 4 ). [7] [8] Генератор случайных чисел CRAY RANF представляет собой ГСЧ Лемера с модулем степени двойки m = 2 48 и a = 44 485 709 377 909. [9] Научная библиотека GNU включает в себя несколько генераторов случайных чисел формы Лемера, включая MINSTD, RANF и печально известный генератор случайных чисел IBM RANDU . [9]
Чаще всего модуль выбирается как простое число, что делает выбор взаимно простого начального числа тривиальным ( подойдет любое 0 < X 0 < m ). Это обеспечивает наилучшее качество вывода, но вносит некоторую сложность реализации, и диапазон вывода вряд ли будет соответствовать желаемому приложению; преобразование в нужный диапазон требует дополнительного умножения.
Использование модуля m , который является степенью двойки , обеспечивает особенно удобную компьютерную реализацию, но имеет свою цену: период не превышает m /4, а младшие биты имеют периоды короче этого. Это связано с тем, что младшие k бит сами по себе образуют генератор по модулю 2 k ; биты старшего порядка никогда не влияют на биты младшего порядка. [10] Значения X i всегда нечетны (бит 0 никогда не меняется), биты 2 и 1 чередуются (младшие 3 бита повторяются с периодом 2), младшие 4 бита повторяются с периодом 4 и так далее. Следовательно, приложение, использующее эти случайные числа , должно использовать старшие биты; уменьшение диапазона до меньшего с использованием операции по модулю с четным модулем приведет к катастрофическим результатам. [11]
Чтобы достичь этого периода, множитель должен удовлетворять a ≡ ±3 (mod 8), [12] и начальное число X 0 должно быть нечетным.
Использование составного модуля возможно, но генератор должен быть заполнен значением, взаимно простым с m , иначе период будет значительно уменьшен. Например, модуль F 5 = 2 32 + 1 может показаться привлекательным, поскольку выходные данные можно легко сопоставить с 32-битным словом 0 ≤ X i − 1 < 2 32 . Однако начальное число X 0 = 6700417 (которое делит 2 32 + 1) или любое другое кратное число приведет к выходу с периодом всего 640.
Более популярная реализация для больших периодов — комбинированный линейный конгруэнтный генератор ; объединение (например, путем суммирования их выходов) нескольких генераторов эквивалентно выходу одного генератора, модуль которого является произведением модулей составляющих генераторов. [13] и чей период является наименьшим общим кратным составляющих периодов. Хотя периоды имеют общий делитель 2, модули можно выбрать так, чтобы это был единственный общий делитель, и результирующий период будет равен ( m 1 - 1)( m 2 - 1)···( m k - 1)/ 2 к -1 . [2] : 744 Одним из примеров является генератор Вихмана-Хилла .
Хотя ГСЧ Лемера можно рассматривать как частный случай линейного конгруэнтного генератора с c = 0 , это особый случай, который подразумевает определенные ограничения и свойства. В частности, для ГСЧ Лемера начальное начальное число X 0 должно быть взаимно простым с модулем m , что вообще не требуется для LCG. Выбор модуля m и множителя a также является более ограничительным для ГСЧ Лемера. В отличие от LCG, максимальный период ГСЧ Лемера равен m − 1, и он таков, когда m простое число, а a является примитивным корнем по модулю m .
С другой стороны, дискретные логарифмы (по основанию a или любого примитивного корня по модулю m ) X k in представляют собой линейную конгруэнтную последовательность по модулю Эйлера .
Простой модуль требует вычисления произведения двойной ширины и явного шага приведения. Если используется модуль чуть меньше степени 2 ( простые числа Мерсенна 2 31 − 1 и 2 61 − 1 популярны, а также 2 32 − 5 и 2 64 − 59), приведение по модулю m = 2 e − d может реализовать дешевле, чем обычное деление двойной ширины с использованием тождества 2 e ≡ d (mod m ) .
Базовый шаг сокращения делит произведение на две e -битные части, умножает старшую часть на d и складывает их: ( ax mod 2 e ) + d ⌊ ax /2 e ⌋ . За этим можно последовать вычитанием m до тех пор, пока результат не окажется в пределах диапазона. Количество вычитаний ограничено ad / m , которое можно легко ограничить одним, если d мало и выбрано a < m / d . (Это условие также гарантирует, что d ⌊ ax /2 e ⌋ является произведением одинарной ширины; если оно нарушается, необходимо вычислить произведение двойной ширины.)
Когда модуль является простым числом Мерсенна ( d = 1), процедура особенно проста. Не только умножение на d тривиально, но и условное вычитание можно заменить безусловным сдвигом и сложением. Чтобы убедиться в этом, обратите внимание, что алгоритм гарантирует, что x ≢ 0 (mod m ) , а это означает, что x = 0 и x = m невозможны. Это позволяет избежать необходимости рассматривать эквивалентные e -битные представления состояния; только значения, у которых старшие биты не равны нулю, требуют сокращения.
Младшие биты e произведения ax не могут представлять значение, большее m , а старшие биты никогда не будут содержать значение, большее, чем a − 1 ≤ m − 2. Таким образом, первый шаг сокращения дает значение не более m + a − 1. ≤ 2 m − 2 = 2 e +1 − 4. Это ( e + 1)-битное число, которое может быть больше m (т. е. может иметь установленный бит e ), но старшая половина не превышает 1, и если это так, младшие биты e будут строго меньше m . Таким образом, независимо от того, равен ли старший бит 1 или 0, второй шаг сокращения (сложение половин) никогда не приведет к переполнению e битов, и сумма будет желаемым значением.
Если d > 1, условного вычитания также можно избежать, но процедура более сложная. Основная проблема модуля, такого как 2 32 − 5, заключается в обеспечении того, чтобы мы создавали только одно представление для таких значений, как 1 ≡ 2 32 − 4. Решение состоит в том, чтобы временно добавить d , чтобы диапазон возможных значений был от d до 2. e - 1 и уменьшайте значения, превышающие e бит, таким образом, чтобы никогда не генерировались представления меньше d . Наконец, вычитание временного смещения дает желаемое значение.
Начнем с предположения, что у нас есть частично уменьшенное значение y, ограниченное так, что 0 ≤ y < 2 m = 2 e +1 − 2 d . В этом случае один шаг вычитания смещения даст 0 ≤ y ′ = (( y + d ) mod 2 e ) + d ⌊ ( y + d )/2 e ⌋ − d < m . Чтобы убедиться в этом, рассмотрим два случая:
(В частности, в случае генератора Лемера нулевое состояние или его образ y = m никогда не возникнет, поэтому смещение d - 1 будет работать точно так же, если это более удобно. Это уменьшает смещение до 0 в Простой случай Мерсенна, когда d = 1.)
Уменьшение оси большего продукта до значения менее 2 m = 2 e +1 − 2 d может быть выполнено за один или несколько шагов уменьшения без смещения.
Если ad ≤ m , то достаточно одного дополнительного шага сокращения. Поскольку x < m , ax < am ≤ ( a − 1)2 e , и один шаг сокращения преобразует это не более чем в 2 e − 1 + ( a − 1) d = m + ad − 1. Это находится в пределах 2 м , если ad − 1 < m , что является исходным предположением.
Если ad > m , то на первом этапе сокращения возможно получить сумму, большую, чем 2 m = 2 e +1 - 2 d , что слишком велико для последнего шага сокращения. (Как упоминалось выше, также требуется умножение на d , чтобы получить произведение размером больше e бит.) Однако, пока d 2 < 2 e , первое сокращение даст значение в диапазоне, необходимом для предыдущего случая двух шаги по сокращению, которые необходимо применить.
Если продукт двойной ширины недоступен, для вычисления ax mod m можно использовать метод Шраге [14] [15] , также называемый методом приближенного факторинга [16] , но за это приходится платить:
Хотя этот метод популярен для переносимых реализаций в языках высокого уровня , в которых отсутствуют операции двойной ширины, [2] : 744 на современных компьютерах деление на константу обычно реализуется с использованием умножения двойной ширины, поэтому этого метода следует избегать, если важна эффективность. беспокойство. Даже в языках высокого уровня, если множитель 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 ) .
Это равенство имеет место, поскольку
поэтому, если мы факторизуем x = ( x mod q ) + q ⌊ x / q ⌋ , мы получим:
Причина, по которой он не переполняется, заключается в том, что оба члена меньше 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] одним условным добавлением. [17]
Этот метод можно расширить, чтобы разрешить отрицательное значение r (− q ≤ r < 0), изменив окончательное сокращение на условное вычитание.
Этот метод также можно расширить, чтобы обеспечить большее значение a, применяя его рекурсивно. [16] : 102 Из двух членов, вычтенных для получения окончательного результата, только второй ( r ⌊ x / q ⌋ ) рискует переполниться. Но это само по себе является модульным умножением на константу времени компиляции r и может быть реализовано тем же методом. Поскольку каждый шаг в среднем уменьшает вдвое размер множителя (0 ≤ r < a , среднее значение ( a −1)/2), это может потребовать один шаг на бит и быть крайне неэффективным. Однако каждый шаг также делит x на постоянно возрастающее частное q = ⌊ m / a ⌋ , и быстро достигается точка, в которой аргумент равен 0, и рекурсия может быть прекращена.
Используя код C , ГСЧ Парка-Миллера можно записать следующим образом:
uint32_t lcg_parkmiller ( uint32_t * state ) { return * state = ( uint64_t ) * state * 48271 % 0x7fffffff ; }
Эту функцию можно вызывать повторно для генерации псевдослучайных чисел, если вызывающая сторона старается инициализировать состояние любым числом, большим нуля и меньшим модуля. В этой реализации требуется 64-битная арифметика; в противном случае произведение двух 32-битных целых чисел может переполниться.
Чтобы избежать 64-битного деления, выполните сокращение вручную:
uint32_t lcg_parkmiller ( uint32_t * состояние ) { uint64_t product = ( uint64_t ) * состояние * 48271 ; uint32_t x = ( продукт & 0x7fffffff ) + ( продукт >> 31 ); х = ( х и 0x7ffffff ) + ( х >> 31 ); вернуть * состояние = х ; }
Чтобы использовать только 32-битную арифметику, используйте метод Шраге:
uint32_t lcg_parkmiller ( uint32_t * state ) { // Предварительно вычисленные параметры для метода Шраге const uint32_t M = 0x7ffffff ; const uint32_t A = 48271 ; const uint32_t Q = M / A ; // 44488 const uint32_t R = M % A ; // 3399 uint32_t div = * состояние / Q ; // максимум: M/Q = A = 48,271 uint32_t rem = * состояние % Q ; // максимум: Q - 1 = 44 487 int32_t s = rem * A ; // максимум: 44 487 * 48 271 = 2 147 431 977 = 0x7fff3629 int32_t t = div * R ; // максимум: 48 271 * 3 399 = 164 073 129 int32_t result = s - t ; если ( результат < 0 ) результат += M ; возврат * состояние = результат ; }
или используйте два умножения 16×16 бит:
uint32_t lcg_parkmiller ( uint32_t * состояние ) { const uint32_t A = 48271 ; uint32_t low = ( * состояние & 0x7fff ) * A ; // максимум: 32 767 * 48 271 = 1 581 695 857 = 0x5e46c371 uint32_t high = ( * состояние >> 15 ) * A ; // максимум: 65 535 * 48 271 = 3 163 439 985 = 0xbc8e4371 uint32_t x = low + (( high & 0xffff ) << 15 ) + ( high >> 16 ); // максимум: 0x5e46c371 + 0x7fff8000 + 0xbc8e = 0xde46ffff х = ( х и 0x7ffffff ) + ( х >> 31 ); вернуть * состояние = х ; }
Другой популярный генератор Лемера использует простой модуль 2 32 −5:
uint32_t lcg_rand ( uint32_t * state ) { return * state = ( uint64_t ) * state * 279470273u % 0xfffffffb ; }
Это также можно записать без 64-битного деления:
uint32_t lcg_rand ( uint32_t * состояние ) { uint64_t product = ( uint64_t ) * состояние * 279470273u ; uint32_t х ; // Не требуется, поскольку 5 * 279470273 = 0x5349e3c5 умещается в 32 бита. // продукт = (продукт & 0xffffffff) + 5 * (продукт >> 32); // Он понадобится множителю больше 0x33333333 = 858 993 459.// Результат умножения умещается в 32 бита, но сумма может составлять 33 бита. продукт = ( продукт & 0xffffffff ) + 5 * ( uint32_t ) ( продукт >> 32 ); продукт += 4 ; // Эта сумма гарантированно будет 32 бита. x = ( uint32_t ) произведение + 5 * ( uint32_t )( произведение >> 32 ); вернуть * состояние = х - 4 ; }
Многие другие генераторы Лемера обладают хорошими свойствами. Следующий генератор Лемера по модулю 2 128 требует 128-битной поддержки со стороны компилятора и использует множитель, вычисленный Л'Экуайером. [18] Его период равен 2 126 :
статическое беззнаковое состояние __int128 ; /* Состояние должно быть заполнено нечетным значением. */ void семя ( беззнаковое семя __int128 ) { state = семя << 1 | 1 ; } uint64_t next ( void ) { // GCC не может записывать 128-битные литералы, поэтому мы используем выражение const unsigned __int128 mult = ( unsigned __int128 ) 0x12e15e35b500f16e << 64 | 0x2e714eb2b37916a5 ; состояние *= мульт ; состояние возврата >> 64 ; }
Генератор вычисляет нечетное 128-битное значение и возвращает его старшие 64 бита.
Этот генератор проходит BigCrush из TestU01 , но не проходит тест TMFn из PractRand. Этот тест был разработан, чтобы точно выявить дефект генератора этого типа: поскольку модуль представляет собой степень 2, период младшего бита на выходе составляет всего 2 62 , а не 2 126 . Линейные конгруэнтные генераторы с модулем степени 2 ведут себя аналогично.
Следующая основная процедура повышает скорость приведенного выше кода для целочисленных рабочих нагрузок (если компилятор позволяет оптимизировать объявление константы вне цикла вычислений):
uint64_t next ( void ) { uint64_t result = состояние >> 64 ; // GCC не может записывать 128-битные литералы, поэтому мы используем выражение const unsigned __int128 mult = ( unsigned __int128 ) 0x12e15e35b500f16e << 64 | 0x2e714eb2b37916a5 ; состояние *= мульт ; вернуть результат ; }
Однако, поскольку умножение отложено, оно не подходит для хеширования, поскольку первый вызов просто возвращает старшие 64 бита начального состояния.
ZX81 использует p=65537 и a=75 [...](Обратите внимание, что в руководстве ZX81 неверно указано, что 65537 — это простое число Мерсенна, равное 2 16 − 1. В руководстве ZX Spectrum это исправлено и правильно указано, что это простое число Ферма, равное 2 16 + 1.)
ZX Spectrum использует p=65537 и a=75 и хранит в памяти некоторое количество би-1.
Жребий определяется с использованием модульной арифметики, например,
, ... Функция CRAY RANF выбрасывает только три из шести возможных результатов (три стороны которых зависят от начального числа)!
lrand48() % 6 + 1