stringtranslate.com

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

Генератор случайных чисел Лемера [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  Одним из примеров является генератор Вихмана-Хилла .

Отношение к LCG

Хотя ГСЧ Лемера можно рассматривать как частный случай линейного конгруэнтного генератора с 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 ed может реализовать дешевле, чем обычное деление двойной ширины с использованием тождества 2 ed (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 edm . Чтобы убедиться в этом, рассмотрим два случая:

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  -  dy  - 2 e  +  dy  -  m  <  m , по желанию. Поскольку умноженная старшая часть равна d , сумма равна как минимум d , и вычитание смещения никогда не приводит к опустошению.

(В частности, в случае генератора Лемера нулевое состояние или его образ 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 = ( mr )/ a . Затем на каждой итерации вычисляйте axa ( x mod q ) − r x / q (mod m ) .

Это равенство имеет место, поскольку

поэтому, если мы факторизуем x = ( x mod q ) + q x / q , мы получим:

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

Если a выбрано так, что r  ≤  q (и, следовательно, r / q  ≤ 1), то второй член также меньше m : r x / q rx / q = x ( r / q ) ≤ x (1 ) знак равно Икс < м . Таким образом, разница лежит в диапазоне [1− mm −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, и рекурсия может быть прекращена.

Пример кода 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 );                х = ( х и 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 бита начального состояния.

Рекомендации

  1. ^ WH Пейн; Дж. Р. Рабунг; Т.П. Богё (1969). «Кодирование генератора псевдослучайных чисел Лемера» (PDF) . Коммуникации АКМ . 12 (2): 85–86. дои : 10.1145/362848.362860. S2CID  2749316.[1]
  2. ^ abc L'Ecuyer, Пьер (июнь 1988 г.). «Эффективные и портативные комбинированные генераторы случайных чисел» (PDF) . Коммуникации АКМ . 31 (6): 742–774. дои : 10.1145/62959.62969. S2CID  9593394.
  3. ^ Парк, Стивен К.; Миллер, Кейт В. (1988). «Генератор случайных чисел: хорошие найти сложно» (PDF) . Коммуникации АКМ . 31 (10): 1192–1201. дои : 10.1145/63039.63042. S2CID  207575300.
  4. ^ Марсалья, Джордж (1993). «Техническая переписка: Замечания по выбору и реализации генераторов случайных чисел» (PDF) . Коммуникации АКМ . 36 (7): 105–108. дои : 10.1145/159544.376068. S2CID  26156905.
  5. ^ Салливан, Стивен (1993). «Техническое соответствие: еще один тест на случайность» (PDF) . Коммуникации АКМ . 36 (7): 108. дои : 10.1145/159544.376068. S2CID  26156905.
  6. ^ Парк, Стивен К.; Миллер, Кейт В.; Стокмейер, Пол К. (1988). «Техническая переписка: Ответ» (PDF) . Коммуникации АКМ . 36 (7): 108–110. дои : 10.1145/159544.376068. S2CID  26156905.
  7. ^ Викерс, Стив (1981). «Глава 5. Функции». Базовое программирование ZX81 (2-е изд.). Синклер Рисерч Лтд . Проверено 21 апреля 2024 г. ZX81 использует p=65537 и a=75 [...](Обратите внимание, что в руководстве ZX81 неверно указано, что 65537 — это простое число Мерсенна, равное 2 16  − 1. В руководстве ZX Spectrum это исправлено и правильно указано, что это простое число Ферма, равное 2 16  + 1.)
  8. ^ Викерс, Стив (1983). «Глава 11. Случайные числа». Базовое программирование Sinclair ZX Spectrum (2-е изд.). Sinclair Research Ltd., стр. 73–75 . Проверено 26 мая 2022 г. ZX Spectrum использует p=65537 и a=75 и хранит в памяти некоторое количество би-1.
  9. ^ ab Научная библиотека GNU: Другие генераторы случайных чисел.
  10. ^ Кнут, Дональд (1981). Получисловые алгоритмы . Искусство компьютерного программирования . Том. 2 (2-е изд.). Ридинг, Массачусетс: Addison-Wesley Professional. стр. 12–14.
  11. ^ Бике, Стивен; Розенберг, Роберт (май 2009 г.). Быстрая генерация высококачественных псевдослучайных чисел и перестановок с использованием MPI и OpenMP на Cray XD1. Группа пользователей Cray, 2009. Жребий определяется с использованием модульной арифметики, например, , ... Функция CRAY RANF выбрасывает только три из шести возможных результатов (три стороны которых зависят от начального числа)!lrand48() % 6 + 1
  12. ^ Гринбергер, Мартин (апрель 1961 г.). «Заметки о новом генераторе псевдослучайных чисел». Журнал АКМ . 8 (2): 163–167. дои : 10.1145/321062.321065 . S2CID  17196519.
  13. ^ Л'Экуйер, Пьер; Тэдзука, Шу (октябрь 1991 г.). «Структурные свойства двух классов комбинированных генераторов случайных чисел». Математика вычислений . 57 (196): 735–746. дои : 10.2307/2938714. JSTOR  2938714.
  14. ^ Шраге, Линус (июнь 1979 г.). «Более портативный генератор случайных чисел на Фортране» (PDF) . Транзакции ACM в математическом программном обеспечении . 5 (2): 132–138. CiteSeerX 10.1.1.470.6958 . дои : 10.1145/355826.355828. S2CID  14090729. 
  15. ^ Джайн, Радж (9 июля 2010 г.). «Анализ производительности компьютерных систем, глава 26: Генерация случайных чисел» (PDF) . стр. 19–22 . Проверено 31 октября 2017 г.
  16. ^ аб Л'Экуйер, Пьер; Коте, Серж (март 1991 г.). «Реализация пакета случайных чисел с возможностью разделения». Транзакции ACM в математическом программном обеспечении . 17 (1): 98–111. дои : 10.1145/103147.103158 . S2CID  14385204. Здесь исследуются несколько различных реализаций модульного умножения на константу.
  17. Фенерти, Пол (11 сентября 2006 г.). «Метод Шраге» . Проверено 31 октября 2017 г.
  18. ^ Л'Экуайер, Пьер (январь 1999 г.). «Таблицы линейных конгруэнтных генераторов разных размеров и хорошей структуры решетки» (PDF) . Математика вычислений . 68 (225): 249–260. CiteSeerX 10.1.1.34.1024 . дои : 10.1090/s0025-5718-99-00996-5. 

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