레머 난수 발생기

Lehmer random number generator

Lehmer[1] 난수 발생기(D. H. Lehmer의 이름)는 때로는 Park-Miller 난수 발생기(Stephen K의 이름)라고도 한다.Park와 Keith W. Miller)는 정수 modulo n의 곱셈 그룹에서 작동하는 선형융합발전기(LCG)의 일종이다.일반적인 공식은 다음과 같다.

여기서 계량 m은 소수 또는 소수 검정력이고, 승수 a는 높은 승수 순서 modulo m(예: 원시 뿌리 modulo n)의 요소이며, 시드 X는0 m과 동일하다.

다른 명칭은 MLCG([2]승수 선형 결합 발전기)와 MCG(승수 결합 발전기)가 있다.

공통 사용 중인 매개 변수

1988년, 박과 밀러는[3] 특정한31 매개변수 m = 2 - 1 = 2,147,483,647 (메르센 프라임M31)과 a = 7 = 165,807 (원시 루트 모듈로 M31)을 가진 Lehmer RNG를 제안했는데, 현재 MINSTD로 알려져 있다.MINSTD는 이후 마르사글리아와 설리반(1993)으로부터 비판을 받았지만,[4][5] 오늘날(특히, CarbonLib과 C++11의)에도 여전히 사용되고 있다.minstd_rand0)박, 밀러, 스톡마이어는 비판(1993)에 대해 다음과 같이 대답했다.[6]

그 지역의 역동적인 특성을 감안할 때, 비전문가가 어떤 발전기를 사용할지에 대한 결정을 내리기는 어렵다."내가 이해하고 이행하고 항만할 수 있는 것을 주옵소서...첨단일 필요는 없고, 합리적으로 좋고 효율적인지 확인만 하면 돼."우리의 기사 및 관련 최소 표준 발전기는 이 요청에 대한 대응 시도였습니다.5년 후, 우리는 16807 대신에 승수 a = 48271의 사용을 제안하는 것 외에는 우리의 반응을 바꿀 필요가 없다고 본다.

이 수정된 상수는 C++11에서 사용된다.minstd_rand난수 생성기

싱클레어 ZX81과 그 계승자는 매개변수 m = 216+1 = 65,537(페르마트 프라임4 F) 및 a = 75(원시 루트 모듈로4 F)를 가진 Lehmer RNG를 사용한다.[7]CRAY 난수 발생기 RANF는 2의 계량 m = 2와48 a = 44,485,709,377,909의 전력을 가진 Lehmer RNG이다.[8]GNU 과학 라이브러리에는 MINSTD, RANF, 그리고 악명 높은 IBM 난수 생성기 RANDU를 포함하여 Lehmer 형식의 난수 생성기 여러 개가 포함되어 있다.[8]

계량선택

가장 일반적으로는 계수를 프라임 숫자로 선택하는데, 복음 씨앗의 선택은 사소한 것으로 만든다(0< X < m이면0 된다).이것은 최상의 품질의 출력을 생산하지만, 일부 구현 복잡성을 야기하고 출력 범위가 원하는 애플리케이션과 일치하지 않을 가능성이 높다. 원하는 범위로 변환하려면 추가적인 곱셈이 필요하다.

두 개의 힘인 계량 m을 사용하면 특히 편리한 컴퓨터 구현을 할 수 있지만, 비용이 많이 들지만, 기간은 최대 m/4이고, 낮은 비트는 그보다 짧은 기간을 갖는다.이것은 낮은 k 비트가 스스로 모듈로-2k 발전기를 형성하기 때문이다; 높은 순서의 비트는 결코 낮은 순서의 비트에 영향을 미치지 않는다.[9]값i X는 항상 홀수(비트 0은 절대 변하지 않음), 비트 2와 1을 교대로(낮은 3비트는 주기 2로 반복), 낮은 4비트는 주기 4로 반복하는 등이다.따라서 이러한 무작위 숫자를 사용하는 애플리케이션은 가장 중요한 비트를 사용해야 한다. 짝수 계수를 사용하는 모듈로 연산을 사용하여 더 작은 범위로 축소하면 비참한 결과를 초래할 것이다.[10]

이 기간을 달성하려면 승수는 ≡ ±3(모드 8)[11]을 만족해야 하며 시드 X는0 홀수여야 한다.

복합계수를 사용할 수 있지만, 발전기는 m에 대한 값 조합으로 시드해야 하며, 그렇지 않으면 기간이 크게 단축된다.예를 들어, 출력은 32비트 단어 0 x X-1i < 2에32 쉽게 매핑될 수 있기 때문에5 F = 232+1의 계수는 매력적으로 보일 수 있다.그러나 X0 = 6700417(232+1) 또는 임의의 배수인 경우 주기가 640에 불과한 출력이 발생한다.

대규모 기간 동안 더 인기 있는 구현은 조합된 선형 결합 발전기이다. (예를 들어 출력을 합산하여) 여러 발전기가 조합된 것은 계량기가 구성품 발전기 모듈리의 산물이고 구성품 기간 중 최소 공통 배수인 단일 발전기의 출력과 동일하다.[12]기간은 공통점 2를 공유하지만, 모듈리는 그것만이 공통점이고 그 결과 기간은 (m-11)(m-12)(m-12)(m-1k)/2k−1)이다.[2]: 744 이것의 한 예는 위크만-힐 발전기다.

LCG와의 관계

Lehmer RNG는 c=0으로 선형 결합 발전기의 특정 사례로 볼 수 있지만, 일정한 제한과 속성을 내포하는 특수한 경우다.특히 Lehmer RNG의 경우, 초기 시드 X는0 일반적으로 LCG에 필요하지 않은 계량 m에 복사되어야 한다.계량 m과 승수 a의 선택은 또한 레머 RNG에 더 제한적이다. LCG와 대조적으로, 레머 RNG의 최대 주기는 m-1이며, m이 프라임이고 a가 원시 루트 모듈로 m이다.

반면, 에서 X의k 이산형 로그(또는 원시 루트 모듈로 m)는 오일러 에 선형 유사 시퀀스를 나타낸다.

실행

프라임 계수는 이중 폭의 제품과 명시적인 감소 단계를 계산해야 한다.전력 2보다 조금 적은 계수를 사용할 경우(메르센 소수 2-1과31 2-1은64 2-596132, 2-59) 감량모듈로 me = 2 - d는 아이덴티티 2e ≡ d(mod m)를 사용하여 일반 이중폭 분할보다 저렴하게 구현할 수 있다.

기본 감속 단계는 제품을 두 개의 e비트 부분으로 나누고 높은 부분을 d로 곱한 후 다음과 같은 (ax mod 2e) + d ⌊ax/2e⌋을 더한다.결과가 범위 내에 있을 때까지 m을 빼면 다음과 같은 결과를 얻을 수 있다.소급 횟수는 ad/m으로 제한되는데, d가 작아서 < m/d를 선택하면 쉽게 1로 제한될 수 있다.(이 조건은 d axax/2e⌋가 단폭 제품임을 보장하기도 하며, 이를 어길 경우 반드시 이중폭 제품을 계산해야 한다.)

계수가 메르센 전성기(d = 1)일 때는 시술이 특히 간단하다.d에 의한 곱셈은 사소한 것일 뿐만 아니라 조건부 뺄셈은 무조건적인 이동과 덧셈으로 대체할 수 있다.이를 확인하기 위해 알고리즘은 x 0 0(mod m)을 보장하며, x = 0과 x = m 모두 불가능하다는 것을 의미한다.따라서 하이비트가 0이 아닌 경우에만 해당 상태의 동등한 전자비트 표현을 고려할 필요가 없다.

제품 도끼의 낮은 e비트는 m보다 큰 값을 나타낼 수 없으며, 높은 비트는 결코 a-1 m m-2보다 큰 값을 보유하지 않는다.따라서 첫 번째 감소 단계는 최대 m+a-1 ≤ 2m-2 = 2-4의e+1 값을 생성한다.이 숫자는 m보다 클 수 있는 e+1비트 숫자지만(즉, 비트 e가 설정되어 있을 수 있음), 하이 반은 최대 1이며, 로우 e 비트는 m보다 확실히 작다.따라서 높은 비트가 1이든 0이든 두 번째 감소 단계(반쪽 추가)는 결코 e비트가 넘치지 않고 합은 원하는 값이 될 것이다.

d > 1이면 조건부 뺄셈도 피할 수 있지만 절차가 더 복잡하다.2-5와32 같은 계수의 근본적인 어려움은 우리가 1 232 2-4와 같은 값에 대해 하나의 표현만을 생산하도록 하는 것이다.해결책은 d에서 2-1까지e 가능한 값의 범위가 되도록 일시적으로 d를 추가하고 d보다 작은 표현을 결코 생성하지 않는 방식으로 e비트보다 큰 값을 줄이는 것이다.마지막으로 임시 오프셋을 빼면 원하는 값이 생성된다.

먼저 0 ≤ y < 2m = 2-2de+1 정도로 경계된 y의 부분적인 감소된 값이 있다고 가정한다.이 경우 단일 오프셋 감산 단계는 0 ≤ y y = (y+d) mode 2) + d ⌊(y+d)/22e - d < m을 생성한다.이를 확인하려면 다음 두 가지 경우를 고려하십시오.

0 ≤ y < m = 2e − d
이 경우 y + d < 2와e y′ = y < m, 원하는 대로.
m ≤ y < 2m
이 경우 2e ≤ y + d < 2는e+1 e+1비트 번호로 ⌊(y+d)/22e = 1이다.따라서 y′ = (y+d) - 2e +d - d = y - 2e + d = y - m < m.곱한 높은 부분은 d이기 때문에 합은 최소한 d이고 오프셋을 빼면 결코 과소 유동을 유발하지 않는다.

(특히 Lehmer 발전기의 경우, 영점 상태나 영점 y = m은 절대 발생하지 않으므로, 그것이 더 편리하다면, d-1의 오프셋은 똑같이 작동할 것이다.이렇게 하면 메르센 프라임 케이스에서 오프셋이 0으로 감소하는데, 이때 d = 1.)

더 큰 제품 도끼를 2m = 2e+1 - 2d 미만으로 줄이는 것은 오프셋 없이 하나 이상의 축소 단계를 통해 수행할 수 있다.

ad ad m일 경우, 하나의 추가 감소 단계로 충분하다.x < m, 도끼 < am am a (a-1)2e, 그리고 하나의 감소 단계는 이것을 최대 2e - 1 + (a-1)d = m + ad - 1로 변환한다.이는 초기 가정인 광고 - 1 < m일 경우 2m 한계 이내다.

add > m인 경우, 1차 감량 단계에서 2m = 2-2de+1 이상의 합을 낼 수 있어 최종 감량 단계에는 너무 크다.(또한 위에 언급한 바와 같이 e비트보다 큰 제품을 생산하기 위해서는 d에 의한 곱셈이 필요하다.)단, d2 < 2가e 되는 한, 첫 번째 감소는 앞의 두 감소단계의 적용에 필요한 범위의 값을 산출할 것이다.

슈라지의 방법

이중 폭의 제품을 사용할 수 없는 경우, 대략적인 인수법이라고도 불리는 [13][14]Schrage의 방법을 도끼 mod m을 계산하는 데 사용할 수 있지만,[15] 이것은 다음과 같은 원가로 이루어진다.

  • 계수는 부호화된 정수로 표현 가능해야 하며, 산술 연산은 ±m 범위를 허용해야 한다.
  • 승수 a의 선택은 제한된다.우리는 m m mod a ≤m/a⌋를 요구하는데, 일반적으로 ≤m을 선택함으로써 달성된다.
  • 1회 반복(잔차 포함)이 필요하다.

이 기법은 이중 폭 연산이 부족한 고급 언어의 휴대용 구현에 인기가 있지만,[2]: 744 상수에 의한 현대의 컴퓨터 분할에서는 대개 이중 폭 곱셈을 사용하여 구현되므로 효율성이 우려되는 경우에는 이 기법을 피해야 한다.고수준 언어에서도 승수 a를 √m으로 제한하면 2개의 단폭 곱셈을 사용하여 이중폭 제품 도끼를 계산할 수 있고, 위에서 설명한 기법을 사용하여 줄일 수 있다.

Schrage의 방법을 사용하려면 첫 번째 인자 m = qa + r, 즉, 보조 상수 r = m mod a, q = ⌊m/a⌋ = (m-r)/a를 사전 계산한다.그런 다음, 각 반복은 a(x mod q) - r ⌊x/q⌋ (mod m)를 계산한다.

이 평등은 때문이다.

x = (x mod q) + q qx/q⌋을 인자화하면 다음과 같은 결과를 얻을 수 있다.

그것이 넘치지 않는 이유는 두 용어 모두 m보다 작기 때문이다.x mod q < q ≤ m/a이므로, 첫 번째 항은 am/a = m보다 엄격히 작으며, 단폭 제품으로 계산할 수 있다.

r ≤ q(따라서 r/q ≤ 1)가 되도록 a를 선택한 경우, 두 번째 항도 m: r ⌊x/q⌋ ≤ rx/q = x(r/q) x(1) = x < m보다 작다.따라서 그 차이는 [1-m, m-1] 범위에 있으며, 단일 조건부 추가로 [0, m-1]로 줄일 수 있다.[16]

이 기법은 음의 r(-q ≤ r < 0)을 허용하도록 확장될 수 있으며, 최종 감소를 조건부 빼기로 변경할 수 있다.

또한 이 기법은 반복적으로 적용하여 더 큰 a를 허용하도록 확장할 수 있다.[15]: 102 최종 결과를 도출하기 위해 감산한 두 용어 중 두 번째(r ⌊x/q⌋) 위험만 넘친다.그러나 이것은 그 자체로 컴파일 시간 상수 r에 의한 모듈형 곱셈이며, 동일한 기법에 의해 구현될 수도 있다.각 단계는 평균적으로 승수의 절반 크기(0 r r < a, 평균 값 (a-1)/2)로 나타나기 때문에, 이것은 비트당 한 단계가 필요한 것으로 보이며, 대단히 비효율적이다.그러나 각 단계 역시 x를 계속 상승하는 지수 q = ⌊m/a⌋로 나누며, 빠르게 인수가 0인 지점에 도달하여 재귀가 종료될 수 있다.

샘플 C99 코드

C코드를 이용하여 Park-Miller RNG를 다음과 같이 작성할 수 있다.

uint32_t lcg_parkmiller(uint32_t *주) {  돌아오다 *주 = (uint64_t)*주 * 48271 % 0x7fffffff; }

이 함수는 발신자가 0보다 크고 계수보다 작은 숫자로 상태를 초기화하도록 주의하기만 하면 가성수를 생성하기 위해 반복적으로 호출할 수 있다.이 구현에서는 64비트 산술화가 필요하다. 그렇지 않으면 32비트 정수 두 개의 산물이 넘칠 수 있다.

64비트 분할을 방지하려면 손으로 다음을 줄이십시오.

uint32_t lcg_parkmiller(uint32_t *주) {  uint64_t 상품 = (uint64_t)*주 * 48271;  uint32_t x = (상품 & 0x7fffffff) + (상품 >> 31);   x = (x & 0x7fffffff) + (x >> 31);  돌아오다 *주 = x; }

32비트 산술만 사용하려면 Schrage의 방법을 사용하십시오.

uint32_t lcg_parkmiller(uint32_t *주) {  // Schrage의 방법에 대한 사전 계산된 파라미터  경시하다 uint32_t M = 0x7fffffff;  경시하다 uint32_t A = 48271;  경시하다 uint32_t Q = M / A;    // 44488  경시하다 uint32_t R = M % A;    //  3399   uint32_t 칸막이하다 = *주 / Q; // 최대: M / Q = A = 48,271  uint32_t 렘 = *주 % Q; // 최대: Q - 1 = 44,487   int32_t s = 렘 * A; // 최대: 44,487 * 48,271 = 2,190,431,977 = 0x7fff3629  int32_t t = 칸막이하다 * R; // 최대: 48,271 * 3,399 = 164,073,1930  int32_t 결과 = s - t;   만일 (결과 < 0)   결과 += M;   돌아오다 *주 = 결과; }

또는 두 개의 16×16비트 곱하기 사용:

uint32_t lcg_parkmiller(uint32_t *주) {  경시하다 uint32_t A = 48271;   uint32_t 낮은  = (*주 & 0x7fff) * A;   // 최대: 32,767 * 48,271 = 1,581,695,857 = 0x5e46c371  uint32_t 높은 = (*주 >> 15)    * A;   // 최대: 65,535 * 48,271 = 3,439,985 = 0xbc8e4371  uint32_t x = 낮은 + ((높은 & 0xffff) << 15) + (높은 >> 16); // 최대: 0x5e46c371 + 0x7fff8000 + 0xbc8e = 0xde46ffff   x = (x & 0x7fffffff) + (x >> 31);  돌아오다 *주 = x; }

또 다른 인기 있는 레머 발전기는 프라임 계량 2-5를32 사용한다.

uint32_t lcg_rand(uint32_t *주) {     돌아오다 *주 = (uint64_t)*주 * 279470273u % 0xfffffffb; }

이것은 또한 64비트 분할 없이 작성될 수 있다.

uint32_t lcg_rand(uint32_t *주) {  uint64_t 상품 = (uint64_t)*주 * 279470273u;  uint32_t x;   // 5 * 279470273 = 0x5349e3c5가 32비트에 적합하므로 필요하지 않음.  // 제품 = (제품 & 0xffffffffffff) + 5 * (제품 >> 32);  // 0x33333333 = 858,993,459보다 큰 승수가 필요할 것이다.   // 곱셈 결과는 32비트에 맞지만 합계는 33비트가 될 수 있다.  상품 = (상품 & 0xffffffff) + 5 * (uint32_t)(상품 >> 32);   상품 += 4;  // 이 합계는 32비트가 보장된다.  x = (uint32_t)상품 + 5 * (uint32_t)(상품 >> 32);  돌아오다 *주 = x - 4; }

많은 다른 레머 발전기들은 좋은 성질을 가지고 있다.다음의 modulo-2128 Lehmer 발생기는 컴파일러의 128비트 지원을 필요로 하며 L'Ecuyer가 계산한 승수를 사용한다.[17]이 기간은126 2:

정태의 서명이 없는 __integraph 주;  /* 상태는 홀수 값으로 시드되어야 한다.*/ 공허하게 하다 씨를 뿌리다(서명이 없는 __integraph 씨를 뿌리다) {  주 = 씨를 뿌리다 << 1   1; }  uint64_t 다음에(공허하게 하다) {  // GCC는 128비트 리터럴을 쓸 수 없기 때문에 우리는 표현을 사용한다.  경시하다 서명이 없는 __integraph 다중의 =   (서명이 없는 __integraph)0x12e15e35b500f16e << 64     0x2e714eb2b37916a5;  주 *= 다중의;  돌아오다 주 >> 64; }

발전기는 홀수 128비트 값을 계산하고 상위 64비트를 반환한다.

이 발전기는 TestU01에서 BigCrush를 통과하지만 LauncleRand에서 TMFn 테스트를 통과하지 못한다.이 시험은 이러한 유형의 발전기의 결함을 정확히 포착하기 위해 설계되었다. 계수는 2의 전력이기 때문에 출력에서 가장 낮은 비트의 주기는 2가126 아니라 2에62 불과하다.계량 2의 검정력을 갖는 선형 응집 발전기는 유사한 동작을 가진다.

다음 핵심 루틴은 정수 워크로드에 대해 위의 코드 속도를 개선한다(컴파일러가 계산 루프를 통해 상수 선언을 최적화할 수 있는 경우). 

uint64_t 다음에(공허하게 하다) {  uint64_t 결과 = 주 >> 64;  // GCC는 128비트 리터럴을 쓸 수 없기 때문에 우리는 표현을 사용한다.  경시하다 서명이 없는 __integraph 다중의 =   (서명이 없는 __integraph)0x12e15e35b500f16e << 64     0x2e714eb2b37916a5;  주 *= 다중의;  돌아오다 결과; }

그러나 곱셈은 지연되기 때문에 첫 번째 호출은 단순히 시드 상태의 상위 64비트를 반환하기 때문에 해싱에는 적합하지 않다.

참조

  1. ^ 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.
  2. ^ 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.
  3. ^ 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.
  4. ^ 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.
  5. ^ Sullivan, Stephen (1993). "Technical correspondence: Another test for randomness" (PDF). Communications of the ACM. 36 (7): 108. doi:10.1145/159544.376068.
  6. ^ 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.
  7. ^ Vickers, Steve (1981). "Chapter 5 - Functions". ZX81 Basic Programming (2nd ed.). Sinclair Research Ltd. The ZX81 uses p=65537 & a=75 [...] (ZX81 설명서에 65537이 2-1과16 같은 메르센 프라임이라고 잘못 기재되어 있다는 점에 유의하십시오.ZX Spectrum 매뉴얼은 216+1과 동일한 Fermat 프라임임을 올바르게 명시하고 수정하였다.)
  8. ^ a b GNU 과학 라이브러리:기타 난수 생성기
  9. ^ Knuth, Donald (1981). Seminumerical Algorithms. The Art of Computer Programming. Vol. 2 (2nd ed.). Reading, MA: Addison-Wesley Professional. pp. 12–14.
  10. ^ 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)!
  11. ^ 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.
  12. ^ 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. JSTOR 2938714.
  13. ^ Schrage, Linus (June 1979). "A More Portable Fortran Random Number Generator" (PDF). ACM Transactions on Mathematical Software. 5 (2): 132–138. CiteSeerX 10.1.1.470.6958. doi:10.1145/355826.355828.
  14. ^ Jain, Raj (9 July 2010). "Computer Systems Performance Analysis Chapter 26: Random-Number Generation" (PDF). pp. 19–22. Retrieved 2017-10-31.
  15. ^ 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. 이것은 상수에 의한 모듈형 곱셈의 여러 가지 다른 구현을 탐구한다.
  16. ^ Fenerty, Paul (11 September 2006). "Schrage's Method". Retrieved 2017-10-31.
  17. ^ 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.

외부 링크