Generátor náhodných čísel Lehmer - Lehmer random number generator

The Generátor náhodných čísel Lehmer^[1] (pojmenoval podle D. H. Lehmer ), někdy označovaný také jako Generátor náhodných čísel Park – Miller (po Stephen K. Park a Keith W. Miller), je typ lineární shodný generátor (LCG), která působí v multiplikativní skupina celých čísel modulo n. Obecný vzorec je:

${ displaystyle X_ {k + 1} = a cdot X_ {k} ~~ { bmod {~}} ~ m}$

kde modul m je prvočíslo nebo a mocnina prvočísla, multiplikátor A je prvek vysoké multiplikativní pořadí modulo m (např primitivní kořenový modul č ) a semeno X₀ je coprime na m.

Jiná jména jsou multiplikativní lineární kongruenční generátor (MLCG)^[2] a multiplikativní kongruenční generátor (MCG).

Parametry běžně používané

V roce 1988 Park a Miller^[3] navrhl Lehmer RNG se zvláštními parametry m = 2³¹ - 1 = 2 147 483 647 (a Mersenne prime M₃₁) a A = 7⁵ = 16 807 (primitivní kořenový modul M₃₁), nyní známý jako MINSTD. Ačkoli MINSTD byl později kritizován Marsaglia a Sullivan (1993),^[4][5] dodnes se používá (zejména v CarbonLib a C ++ 11 je minstd_rand0). Park, Miller a Stockmeyer reagovali na kritiku (1993),^[6] rčení:

Vzhledem k dynamické povaze oblasti je pro nespecialisty obtížné rozhodovat o tom, jaký generátor použít. „Dej mi něco, čemu rozumím, implementuji a přenesu ... nemusí to být nejmodernější, jen se ujisti, že je to přiměřeně dobré a efektivní.“ Náš článek a přidružený generátor minimálních standardů byl pokusem odpovědět na tento požadavek. O pět let později nevidíme potřebu měnit naši reakci jinak, než navrhovat použití multiplikátoru A = 48271 místo 16807.

Tato revidovaná konstanta se používá v C ++ 11 je minstd_rand generátor náhodných čísel.

The Sinclair ZX81 a jeho nástupci používají Lehmer RNG s parametry m = 2¹⁶+1 = 65 537 (a Fermat prime F₄) a A = 75 (primitivní kořenový modul F₄).^[7]The CRAY generátor náhodných čísel RANF je Lehmer RNG s modulem síly dvou m = 2⁴⁸ a A = 44,485,709,377,909.^[8] The Vědecká knihovna GNU zahrnuje několik generátorů náhodných čísel ve formě Lehmer, včetně MINSTD, RANF a nechvalně známých IBM generátor náhodných čísel RANDU.^[8]

Volba modulu

Nejčastěji je modul zvolen jako prvočíslo, takže volba coprime seed triviální (libovolná 0 < X₀ < m udělám). To produkuje výstup nejvyšší kvality, ale přináší určitou složitost implementace a je nepravděpodobné, že by rozsah výstupu odpovídal požadované aplikaci; převod do požadovaného rozsahu vyžaduje další násobení.

Pomocí modulu m což je síla dvou umožňuje obzvláště pohodlnou počítačovou implementaci, ale stojí za to: období je nanejvýš m/ 4 a nízké bity mají období kratší než toto. Je to proto, že nízké k bity tvoří modulo-2^k generátor sami; bity vyššího řádu nikdy neovlivní bity nižšího řádu.^[9] Hodnoty X_i jsou vždy liché (bit 0 se nikdy nezmění), bity 2 a 1 se střídají (dolní 3 bity se opakují s periodou 2), dolní 4 bity se opakují s periodou 4 atd. Proto aplikace pomocí těchto náhodných čísel musí používat nejvýznamnější bity; snížení na menší rozsah pomocí operace modulo s rovnoměrným modulem přinese katastrofální výsledky.^[10]

K dosažení tohoto období musí multiplikátor uspokojit A ≡ ± 3 (mod 8)^[11] a semeno X₀ musí být liché.

Použití složeného modulu je možné, ale generátor musí být nasazen s hodnotou coprime na m, nebo se období výrazně zkrátí. Například modul o F₅ = 2³²+1 se může zdát atraktivní, protože výstupy lze snadno mapovat na 32bitové slovo 0 ≤ X_i−1 < 2³². Avšak semeno X₀ = 6700417 (což dělí 2³²+1) nebo jakýkoli násobek by vedl k výstupu s periodou pouze 640.

Více populární implementace pro velká období je a kombinovaný lineární shodný generátor; kombinace (např. sečtením jejich výstupů) několika generátorů je ekvivalentní výstupu jednoho generátoru, jehož modul je produktem modulů generátorů komponent.^[12] a jehož období je nejmenší společný násobek dílčích období. Ačkoli období budou sdílet a společný dělitel ze 2 lze zvolit moduly, takže je to jediný společný dělitel a výsledná perioda je (m₁−1)(m₂−1)(m₂···(m_k−1)/2^k−1.^[2]:744 Jedním z příkladů je Wichmann – Hill generátor.

Vztah k LCG

Zatímco Lehmer RNG lze považovat za konkrétní případ lineární shodný generátor s C=0, jedná se o speciální případ, který implikuje určitá omezení a vlastnosti. Zejména pro Lehmer RNG počáteční semeno X₀ musí být coprime do modulu m to se u LCG obecně nevyžaduje. Volba modulu m a multiplikátor A je také přísnější pro Lehmer RNG. Na rozdíl od LCG se maximální doba Lehmer RNG rovná m−1 a je to tak, když m je hlavní a A je primitivní kořenový modul m.

Na druhou stranu diskrétní logaritmy (na základnu A nebo jakýkoli primitivní kořenový modul m) z X_k v ${ displaystyle mathbb {Z} _ {m}}$ představují lineární kongruentní sekvenci modulo Eulerův totient ${ displaystyle varphi (m)}$.

Implementace

Prime modulus vyžaduje výpočet produktu s dvojnásobnou šířkou a explicitní redukční krok. Pokud je použit modul menší než síla 2 ( Mersenne připraví 2³¹-1 a 2⁶¹−1 jsou populární, stejně jako 2³²-5 a 2⁶⁴-59), redukční modulo m = 2^E − d lze implementovat levněji než obecné dělení s dvojnásobnou šířkou pomocí identity 2^E ≡ d (mod m).

Základní krok redukce rozděluje produkt na dva E-bitové části, vynásobí horní část da přidává je: (sekera mod 2^E) + d ⌊sekera/2^E⌋. Poté může následovat odečtení m dokud není výsledek v rozsahu. Počet odečtení je omezen na inzerát/m, který lze snadno omezit na jeden, pokud d je malý a A < m/d je vybrán. (Tato podmínka to také zajišťuje d ⌊sekera/2^E⌋ je produkt s jednou šířkou; pokud je porušen, musí být vypočítán produkt s dvojí šířkou.)

Když je modul Mersennovy prime (d = 1), postup je obzvláště jednoduchý. Násobení je nejen d triviální, ale podmiňovací způsob odčítání lze nahradit bezpodmínečným posunem a sčítáním. Chcete-li to vidět, všimněte si, že to algoritmus zaručuje X ≢ 0 (mod m), což znamená, že x = 0 a x =m jsou oba nemožné. Tím se vyhnete nutnosti uvažovat o ekvivalentu E-bitové reprezentace státu; pouze hodnoty, kde vysoké bity jsou nenulové, je třeba snížit.

Nízká E kousky produktu sekera nemůže představovat hodnotu větší než ma vysoké bity nikdy neudrží hodnotu větší než A−1 ≤ m − 2. První redukční krok tedy produkuje hodnotu maximálně m+A−1 ≤ 2m−2 = 2^E+1-4. Tohle je E+ 1-bitové číslo, které může být větší než m (tj. může mít bit E set), ale horní polovina je maximálně 1, a pokud ano, tak nízká E bity budou přísně menší než m. Ať už je bit vysoké hodnoty 1 nebo 0, druhý redukční krok (přidání polovin) tedy nikdy nepřeteče E bitů a součet bude požadovanou hodnotou.

Li d > 1 lze také zabránit podmíněnému odčítání, ale postup je složitější. Základní výzva modulu jako 2³²−5 spočívá v zajištění, že pro hodnoty, jako je 1 ≡ 2, vytvoříme pouze jednu reprezentaci³²-4. Řešením je dočasně přidat d takže rozsah možných hodnot je d přes 2^E-1, a snížit hodnoty větší než E bity způsobem, který nikdy negeneruje reprezentace menší než d. Nakonec odečtení dočasného posunutí vytvoří požadovanou hodnotu.

Začněte tím, že předpokládáme, že máme částečně sníženou hodnotu y který je ohraničen tak 0 ≤y < 2m = 2^E+1−2d. V tomto případě bude jediný krok odčítání offsetu produkovat 0 ≤y′ = ((y+d) mod 2^E) + d ⌊(y+d)/2^E⌋ − d < m. Chcete-li to vidět, zvažte dva případy:

0 ≤ y < m = 2^E − d

V tomto případě, y + d < 2^E a y′ = y < m, podle přání.

m ≤ y < 2m

V tomto případě 2^E ≤ y + d < 2^E+1 je E+ 1-bitové číslo a ⌊ (y+d)/2^E⌋ = 1. Tedy y′ = (y+d) − 2^E +d − d = y − 2^E + d = y − m < m, podle přání. Protože znásobená vysoká část je d, součet je minimálně d a odečtení offsetu nikdy nezpůsobí podtečení.

(V případě generátoru Lehmer konkrétně nulový stav nebo jeho obraz y = m nikdy nenastane, takže offset of d−1 bude fungovat stejně, pokud je to pohodlnější. To snižuje offset na 0 v případě Mersenne prime, když d = 1.)

Zmenšení většího produktu sekera na méně než 2m = 2^E+1 − 2d lze provést jedním nebo více redukčními kroky bez odsazení.

Li inzerát ≤ m, pak postačí jeden další redukční krok. Od té doby X < m, sekera < dopoledne ≤ (A−1)2^Ea jeden redukční krok to převede na maximálně 2^E − 1 + (A−1)d = m + inzerát - 1. To je v limitu 2m -li inzerát − 1 < m, což je počáteční předpoklad.

Li inzerát > m, pak je možné, aby první redukční krok vytvořil součet větší než 2m = 2^E+1−2d, což je pro konečný krok redukce příliš velké. (Vyžaduje také násobení pomocí d vyrábět produkt větší než E bity, jak je uvedeno výše.) Nicméně, pokud d² < 2^E, první redukce vyprodukuje hodnotu v rozsahu požadovaném pro předchozí případ použití dvou redukčních kroků.

Schrageova metoda

Pokud produkt s dvojitou šířkou není k dispozici, Schrageova metoda,^[13][14] také se nazývá metoda přibližné faktorizace,^[15] lze použít k výpočtu sekera mod m, ale to stojí za cenu:

• Modul musí být reprezentovatelný v a celé číslo se znaménkem; aritmetické operace musí umožňovat rozsah ±m.
• Volba multiplikátoru A je omezen. To požadujeme m mod A ≤ ⌊m/A⌋, běžně se dosahuje výběrem A ≤ √m.
• Je požadováno jedno rozdělení (se zbytkem) na jednu iteraci.

Zatímco tato technika je populární pro přenosné implementace v jazyky na vysoké úrovni které nemají operace s dvojnásobnou šířkou,^[2]:744 na moderních počítačích dělení konstantou je obvykle implementováno pomocí násobení s dvojnásobnou šířkou, takže této technice je třeba se vyhnout, pokud jde o efektivitu. I v jazycích vysoké úrovně, pokud je multiplikátor A je omezeno na √m, pak produkt s dvojitou šířkou sekera lze vypočítat pomocí dvou násobení jedné šířky a snížit pomocí výše popsaných technik.

Chcete-li použít Schrageovu metodu, první faktor m = qa + r, tj. předem vypočítat pomocné konstanty r = m mod A a q = ⌊m/A⌋ = (m−r)/A. Poté vypočítejte každou iteraci sekera ≡ A(X mod q) − r ⌊X/q⌋ (mod m).

Tato rovnost platí, protože

${ displaystyle { begin {aligned} aq & = a cdot (mr) / a & = (a / a) cdot (mr) & = (mr) & equiv -r { pmod {m}} konec {zarovnáno}}}$

takže pokud zohledníme X = (X mod q) + q⌊X/q⌋, dostaneme:

${ displaystyle { begin {zarovnáno} 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 {zarovnáno}}}$

Důvod, proč nepřeteče, je ten, že oba termíny jsou menší než m. Od té doby X modq < q ≤ m/A, první termín je přísně menší než dopoledne/A = m a lze jej vypočítat pomocí produktu s jednou šířkou.

Li A je vybrán tak, aby r ≤ q (a tudíž r/q ≤ 1), pak je druhý člen také menší než m: r ⌊X/q⌋ ≤ rx/q = X(r/q) ≤ X(1) = X < m. Rozdíl tedy spočívá v rozsahu [1−m, m−1] a lze jej snížit na [0,m−1] s jediným podmíněným přidáním.^[16]

Tuto techniku ​​lze rozšířit tak, aby umožňovala zápor r (−q ≤ r <0), změna konečné redukce na podmíněné odečtení.

Technika může být také rozšířena, aby umožnila větší A rekurzivním použitím.^[15]:102 Ze dvou termínů odečtených k dosažení konečného výsledku, pouze druhý (r ⌊X/qRisks) rizika přetečení. Ale toto je samo o sobě modulární násobení a časová konstanta kompilace r, a mohou být implementovány stejnou technikou. Protože každý krok v průměru zmenší velikost multiplikátoru na polovinu (0 ≤r < A, průměrná hodnota (A-1) / 2), zdálo by se, že to vyžaduje jeden krok na bit a bylo by to efektivně neúčinné. Každý krok se však také dělí X stále rostoucím podílem q = ⌊m/A⌋, a rychle je dosaženo bodu, kde je argument 0 a rekurze může být ukončena.

Ukázkový kód C99

Použitím C kódu lze Park-Miller RNG zapsat následovně:

uint32_t lcg_parkmiller(uint32_t *Stát){	vrátit se *Stát = (uint64_t)*Stát * 48271 % 0x7fffffff;}

Tuto funkci lze volat opakovaně za účelem generování pseudonáhodných čísel, pokud je volající opatrný při inicializaci stavu na libovolné číslo větší než nula a menší než modul. V této implementaci je vyžadována 64bitová aritmetika; jinak může produkt dvou 32bitových celých čísel přetéct.

Abyste se vyhnuli 64bitovému dělení, proveďte redukci ručně:

uint32_t lcg_parkmiller(uint32_t *Stát){	uint64_t produkt = (uint64_t)*Stát * 48271;	uint32_t X = (produkt & 0x7fffffff) + (produkt >> 31);	X = (X & 0x7fffffff) + (X >> 31);	vrátit se *Stát = X;}

Chcete-li použít pouze 32bitovou aritmetiku, použijte Schrageovu metodu:

uint32_t lcg_parkmiller(uint32_t *Stát){	// Předpočítané parametry pro Schrageovu metodu	konst uint32_t M = 0x7fffffff;	konst uint32_t A = 48271;	konst uint32_t Q = M / A;    // 44488	konst uint32_t R = M % A;    //  3399	uint32_t div = *Stát / Q;	// max: M / Q = A = 48 271	uint32_t rem = *Stát % 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 výsledek = s - t;	-li (výsledek < 0)		výsledek += M;	vrátit se *Stát = výsledek;}

nebo použijte dva násobky 16 × 16 bitů:

uint32_t lcg_parkmiller(uint32_t *Stát){	konst uint32_t A = 48271;	uint32_t nízký  = (*Stát & 0x7fff) * A;			// max: 32 767 * 48 271 = 1 561 695 857 = 0x5e46c371	uint32_t vysoký = (*Stát >> 15)    * A;			// max: 65 535 * 48 271 = 3 163 349 985 = 0xbc8e4371	uint32_t X = nízký + ((vysoký & 0xffff) << 15) + (vysoký >> 16);	// max: 0x5e46c371 + 0x7fff8000 + 0xbc8e = 0xde46ffff	X = (X & 0x7fffffff) + (X >> 31);	vrátit se *Stát = X;}

Další populární generátor Lehmer používá primární modul 2³²−5:

uint32_t lcg_rand(uint32_t *Stát){    vrátit se *Stát = (uint64_t)*Stát * 279470273u % 0xfffffffb;}

To lze také napsat bez 64bitového dělení:

uint32_t lcg_rand(uint32_t *Stát){	uint64_t produkt = (uint64_t)*Stát * 279470273u;	uint32_t X;	// Není vyžadováno, protože 5 * 279470273 = 0x5349e3c5 se vejde do 32 bitů.	// product = (product & 0xffffffff) + 5 * (product >> 32);	// Násobitel větší než 0x33333333 = 858 993 459 by to potřeboval.	// Výsledek násobení se vejde do 32 bitů, ale součet může být 33 bitů.	produkt = (produkt & 0xffffffff) + 5 * (uint32_t)(produkt >> 32);	produkt += 4;	// Tento součet je zaručeně 32 bitů.	X = (uint32_t)produkt + 5 * (uint32_t)(produkt >> 32);	vrátit se *Stát = X - 4;}

Mnoho dalších generátorů Lehmer má dobré vlastnosti. Následující modulo-2¹²⁸ Lehmerův generátor vyžaduje 128bitovou podporu kompilátoru a používá multiplikátor vypočítaný společností L'Ecuyer.^[17] Má období 2¹²⁶:

statický nepodepsaný __int128 Stát;/ * Stát musí být naočkován lichou hodnotou. * /prázdnota semínko(nepodepsaný __int128 semínko){	Stát = semínko << 1 | 1;}uint64_t další(prázdnota){	// GCC nemůže psát 128bitové literály, proto používáme výraz	konst nepodepsaný __int128 mult =		(nepodepsaný __int128)0x12e15e35b500f16e << 64 |		0x2e714eb2b37916a5;	Stát *= mult;	vrátit se Stát >> 64;}

Generátor počítá lichou 128bitovou hodnotu a vrací svých horních 64 bitů.

Tento generátor předává BigCrush z TestU01, ale selže test TMFn z PractRand. Tento test byl navržen tak, aby přesně zachytil vadu tohoto typu generátoru: protože modul je výkon 2, doba nejnižšího bitu na výstupu je pouze 2⁶², spíše než 2¹²⁶. Lineární kongruentní generátory s modulem síly 2 mají podobné chování.

Následující základní rutina vylepšuje rychlost výše uvedeného kódu pro celočíselná pracovní vytížení (pokud je povoleno optimalizovat konstantní deklaraci z smyčky výpočtu kompilátorem):

uint64_t další(prázdnota){	uint64_t výsledek = Stát >> 64;	// GCC nemůže psát 128bitové literály, proto používáme výraz	konst nepodepsaný __int128 mult =		(nepodepsaný __int128)0x12e15e35b500f16e << 64 |		0x2e714eb2b37916a5;	Stát *= mult;	vrátit se výsledek;}

Protože je však násobení odloženo, není vhodné pro hašování, protože první volání jednoduše vrátí horních 64 bitů počátečního stavu.

Reference

• Lehmer, D. H. (1949). "Matematické metody ve velkých výpočetních jednotkách". Sborník druhého sympozia o digitálních počítacích strojích velkého rozsahu. str.141 –146. PAN  0044899. (deníková verze: Annals of the Computation Laboratory of Harvard University, Sv. 26 (1951)).
• Steve Park, Generátory náhodných čísel

externí odkazy

• Připravuje jen méně než mocninu dvou může být užitečné pro výběr modulů. Část Prime Stránky.