Lehmer slumptalsgenerator

Lehmer slumptalsgenerator (uppkallad efter DH Lehmer ), ibland även kallad Park–Miller slumptalsgenerator (efter Stephen K. Park och Keith W. Miller), är en typ av linjär kongruentialgenerator (LCG) som verkar i multiplikativ grupp av heltal modulo n . Den allmänna formeln är

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

₀ där modulen m är ett primtal eller en potens av ett primtal , multiplikatorn a är ett element av hög multiplikativ ordning modulo m (t.ex. en primitiv rot modulo n ), och fröet X är coprime till m .

Andra namn är multiplikativ linjär kongruentialgenerator (MLCG) och multiplikativ kongruentialgenerator (MCG) .

Parametrar i vanligt bruk

1988 föreslog Park och Miller en Lehmer RNG med speciella parametrar m = 2 ³¹ − 1 = 2 147 483 647 (en Mersenne prime M ₃₁ ) och a = 7 ⁵ = 16 807 (en primitiv rot modulo M ₃₁ ), nu känd som MINSTD . Även om MINSTD senare kritiserades av Marsaglia och Sullivan (1993), används den fortfarande idag (särskilt i CarbonLib och C++11s minstd_rand0 ) . Park, Miller och Stockmeyer svarade på kritiken (1993) och sa:

Med tanke på områdets dynamiska natur är det svårt för icke-specialister att fatta beslut om vilken generator som ska användas. "Ge mig något jag kan förstå, implementera och porta... det behöver inte vara toppmodernt, bara se till att det är någorlunda bra och effektivt." Vår artikel och den tillhörande minimala standardgeneratorn var ett försök att svara på denna begäran. Fem år senare ser vi inget behov av att ändra vårt svar annat än att föreslå användningen av multiplikatorn a = 48271 i stället för 16807.

Denna reviderade konstant används i C++11s minstd_rand slumptalsgenerator .

Sinclair ZX81 och dess efterföljare använder Lehmer RNG med parametrarna m = 2 ¹⁶ + 1 = 65 537 (en Fermat prime F ₄ ) och a = 75 (en primitiv rotmodulo F ₄ ). CRAY - slumptalsgeneratorn RANF är en Lehmer RNG med tvåeffektsmodulen m = 2 ⁴⁸ och a = 44 485 709 377 909. GNU Scientific Library innehåller flera slumptalsgeneratorer av Lehmer-formen, inklusive MINSTD, RANF och den ökända IBM slumptalsgeneratorn RANDU .

Val av modul

₀ Vanligast är att modulen väljs som ett primtal, vilket gör valet av ett coprime-frö trivialt (vilket som helst 0 < X < m räcker). Detta ger den bästa kvaliteten på resultatet, men introducerar en viss implementeringskomplexitet, och det är osannolikt att utdataområdet matchar den önskade applikationen; omvandling till önskat intervall kräver ytterligare en multiplikation.

Att använda en modul m som är en potens av två ger en särskilt bekväm datorimplementering, men kostar pengar: perioden är som mest m /4 och de lägre bitarna har perioder kortare än så. Detta beror på att de lägsta k bitarna bildar en modulo-2k ^generator helt av sig själva; de högre ordningens bitar påverkar aldrig lägre ordningens bitar. Värdena X _i är alltid udda (bit 0 ändras aldrig), bitar 2 och 1 växlar (de lägre 3 bitarna upprepas med en period på 2), de lägre 4 bitarna upprepas med en period på 4, och så vidare. Därför måste applikationen som använder dessa slumptal använda de mest signifikanta bitarna; att reducera till ett mindre område med en modulo-operation med en jämn modul kommer att ge katastrofala resultat.

₀ För att uppnå denna period måste multiplikatorn uppfylla en ≡ ±3 (mod 8), och frö- X måste vara udda.

₀ Det är möjligt att använda en sammansatt modul, men generatorn måste seedas med ett värde coprime till m , annars kommer perioden att reduceras kraftigt. Till exempel kan en modul på F ₅ = 2 ³² + 1 verka attraktiv, eftersom utdata enkelt kan mappas till ett 32-bitars ord 0 ≤ X _i − 1 < 2 ³² . Men ett frö på X = 6700417 (som delar 2 ³² + 1) eller vilken multipel som helst skulle leda till en utdata med en period på endast 640.

En mer populär implementering för stora perioder är en kombinerad linjär kongruentialgenerator ; att kombinera (t.ex. genom att summera deras utgångar) flera generatorer är ekvivalent med utsignalen från en enda generator vars modul är produkten av komponentgeneratorernas moduler. och vars period är den minsta gemensamma multipeln av komponentperioderna. Även om perioderna kommer att dela en gemensam divisor på 2, kan modulerna väljas så att det är den enda gemensamma divisorn och den resulterande perioden är ( m ₁ − 1)( m ₂ − 1)···( m _k − 1)/ 2 ^{k −1} . Ett exempel på detta är Wichmann–Hill -generatorn.

Relation till LCG

₀ Medan Lehmer RNG kan ses som ett särskilt fall av den linjära kongruentialgeneratorn med c = 0 , är det ett specialfall som innebär vissa restriktioner och egenskaper. I synnerhet för Lehmer RNG måste det initiala fröet X vara coprime till modulen m , vilket inte krävs för LCGs i allmänhet. Valet av modulen m och multiplikatorn a är också mer restriktivt för Lehmer RNG. I motsats till LCG är den maximala perioden för Lehmer RNG lika med m − 1, och den är sådan när m är primtal och a är en primitiv rot modulo m .

Å andra sidan representerar de diskreta logaritmerna (för att basera en eller valfri primitiv rot modulo m ) av X _k i $\mathbb {Z} _{m}$ en linjär kongruentiell sekvens modulo Euler totient $\varphi (m)$ .

Genomförande

En primmodul kräver beräkning av en dubbelbreddsprodukt och ett explicit reduktionssteg. Om en modul bara mindre än en potens av 2 används ( Mersennes primtal 2 ³¹ − 1 och 2 ⁶¹ − 1 är populära, liksom 2 ³² − 5 och 2 ^{64 − 59),} kan reduktionsmodulo m = 2 ^e − d implementeras billigare än en allmän dubbelbreddsindelning med identiteten 2 ^e ≡ d (mod m ) .

Det grundläggande reduktionssteget delar upp produkten i två e -bitdelar, multiplicerar den höga delen med d och adderar dem: ( ax mod 2 ^e ) + d ⌊ ax /2 ^e ⌋ . Detta kan följas genom att subtrahera m tills resultatet är inom intervallet. Antalet subtraktioner är begränsat till ad / m , vilket enkelt kan begränsas till en om d är litet och a < m / d är valt. (Detta villkor säkerställer också att d ⌊ ax /2 ^e ⌋ är en enkelbreddsprodukt; om den överträds måste en dubbelbreddsprodukt beräknas.)

När modulen är ett Mersenne-primtal ( d = 1), är proceduren särskilt enkel. Multiplikation med d är inte bara trivial, utan den villkorliga subtraktionen kan ersättas med en ovillkorlig förskjutning och addition. För att se detta, notera att algoritmen garanterar att x ≢ 0 (mod m ) , vilket betyder att x = 0 och x = m båda är omöjliga. Detta undviker behovet av att överväga likvärdiga e -bitars representationer av tillståndet; endast värden där de höga bitarna inte är noll behöver reduktion.

Produktaxelns låga e- bitar kan inte representera ett värde som är större än m , och de höga bitarna kommer aldrig att hålla ett värde större än a − 1 ≤ m − 2. Det första reduktionssteget ger alltså högst ett värde m + a − 1 ≤ 2 m − 2 = 2 ^{e +1} − 4. Detta är ett ( e + 1)-bittal, som kan vara större än m (dvs. kan ha bit e set), men den övre halvan är högst 1, och om så är fallet kommer de låga e- bitarna att vara strikt mindre än m . Så oavsett om den höga biten är 1 eller 0, kommer ett andra reduktionssteg (tillägg av halvorna) aldrig att svämma över e bitar, och summan kommer att vara det önskade värdet.

Om d > 1 kan villkorlig subtraktion också undvikas, men proceduren är mer intrikat. Den grundläggande utmaningen med en modul som 2 ³² − 5 ligger i att säkerställa att vi endast producerar en representation för värden som 1 ≡ 2 ³² − 4. Lösningen är att tillfälligt addera d , så att intervallet av möjliga värden är d till 2 ^e − 1, och reducera värden större än e bitar på ett sätt som aldrig genererar representationer mindre än d . Att slutligen subtrahera den temporära offset ger det önskade värdet.

Börja med att anta att vi har ett delvis reducerat värde y avgränsat så att 0 ≤ y < 2 m = 2 ^{e +1} − 2 d . I detta fall kommer ett enda offset-subtraktionssteg att producera 0 ≤ y ′ = (( y + d ) mod 2 ^e ) + d ⌊ ( y + d )/2 ^e ⌋ − d < m . För att se detta, överväg två fall:

0 ≤ y < m = 2 ^e − d: I detta fall, y + d < 2 ^e och y ′ = y < m , enligt önskemål.
m ≤ y < 2 m: I det här fallet är 2 ^e ≤ y + d < 2 ^{e +1} ett ( e + 1)-bittal och ⌊ ( y + d )/2 ^e ⌋ = 1. Alltså y ′ = ( y + d ) − 2 ^e + d − d = y − 2 ^e + d = y − m < m , enligt önskemål. Eftersom den multiplicerade höga delen är d , är summan minst d , och subtrahering av offset orsakar aldrig underflöde.

(För fallet med en Lehmer-generator specifikt, kommer ett nolltillstånd eller dess bild y = m aldrig att inträffa, så en offset på d − 1 kommer att fungera precis likadant, om det är mer praktiskt. Detta minskar förskjutningen till 0 i Mersenne primtal, när d = 1.)

Att reducera en större produktyxa till mindre än 2 m = 2 ^{e +1} − 2 d kan göras genom ett eller flera reduktionssteg utan förskjutning.

Om ad ≤ m räcker ett ytterligare reduktionssteg. Eftersom x < m , ax < am ≤ ( a − 1)2 ^e , och ett reduktionssteg omvandlar detta till högst 2 ^e − 1 + ( a − 1) d = m + ad − 1. Detta är inom gränsen för 2 m om ad − 1 < m , vilket är det initiala antagandet.

Om ad > m , så är det möjligt för det första reduktionssteget att producera en summa som är större än 2 m = 2 ^{e +1} − 2 d , vilket är för stort för det sista reduktionssteget. (Det kräver också multiplikation med d för att producera en produkt som är större än e bitar, som nämnts ovan.) Så länge som d ² < 2 ^e , kommer emellertid den första reduktionen att producera ett värde inom det område som krävs för det föregående fallet med två minskningssteg att tillämpa.

Schrages metod

Om en produkt med dubbel bredd inte är tillgänglig kan Schrages metod , även kallad den approximativa faktoriseringsmetoden, användas för att beräkna ax mod m , men detta kommer till kostnaden:

Modulen måste kunna representeras i ett heltal med tecken ; aritmetiska operationer måste tillåta ett intervall på ± m .
Valet av multiplikator a är begränsat. Vi kräver att m mod a ≤ ⌊ m / a ⌋ , vanligtvis uppnås genom att välja a ≤ √ m .
En division (med resten) per iteration krävs.

Även om den här tekniken är populär för bärbara implementeringar på högnivåspråk som saknar dubbelbreddsoperationer, implementeras division med en konstant på moderna datorer vanligtvis med dubbelbreddsmultiplikation, så denna teknik bör undvikas om effektivitet är ett problem. Även i högnivåspråk, om multiplikatorn a är begränsad till √ m , kan den dubbelbredda produktaxeln beräknas med två enkelbreddsmultiplikationer och reduceras med de tekniker som beskrivs ovan.

För att använda Schrages metod, först faktorn m = qa + r , dvs förberäkna hjälpkonstanterna r = m mod a och q = ⌊ m / a ⌋ = ( m − r )/ a . Beräkna sedan, varje iteration, ax ≡ a ( x mod q ) − r ⌊ x / q ⌋ (mod m ) .

Denna jämlikhet gäller eftersom

{\begin{aligned}aq&=a\ cdot (mr)/a\\&=(a/a)\cdot (mr)\\&=(mr)\\&\equiv -r{\pmod {m}}\\\end{aligned}}

så om vi faktoriserar x = ( x mod q ) + q ⌊ x / q ⌋ , får vi:

{\begin{aligned}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{aligned}}

Anledningen till att det inte svämmar över är att båda termerna är mindre än m . Eftersom x mod q < q ≤ m / a , är den första termen strikt mindre än am / a = m och kan beräknas med en enkelbreddsprodukt.

Om a väljs så att r ≤ q (och därmed r / q ≤ 1), så är den andra termen också mindre än m : r ⌊ x / q ⌋ ≤ rx / q = x ( r / q ) ≤ x (1 ) = x < m . Skillnaden ligger alltså i området [1− m , m −1] och kan reduceras till [0, m −1] med en enda villkorlig addering.

Denna teknik kan utökas för att tillåta ett negativt r (− q ≤ r < 0), vilket ändrar den slutliga reduktionen till ett villkorligt subtrahera.

Tekniken kan också utökas för att tillåta större a genom att applicera den rekursivt. Av de två termer som subtraheras för att producera det slutliga resultatet är det bara den andra ( r ⌊ x / q ⌋ ) som riskerar att svämma över. Men detta är i sig en modulär multiplikation med en kompileringstidskonstant r och kan implementeras med samma teknik. Eftersom varje steg i genomsnitt halverar storleken på multiplikatorn (0 ≤ r < a , medelvärde ( a −1)/2), verkar detta kräva ett steg per bit och vara spektakulärt ineffektivt. Men varje steg delar också x med en ständigt ökande kvot q = ⌊ m / a ⌋ , och snabbt nås en punkt där argumentet är 0 och rekursionen kan avslutas.

Exempel på C99-kod

Med C- kod kan Park-Miller RNG skrivas enligt följande:

  

	       
 uint32_t  lcg_parkmiller  (  uint32_t  *  state  )  {  return  *  state  =  (  uint64_t  )  *  state  *  48271  %  0x7ffffffff  ;  }

Denna funktion kan anropas upprepade gånger för att generera pseudoslumptal, så länge som den som ringer är noga med att initiera tillståndet till valfritt tal som är större än noll och mindre än modulen. I denna implementering krävs 64-bitars aritmetik; annars kan produkten av två 32-bitars heltal svämma över.

För att undvika 64-bitarsdelningen, gör reduktionen för hand:

  

	     
	         

	        
	   
 uint32_t  lcg_parkmiller  (  uint32_t  *  state  )  {  uint64_t  product  =  (  uint64_t  )  *  state  *  48271  ;  uint32_t  x  =  (  produkt  &  0x7fffffff  )  +  (  produkt  >>  31  );  x  =  (  x  &  0x7ffffffff  )  +  (  x  >>  31  );  return  *  state  =  x  ;  }

För att bara använda 32-bitars aritmetik, använd Schrages metod:

  

	
	    
	    
	          
	          

	     	
	     	

	     	
	     	
	     

	   0
		  

	   
 uint32_t  lcg_parkmiller  (  uint32_t  *  state  )  {  // Förberäknade parametrar för Schrages metod  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  =  *  tillstånd  /  Q  ;  // max: M / Q = A = 48 271  uint32_t  rem  =  *  tillstånd  %  Q  ;  // max: Q - 1 = 44,487  int32_t  s  =  rem  *  A  ;  // max: 44 487 * 48 271 = 2 147 431 977 = 0x7fff3629  int32_t  t  =  div  *  R  ;  // max: 48 271 * 3 399 = 164 073 129  int32_t  resultat  =  s  -  t  ;  if  (  resultat  <  )  resultat  +=  M  ;  return  *  state  =  resultat  ;  }

eller använd två 16×16-bitars multiplikationer:

  

	    

	        			
	          			
	             	

	        
	   
 uint32_t  lcg_parkmiller  (  uint32_t  *  state  )  {  const  uint32_t  A  =  48271  ;  uint32_t  low  =  (  *  tillstånd  &  0x7fff  )  *  A  ;  // max: 32 767 * 48 271 = 1 581 695 857 = 0x5e46c371  uint32_t  hög  =  (  *  tillstånd  >>  15  )  *  A  ;  // max: 65 535 * 48 271 = 3 163 439 985 = 0xbc8e4371  uint32_t  x  =  låg  +  ((  hög  &  0xffff  )  <<  15  )  +  (  hög  >>  16  );  // max: 0x5e46c371 + 0x7fff8000 + 0xbc8e = 0xde46ffff  x  =  (  x  &  0x7fffffff  )  +  (  x  >>  31  );  return  *  state  =  x  ;  }

En annan populär Lehmer-generator använder primmodulen 2 ³² −5:

  

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

Detta kan också skrivas utan en 64-bitars division:

  

	     
	 

	
	
	

	
	          

	  
	
	        
	     
 uint32_t  lcg_rand  (  uint32_t  *  state  )  {  uint64_t  product  =  (  uint64_t  )  *  state  *  279470273u  ;  uint32_t  x  ;  // Krävs inte eftersom 5 * 279470273 = 0x5349e3c5 passar i 32 bitar.  // produkt = (produkt & 0xffffffff) + 5 * (produkt >> 32);  // En multiplikator större än 0x33333333 = 858 993 459 skulle behöva den.  // Multipliceringsresultatet passar in i 32 bitar, men summan kan vara 33 bitar.  produkt  =  (  produkt  &  0xffffffff  )  +  5  *  (  uint32_t  )(  produkt  >>  32  );  produkt  +=  4  ;  // Denna summa är garanterat 32 bitar.  x  =  (  uint32_t  )  produkt  +  5  *  (  uint32_t  )(  produkt  >>  32  );  return  *  state  =  x  -  4  ;  }

Många andra Lehmer-generatorer har bra egenskaper. Följande modulo-2 ¹²⁸ Lehmer-generator kräver 128-bitars stöd från kompilatorn och använder en multiplikator som beräknas av L'Ecuyer. Den har en period på 2 ¹²⁶ :

   


   

	      


 

	
	    
		    
		
	  
	   
 statisk  osignerad  __int128-  tillstånd  ;  /* Tillståndet måste seedas med ett udda värde. */   void  seed  (  unsigned  __int128  seed  )  {  state  =  seed  <<  1  |  1  ;  }  uint64_t  next  (  void  )  {  // GCC kan inte skriva 128-bitars bokstaver, så vi använder ett uttryck  const  unsigned  __int128  mult  =  (  unsigned  __int128  )  0x12e15e35b500f16e  <<  64  |  0x2e714eb2b37916a5  ;  state  *=  mult  ;  returtillstånd  >>  64  ;  _  }

Generatorn beräknar ett udda 128-bitars värde och returnerar dess övre 64 bitar.

Denna generator klarar BigCrush från TestU01 , men klarar inte TMFn-testet från PractRand . Det testet har utformats för att fånga exakt defekten hos denna typ av generator: eftersom modulen är en potens av 2, är perioden för den lägsta biten i utgången endast 2 62 , snarare ^än 2 ¹²⁶ . Linjära kongruentialgeneratorer med en effekt-av-2-modul har ett liknande beteende.

Följande kärnrutin förbättrar hastigheten på ovanstående kod för heltalsarbetsbelastningar (om den konstanta deklarationen tillåts optimeras ur en beräkningsslinga av kompilatorn):

 

	     
	
	    
		    
		
	  
	 
 uint64_t  next  (  void  )  {  uint64_t  result  =  state  >>  64  ;  // GCC kan inte skriva 128-bitars bokstaver, så vi använder ett uttryck  const  unsigned  __int128  mult  =  (  unsigned  __int128  )  0x12e15e35b500f16e  <<  64  |  0x2e714eb2b37916a5  ;  state  *=  mult  ;  returnera  resultat  ;  }

Men eftersom multiplikationen skjuts upp är den inte lämplig för hashning, eftersom det första anropet helt enkelt returnerar de övre 64 bitarna i frötillståndet.

Lehmer, DH (1949). "Matematiska metoder i storskaliga beräkningsenheter". Handlingar från ett andra symposium om storskaliga digitala beräkningsmaskiner . s. 141 –146. MR 0044899 . (tidskriftsversion: Annals of the Computation Laboratory of Harvard University , Vol. 26 (1951)).
Steve Park, slumptalsgeneratorer

externa länkar

Primer bara mindre än en potens av två kan vara användbara för att välja moduli. En del av Prime Pages .