Sign in Join
Wiki copytemporary

Marsaglia's pseudo random number generators: what changed

Revision 2 to 3. Removed lines have a minus sign; added lines have a plus sign.

OldNewLine
11![](/wiki-images/8/84/Fpc_source_logo.png)
22
33## Marsaglia's pseudo random number generators
44
55Although Free Pascal has a reasonably good pseudo random number generator (PRNG), a [[A simple implementation of the Mersenne twister|Mersenne Twister]], it is rather slow.
66
77In 1999 (and 2003), Professor George Marsaglia described a set of PRNGs with good properties that are much faster and often just as good. He was kind enough to provide sourcecode - in C - and a description of each. The XorShift PRNG from the 2003 paper is added as well. Some of these PRNGs, either stand-alone or in combination, have over time been proven to be equal or superior to either Delphi's LNG or even Free Pascal's Mersenne twister. Notably SWR, XorShift and combined e.g. SWB+KISS which has a period of 2^7700 and passes all tests for randomness.
88
99Here we present a class that encapsulates these PRNGs for 32 bit random values with a short explanation. The class contains the following PRNGs with Marsaglia's comments from his usenet post. It is extremely fast.
1010
1111Marsaglia describes them as follows:
1212
1313### Multiply with carry
1414
15-`     The MWC generator concatenates two 16-bit multiply-`\
16-`     with-carry generators, x(n)=36969x(n-1)+carry,`\
17-`     y(n)=18000y(n-1)+carry mod 2^16, has period about`\
18-`     2^60 and seems to pass all tests of randomness. A`\
19-`     favorite stand-alone generator---faster than KISS,`\
20-`     which contains it.  `
15+```text
16+ The MWC generator concatenates two 16-bit multiply-
17+ with-carry generators, x(n)=36969x(n-1)+carry,
18+ y(n)=18000y(n-1)+carry mod 2^16, has period about
19+ 2^60 and seems to pass all tests of randomness. A
20+ favorite stand-alone generator---faster than KISS,
21+ which contains it.
22+```
2123
2224### fibonaci
2325
24-`     FIB is the classical Fibonacci sequence`\
25-`     x(n)=x(n-1)+x(n-2),but taken modulo 2^32.`\
26-`     Its period is 3*2^31 if one of its two seeds is odd`\
27-`     and not 1 mod 8. It has little worth as a RNG by`\
28-`     itself, but provides a simple and fast component for`\
29-`     use in combination generators.`\
30-`     The classical Fibonacci sequence mod 2^32 from FIB`\
31-`     fails several tests. It is not suitable for use by`\
32-`     itself, but is quite suitable for combining with`\
33-`     other generators.`
26+```text
27+ FIB is the classical Fibonacci sequence
28+ x(n)=x(n-1)+x(n-2),but taken modulo 2^32.
29+ Its period is 3*2^31 if one of its two seeds is odd
30+ and not 1 mod 8. It has little worth as a RNG by
31+ itself, but provides a simple and fast component for
32+ use in combination generators.
33+ The classical Fibonacci sequence mod 2^32 from FIB
34+ fails several tests. It is not suitable for use by
35+ itself, but is quite suitable for combining with
36+ other generators.
37+```
3438
3539### Congruential generator
3640
37-`     CONG is a congruential generator with the widely used 69069`\
38-`     multiplier: x(n)=69069x(n-1)+1234567. It has period`\
39-`     2^32. The leading half of its 32 bits seem to pass`\
40-`     tests, but bits in the last half are too regular.`\
41-`     The last half of the bits of CONG are too regular,`\
42-`     and it fails tests for which those bits play a`\
43-`     significant role. CONG+FIB will also have too much`\
44-`     regularity in trailing bits, as each does. But keep`\
45-`     in mind that it is a rare application for which`\
46-`     the trailing bits play a significant role. CONG`\
47-`     is one of the most widely used generators of the`\
48-`     last 30 years, as it was the system generator for`\
49-`     VAX and was incorporated in several popular`\
50-`     software packages, all seemingly without complaint.`\
51-`     (note: this is about the same as Delphi uses)`
41+``` pascal
42+ CONG is a congruential generator with the widely used 69069
43+ multiplier: x(n)=69069x(n-1)+1234567. It has period
44+ 2^32. The leading half of its 32 bits seem to pass
45+ tests, but bits in the last half are too regular.
46+ The last half of the bits of CONG are too regular,
47+ and it fails tests for which those bits play a
48+ significant role. CONG+FIB will also have too much
49+ regularity in trailing bits, as each does. But keep
50+ in mind that it is a rare application for which
51+ the trailing bits play a significant role. CONG
52+ is one of the most widely used generators of the
53+ last 30 years, as it was the system generator for
54+ VAX and was incorporated in several popular
55+ software packages, all seemingly without complaint.
56+ (note: this is about the same as Delphi uses)
57+```
5258
5359### four lagged Fibonaci
5460
55-`     LFIB4 is an extension of what (Marsaglia) have previously`\
56-`     defined as a lagged Fibonacci generator:`\
57-`     x(n)=x(n-r) op x(n-s), with the x's in a finite`\
58-`     set over which there is a binary operation op, such`\
59-`     as +,- on integers mod 2^32, * on odd such integers,`\
60-`     exclusive-or(xor) on binary vectors. Except for`\
61-`     those using multiplication, lagged Fibonacci`\
62-`     generators fail various tests of randomness, unless`\
63-`     the lags are very long. (See SWB below).`\
64-`     To see if more than two lags would serve to overcome`\
65-`     the problems of 2-lag generators using +,- or xor, I`\
66-`     have developed the 4-lag generator LFIB4 using`\
67-`     addition: x(n)=x(n-256)+x(n-179)+x(n-119)+x(n-55)`\
68-`     mod 2^32. Its period is 2^31*(2^256-1), about 2^287,`\
69-`     and it seems to pass all tests---in particular,`\
70-`     those of the kind for which 2-lag generators using`\
71-`     +,-,xor seem to fail. For even more confidence in`\
72-`     its suitability, LFIB4 can be combined with KISS,`\
73-`     with a resulting period of about 2^410: just use`\
74-`     (KISS+LFIB4) in any Pascal expression.`
61+```text
62+ LFIB4 is an extension of what (Marsaglia) have previously
63+ defined as a lagged Fibonacci generator:
64+ x(n)=x(n-r) op x(n-s), with the x's in a finite
65+ set over which there is a binary operation op, such
66+ as +,- on integers mod 2^32, * on odd such integers,
67+ exclusive-or(xor) on binary vectors. Except for
68+ those using multiplication, lagged Fibonacci
69+ generators fail various tests of randomness, unless
70+ the lags are very long. (See SWB below).
71+ To see if more than two lags would serve to overcome
72+ the problems of 2-lag generators using +,- or xor, I
73+ have developed the 4-lag generator LFIB4 using
74+ addition: x(n)=x(n-256)+x(n-179)+x(n-119)+x(n-55)
75+ mod 2^32. Its period is 2^31*(2^256-1), about 2^287,
76+ and it seems to pass all tests---in particular,
77+ those of the kind for which 2-lag generators using
78+ +,-,xor seem to fail. For even more confidence in
79+ its suitability, LFIB4 can be combined with KISS,
80+ with a resulting period of about 2^410: just use
81+ (KISS+LFIB4) in any Pascal expression.
82+```
7583
7684### 3-shift-register generator
7785
78-`     SHR3 is a 3-shift-register generator with period`\
79-`     2^32-1. It uses y(n)=y(n-1)(I+L^17)(I+R^13)(I+L^5),`\
80-`     with the y's viewed as binary vectors, L the 32x32`\
81-`     binary matrix that shifts a vector left 1, and R its`\
82-`     transpose. SHR3 seems to pass all except those`\
83-`     related to the binary rank test, since 32 successive`\
84-`     values, as binary vectors, must be linearly`\
85-`     independent, while 32 successive truly random 32-bit`\
86-`     integers, viewed as binary vectors, will be linearly`\
87-`     independent only about 29% of the time.`\
88-`     (note: this is an ancestor of XORSHIFT (2003) and very fast)`
86+``` pascal
87+ SHR3 is a 3-shift-register generator with period
88+ 2^32-1. It uses y(n)=y(n-1)(I+L^17)(I+R^13)(I+L^5),
89+ with the y's viewed as binary vectors, L the 32x32
90+ binary matrix that shifts a vector left 1, and R its
91+ transpose. SHR3 seems to pass all except those
92+ related to the binary rank test, since 32 successive
93+ values, as binary vectors, must be linearly
94+ independent, while 32 successive truly random 32-bit
95+ integers, viewed as binary vectors, will be linearly
96+ independent only about 29% of the time.
97+ (note: this is an ancestor of XORSHIFT (2003) and very fast)
98+```
8999
90100### Subtract with borrow
91101
92-`     SWB is a subtract-with-borrow generator that I`\
93-`     developed to give a simple method for producing`\
94-`     extremely long periods:`\
95-`     x(n)=x(n-222)-x(n-237)- borrow mod 2^32.`\
96-`     The 'borrow' is 0, or set to 1 if computing x(n-1)`\
97-`     caused overflow in 32-bit integer arithmetic. This`\
98-`     generator has a very long period, 2^7098(2^480-1),`\
99-`     about 2^7578. It seems to pass all tests of`\
100-`     randomness, except for the Birthday Spacings test,`\
101-`     which it fails badly, as do all lagged Fibonacci`\
102-`     generators using +,- or xor. I would suggest`\
103-`     combining SWB with KISS, MWC, SHR3, or CONG.`\
104-`     KISS+SWB has period >2^7700 and is highly`\
105-`     recommended.`\
106-`     Subtract-with-borrow has the same local behaviour`\
107-`     as lagged Fibonacci using +,-,xor---the borrow`\
108-`     merely provides a much longer period.`\
109-`     SWB fails the birthday spacings test, as do all`\
110-`     lagged Fibonacci and other generators that merely`\
111-`     combine two previous values by means of =,- or xor.`\
112-`     Those failures are for a particular case: m=512`\
113-`     birthdays in a year of n=2^24 days. There are`\
114-`     choices of m and n for which lags >1000 will also`\
115-`     fail the test. A reasonable precaution is to always`\
116-`     combine a 2-lag Fibonacci or SWB generator with`\
117-`     another kind of generator, unless the generator uses`\
118-`     *, for which a very satisfactory sequence of odd`\
119-`     32-bit integers results.`
102+``` pascal
103+ SWB is a subtract-with-borrow generator that I
104+ developed to give a simple method for producing
105+ extremely long periods:
106+ x(n)=x(n-222)-x(n-237)- borrow mod 2^32.
107+ The 'borrow' is 0, or set to 1 if computing x(n-1)
108+ caused overflow in 32-bit integer arithmetic. This
109+ generator has a very long period, 2^7098(2^480-1),
110+ about 2^7578. It seems to pass all tests of
111+ randomness, except for the Birthday Spacings test,
112+ which it fails badly, as do all lagged Fibonacci
113+ generators using +,- or xor. I would suggest
114+ combining SWB with KISS, MWC, SHR3, or CONG.
115+ KISS+SWB has period >2^7700 and is highly
116+ recommended.
117+ Subtract-with-borrow has the same local behaviour
118+ as lagged Fibonacci using +,-,xor---the borrow
119+ merely provides a much longer period.
120+ SWB fails the birthday spacings test, as do all
121+ lagged Fibonacci and other generators that merely
122+ combine two previous values by means of =,- or xor.
123+ Those failures are for a particular case: m=512
124+ birthdays in a year of n=2^24 days. There are
125+ choices of m and n for which lags >1000 will also
126+ fail the test. A reasonable precaution is to always
127+ combine a 2-lag Fibonacci or SWB generator with
128+ another kind of generator, unless the generator uses
129+ *, for which a very satisfactory sequence of odd
130+ 32-bit integers results.
131+```
120132
121133### Keep it simple sir
122134
123-`     The KISS generator, (Keep It Simple Stupid), is`\
124-`     designed to combine the two multiply-with-carry`\
125-`     generators in MWC with the 3-shift register SHR3 and`\
126-`     the congruential generator CONG, using addition and`\
127-`     exclusive-or. Period about 2^123.`\
128-`     It is one of (Marsaglia's) favorite generators.`
135+```text
136+ The KISS generator, (Keep It Simple Stupid), is
137+ designed to combine the two multiply-with-carry
138+ generators in MWC with the 3-shift register SHR3 and
139+ the congruential generator CONG, using addition and
140+ exclusive-or. Period about 2^123.
141+ It is one of (Marsaglia's) favorite generators.
142+```
129143
130144### XorShift
131145
132146This was not in the original post, it was only published in 2003, but since it is a similar algorithm and it is of importance for its speed and properties.\
133147I feel it is warranted to include it here.\
134148It is similar to the SHR3 PRNG and is described in detail in the pdf you find in the link below.
135149
136150### Implementation
137151
138152The class is a fully static sealed class so you do not have to instantiate it or derive from it. If you use just one of the PRNGs it is easy to factor it out. Note that several Prng's maintain state:
139153
140154Here I have implemented it as a static class, but it is also possible to write them as a function with writable typed constants. I would recommend the static class approach, though, since calling the class constructor can reset the PRNG, Something that is not (easily) possible with the functional approach unless you would use global variables or double the state fields. Note that if you assign the output to a signed integer, you will get signed result, if needed.
141155
142156``` pascal
143157unit uMarsaglia;
144158{ author: Thaddy de Koning, based on the original C sourcecode in the original post by George Marsaglia }
145159{ http://www.ciphersbyritter.com/NEWS4/RANDC.HTM#369B5E30.65A55FD1@stat.fsu.edu }
146160{$mode objfpc}
147161interface
148162
149163type
150164 { It is not meant to be extended, since it is based on a single article so I sealed the class }
151165 TMarsagliaPrngs = class sealed
152166 strict private
153167 {static variables:}
154168 class var T: array [0..255] of dword;
155169 class var a,b,v,w,x,y,z,bro,jcong,jsr:dword;
156170 class var c:byte;
157171 class function wnew:dword;inline;static;
158172 class function znew:dword;inline;static;
159173 public
160174 { You can also call the create explicitly to reset the static
161175 variables.}
162176 class constructor create;
163177 class procedure settable(const i1,i2,i3,i4,i5,i6:dword);static;
164178 class procedure Randomize;inline;static;
165179 class function MWC:dword;inline;static;
166180 class function SHR3:dword;inline;static;
167181 class function CONG:dword;inline;static;
168182 class function FIB:dword;inline;static;
169183 class function KISS:dword;inline;static;
170184 class function LFIB4:dword;inline;static;
171185 class function SWB:dword;inline;static;
172186 // From XSorShift Prng's, 2003
173187 class function XOS:dword;inline;static;
174188 end;
175189
176190 { a shorthand alias: }
177191 Prng = type TMarsagliaPrngs;
178192
179193 { Because many simulations call for uniform
180194 random variables in 0<x<1 or -1<x<1, here's a utility function.
181195 Using Signed = false will provide a random single in range (0,1),
182196 while Signed true will provide one in range (-1,1).
183197 Simply feed it the output of one of the Prng's}
184198 function ScaleToFloat(const value:dword; Signed:Boolean=true):single;inline;
185199
186200
187201implementation
188202
189203 class function TMarsagliaPrngs.znew:dword;inline;static;
190204 begin
191205 z:=36969*(z and 65535)+(z >> 16);
192206 Result:= z;
193207 end;
194208
195209 class function TMarsagliaPrngs.wnew:dword;inline;static;
196210 begin
197211 w:=18000*(w and 65535)+(w >> 16);
198212 Result := w;
199213 end;
200214
201215
202216 class constructor TMarsagliaPrngs.create;
203217 begin
204218 a:=224466889; b:=7584631;c:=0;v:=2463534242;w:=521288629;x:=0;y:=0;
205219 z:=362436069;jsr:=123456789; jcong:=380116160;
206220 settable(12345,65435,34221,12345,9983651,95746118);
207221 end;
208222
209223 { Use random seeds to reset z,w,jsr,jcong,a,b, and the table t[256]}
210224 class procedure TMarsagliaPrngs.settable(const i1,i2,i3,i4,i5,i6:dword);inline;static;
211225 var
212226 i:integer;
213227 begin
214228 z:=i1; w:=i2; jsr:=i3; jcong:=i4; a:=i5; b:=i6;
215229 { Example procedure to set the table using KISS:}
216230 for i:=0 to 255 do t[i]:=KISS;
217231 end;
218232
219233 class function TMarsagliaPrngs.MWC:dword;inline;static;
220234 begin
221235 Result:= (znew << 16)+wnew;
222236 end;
223237
224238 class function TMarsagliaPrngs.SHR3:dword;inline;static;
225239 begin
226240 jsr:=jsr xor(jsr << 17); jsr:= jsr xor(jsr >> 13); jsr:=jsr xor(jsr << 5);
227241 result := jsr;
228242 end;
229243
230244 class function TMarsagliaPrngs.CONG:dword;inline;static;
231245 begin
232246 jcong:=69069*jcong+1234567;
233247 Result := jcong;
234248 end;
235249
236250 class function TMarsagliaPrngs.FIB:dword;inline;static;
237251 begin
238252 b:=a+b;a:=b-a;
239253 Result := a;
240254 end;
241255
242256 class function TMarsagliaPrngs.KISS:dword;inline;static;
243257 begin
244258 Result:= (MWC xor CONG)+SHR3;
245259 end;
246260
247261 class function TMarsagliaPrngs.LFIB4:dword;inline;static;
248262 begin
249263 inc(c);
250264 t[c]:=t[c]+t[(c+58) and 255]+t[(c+119) and 255]+t[(c+178) and 255];
251265 Result:=t[c];
252266 end;
253267
254268 class function TMarsagliaPrngs.SWB:dword;inline;static;
255269 begin
256270 inc(c);bro:=ord(x<y);
257271 x:=t[(c+34) and 255];
258272 y:=t[(c+19) and 255]+bro;
259273 t[c] := x-y;
260274 Result:= t[c];
261275 end;
262276
263277 class function TMarsagliaPrngs.XOS:dword;inline;static;
264278 var tmp:dword;
265279 begin
266280 tmp:=(x xor (x<<15));
267281 x:=y;
268282 y:=z;
269283 z:=w;
270284 w:=(w xor (w>>21)) xor (tmp xor(tmp>>4));
271285 Result := w;
272286 end;
273287
274288 class procedure TMarsagliaPrngs.Randomize;inline;static;
275289 begin
276290 {$WARNING This is not implemented here yet, so use system.random}
277291 system.randomize;
278292 for i:=0 to 255 do t[i]:=system.random(256);
279293 end;
280294
281295 function ScaleToFloat(const value:dword; Signed:Boolean=true):single;inline;
282296 begin
283297 case Signed of
284298 False:Result:= dword(value) * 2.328306e-10;
285299 True :Result:= longint(value) * 4.656613e-10;
286300 end;
287301 end;
288302
289303 end.
290304```
291305
292306### Example
293307
294308The original C code contained an example, here translated to Object Pascal. It is just a proof that we get the same results as intended.
295309
296310``` pascal
297311{$mode objfpc}{$I-}
298312{ It is basically Marsaglia's C test program in Pascal }
299313uses
300314 uMarsaglia;
301315var
302316 i:integer;
303317 k:dword;
304318begin
305319 // Set the table to known values, so we can verify the results.
306320 // The one million iterations per prng are not a mistake.
307321 // The whole program will finish in less than .5 seconds even on a Raspberry pi one.
308322 // Order is important: some seeds are taken from other prng's
309323 Prng.settable(12345,65435,34221,12345,9983651,95746118);
310324 for i:=1 to 1000000 do k:=Prng.LFIB4;writeln('LFIB4 ','Pass: ',k-1064612766=0);
311325 for i:=1 to 1000000 do k:=Prng.SWB; writeln('SWB ','Pass: ',k- 627749721=0);
312326 for i:=1 to 1000000 do k:=Prng.KISS; writeln('KISS ','Pass: ',k-1372460312=0);
313327 for i:=1 to 1000000 do k:=Prng.CONG; writeln('CONG ','Pass: ',k-1529210297=0);
314328 for i:=1 to 1000000 do k:=Prng.SHR3; writeln('SHR3 ','Pass: ',k-2642725982=0);
315329 for i:=1 to 1000000 do k:=Prng.MWC; writeln('MWC ','Pass: ',k- 904977562=0);
316330 for i:=1 to 1000000 do k:=Prng.FIB; writeln('FIB ','Pass: ',k-3519793928=0);
317331 for i:=1 to 1000000 do k:=Prng.XOS; writeln('XOS ','Pass: ',k-1110212780=0);
318332 // Scale to float demo
319333 for i := 0 to 9 do writeln(ScaleToFloat(Prng.SHR3):2:10);
320334 writeln('..........................');
321335 for i := 0 to 9 do writeln(ScaleToFloat(Prng.SHR3, false):2:10);
322336end.
323337```
324338
325339### Tests
326340
327341The above routines are tested with FIPS-140-2 and 1000 series, which is just ***indicative***. A full Dieharder test for any of them can be found on the web.
328342
329343Note there can always be a "failure" since statistics predict that a failure can always happen by chance. Any of the algorithms with few failures is usable. Particulary the fast KISS+SWB that passes Dieharder and XOS that is blazing fast and as good as FPC's Mersenne twister. ote you can also test against Crush/BigCrush, but that test suite is not suitable for 32 bit Prng's as per their documentation.
330344
331345![](/wiki-images/5/52/Selection_018.png)
332346
333347### Can we have one of these as the standard random, plz?
334348
335349**No you can't.**
336350
337351Programming languages tend not to change their PRNG and Free Pascal has a good one. It is important for a programming language to have a predictable PRNG. Many software relies on that. It would break backwards compatibility if the PRNG in a programming language would change. This is not a minor, but a very major issue. E.g. many existing datasets would become untestible and unverifiable. Even if some of the prng's that are presented here are superior in speed and properties, because of the above I hope the default Random will never be replaced. And you should hope so too..
338352
339353But if you need higher speed or better properties, just use one or a combination of these...
340354
341355### References
342356
3433571. [The original post by George Marsaglia](http://www.ciphersbyritter.com/NEWS4/RANDC.HTM#369B5E30.65A55FD1@stat.fsu.edu)\
3443582. [George Marsaglia, XorShift RNG's, Journal of Statistical Software volume 8 issue 14, July 2003](http://www.jstatsoft.org/v08/i14/paper)
3453593. [Diehard tests](https://en.wikipedia.org/wiki/Diehard_tests)
3463604. [Dieharder tests homepage, diegarder available for most linux distro's](http://www.phy.duke.edu/~rgb/General/dieharder.php)
3473615. [FIPS-140-2](https://en.wikipedia.org/wiki/FIPS_140-2)
3483626. [rngtest 2, available for most linux distro's](https://wiki.archlinux.org/index.php/Rng-tools)
349363
350364## See also
351365
352366- [[Dev random]]
353367- [[Functions for descriptive statistics]]
354368- [[Generating Random Numbers]]
355369
356370### TODO
357371
358372Provide X-platform randomize (other than FPC's random) based on timer, Windows RtlGenRandom and Unix /dev/random, /dev/urandom to seed the prng's
359373
360374---
361375*From the Free Pascal wiki, page "Marsaglia's pseudo random number generators", as it stood on 2024-03-03 after a change by Thaddy, out of the wiki's dump of 2024-05-26. Copyright stays with the author who wrote it, according to the old wiki's [[Licensing]] page. If you have concerns about this material, [message the admin](/wiki-contact).*
362376Category: Code
363377Category: Mathematics
364378Category: Modelling and Simulation
365379Category: Statistical algorithms

Back to history