Generator liczb losowych Lehmera

Generator liczb losowych Lehmera (nazwany na cześć DH Lehmera ), czasami nazywany także generatorem liczb losowych Parka-Millera (od nazwiska Stephena K. Parka i Keitha W. Millera ), jest rodzajem liniowego generatora kongruencji (LCG), który działa w multiplikatywna grupa liczb całkowitych modulo n . Ogólna formuła to

{\ Displaystyle X_ {k + 1} = a \ cdot X_ {k} {\ bmod {m.}}}

₀ gdzie moduł m jest liczbą pierwszą lub potęgą liczby pierwszej , mnożnik a jest elementem wysokiego rzędu modulo m (np. prymitywny pierwiastek modulo n ), a ziarno X jest względnie pierwsze do m .

Inne nazwy to multiplikatywny generator kongruencji liniowej (MLCG) i multiplikatywny generator kongruencji (MCG) .

Parametry w powszechnym użyciu

W 1988 Park i Miller zasugerowali RNG Lehmera o określonych parametrach m = 2 ³¹ - 1 = 2 147 483 647 ( liczba pierwsza Mersenne'a M ₃₁ ) i a = 7 ⁵ = 16 807 (prymitywny pierwiastek modulo M ₃₁ ), obecnie znany jako MINSTD . Chociaż MINSTD został później skrytykowany przez Marsaglię i Sullivana ( 1993), nadal jest używany (w szczególności w minstd_rand0 w CarbonLib i C ++ 11 ). Park, Miller i Stockmeyer odpowiedzieli na krytykę (1993), mówiąc:

Biorąc pod uwagę dynamiczny charakter tego obszaru, osobom niebędącym specjalistami trudno jest podjąć decyzję o tym, jakiego generatora użyć. „Daj mi coś, co mogę zrozumieć, wdrożyć i przenieść… to nie musi być najnowocześniejsze, po prostu upewnij się, że jest w miarę dobre i wydajne”. Nasz artykuł i powiązany z nim generator minimalnych standardów był próbą odpowiedzi na tę prośbę. Pięć lat później nie widzimy potrzeby zmiany naszej odpowiedzi poza zasugerowaniem użycia mnożnika a = 48271 zamiast 16807.

Ta poprawiona stała jest używana w generatorze liczb losowych minstd_rand języka C++ 11 .

Sinclair ZX81 i jego następcy używają RNG Lehmera z parametrami m = 2 ¹⁶ + 1 = 65,537 ( liczba pierwsza Fermata F ₄ ) i a = 75 (pierwotny pierwiastek modulo F ₄ ). Generator liczb losowych CRAY RANF to RNG Lehmera o module potęgi dwóch m = 2 ⁴⁸ i a = 44 485 709 377 909. Biblioteka naukowa GNU zawiera kilka generatorów liczb losowych w postaci Lehmera, w tym MINSTD, RANF i niesławny generator liczb losowych IBM RANDU .

Wybór modułu

₀ Najczęściej moduł jest wybierany jako liczba pierwsza, przez co wybór ziarna względnie pierwszego jest trywialny (wystarczy dowolne 0 < X < m ) . Daje to dane wyjściowe najlepszej jakości, ale wprowadza pewną złożoność implementacji, a zakres danych wyjściowych raczej nie pasuje do pożądanej aplikacji; przeliczenie na żądany zakres wymaga dodatkowego pomnożenia.

Użycie modułu m , który jest potęgą dwójki , zapewnia szczególnie wygodną implementację komputerową, ale wiąże się z kosztami: okres wynosi co najwyżej m /4, a niższe bity mają okresy krótsze. Dzieje się tak, ponieważ najniższe k bitów same tworzą generator modulo-2 ^{k ;} bity wyższego rzędu nigdy nie wpływają na bity niższego rzędu. Wartości Xi _{_} są zawsze nieparzyste (bit 0 nigdy się nie zmienia), bity 2 i 1 występują naprzemiennie (3 niższe bity powtarzają się z okresem 2), dolne 4 bity powtarzają się z okresem 4 i tak dalej. Dlatego aplikacja korzystająca z tych liczb losowych musi używać najbardziej znaczących bitów; redukcja do mniejszego zakresu za pomocą operacji modulo z parzystym modułem da katastrofalne rezultaty.

₀ Aby osiągnąć ten okres, mnożnik musi spełniać a ≡ ± 3 (mod 8), a ziarno X musi być nieparzyste.

₀ Użycie modułu złożonego jest możliwe, ale generator musi być zaszczepiony wartością względnie pierwszą do m , w przeciwnym razie okres zostanie znacznie skrócony. Na przykład moduł F ₅ = 2 ³² + 1 może wydawać się atrakcyjny, ponieważ dane wyjściowe można łatwo odwzorować na 32-bitowe słowo 0 ≤ X _i − 1 < 2 ³² . Jednak ziarno X = 6700417 (które dzieli 2 ³² + 1) lub dowolna wielokrotność prowadziłoby do wyniku o okresie tylko 640.

Bardziej popularną implementacją dla dużych okresów jest połączony liniowy generator kongruencji ; łączenie (np. sumowanie ich wyjść) kilku generatorów jest równoważne z wyjściem jednego generatora, którego moduł jest iloczynem modułów generatorów składowych. i którego okres jest najmniejszą wspólną wielokrotnością okresów składowych. Chociaż okresy będą miały wspólny dzielnik 2, moduły można wybrać tak, aby był to jedyny wspólny dzielnik, a wynikowy okres to ( m ₁ − 1)( m ₂ − 1)···( m _k − 1)/2 ^{k −1} . Jednym z przykładów jest Wichmanna-Hilla .

Stosunek do LCG

₀ Chociaż RNG Lehmera można postrzegać jako szczególny przypadek liniowego generatora kongruencji z c = 0 , jest to szczególny przypadek, który implikuje pewne ograniczenia i właściwości. W szczególności dla Lehmer RNG początkowe ziarno X musi być względnie pierwsze w stosunku do modułu m , co nie jest ogólnie wymagane dla LCG. Wybór modułu m i mnożnika a jest również bardziej restrykcyjny dla Lehmer RNG. W przeciwieństwie do LCG, maksymalny okres Lehmer RNG wynosi m − 1, i jest tak, gdy m jest liczbą pierwszą, a a jest pierwiastkiem pierwotnym modulo m .

Z drugiej strony, $dyskretne$ logarytmy ( ${\ Displaystyle \ varphi (m)}$ podstawie a lub dowolnym pierwotnym pierwiastku modulo m ) z X _k in liniową sekwencję kongruencji modulo totient Eulera .

Realizacja

Podstawowy moduł wymaga obliczenia iloczynu o podwójnej szerokości i wyraźnego kroku redukcji. Jeśli używany jest moduł nieco mniejszy niż potęga 2 ( liczby pierwsze Mersenne'a 2 ³¹ - 1 i 2 ⁶¹ - 1, podobnie jak 2 ³² - 5 i 2 ⁶⁴ - 59), redukcja modulo m = 2 ^e - d może być zaimplementowany taniej niż ogólny podział o podwójnej szerokości przy użyciu tożsamości 2 ^e ≡ d (mod m ) .

Podstawowy krok redukcji dzieli iloczyn na dwie e -bitowe części, mnoży wyższą część przez d i dodaje je: ( ax mod 2 ^e ) + d ⌊ ax /2 ^e ⌋ . Po tym można odjąć m , aż wynik znajdzie się w zakresie. Liczba odejmowań jest ograniczona do ad / m , co można łatwo ograniczyć do jednego , jeśli d jest małe i a < m / d jest wybrany. (Ten warunek zapewnia również, że d ⌊ ax /2 ^e ⌋ jest iloczynem o pojedynczej szerokości; jeśli zostanie naruszony, należy obliczyć iloczyn o podwójnej szerokości).

Gdy moduł jest liczbą pierwszą Mersenne'a ( d = 1), procedura jest szczególnie prosta. Nie tylko mnożenie przez d jest trywialne, ale warunkowe odejmowanie można zastąpić bezwarunkowym przesunięciem i dodawaniem. Aby to zobaczyć, zauważ, że algorytm gwarantuje, że x ≢ 0 (mod m ) , co oznacza, że x = 0 i x = m są niemożliwe. Pozwala to uniknąć konieczności rozważenia równoważnego e -bitowe reprezentacje stanu; tylko wartości, w których wysokie bity są niezerowe, wymagają redukcji.

Niższe e bity iloczynu ax nie mogą reprezentować wartości większej niż m , a wysokie bity nigdy nie będą miały wartości większej niż a − 1 ≤ m − 2. Zatem pierwszy krok redukcji daje wartość co najwyżej m + a − 1 ≤ 2 m − 2 = 2 ^{e +1} − 4. Jest to liczba bitowa ( e + 1), która może być większa niż m (tj. może mieć ustawiony bit e ), ale górna połowa to co najwyżej 1, oraz jeśli tak, niskie e bitów będzie ściśle mniejsza niż m . Zatem niezależnie od tego, czy wysoki bit ma wartość 1, czy 0, drugi krok redukcji (dodanie połówek) nigdy nie spowoduje przepełnienia e bitów, a suma będzie wymaganą wartością.

Jeśli d > 1, można również uniknąć odejmowania warunkowego, ale procedura jest bardziej skomplikowana. Podstawowe wyzwanie modułu takiego jak 2 ³² − 5 polega na zapewnieniu, że tworzymy tylko jedną reprezentację dla wartości takich jak 1 ≡ 2 ³² − 4. Rozwiązaniem jest tymczasowe dodanie d , tak aby zakres możliwych wartości wynosił od d do 2 ^e - 1 i zmniejsz wartości większe niż e bitów w sposób, który nigdy nie generuje reprezentacji mniejszych niż d . Ostatecznie odjęcie tymczasowego przesunięcia daje pożądaną wartość.

Zacznijmy od założenia, że mamy częściowo zredukowaną wartość y ograniczoną tak, że 0 ≤ y < 2 m = 2 ^{e +1} − 2 d . W tym przypadku pojedynczy krok odejmowania przesunięcia da 0 ≤ y ′ = (( y + d ) mod 2 ^e ) + d ⌊ ( y + d )/2 ^e ⌋ − d < m . Aby to zobaczyć, rozważ dwa przypadki:

0 ≤ < m = 2 ^mi - re: W tym przypadku y + re < 2 ^mi y y ′ = y < m , zgodnie z życzeniem.
m ≤ y < 2 m: W tym przypadku 2 ^e ≤ y + d < 2 ^{e +1} to ( e + 1)-bitowa liczba i ⌊ ( y + d )/2 ^mi ⌋ = 1. Zatem y ′ = ( y + re ) - 2 ^mi + re - re = y - 2 ^mi + re = y - m < m , zgodnie z życzeniem. Ponieważ pomnożona wysoka część to d , suma wynosi co najmniej d , a odjęcie przesunięcia nigdy nie powoduje niedomiaru.

(W szczególności w przypadku generatora Lehmera stan zerowy lub jego obraz y = m nigdy nie wystąpią, więc przesunięcie d - 1 będzie działać tak samo, jeśli jest to wygodniejsze. Zmniejsza to przesunięcie do 0 w Przypadek pierwszy Mersenne'a, gdy d = 1.)

Zmniejszenie większej osi produktu do mniej niż 2 m = 2 ^{e +1} - 2 d można wykonać za pomocą jednego lub więcej kroków redukcji bez przesunięcia.

Jeśli ad ≤ m wystarczy jeden dodatkowy krok redukcji. Ponieważ x < m , ax < am ≤ ( a − 1)2 ^e , a jeden krok redukcji przekształca to najwyżej w 2 ^e − 1 + ( a − 1) d = m + ad − 1. To mieści się w granicach 2 m jeśli ad − 1 < m , co jest początkowym założeniem.

Jeśli ad > m , to możliwe jest, że pierwszy krok redukcji da sumę większą niż 2 m = 2 ^{e +1} - 2 d , co jest zbyt duże dla końcowego etapu redukcji. (Wymaga to również pomnożenia przez d , aby otrzymać iloczyn większy niż e bitów, jak wspomniano powyżej). Jednak dopóki d ² < 2 ^e , pierwsza redukcja da wartość w zakresie wymaganym dla poprzedniego przypadku dwóch kroki redukcji, które należy zastosować.

metoda Schrage'a

Jeśli iloczyn o podwójnej szerokości nie jest dostępny, do obliczenia ax mod m można zastosować metodę Schrage'a , zwaną również metodą faktoryzacji przybliżonej , ale wiąże się to z kosztami:

Moduł musi być reprezentowalny w postaci liczby całkowitej ze znakiem ; operacje arytmetyczne muszą dopuszczać zakres ± m .
Wybór mnożnika a jest ograniczony. Wymagamy, aby m mod a ≤ ⌊ m / a ⌋ , zwykle osiągane przez wybór a ≤ √ m .
Wymagane jest jedno dzielenie (z resztą) na iterację.

Chociaż ta technika jest popularna w przenośnych implementacjach w językach wysokiego poziomu , w których brakuje operacji podwójnej szerokości, na nowoczesnych komputerach dzielenie przez stałą jest zwykle implementowane przy użyciu mnożenia podwójnej szerokości, więc tej techniki należy unikać, jeśli zależy Ci na wydajności. Nawet w językach wysokiego poziomu, jeśli mnożnik a jest ograniczony do √ m , to ax iloczynu podwójnej szerokości można obliczyć za pomocą dwóch mnożeń o pojedynczej szerokości i zmniejszyć za pomocą technik opisanych powyżej.

Aby użyć metody Schrage'a, pierwszy czynnik m = qa + r , tj. wstępnie obliczyć stałe pomocnicze r = m mod a i q = ⌊ m / a ⌋ = ( m − r )/ a . Następnie w każdej iteracji oblicz ax ≡ a ( x mod q ) − r ⌊ x / q ⌋ (mod m ) .

Ta równość zachodzi, ponieważ

{\ Displaystyle {\ rozpocząć {wyrównane} aq & = a \ cdot (mr)/a\\&=(a/a)\cdot (mr)\\&=(mr)\\&\equiv -r{\pmod {m}}\\\end{wyrównane}}}

więc jeśli rozłożymy x = ( x mod q ) + q ⌊ x / q ⌋ , otrzymamy:

{\ Displaystyle {\ rozpocząć {wyrównane} ax & = a (x {\ bmod {q}}) + aq \ lfloor x/q \ rfloor \\& = a (x {\ bmod {q}}) + (mr) \lfloor x/q\rfloor \\&\equiv a(x{\bmod {q}})-r\lfloor x/q\rfloor {\pmod {m}}\\\end{wyrównane}}}

Powodem, dla którego nie jest przepełniony, jest to, że oba wyrazy są mniejsze niż m . Ponieważ x mod q < q ≤ m / a , pierwszy wyraz jest ściśle mniejszy niż am / a = m i można go obliczyć za pomocą iloczynu o pojedynczej szerokości.

Jeśli a jest wybrane tak, że r ≤ q (a zatem r / q ≤ 1), to drugi wyraz jest również mniejszy niż m : r ⌊ x / q ⌋ ≤ rx / q = x ( r / q ) ≤ x (1 ) = x < m . Zatem różnica mieści się w przedziale [1− m , m −1] i można ją sprowadzić do [0, m −1] z pojedynczym dodaniem warunkowym.

Technikę tę można rozszerzyć, aby umożliwić ujemne r (- q ≤ r <0), zmieniając ostateczną redukcję na odejmowanie warunkowe.

Technikę można również rozszerzyć, aby umożliwić większe a , stosując ją rekurencyjnie. Z dwóch terminów odjętych w celu uzyskania końcowego wyniku tylko drugi ( r ⌊ x / q ⌋ ) grozi przepełnieniem. Ale to samo w sobie jest modułowym mnożeniem przez stałą czasu kompilacji r i może być zaimplementowane tą samą techniką. Ponieważ każdy krok zmniejsza średnio o połowę wielkość mnożnika (0 ≤ r < a , średnia wartość ( a −1)/2), wymagałoby to jednego kroku na bit i byłoby spektakularnie nieefektywne. Jednak każdy krok dzieli również x przez stale rosnący iloraz q = ⌊ m / a ⌋ i szybko osiąga się punkt, w którym argument wynosi 0, a rekurencja może zostać zakończona.

Przykładowy kod C99

Używając kodu C , Park-Miller RNG można zapisać w następujący sposób:

  

	       
 uint32_t  lcg_parkmiller  (  uint32_t  *  stan  )  {  powrót  *  stan  =  (  uint64_t  )  *  stan  *  48271  %  0x7fffffff  ;  }

Ta funkcja może być wywoływana wielokrotnie w celu generowania liczb pseudolosowych, o ile osoba wywołująca jest ostrożna, aby zainicjować stan na dowolną liczbę większą od zera i mniejszą niż moduł. W tej implementacji wymagana jest arytmetyka 64-bitowa; w przeciwnym razie iloczyn dwóch 32-bitowych liczb całkowitych może się przepełnić.

Aby uniknąć podziału 64-bitowego, wykonaj redukcję ręcznie:

  

	     
	         

	        uint32_t  lcg_parkmiller  (  uint32_t  *  stan  )  {  uint64_t  produkt  =  (  uint64_t  )  *  stan  *  48271  ;  uint32_t  x  =  (  produkt  &  0x7fffffff  )  +  (  produkt  >>  31  );  x  =  (  x  &  0x7fffffff  )  +  (  x  >>  
	   
 31  );  powrót  *  stan  =  x  ;  }

Aby użyć tylko arytmetyki 32-bitowej, użyj metody Schrage:

  

	
	    
	    
	          
	          

	  uint32_t  lcg_parkmiller  (  uint32_t  *  stan  )  {  // Wstępnie obliczone parametry metody Schrage  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     	
	     	

	     	
	     	
	      =  *  stan  /  Q  ;  // max: M / Q = A = 48,271  uint32_t  rem  =  *  state  %  Q  ;  // maks.: Q - 1 = 44 487  int32_t  s  =  rem  *  A  ;  // maks.: 44 487 * 48 271 = 2 147 431 977 = 0x7fff3629  int32_t  t  =  div  *  R  ;  // maks.: 48 271 * 3 399 = 164 073 129  int32_t  wynik  =  s  -  t  ; 

	   0
		  

	   
 jeśli  (  wynik  <  )  wynik  +=  M  ;  powrót  *  stan  =  wynik  ;  }

lub użyj dwóch mnożników 16 × 16-bitowych:

  

	    

	        			
	           uint32_t  lcg_parkmiller  (  stan  uint32_t  *  )  {  const  uint32_t  A  =  48271  ;  uint32_t  low  =  (  *  stan  &  0x7fff  )  *  A  ;  // maks.: 32 767 * 48 271 = 1 581 695 857 = 0x5e46c371  uint32_t  high  =  (  *  stan  >>  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  );  powrót  *  stan  =  x  ; 
 }

Inny popularny generator Lehmera wykorzystuje główny moduł 2 ³² −5:

  

           
 uint32_t  lcg_rand  (  uint32_t  *  stan  )  {  return  *  stan  =  (  uint64_t  )  *  stan  *  279470273u  %  0xfffffffb  ;  }

Można to również zapisać bez dzielenia 64-bitowego:

  

	     
	 

	
	
	

	
	 uint32_t  lcg_rand  (  uint32_t  *  stan  )  {  uint64_t  produkt  =  (  uint64_t  )  *  stan  *  279470273u  ;  uint32_t  x  ;  // Niewymagane, ponieważ 5 * 279470273 = 0x5349e3c5 mieści się w 32 bitach.  // produkt = (produkt & 0xffffffff) + 5 * (produkt >> 32);  // Potrzebny byłby mnożnik większy niż 0x33333333 = 858 993 459.  // Wynik mnożenia mieści się w 32 bitach, ale suma może wynosić 33 bity.  produkt           

	  
	
	        
	     =  (  produkt  &  0xffffffff  )  +  5  *  (  uint32_t  )(  produkt  >>  32  );  iloczyn  +=  4  ;  // Ta suma to gwarantowane 32 bity.  x  =  (  uint32_t  )  produkt  +  5  *  (  uint32_t  )(  produkt  >>  32  );  powrót  *  stan  =  x  -  
 4  ;  }

Wiele innych generatorów Lehmer ma dobre właściwości. Poniższy generator modulo-2 ¹²⁸ Lehmer wymaga obsługi 128-bitowej przez kompilator i wykorzystuje mnożnik obliczony przez L'Ecuyer. Ma okres 2 ¹²⁶ :

   


   

	      


 

	
	    
		 statyczny  stan  bez znaku  __int128  ;  /* Stan musi być zaszczepiony nieparzystą wartością. */   void  seed  (  unsigned  __int128  seed  )  {  state  =  seed  <<  1  |  1  ;  }  uint64_t  next  (  void  )  {  // GCC nie może zapisywać literałów 128-bitowych, więc używamy wyrażenia  const  unsigned  __int128  mult  =  (  unsigned     
		
	  
	   
 __int128  )  0x12e15e35b500f16e  <<  64  |  0x2e714eb2b37916a5  ;  stan  *=  mnożenie  ;  stan  powrotu  >>  64  ;  }

Generator oblicza nieparzystą wartość 128-bitową i zwraca jej górne 64 bity.

Ten generator przechodzi test BigCrush z TestU01 , ale nie przechodzi testu TMFn z PractRand . Ten test został zaprojektowany, aby dokładnie wykryć wadę tego typu generatora: ponieważ moduł jest potęgą 2, okres najniższego bitu na wyjściu wynosi tylko 2 62 , a ^nie 2 ¹²⁶ . Liniowe generatory kongruencji z modułem potęgi 2 mają podobne zachowanie.

Następująca podstawowa procedura poprawia szybkość powyższego kodu dla obciążeń całkowitych (jeśli kompilator zezwala na optymalizację deklaracji stałej poza pętlą obliczeniową):

 

	     
	
	    
		    
		
	  
	 
 uint64_t  następny  (  pusty  )  {  uint64_t  wynik  =  stan  >>  64  ;  // GCC nie może pisać literałów 128-bitowych, więc używamy wyrażenia  const  unsigned  __int128  mult  =  (  unsigned  __int128  )  0x12e15e35b500f16e  <<  64  |  0x2e714eb2b37916a5  ;  stan  *=  mnożenie  ;  zwróć  wynik  ;  }

Jednakże, ponieważ mnożenie jest odroczone, nie nadaje się do mieszania, ponieważ pierwsze wywołanie po prostu zwraca górne 64 bity stanu początkowego.

Lehmer, DH (1949). „Metody matematyczne w wielkoskalowych jednostkach obliczeniowych”. Materiały z drugiego sympozjum na temat wielkoskalowych cyfrowych maszyn liczących . s. 141 –146. MR 0044899 . (wersja dziennikowa: Annals of the Computation Laboratory of Harvard University , Vol. 26 (1951)).
Steve Park, Generatory liczb losowych

Linki zewnętrzne

Liczby pierwsze mniejsze niż potęga dwójki mogą być przydatne przy wyborze modułów. Część Prime Pages .