View source for Multiply-with-carry pseudorandom number generator

{{Short description|Method for generating sequences of random integers}}
In [[computer science]], '''multiply-with-carry''' (MWC) is a method invented by [[George Marsaglia]]{{r|Marsaglia95}} for generating sequences of random integers based on an initial set from two to many thousands of randomly chosen seed values. The main advantages of the MWC method are that it invokes simple computer integer arithmetic and leads to very fast generation of sequences of random numbers with immense periods, ranging from around <math>2^{60}</math> to <math>2^{2000000}</math>.

As with all [[pseudorandom number generator]]s, the resulting sequences are functions of the supplied seed values.

==General theory==
An MWC generator is a special form of [[Lehmer random number generator]]  <math>x_n = bx_{n-1} \bmod p</math> which allows efficient implementation of a prime modulus <math>p</math> much larger than the machine word size.

Normal Lehmer generator implementations choose a modulus close to the machine word size.  An MWC generator instead maintains its state in base <math>b</math>, so multiplying by <math>b</math> is done implicitly by shifting one word. The base <math>b</math> is typically chosen to equal the computer's word size, as this makes arithmetic modulo <math>b</math> trivial.  This may vary from <math>b = 2^8</math> for a [[microcontroller]] to <math>b = 2^{64}</math>.  (This article uses <math>b = 2^{32}</math> for examples.)

The initial state ("seed") values are arbitrary, except that they must not be all zero, nor all at the maximum permitted values (<math>x_0 = b-1</math> and <math>c_0 = a-1</math> ).  (This is commonly done by choosing  <math>c_0</math> between 1 and <math>a-2</math>.). The MWC sequence is then a sequence of pairs <math>x_n, c_n</math> determined by

: <math>x_n=(ax_{n-1}+c_{n-1})\,\bmod\,b,\ c_n=\left\lfloor\frac{ax_{n-1}+c_{n-1}}{b}\right\rfloor</math>

This is called a lag-1 MWC sequence. Sometimes an odd base is preferred, in which case <math>b = 2^k-1</math> can be used, which is almost as simple to implement. A lag-<math>r</math> sequence is a generalization of the lag-1 sequence allowing longer periods{{r|Marsaglia03}}. The lag-<math>r</math> MWC sequence is then a sequence of pairs
<math>x_n, c_n</math> (for <math>n > r</math>) determined by

: <math>x_n=(ax_{n-r}+c_{n-1})\,\bmod\,b,\ c_n=\left\lfloor\frac{ax_{n-r}+c_{n-1}}{b}\right\rfloor</math>

and the MWC generator output is the sequence of <math>x</math>'s,
: <math>x_r, x_{r+1}, x_{r+2}, ...</math>

In this case, the initial state ("seed") values must not be all zero nor <math>x_0 = b-1</math> and <math>c_{r-1} = a-1</math>.

The MWC multiplier <math>a</math> and lag <math>r</math> determine the modulus <math>p = ab^r-1</math>. In practice, <math>a</math> is chosen so the modulus is prime and the sequence has long period. If the modulus is prime, the period of a lag-<math>r</math> MWC generator is the order of <math>b</math> in the multiplicative group of numbers modulo <math>p</math>.  While it is theoretically possible to choose a non-prime modulus, a prime modulus eliminates the possibility of the initial seed sharing a common divisor with the modulus, which would reduce the generator's period.

Because 2 is a [[quadratic residue]] of numbers of the form <math>8k \pm 1</math>, <math>b = 2^k</math> cannot be a primitive root of <math>p = ab^r-1</math>. Therefore, MWC generators with base <math>2^k</math> have their parameters chosen so their period is (''ab<sup>r</sup>''−1)/2.  This is one of the difficulties that use of ''b'' = 2<sup>''k''</sup>&nbsp;−&nbsp;1 overcomes.

The basic form of an MWC generator has parameters ''a'', ''b'' and ''r'', and ''r''+1 words of state.  The state consists of ''r'' residues modulo ''b''
: 0 ≤ ''x''<sub>0</sub>, ''x''<sub>1</sub>, ''x''<sub>2</sub>,..., ''x''<sub>''r''−1</sub> < b,
and a carry ''c''<sub>''r''−1</sub>&nbsp;&lt;&nbsp;''a''.

Although the theory of MWC generators permits ''a''&nbsp;&gt;&nbsp;''b'', ''a'' is almost always chosen smaller for convenience of implementation.

The state transformation function of an MWC generator is one step of [[Montgomery reduction]] modulo ''p''.  The state is a large integer with most significant word ''c''<sub>''n''−1</sub> and least significant word ''x''<sub>''n''−''r''</sub>.  Each step, ''x''<sub>''n''−''r''</sub>·(''ab<sup>r</sup>''−1) is added to this integer.  This is done in two parts:  −1·''x''<sub>''n''−''r''</sub> is added to ''x''<sub>''n''−''r''</sub>, resulting in a least significant word of zero.  And second, ''a''·''x''<sub>''n''−''r''</sub> is added to the carry.  This makes the integer one word longer, producing two new most significant words ''x<sub>n</sub>'' and ''c<sub>n</sub>''.

So far, this has simply added a multiple of ''p'' to the state, resulting in a different representative of the same [[residue class]] modulo ''p''.  But finally, the state is shifted down one word, dividing by ''b''.  This discards the least significant word of zero (which, in practice, is never computed at all) and effectively multiplies the state by ''b''<sup>−1</sup> (mod ''p'').

Thus, a multiply-with-carry generator ''is'' a Lehmer generator with modulus ''p'' and multiplier ''b''<sup>−1</sup> (mod ''p'').  This is the same as a generator with multiplier ''b'', but producing output in reverse order, which does not affect the quality of the resultant pseudorandom numbers.

Couture and L'Ecuyer{{r|Couture97}} have proved the surprising result that the lattice associated with a multiply-with-carry generator is very close to the lattice associated with the Lehmer generator it simulates. Thus, the mathematical techniques developed for Lehmer generators (such as the [[spectral test]]) can be applied to multiply-with-carry generators.

==Efficiency==
A [[linear congruential generator]] with base ''b'' = 2<sup>32</sup> is implemented as
:<math>x_{n+1}=(ax_n+c)\ \bmod\,2^{32},</math>
where ''c'' is a constant. 
If ''a'' ≡ 1 (mod 4) and ''c'' is odd, the resulting base-2<sup>32</sup> congruential sequence will have period 2<sup>32</sup>.<ref>{{Cite journal
 |last1=Hull |first1=T. E.  |last2=Dobell |first2=A. R.
 |date=July 1962
 |title=Random Number Generators
 |url=http://chagall.med.cornell.edu/BioinfoCourse/PDFs/Lecture4/random_number_generator.pdf
 |journal=SIAM Review |volume=4 |issue=3 |pages=230–254 
 |doi=10.1137/1004061
|bibcode=1962SIAMR...4..230H |hdl=1828/3142 |hdl-access=free}}</ref>

This can be computed using only the low 32 bits of the product of ''a'' and the current ''x''.
However, many [[microprocessor]]s can compute a full 64-bit product in almost the same time as the low 32 bits.  Indeed, many compute the 64-bit product and then ignore the high half.

A lag-1 multiply-with-carry generator allows us to make the period nearly 2<sup>63</sup> by using almost the same computer operations, except that the top half of the 64-bit product is not ignored after the product is computed.
Instead, a 64-bit sum is computed, and the top half is used as a ''new'' carry value ''c'' rather than the fixed additive constant of the standard congruential sequence:
Compute ''ax''+''c'' in 64 bits, then use the top half as the new ''c'', and bottom half as the new ''x''.

==Choice of multiplier==
With multiplier ''a'' specified, 
each pair of input values ''x'', ''c'' is converted to a new pair,

:<math>x\leftarrow (ax+c)\,\bmod\,2^{32},\ \ c\leftarrow \left\lfloor\frac{ax+c}{2^{32}}\right\rfloor.</math>

If ''x'' and ''c'' are not both zero, then the period of the resulting multiply-with-carry sequence will be the order of ''b'' = 2<sup>32</sup>
in the [[Multiplicative group of integers modulo n|multiplicative group of residues modulo
''ab<sup>r</sup>''&nbsp;−&nbsp;1]], that is, the smallest ''n'' such that
''b<sup>n</sup>'' ≡ 1 (mod ''ab<sup>r</sup>''&nbsp;−&nbsp;1).

If ''p'' = ''ab<sup>r</sup>''&nbsp;−&nbsp;1 is prime, then [[Fermat's little theorem]] guarantees that the order of any element must divide ''p''&nbsp;−&nbsp;1 = ''ab<sup>r</sup>''&nbsp;−&nbsp;2, so one way to ensure a large order is to choose ''a'' such that ''p'' is a "[[safe prime]]",
that is both ''p'' and (''p''&nbsp;−&nbsp;1)/2 = ''ab<sup>r</sup>''/2&nbsp;−&nbsp;1 are prime.  In such a case, for ''b''= 2<sup>32</sup> and ''r'' = 1, the period will be ''ab<sup>r</sup>''/2&nbsp;−&nbsp;1, approaching 2<sup>63</sup>, which in practice may be an acceptably large subset of the number of possible 32-bit pairs (''x'', ''c'').

More specifically, in such a case, the order of ''any'' element divides ''p''&nbsp;−&nbsp;1, and there are only four possible divisors: 1, 2, ''ab<sup>r</sup>''/2&nbsp;−&nbsp;1, or ''ab<sup>r</sup>''&nbsp;−&nbsp;2.  The first two apply only to the elements 1 and −1, and quadratic reciprocity arguments show that the fourth option cannot apply to ''b'', so only the third option remains.

Following are some maximal values of ''a'' for computer applications which satisfy the above safe prime condition, for lag-1 generators:
<!--Pari/GP script: maxsafeprime(bb,ab)=local(b=2^bb,am=2^ab);forstep(a=am-1,am>>1,-1,p=a*b-1; if(gcd(p,16294579238595022365)==1&&gcd(p>>1,16294579238595022365)==1&&gcd(p,7145393598349078859)==1&& gcd(p>>1,7145393598349078859)==1&&#factor(p,10^5)[,1]==1&&#factor(p>>1,10^5)[,1]==1&&ispseudoprime(p,1)&&isprime(p>>1)&&isprime(p), print(a-am" "a" "p>>1);break)) -->

{| class="wikitable" border="1"
|-
! Bits in ''a''
! ''b''
! Maximum ''a'' such that ''ab''−1 is a safe prime
! Period
|-
| 15
| 2<sup>16</sup>
| 2<sup>15</sup>−50 = 32,718
| 1,072,103,423
|-
| 16
| 2<sup>16</sup>
| 2<sup>16</sup>−352 = 65,184
| 2,135,949,311
|-
| 31
| 2<sup>32</sup>
| 2<sup>31</sup>−563 = 2,147,483,085
| 4,611,684,809,394,094,079
|-
| 32
| 2<sup>32</sup>
| 2<sup>32</sup>−178 = 4,294,967,118
| 9,223,371,654,602,686,463
|-
| 64
| 2<sup>64</sup>
| 2<sup>64</sup>−742 = 18,446,744,073,709,550,874
| 2<sup>63</sup>(2<sup>64</sup>−742)−1 ≈ {{scinote|{{#expr:2^63*(2^64-10408)}}|4}}
|-
| 128
| 2<sup>128</sup>
| 2<sup>128</sup>−10,408
| 2<sup>127</sup>(2<sup>128</sup>−10,408)−1 ≈ {{scinote|{{#expr:2^127*(2^128-10408)}}|4}}
|-
| 256
| 2<sup>256</sup>
| 2<sup>256</sup>−9166
| 2<sup>255</sup>(2<sup>256</sup>−9166)−1 ≈ {{scinote|{{#expr:2^255*(2^256-9166)}}|4}}
|-
| 512
| 2<sup>512</sup>
| 2<sup>512</sup>−150,736
| 2<sup>511</sup>(2<sup>512</sup>−150,736)−1 ≈ {{scinote|{{#expr:2^511*(2^512-150736)}}|4}}
|}
<!-- Some more, if anyone wants:
? maxsafeprime(64,32) \\ lag-2
-4222 4294963074 39614042316055429195908513791
? maxsafeprime(128,32) \\ lag-4
-5521 4294961775 730749879315977573851214025712537694431818547199
? maxsafeprime(128,64) \\ lag-2
-4057 18446744073709547559 3138550867693339691655113412480160072595786546890200317951
? maxsafeprime(256,32) \\ lag-8
-22807 4294944489 248660297769803703342456611381682873477670870634403724088059964492020514268892838756351
? maxsafeprime(256,64) \\ lag-4
-21466 18446744073709530150 1067993517960453798401017068970050576113970162931011395826775387590936027392071118079860618035199
? maxsafeprime(512,32) \\ lag-16
-67063 4294900233 28792598701164853979792802082771056557436637982777253258006290527115210359966947958901407759016872502237208239751280948631518700106842009223244338317635401163997183
? maxsafeprime(512,64) \\ lag-8
-26086 18446744073709525530 123665200736552092152212430268529625273557622951754879243663130059861786904854023122189956281348304067952364722086662308590156190815814319471879561331834772624219399719485439
? maxsafeprime(1024,32) \\ lag-32
-342523 4294624773 386020873561586592074312812522001851875274492595466457251272730611699635012941082972544756647960056099096146138538415983914543911174668233250385795466842403122622602773225563630137131067198828737446415628624118888350116673947397380305414998889062040317435675434112752838821870310259311800378970489743273380296192425983
? maxsafeprime(1024,64) \\ lag-16
-152947 18446744073709398669 1658079259093474837955046990989813605647917294419555600143486010512540342085706465869156251304268794682832026731777412049076520831664483243035188343024949666509083731849181029064715081772539350969645500928870083606320829031526419342763296573142748690249657038997521828687436494817689824137958118622243245482394468296933351882751
? maxsafeprime(2048,32) \\ lag-64
-1749331 4293217965 69371975520183243822837632971278774856132795479963488742706318105552227461308599472167060693836765045724389558824030269595037151601383880507414964542078695818449372423504269991748368637137720087498466622142976226949894222824351525908961377017655342894351812985614998321172518146432497632380597797298854096266029258649761736536128160679625525193903755023032720248557704312209315766210756914099716950140053001234861530530470624087454335210719938647927313602780385238554702492456709577138556090790697041120625025891146398536158939286917009857603714093126373277598762864327152198316564372746624350294958143976955590611549311467519
time = 3min, 40,326 ms.
? maxsafeprime(2048,64) \\ lag-32
-332116 18446744073709219500 298071770112990595075757018580978828533608570182393883969695435903079275381608290209678233651762705213235666843932624307449818244274650520216029820130416772033982119613360882096743826673135638907502306566491761156503684445718874265540350778839560442913813248433348212305532331813493525835522564261912018620160726293562134908776979459110820699918337349980978809892711019892989688428080233423943299299553034323044182299930806193772627242438760165350131715351846812909177033297376337943631962139503149756400280910165255374105421571036754611380470994838686293814993994238634318615715024647997514398437529479730014419868867795105397866495999
time = 2min, 47,322 ms.
? maxsafeprime(4096,32) \\ lag-128
-3674707 4291292589 [redacted, 2.240889e1242]
time = 1h, 30min, 1,146 ms.
? maxsafeprime(4096,64) \\lag-64
-5098417 18446744073704453199 [redacted, 9.632787e1251]
time = 2h, 18min, 8,085 ms.
maxsafeprime(8192,32) \\ lag-256
-31313656 4263653640 [redacted, 2.325286e2475]
time = 59h, 20min, 50,916 ms.
-->

While a safe prime ensures that almost ''any'' element of the group has a large order, the period of the generator is specifically the order of ''b''.  For small moduli, more computationally expensive methods<!--The expensive part is finding the factors of ab^r−1, required for the test--> can be used to find multipliers ''a'' where the period is ''ab''/2&nbsp;−&nbsp;1.
The following are again maximum values of ''a'' of various sizes.
<!-- Pari/GP script: maxord(bbits,abits,r)=local(b=2^bbits,am=2^abits);forstep(a=am-1,am>>1,-1,p=a*b^r-1; if(isprime(p),o=znorder(Mod(b,p));if(p==o*2+1,print(a-am" "a" "o);break))) -->

{| class="wikitable" border="1"
|-
! Bits in ''a''
! ''b<sup>r</sup>''
! Maximum ''a'' such that<br/>''b'' has order ''ab<sup>r</sup>''/2&minus;1
! Period
|-
| 8
| 2<sup>8</sup>
| 2<sup>8</sup>−7 = 249
| 31,871
|-
| 8
| (2<sup>8</sup>)<sup>2</sup> = 2<sup>16</sup>
| 2<sup>8</sup>−32 = 224
| 7,340,031
|-
| 15
| 2<sup>16</sup>
| 2<sup>15</sup>−29 = 32,739
| 1,072,791,551
|-
| 16
| 2<sup>16</sup>
| 2<sup>16</sup>−22 = 65,514
| 2,146,762,751
|-
| 8
| (2<sup>8</sup>)<sup>4</sup> = 2<sup>32</sup>
| 2<sup>8</sup>−64 = 192
| 412,316,860,415
|-
| 15
| (2<sup>16</sup>)<sup>2</sup> = 2<sup>32</sup>
| 2<sup>15</sup>−26 = 32,742
| 70,312,909,602,815
|-
| 16
| (2<sup>16</sup>)<sup>2</sup> = 2<sup>32</sup>
| 2<sup>16</sup>−2 = 65,534
| 140,733,193,388,031
|-
| 31
| 2<sup>32</sup>
| 2<sup>31</sup>−68 = 2,147,483,580
| 4,611,685,872,398,499,839
|-
| 32
| 2<sup>32</sup>
| 2<sup>32</sup>−76 = 4,294,967,220
| 9,223,371,873,646,018,559
|-
| 8
| (2<sup>8</sup>)<sup>8</sup> = 2<sup>64</sup>
| 2<sup>8</sup>−41 = 215
| 2<sup>63</sup>(2<sup>8</sup>−41)−1 ≈ {{scinote|{{#expr:2^63*(2^8-41)}}|4}}
|-
| 15
| (2<sup>16</sup>)<sup>4</sup> = 2<sup>64</sup>
| 2<sup>15</sup>−50 = 32,718
| 2<sup>63</sup>(2<sup>15</sup>−50)−1 ≈ {{scinote|{{#expr:2^63*(2^15-50)}}|4}}
|-
| 16
| (2<sup>16</sup>)<sup>4</sup> = 2<sup>64</sup>
| 2<sup>16</sup>−56 = 65,480
| 2<sup>63</sup>(2<sup>16</sup>−56)−1 ≈ {{scinote|{{#expr:2^63*(2^16-56)}}|4}}
|-
| 31
| (2<sup>32</sup>)<sup>2</sup> = 2<sup>64</sup>
| 2<sup>31</sup>−38 = 2,147,483,610
| 2<sup>63</sup>(2<sup>31</sup>−38)−1 ≈ {{scinote|{{#expr:2^63*(2^31-38)}}|4}}
|-
| 32
| (2<sup>32</sup>)<sup>2</sup> = 2<sup>64</sup>
| 2<sup>32</sup>−43 = 4,294,967,253
| 2<sup>63</sup>(2<sup>32</sup>−43)−1 ≈ {{scinote|{{#expr:2^63*(2^32-43)}}|4}}
|-
| 63
| 2<sup>64</sup>
| 2<sup>63</sup>−140 = 9,223,372,036,854,775,668
| 2<sup>63</sup>(2<sup>63</sup>−140)−1 ≈ {{scinote|{{#expr:2^63*(2^63-140)}}|4}}
|-
| 64
| 2<sup>64</sup>
| 2<sup>64</sup>−116 = 18,446,744,073,709,551,500
| 2<sup>63</sup>(2<sup>64</sup>−116)−1 ≈ {{scinote|{{#expr:2^63*(2^64-116)}}|4}}
|}
<!-- Some more, if anyone wants:
? maxord(8,8,3)
-1 255 2139095039
? maxord(8,8,5)
-1 255 140187732541439
? maxord(8,8,6)
-22 234 32932572275146751
? maxord(16,15,3)
-53 32715 4604226931544555519
? maxord(16,16,3)
-31 65505 9219009174715760639
? maxord(8,8,12)
-4 252 9982748476797306536786537742335
? maxord(16,15,6)
-51 32717 1296053896489593166523988711571455
? maxord(16,16,6)
-83 65453 2592860460523071844255116090671103
? maxord(32,31,3)
-14 2147483634 85070591175637478265993288703134400511
? maxord(32,32,3)
-22 4294967274 170141182588959444074779590186900652031
? maxord(8,8,16)
-26 230 39132472195907923298288079854653344317439
? maxord(16,15,8)
-11 32757 5573314746614590623834881007821215651332095
? maxord(16,16,8)
-121 65415 11129785516066594793728324972574558776197119
? maxord(32,31,4)
-11 2147483637 365375407461172711485759659130510414039257972735
? maxord(32,32,4)
-76 4294967220 730750805734721516106180804749906427420774236159
? maxord(64,63,2)
-3 9223372036854775805 1569275433846670190448523805420508908830526949968356311039
? maxord(64,64,2)
-142 18446744073709551474 3138550867693340357757846660217202302151580594576474243071
? maxord(128,127,1)
-213 170141183460469231731687303715884105515 28948022309329048855892746252171976927077424089330194651015000310494967889919
? maxord(128,128,1)
-173 340282366920938463463374607431768211283 57896044618658097711785492504343953897200567594159104930146888461108614529023
? maxord(8,8,32)
-107 149 8626510648180056559056038383147249135068613857590222020939590008589528158175231
time = 3,831 ms.
? maxord(16,15,16)
-306 32462 1879421400410879167919980657676013432366425121108011994924436046032438007185801215
? maxord(16,16,16)
-86 65450 3789296120291172495236360484409311784498260248183087458191249436658957167466905599
? maxord(32,31,8)
-38 2147483610 124330806902396965029837848988856494860043788487207220712945277473595778102783880724479
? maxord(32,32,8)
-107 4294967189 248661612010016546881274668912362722085425005248729679108461812335599004082914252029951
? maxord(64,63,4)
-239 9223372036854775569 533996758980227506761600762683102675533943422050958204020997234411507528729978848560121529761791
? maxord(64,64,4)
-86 18446744073709551530 1067993517960455036218451015880179654087799905804746347282911463288565609193644162671216468951039
-->

==MWC generators as repeating decimals==
The output of a multiply-with-carry generator is equivalent to the [[radix]]-''b'' expansion of a fraction with denominator ''p'' = ''ab<sup>r</sup>''&nbsp;−&nbsp;1.  Here is an example for the simple case of ''b'' = 10 and ''r'' = 1, so the result is a [[repeating decimal]].

Starting with <math display=inline>c_0=1,x_0=0</math>, the MWC sequence

:<math>x_n=(7x_{n-1}+c_{n-1})\,\bmod\,10,\ c_n=\left\lfloor\frac{7x_{n-1}+c_{n-1}}{10}\right\rfloor,</math>

produces this sequence of states:
:10,01,07,49,67,55,40,04,28,58,61,13,22,16,43,25,37,52,19,64,34,31, 10,01,07,...

with period 22.  Consider just the sequence of ''x<sub>i</sub>'':
:0,1,7,9,7,5,0,4,8,8,1,3,2,6,3,5,7,2,9,4,4,1, 0,1,7,9,7,5,0,...

Notice that if those repeated segments of ''x'' values are put ''in reverse order'':
:<math>1449275\cdots 9710\,1449275\cdots 9710\,144\cdots</math>
we get the expansion
''j''/(''ab''−1) with ''a''=7, ''b''=10, ''j=10'':

:<math>\frac{10}{69}=0.1449275362318840579710\,14492753623\ldots</math>

This is true in general: The sequence of ''x''s produced by a lag-''r'' MWC generator:

:<math> x_n=(ax_{n-r}+c_{n-1})\bmod\,b\,,\ \ c_n=\left\lfloor\frac{ax_{n-r}+c_{n-1}}{b}\right\rfloor,</math>

when put in reverse order, will be the radix-''b'' expansion of a fraction ''j''/(''ab<sup>r</sup>''&nbsp;&minus;&nbsp;1) for some 0 < ''j'' < ''ab<sup>r</sup>''.

==Equivalence to linear congruential generator==
Continuing the above example, if we start with <math display=inline>x_0=34</math>, and generate the ordinary congruential sequence
:<math>x_n=7x_{n-1}\,\bmod\,69</math>,
we get the period 22 sequence
:31,10,1,7,49,67,55,40,4,28,58,61,13,22,16,43,25,37,52,19,64,34, 31,10,1,7,...
and that sequence, reduced mod 10, is
:1,0,1,7,9,7,5,0,4,8,8,1,3,2,6,3,5,7,2,9,4,4, 1,0,1,7,9,7,5,0,...
the same sequence of ''x'''s resulting from the MWC sequence.

This is true in general, (but apparently only for lag-1 MWC sequences{{citation needed|date=April 2018}}):

Given initial values <math>x_0,c_0</math>, the sequence <math>x_1,x_2,\ldots</math> resulting from the lag-1 MWC sequence

:<math>x_n=(ax_{n-1}+c_{n-1})\,\bmod b\,,\ \ c_n=\left\lfloor\frac{ax_{n-1}+c_{n-1}}{b}\right\rfloor</math>

is exactly the [[Lehmer random number generator]] output sequence ''y<sub>n</sub>'' = ''ay''<sub>''n''&nbsp;&minus;&nbsp;1</sub> mod (''ab''&nbsp;&minus;&nbsp;1),
reduced modulo ''b''.

Choosing a different initial value ''y''<sub>0</sub> merely rotates the cycle of ''x'''s.

==Complementary-multiply-with-carry generators==
<!--CMWC originates with Couture & L'Ecuyer, but was adopted (with b=2^k−1) by Marsaglia -->
Establishing the period of a lag-''r'' MWC generator usually entails choosing multiplier ''a'' so that ''p'' = ''ab<sup>r</sup>''&nbsp;−&nbsp;1 is prime.  Then ''p''&nbsp;−&nbsp;1 will have to be factored in order to find the order of ''b'' mod ''p''.  If ''p'' is a [[safe prime]], then this is simple, and the order of ''b'' will be either ''p''&nbsp;−&nbsp;1 or (''p''&nbsp;−&nbsp;1)/2.  In other cases, ''p''&nbsp;&minus;&nbsp;1 may be difficult to factor.

However, the algorithm also permits a ''negative'' multiplier.  This leads to a slight modification of the MWC procedure, and produces a modulus ''p'' = {{abs|−''ab<sup>r</sup>''&nbsp;−&nbsp;1}} = ''ab<sup>r</sup>''&nbsp;+&nbsp;1.  This makes ''p''&nbsp;−&nbsp;1 = ''ab<sup>r</sup>'' easy to factor, making it practical to establish the period of very large generators.

The modified procedure is called complementary-multiply-with-carry (CMWC), and the setup is the same as that for lag-''r'' MWC: multiplier ''a'', base ''b'', and ''r''&nbsp;+&nbsp;1 seeds,

: ''x''<sub>0</sub>, ''x''<sub>1</sub>, ''x''<sub>2</sub>, ..., ''x''<sub>''r''−1</sub>, and ''c''<sub>''r''−1</sub>.

The modification is to the generation of a new pair (''x'', ''c'').  Rearranging the computation to avoid negative numbers, the new ''x'' value is complemented by subtracting it from ''b''&nbsp;−&nbsp;1:
:<math>x_n=(b-1)-(ax_{n-r}+c_{n-1})\,\bmod\,b,\ c_n=\left\lfloor\frac{ax_{n-r}+c_{n-1}}{b}\right\rfloor.</math>

The resulting sequence of ''x''{{'}}s produced by the CMWC RNG will have period the order of ''b'' in the multiplicative group of residues modulo ''ab<sup>r</sup>''+1, and the output ''x''{{'}}s, in reverse order, will form the base ''b'' expansion of ''j''/(''ab<sup>r</sup>''+1) for some 0&nbsp;<&nbsp;''j''&nbsp;<&nbsp;''ab<sup>r</sup>''.

Use of lag-''r'' CMWC makes it much easier to find periods for ''r''{{'}}s as large as 512, 1024, 2048, etc.
(Making ''r'' a power of 2 makes it slightly simpler to access elements in the array containing the ''r'' most recent ''x''{{'}}s.)

Another advantage of this modified procedure is that the period is a multiple of ''b'', so the output is exactly [[equidistributed]] mod ''b''.{{r|Couture97}}  (The ordinary MWC, over its full period, produces each possible output an equal number of times ''except'' that zero is produced one time less, a bias which is negligible if the period is long enough.)

One ''disadvantage'' of the CMWC construction is that, with a power-of-two base, the maximum achievable period is less than for a similar-sized MWC generator; you lose several bits.  Thus, an MWC generator is usually preferable for small lags.  This can remedied by using ''b'' = 2<sup>''k''</sup>−1, or choosing a lag one word longer to compensate.

Some examples:
With ''b'' = 2<sup>32</sup>, and ''a'' = 109111 or 108798 or 108517, the period of the lag-1024 CMWC
: <math>x_n=(b-1)-(ax_{n-1024}+c_{n-1})\,\bmod\,b,\ c_n=\left\lfloor\frac{ax_{n-1024}+c_{n-1}}{b}\right\rfloor.</math>
will be ''a''·2<sup>32762</sup> = ''ab''<sup>1024</sup>/64, about 10<sup>9867</sup>.

With ''b'' = 2<sup>32</sup> and ''a'' = 3636507990, ''p'' = ''ab''<sup>1359</sup>&nbsp;−&nbsp;1 is a safe prime, so the MWC sequence based on that ''a'' has period 3636507990·2<sup>43487</sup> ≈ 10<sup>13101</sup>.

With ''b'' = 2<sup>32</sup>, a CMWC RNG with near record period may be based on the prime ''p'' = 15455296''b''<sup>42658</sup>&nbsp;+&nbsp;1.  The order of ''b'' for that prime is 241489·2<sup>1365056</sup> ≈ 10<sup>410928</sup>.

==More general moduli==
The MWC modulus of ''ab<sup>r</sup>''−1 is chosen to make computation particularly simple, but brings with it some disadvantages, notably that the period is at most half the modulus.  There are several ways to generalize this, at the cost of more multiplications per iteration.

First, it is possible to add additional terms to the product, producing a modulus of the form ''a<sub>r</sub>b<sup>r</sup>''+''a<sub>s</sub>b<sup>s</sup>''−1.  This requires computing ''c<sub>n</sub>b'' + ''x<sub>n</sub>'' =  ''a<sub>r</sub>x''<sub>''n''−''r''</sub> + ''a<sub>s</sub>x''<sub>''n''−''s''</sub>.  (The carry is limited to one word if ''a<sub>r</sub>'' + ''a<sub>s</sub>'' ≤ ''b''.)

However, this does not fix the period issue, which depends on the low bits of the modulus.  Fortunately, the Montgomery reduction algorithm permits other moduli, as long as they are [[relatively prime]] to the base ''b'', and this can be applied to permit a modulus of the form ''a<sub>r</sub>b<sup>r</sup>''−''a''<sub>0</sub>, for a wide range of values ''a''<sub>0</sub>. Goresky and Klapper{{r|Goresky03}} developed the theory of these generalized multiply-with-carry generators, proving, in particular, that choosing a negative ''a''<sub>0</sub> and ''a<sub>r</sub>''–''a''<sub>0</sub> &lt; ''b'' the carry value is always smaller than ''b'', making the implementation efficient. The more general form of the modulus improves also the quality of the generator, albeit one cannot always get full period.

To implement a Goresky-Klapper generator one precomputes ''a''{{su|b=0|p=−1}}&nbsp;(mod&nbsp;''b''), and changes the iteration as follows:{{refn|Note that the paper by Goresky and Klapper{{r|Goresky03}} contains a mistake in equation (4): the last equality is not true—one cannot eliminate the second summand from the computation of the carry.}}
:<math>t = c_{n-1} + \sum_1^r a_i x_{n-i},
\ x_n = a_0^{-1}t \bmod b, \ c_n = \left\lfloor\frac{t - a_0 x_n}{b}\right\rfloor.</math>

In the common case that ''b'' = 2<sup>''k''</sup>, ''a''<sub>0</sub> must be odd for the inverse to exist.

==Implementation==
The following is an implementation of the CMWC algorithm in the [[C (programming language)|C programming language]]. Also, included in the program is a sample initialization function. In this implementation the base is 2<sup>32</sup>&minus;1 and lag ''r''=4096. The period of the resulting generator is about <math>2^{131104}</math>.

<syntaxhighlight lang="c">
// C99 Complementary Multiply With Carry generator
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>

// How many bits in rand()?
// https://stackoverflow.com/a/27593398
#define LOG_1(n) (((n) >= 2) ? 1 : 0)
#define LOG_2(n) (((n) >= 1<<2) ? (2 + LOG_1((n)>>2)) : LOG_1(n))
#define LOG_4(n) (((n) >= 1<<4) ? (4 + LOG_2((n)>>4)) : LOG_2(n))
#define LOG_8(n) (((n) >= 1<<8) ? (8 + LOG_4((n)>>8)) : LOG_4(n))
#define LOG(n)   (((n) >= 1<<16) ? (16 + LOG_8((n)>>16)) : LOG_8(n))
#define BITS_TO_REPRESENT(n) (LOG(n) + !!((n) & ((n) - 1)))
#if ((RAND_MAX | (RAND_MAX >> 1)) != RAND_MAX) 
#error "expected a RAND_MAX that is 2^n - 1!"
#endif
#define RAND_BITS BITS_TO_REPRESENT(RAND_MAX)

// CMWC working parts
#define CMWC_CYCLE 4096         // as Marsaglia recommends
#define CMWC_C_MAX 809430660    // as Marsaglia recommends
struct cmwc_state {
	uint32_t Q[CMWC_CYCLE];
	uint32_t c;	// must be limited with CMWC_C_MAX
	unsigned i;
};

// Collect 32 bits of rand(). You are encouraged to use a better source instead.
uint32_t rand32(void)
{
	uint32_t result = rand();
    for (int bits = RAND_BITS; bits < 32; bits += RAND_BITS)
        result = result << RAND_BITS | rand();
	return result;
}

// Init the state with seed
void initCMWC(struct cmwc_state *state, unsigned int seed)
{
	srand(seed);        
	for (int i = 0; i < CMWC_CYCLE; i++)
		state->Q[i] = rand32();
	do
		state->c = rand32();
	while (state->c >= CMWC_C_MAX);
	state->i = CMWC_CYCLE - 1;
}

// CMWC engine
uint32_t randCMWC(struct cmwc_state *state)  //EDITED parameter *state was missing
{
	uint64_t const a = 18782;	// as Marsaglia recommends
	uint32_t const m = 0xfffffffe;	// as Marsaglia recommends
	uint64_t t;
	uint32_t x;

	state->i = (state->i + 1) & (CMWC_CYCLE - 1);
	t = a * state->Q[state->i] + state->c;
	/* Let c = t / 0xffffffff, x = t mod 0xffffffff */
	state->c = t >> 32;
	x = t + state->c;
	if (x < state->c) {
		x++;
		state->c++;
	}
	return state->Q[state->i] = m - x;
}

int main()
{
	struct cmwc_state cmwc;
	unsigned int seed = time(NULL);

	initCMWC(&cmwc, seed);
	printf("Random CMWC: %u\n", randCMWC(&cmwc));
}
</syntaxhighlight>

The following are implementations of small-state MWC generators with 64-bit output using 128-bit multiplications.
<syntaxhighlight lang="c">
// C99 + __uint128_t MWC, 128 bits of state, period approx. 2^127

/* The state must be neither all zero, nor x = 2^64 - 1, c = MWC_A1 -
   1. The condition 0 < c < MWC_A1 - 1 is thus sufficient. */

uint64_t x, c = 1;

#define MWC_A1 0xff3a275c007b8ee6

uint64_t inline next() {
    const __uint128_t t = MWC_A1 * (__uint128_t)x + c;
    c = t >> 64;
    return x = t;
}</syntaxhighlight>

<syntaxhighlight lang="c">
// C99 + __uint128_t MWC, 256 bits of state, period approx. 2^255

/* The state must be neither all zero, nor x = y = z = 2^64 - 1, c =
   MWC_A3 - 1. The condition 0 < c < MWC_A3 - 1 is thus sufficient. */

uint64_t x, y, z, c = 1;

#define MWC_A3 0xff377e26f82da74a

uint64_t inline next() {
    const __uint128_t t = MWC_A3 * (__uint128_t)x + c;
    x = y;
    y = z;
    c = t >> 64;
    return z = t;
}
</syntaxhighlight>

The following are implementations of small-state Goresky-Klapper's generalized MWC generators with 64-bit output using 128-bit multiplications.
<syntaxhighlight lang="c">
// C99 + __uint128_t Goresky-Klapper GMWC, 128 bits of state, period approx. 2^127

/* The state must be neither all zero, nor x = 2^64 - 1, c = GMWC_A1 +
   GMWC_MINUS_A0. The condition 0 < c < GMWC_A1 + GMWC_MINUS_A0 is thus
   sufficient. */

uint64_t x = 0, c = 1;

#define GMWC_MINUSA0 0x7d084a4d80885f
#define GMWC_A0INV 0x9b1eea3792a42c61
#define GMWC_A1 0xff002aae7d81a646

uint64_t inline next() {
    const __uint128_t t = GMWC_A1 * (__uint128_t)x + c;
    x = GMWC_A0INV * (uint64_t)t;
    c = (t + GMWC_MINUSA0 * (__uint128_t)x) >> 64;
    return x;
}
</syntaxhighlight>

<syntaxhighlight lang="c">
// C99 + __uint128_t Goresky-Klapper GMWC, 256 bits of state, period approx. 2^255

/* The state must be neither all zero, nor x = y = z = 2^64 - 1, c =
   GMWC_A3 + GMWC_MINUS_A0. The condition 0 < c < GMWC_A3 + GMWC_MINUS_A0
   is thus sufficient. */

uint64_t x, y, z, c = 1; /* The state can be seeded with any set of values, not all zeroes. */

#define GMWC_MINUSA0 0x54c3da46afb70f
#define GMWC_A0INV 0xbbf397e9a69da811
#define GMWC_A3 0xff963a86efd088a2

uint64_t inline next() {
    const __uint128_t t = GMWC_A3 * (__uint128_t)x + c;
    x = y;
    y = z;
    z = GMWC_A0INV * (uint64_t)t;
    c = (t + GMWC_MINUSA0 * (__uint128_t)z) >> 64;
    return z;
}
</syntaxhighlight>

==Usage==
Because of simplicity and speed, CMWC is known to be used in game development, particularly in modern [[roguelike]] games. It is informally known as the Mother of All PRNGs, a name originally coined by Marsaglia himself.<ref>{{cite web |title=Marsaglia's Mother of All Random Number Generators |url=http://home.sandiego.edu/~pruski/MotherOfAll.html |website=home.sandiego.edu}}</ref> In libtcod, CMWC4096 replaced MT19937 as the default PRNG.<ref>{{Cite web|url=http://www.roguebasin.com/index.php?title=Random_number_generator#Complementary_multiply-with-carry|title=Random number generator - RogueBasin|website=www.roguebasin.com|access-date=2016-11-30}}</ref>

==See also==
* [[Mersenne Twister]]
* [[List of random number generators]]

==References==
{{Reflist|refs=
<ref name=Marsaglia95>{{cite web
 |first1=George |last1=Marsaglia |authorlink1=George Marsaglia
 |first2=Arif |last2=Zaman
 |date=August 1995
 |title=The Marsaglia Random Number CDROM including the Diehard Battery of Tests of Randomness
 |url=http://ftpmirror.your.org/pub/misc/diehard/
 }}</ref>
<ref name=Couture97>{{cite journal
  | last1 =Couture
  | first1 =Raymond
  | last2 =L'Ecuyer
  | first2 =Pierre
  | title =Distribution properties of Multiply-with-carry random number generators
  | journal =Mathematics of Computation
  | volume =66
  | issue =218
  | pages =591–607
  | date =April 1997
  | doi =10.1090/S0025-5718-97-00827-2
  | citeseerx =10.1.1.154.331
  | url =https://www.ams.org/mcom/1997-66-218/S0025-5718-97-00827-2/S0025-5718-97-00827-2.pdf
  | quote =We shall see that, for the complementary MWC, each bit of the output value is fair, that is, the two binary digits will appear equally often in a full period, a property not shared by MWC generators.
 | bibcode =1997MaCom..66..591C
 }}</ref>
<ref name=Goresky03>{{cite journal
 |title=Efficient multiply-with-carry random number generators with maximal period
 |first1=Mark |last1=Goresky  |authorlink1=Mark Goresky |first2=Andrew |last2=Klapper
 |journal=ACM Transactions on Modeling and Computer Simulation
 |volume=13 |issue=4 |date=October 2003 |pages=310–321
 |doi=10.1145/945511.945514
 |citeseerx=10.1.1.4.9190
 |s2cid=13334372 |url=https://www.math.ias.edu/~goresky/pdf/p1-goresky.pdf  <!--Backup URL http://cs.engr.uky.edu/~klapper/pdf/MWC.pdf -->
 }}</ref>
<1--ref name=Hull62>{{Cite journal
 |last1=Hull |first1=T. E.  |last2=Dobell |first2=A. R.
 |date=July 1962
 |title=Random Number Generators
 |url=http://chagall.med.cornell.edu/BioinfoCourse/PDFs/Lecture4/random_number_generator.pdf
 |journal=SIAM Review |volume=4 |issue=3 |pages=230–254 
 |doi=10.1137/1004061
 |bibcode=1962SIAMR...4..230H |hdl=1828/3142 |hdl-access=free}}</ref-->
<ref name=Marsaglia03>{{cite journal
 |first=George¨
 |last=Marsaglia
 |authorlink=George Marsaglia
 |title=Random number generators
 |journal=Journal of Modern Applied Statistical Methods
 |volume=2
 |issue=1
 |pages=2–13
 |date=May 2003
 |url=http://digitalcommons.wayne.edu/cgi/viewcontent.cgi?article=1725&context=jmasm
 |doi=10.22237/jmasm/1051747320
 |doi-access=free
}}</ref>
}}
*{{cite journal |last=Marsaglia |first=George |authorlink=George Marsaglia |title=Xorshift RNGs |journal=Journal of Statistical Software |date=4 July 2003 |volume=8 |issue=14 |pages=1–6 |doi=10.18637/jss.v008.i14 |doi-access=free}} 
*{{cite journal |first=George |last=Marsaglia |authorlink=George Marsaglia |title=On the randomness of Pi and other decimal expansions |journal=Interstat |date=October 2005 |url=http://interstat.statjournals.net/YEAR/2005/abstracts/0510005.php |citeseerx=10.1.1.694.4783}}
*{{Cite book | last1=Press | first1=William H. | authorlink1=William H. Press | last2=Teukolsky | first2=Saul A. | authorlink2=Saul Teukolsky | last3=Vetterling | first3=William T. | last4=Flannery | first4=Brian P. |authorlink4=Brian P. Flannery | year=2007 | title=Numerical Recipes: The Art of Scientific Computing | edition=3rd | publisher=Cambridge University Press |  location=New York | isbn=978-0-521-88068-8 | chapter=Section 7.1.2.B Multiply With Carry (MWC) | chapter-url=http://apps.nrbook.com/empanel/index.html#pg=347}}

{{Use dmy dates|date=September 2019}}

{{DEFAULTSORT:Multiply-With-Carry}}
[[Category:Pseudorandom number generators]]