Генератор случайных чисел Лехмера [1] (названный после того, как DH Лехмер ), иногда также упоминаются как генератор Парк-Miller случайных чисел (после того, как Стивен К. Парк и Кит У. Миллер), представляет собой тип линейного конгруэнтного генератора (LCG) который действует в мультипликативной группе целых чисел по модулю n . Общая формула:
где модуль т является простым числом или степенью простого числа , множитель является элементом высокого мультипликативного порядка по модулю т (например, первообразный корень ), и семя Х 0 является взаимно прост с т .
Другие названия - мультипликативный линейный конгруэнтный генератор (MLCG) [2] и мультипликативный конгруэнтный генератор (MCG) .
Параметры в общем использовании
В 1988 году Парк и Миллер [3] предложили ГСЧ Лемера с конкретными параметрами m = 2 31 - 1 = 2 147 483 647 ( простое число Мерсенна M 31 ) и a = 7 5 = 16 807 (примитивный корень по модулю M 31 ), теперь известный как РАЗУМ . Хотя MINSTD был позже критикой Marsaglia и Салливан (1993), [4] [5] он все еще используется сегодня (в частности, в CarbonLib и C ++ 11 «ы minstd_rand0
). Парк, Миллер и Стокмейер ответили на критику (1993) [6], сказав:
Учитывая динамичный характер местности, неспециалистам сложно принять решение о том, какой генератор использовать. «Дайте мне то, что я могу понять, реализовать и перенести ... это не обязательно должно быть по последнему слову техники, просто убедитесь, что оно достаточно хорошее и эффективное». Наша статья и связанный с ней генератор минимальных стандартов были попыткой ответить на этот запрос. Пять лет спустя мы не видим необходимости изменять наш ответ, кроме как предложить использовать множитель a = 48271 вместо 16807.
Эта исправленная константа используется в C ++ 11 «с minstd_rand
генератором случайных чисел.
Sinclair ZX81 и его преемники использовать Лехмер RNG с параметрами т = 2 16 +1 = 65,537 (а Ферма простое F 4 ) и а = 75 (первообразный корень по модулю F 4 ). [7] CRAY генератора случайных чисел RANF является Лехмер ГСЧ с мощностью из-два модуля м = 2 48 и = 44,485,709,377,909. [8] ГНУ научная библиотека включает в себя несколько генераторов случайных чисел вида Лемера, в том числе MINSTD, RANF и печально известной IBM генератора случайных чисел Randu . [8]
Выбор модуля
Чаще всего модуль выбирается как простое число, что делает выбор взаимно простого начального числа тривиальным ( подойдет любое 0 < X 0 < m ). Это обеспечивает наилучшее качество вывода, но вносит некоторую сложность в реализацию, и диапазон вывода вряд ли будет соответствовать желаемому приложению; преобразование в желаемый диапазон требует дополнительного умножения.
Используя модуль м , который представляет собой мощность двух марок для особенно удобной компьютерной реализации, но приходит по цене: период составляет не более м / 4, и низкие биты имеют периоды короче , чем это. Это потому , что низкие K битов образуют по модулю 2 K генератора все сами по себе; биты более высокого порядка никогда не влияют на биты более низкого порядка. [9] Значения X i всегда нечетные (бит 0 никогда не меняется), биты 2 и 1 чередуются (младшие 3 бита повторяются с периодом 2), младшие 4 бита повторяются с периодом 4 и так далее. Следовательно, приложение, использующее эти случайные числа, должно использовать самые старшие биты; уменьшение до меньшего диапазона с использованием операции по модулю с четным модулем приведет к катастрофическим результатам. [10]
Для достижения этого периода, множитель должен удовлетворять а ≡ ± 3 ( по модулю 8) [11] и семя X 0 должно быть нечетным.
Использование составного модуля возможно, но генератор должен быть заполнен значением, взаимно простым с m , иначе период будет значительно сокращен. Например, модуль F 5 = 2 32 +1 может показаться привлекательным, поскольку выходные данные могут быть легко отображены в 32-битное слово 0 ≤ X i −1 <2 32 . Однако начальное значение X 0 = 6700417 (которое делит 2 32 +1) или любое кратное значение приведет к выходу с периодом только 640.
Более популярной реализацией для больших периодов является комбинированный линейный конгруэнтный генератор ; объединение (например, суммирования их выходов) нескольких генераторов эквивалентно выходу одного генератора, модуль которого является произведением модулей генераторов компонентов. [12] и период которой является наименьшим общим кратным составляющих периодов. Хотя периоды будут иметь общий делитель 2, модули можно выбрать так, чтобы это был единственный общий делитель, а результирующий период был равен ( m 1 −1) ( m 2 −1) ( m 2 ··· ( m k - 1) / 2 k −1 . [2] : 744 Одним из примеров этого является генератор Вихмана – Хилла .
Отношение к LCG
Хотя ГСЧ Лемера можно рассматривать как частный случай линейного конгруэнтного генератора с c = 0 , это частный случай, который подразумевает определенные ограничения и свойства. В частности, для ГСЧ Лемера начальное начальное число X 0 должно быть взаимно простым с модулем m, что не требуется для ГСЧ в целом. Выбор модуля m и множителя a также более ограничен для ГСЧ Lehmer. В отличие от LCG, максимальный период ГСЧ Лемера равен m -1, и это так, когда m - простое число, а a - первообразный корень по модулю m .
С другой стороны, дискретные логарифмы (для основания a или любого примитивного корня по модулю m ) X k впредставляют собой линейную конгруэнтную последовательность по модулю Эйлера .
Выполнение
Простой модуль требует вычисления произведения двойной ширины и явного шага редукции. Если используется модуль чуть меньше степени 2 ( популярны простые числа Мерсенна 2 31 −1 и 2 61 −1, а также 2 32 −5 и 2 64 −59), уменьшение по модулю m = 2 e - d может быть реализовано дешевле, чем обычное деление двойной ширины с использованием тождества 2 e ≡ d (mod m ) .
Основные стадия восстановления делит продукт на два электронных -битный части, умножают высокую часть по д , и добавляет их: ( ах по модулю 2 е ) + д ⌊ ах / 2 х ⌋ . После этого можно вычесть m, пока результат не окажется в допустимом диапазоне. Количество вычитаний ограничено до ad / m , которое можно легко ограничить до одного, если d мало и выбрано a < m / d . (Это условие также гарантирует, что d ⌊ ax / 2 e ⌋ является произведением одинарной ширины; если оно нарушается, должно быть вычислено произведение двойной ширины.)
Когда модуль представляет собой простое число Мерсенна ( d = 1), процедура особенно проста. Мало того, что умножение на d тривиально, условное вычитание можно заменить безусловным сдвигом и сложением. Чтобы увидеть это, обратите внимание, что алгоритм гарантирует, что x ≢ 0 (mod m) , а это означает, что x = 0 и x = m оба невозможны. Это позволяет избежать необходимости рассматривать эквивалентные e- битовые представления состояния; только значения, у которых старшие биты не равны нулю, нуждаются в уменьшении.
Низкие электронные биты продукта ах не могут представлять значение больше , чем м , а старшие биты никогда не будут держать значение больше -1 ≤ м-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 32 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 . Чтобы убедиться в этом, рассмотрим два случая:
- 0 ≤ y < m = 2 e - d
- В этом случае y + d <2 e и y ′ = y < m , что и нужно.
- м ≤ у <2 м
- В этом случае 2 e ≤ y + d <2 e +1 является e + 1-битным числом, а ⌊ ( y + d ) / 2 e ⌋ = 1. Таким образом, y ′ = ( y + d ) - 2 e + d - d = y - 2 e + d = y - m < m , как и нужно. Поскольку умноженная старшая часть равна d , сумма равна по крайней мере d, и вычитание смещения никогда не вызывает потери значимости.
(В частности, в случае генератора Лемера нулевое состояние или его изображение y = m никогда не произойдет, поэтому смещение d −1 будет работать точно так же, если это более удобно. Это уменьшает смещение до 0 в Простой случай Мерсенна, когда d = 1.)
Уменьшение более крупного продукта ax до менее 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 m, если ad - 1 < m , что является исходным предположением.
Если ad > m , то на первом этапе сокращения возможно получение суммы, превышающей 2 m = 2 e +1 -2 d , что слишком велико для заключительного этапа сокращения. (Это также требует умножения на d для получения произведения, большего, чем e битов, как упомянуто выше.) Однако, пока d 2 <2 e , первое уменьшение даст значение в диапазоне, требуемом для предыдущего случая двух шаги по сокращению для применения.
Метод Шраге
Если двойной ширины продукт не доступен, метод Шраге в , [13] [14] также называемый приближенный метод factoriation, [15] , могут быть использованы для вычисления ах по модулю м , но это происходит за счет:
- Модуль должен быть представлен в виде целого числа со знаком ; арифметические операции должны допускать диапазон ± m .
- Выбор множителя a ограничен. Мы требуем, чтобы m mod a ≤ ⌊ m / a ⌋ , обычно достигается выбором a ≤ √ m .
- Требуется одно деление (с остатком) на итерацию.
Хотя этот метод популярен для переносимых реализаций на языках высокого уровня, в которых отсутствуют операции двойной ширины, [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] с помощью одного условного сложения. [16]
Этот метод может быть расширен, чтобы разрешить отрицательное r (- q ≤ r <0), изменив окончательное сокращение на условное вычитание.
Техника также может быть расширена для увеличения a путем ее рекурсивного применения. [15] : 102 Из двух членов, вычтенных для получения окончательного результата, только второй ( 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 * состояние ) { uint64_t product = ( uint64_t ) * состояние * 48271 ; uint32_t x = ( продукт & 0x7fffffff ) + ( продукт >> 31 );х = ( х & 0x7fffffff ) + ( х >> 31 ); возврат * состояние = х ; }
Чтобы использовать только 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 ; // 3399uint32_t div = * состояние / Q ; // макс: M / Q = A = 48 271 uint32_t rem = * состояние % Q ; // макс: Q - 1 = 44 487int32_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х = ( х & 0x7fffffff ) + ( х >> 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 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 обладают хорошими свойствами. Следующий генератор Лемера по модулю 2 128 требует 128-битной поддержки от компилятора и использует множитель, вычисленный L'Ecuyer. [17] Он имеет период 2 126 :
статическое беззнаковое состояние __int128 ;/ * Состояние должно быть заполнено нечетным значением. * / void seed ( unsigned __int128 seed ) { state = seed << 1 | 1 ; }uint64_t следующий ( аннулируется ) { // GCC не может написать 128-битные литералы, поэтому мы используем выражение сопзИте неподписанный __int128 мульт = ( беззнаковое __int128 ) 0x12e15e35b500f16e << 64 | 0x2e714eb2b37916a5 ; состояние * = мульт ; вернуть состояние >> 64 ; }
Генератор вычисляет нечетное 128-битное значение и возвращает его старшие 64 бита.
Этот генератор проходит BigCrush от TestU01 , но не проходит тест TMFn от PractRand . Этот тест был разработан, чтобы точно выявить дефект генератора этого типа: поскольку модуль имеет степень двойки, период самого младшего бита на выходе составляет всего 2 62 , а не 2 126 . Аналогичное поведение имеют линейные конгруэнтные генераторы с модулем степени двойки.
Следующая основная подпрограмма улучшает скорость вышеуказанного кода для целочисленных рабочих нагрузок (если компилятор разрешает оптимизировать объявление константы вне цикла вычислений):
uint64_t next ( void ) { uint64_t result = state >> 64 ; // GCC не может написать 128-битные литералы, поэтому мы используем выражение сопзИте неподписанный __int128 мульт = ( беззнаковое __int128 ) 0x12e15e35b500f16e << 64 | 0x2e714eb2b37916a5 ; состояние * = мульт ; вернуть результат ; }
However, because the multiplication is deferred, it is not suitable for hashing, since the first call simply returns the upper 64 bits of the seed state.
Рекомендации
- ^ W.H. Payne; J.R. Rabung; T.P. Bogyo (1969). "Coding the Lehmer pseudo-random number generator" (PDF). Communications of the ACM. 12 (2): 85–86. doi:10.1145/362848.362860.
- ^ a b c L'Ecuyer, Pierre (June 1988). "Efficient and Portable Combined Random Number Generators" (PDF). [Communications of the ACM]. 31 (6): 742–774. doi:10.1145/62959.62969.
- ^ Park, Stephen K.; Miller, Keith W. (1988). "Random Number Generators: Good Ones Are Hard To Find" (PDF). Communications of the ACM. 31 (10): 1192–1201. doi:10.1145/63039.63042.
- ^ Marsaglia, George (1993). "Technical correspondence: Remarks on Choosing and Implementing Random Number Generators" (PDF). Communications of the ACM. 36 (7): 105–108. doi:10.1145/159544.376068.
- ^ Sullivan, Stephen (1993). "Technical correspondence: Another test for randomness" (PDF). Communications of the ACM. 36 (7): 108. doi:10.1145/159544.376068.
- ^ Park, Stephen K.; Miller, Keith W.; Stockmeyer, Paul K. (1988). "Technical Correspondence: Response" (PDF). Communications of the ACM. 36 (7): 108–110. doi:10.1145/159544.376068.
- ^ Vickers, Steve (1981). "Chapter 5 - Functions". ZX81 Basic Programming (2nd ed.). Sinclair Research Ltd.
The ZX81 uses p=65537 & a=75 [...]
(Note the ZX81 manual incorrectly states that 65537 is a Mersenne prime that equals 216−1. The ZX Spectrum manual fixed that and correctly states that it is a Fermat prime that equals 216+1.) - ^ a b GNU Scientific Library: Other random number generators
- ^ Knuth, Donald (1981). Seminumerical Algorithms. The Art of Computer Programming. 2 (2nd ed.). Reading, MA: Addison-Wesley Professional. pp. 12–14.
- ^ Bique, Stephen; Rosenberg, Robert (May 2009). Fast Generation of High-Quality Pseudorandom Numbers and Permutations Using MPI and OpenMP on the Cray XD1. Cray User Group 2009.
The die is determined using modular arithmetic, e.g.,
lrand48() % 6 + 1,
... The CRAY RANF function only rolls three of the six possible outcomes (which three sides depends on the seed)! - ^ Greenberger, Martin (April 1961). "Notes on a New Pseudo-Random Number Generator". Journal of the ACM. 8 (2): 163–167. doi:10.1145/321062.321065.
- ^ L'Ecuyer, Pierre; Tezuka, Shu (October 1991). "Structural Properties for Two Classes of Combined Random Number Generators". Mathematics of Computation. 57 (196): 735–746. doi:10.2307/2938714.
- ^ Schrage, Linus (June 1979). "A More Portable Fortran Random Number Generator" (PDF). ACM Transactions on Mathematical Software (TOMS). 5 (2): 132–138. CiteSeerX 10.1.1.470.6958. doi:10.1145/355826.355828.
- ^ Jain, Raj (9 July 2010). "Computer Systems Performance Analysis Chapter 26: Random-Number Generation" (PDF). pp. 19–22. Retrieved 2017-10-31.
- ^ a b L'Ecuyer, Pierre; Côté, Serge (March 1991). "Implementing a random number package with splitting facilities". ACM Transactions on Mathematical Software. 17 (1): 98–111. doi:10.1145/103147.103158. This explores several different implementations of modular multiplication by a constant.
- ^ Fenerty, Paul (11 September 2006). "Schrage's Method". Retrieved 2017-10-31.
- ^ L’Ecuyer, Pierre (January 1999). "Tables of linear congruential generators of different sizes and good lattice structure" (PDF). Mathematics of Computation. 68 (225): 249–260. CiteSeerX 10.1.1.34.1024. doi:10.1090/s0025-5718-99-00996-5.
- Lehmer, D. H. (1949). "Mathematical methods in large-scale computing units". Proceedings of a Second Symposium on Large-Scale Digital Calculating Machinery. pp. 141–146. MR 0044899. (journal version: Annals of the Computation Laboratory of Harvard University, Vol. 26 (1951)).
- Steve Park, Random Number Generators
Внешние ссылки
- Primes just less than a power of two may be useful for choosing moduli. Part of Prime Pages.