Marsaglia's pseudo random number generators: what changed

Revision 0 to 1. Removed lines have a minus sign; added lines have a plus sign.

Title: Marsaglia's pseudo random number generators

New		Line
1	+	<figure>
2	+	<img src="fpc_source_logo.png" title="fpc_source_logo.png" />
3	+	<figcaption>fpc_source_logo.png</figcaption>
4	+	</figure>
5	+
6	+	## Marsaglia's pseudo random number generators
7	+
8	+	Although 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.
9	+
10	+	In 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.
11	+
12	+	Here 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.
13	+
14	+	Marsaglia describes them as follows:
15	+
16	+	### Multiply with carry
17	+
18	+	` The MWC generator concatenates two 16-bit multiply-`\
19	+	` with-carry generators, x(n)=36969x(n-1)+carry,`\
20	+	` y(n)=18000y(n-1)+carry mod 2^16, has period about`\
21	+	` 2^60 and seems to pass all tests of randomness. A`\
22	+	` favorite stand-alone generator---faster than KISS,`\
23	+	` which contains it. `
24	+
25	+	### fibonaci
26	+
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	+
38	+	### Congruential generator
39	+
40	+	` CONG is a congruential generator with the widely used 69069`\
41	+	` multiplier: x(n)=69069x(n-1)+1234567. It has period`\
42	+	` 2^32. The leading half of its 32 bits seem to pass`\
43	+	` tests, but bits in the last half are too regular.`\
44	+	` The last half of the bits of CONG are too regular,`\
45	+	` and it fails tests for which those bits play a`\
46	+	` significant role. CONG+FIB will also have too much`\
47	+	` regularity in trailing bits, as each does. But keep`\
48	+	` in mind that it is a rare application for which`\
49	+	` the trailing bits play a significant role. CONG`\
50	+	` is one of the most widely used generators of the`\
51	+	` last 30 years, as it was the system generator for`\
52	+	` VAX and was incorporated in several popular`\
53	+	` software packages, all seemingly without complaint.`\
54	+	` (note: this is about the same as Delphi uses)`
55	+
56	+	### four lagged Fibonaci
57	+
58	+	` LFIB4 is an extension of what (Marsaglia) have previously`\
59	+	` defined as a lagged Fibonacci generator:`\
60	+	` x(n)=x(n-r) op x(n-s), with the x's in a finite`\
61	+	` set over which there is a binary operation op, such`\
62	+	` as +,- on integers mod 2^32, * on odd such integers,`\
63	+	` exclusive-or(xor) on binary vectors. Except for`\
64	+	` those using multiplication, lagged Fibonacci`\
65	+	` generators fail various tests of randomness, unless`\
66	+	` the lags are very long. (See SWB below).`\
67	+	` To see if more than two lags would serve to overcome`\
68	+	` the problems of 2-lag generators using +,- or xor, I`\
69	+	` have developed the 4-lag generator LFIB4 using`\
70	+	` addition: x(n)=x(n-256)+x(n-179)+x(n-119)+x(n-55)`\
71	+	` mod 2^32. Its period is 2^31*(2^256-1), about 2^287,`\
72	+	` and it seems to pass all tests---in particular,`\
73	+	` those of the kind for which 2-lag generators using`\
74	+	` +,-,xor seem to fail. For even more confidence in`\
75	+	` its suitability, LFIB4 can be combined with KISS,`\
76	+	` with a resulting period of about 2^410: just use`\
77	+	` (KISS+LFIB4) in any Pascal expression.`
78	+
79	+	### 3-shift-register generator
80	+
81	+	` SHR3 is a 3-shift-register generator with period`\
82	+	` 2^32-1. It uses y(n)=y(n-1)(I+L^17)(I+R^13)(I+L^5),`\
83	+	` with the y's viewed as binary vectors, L the 32x32`\
84	+	` binary matrix that shifts a vector left 1, and R its`\
85	+	` transpose. SHR3 seems to pass all except those`\
86	+	` related to the binary rank test, since 32 successive`\
87	+	` values, as binary vectors, must be linearly`\
88	+	` independent, while 32 successive truly random 32-bit`\
89	+	` integers, viewed as binary vectors, will be linearly`\
90	+	` independent only about 29% of the time.`\
91	+	` (note: this is an ancestor of XORSHIFT (2003) and very fast)`
92	+
93	+	### Subtract with borrow
94	+
95	+	` SWB is a subtract-with-borrow generator that I`\
96	+	` developed to give a simple method for producing`\
97	+	` extremely long periods:`\
98	+	` x(n)=x(n-222)-x(n-237)- borrow mod 2^32.`\
99	+	` The 'borrow' is 0, or set to 1 if computing x(n-1)`\
100	+	` caused overflow in 32-bit integer arithmetic. This`\
101	+	` generator has a very long period, 2^7098(2^480-1),`\
102	+	` about 2^7578. It seems to pass all tests of`\
103	+	` randomness, except for the Birthday Spacings test,`\
104	+	` which it fails badly, as do all lagged Fibonacci`\
105	+	` generators using +,- or xor. I would suggest`\
106	+	` combining SWB with KISS, MWC, SHR3, or CONG.`\
107	+	` KISS+SWB has period >2^7700 and is highly`\
108	+	` recommended.`\
109	+	` Subtract-with-borrow has the same local behaviour`\
110	+	` as lagged Fibonacci using +,-,xor---the borrow`\
111	+	` merely provides a much longer period.`\
112	+	` SWB fails the birthday spacings test, as do all`\
113	+	` lagged Fibonacci and other generators that merely`\
114	+	` combine two previous values by means of =,- or xor.`\
115	+	` Those failures are for a particular case: m=512`\
116	+	` birthdays in a year of n=2^24 days. There are`\
117	+	` choices of m and n for which lags >1000 will also`\
118	+	` fail the test. A reasonable precaution is to always`\
119	+	` combine a 2-lag Fibonacci or SWB generator with`\
120	+	` another kind of generator, unless the generator uses`\
121	+	` *, for which a very satisfactory sequence of odd`\
122	+	` 32-bit integers results.`
123	+
124	+	### Keep it simple sir
125	+
126	+	` The KISS generator, (Keep It Simple Stupid), is`\
127	+	` designed to combine the two multiply-with-carry`\
128	+	` generators in MWC with the 3-shift register SHR3 and`\
129	+	` the congruential generator CONG, using addition and`\
130	+	` exclusive-or. Period about 2^123.`\
131	+	` It is one of (Marsaglia's) favorite generators.`
132	+
133	+	### XorShift
134	+
135	+	This 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.\
136	+	I feel it is warranted to include it here.\
137	+	It is similar to the SHR3 PRNG and is described in detail in the pdf you find in the link below.
138	+
139	+	### Implementation
140	+
141	+	The 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:
142	+
143	+	Here 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.
144	+
145	+	``` pascal
146	+	unit uMarsaglia;
147	+	{ author: Thaddy de Koning, based on the original C sourcecode in the original post by George Marsaglia }
148	+	{ http://www.ciphersbyritter.com/NEWS4/RANDC.HTM#369B5E30.65A55FD1@stat.fsu.edu }
149	+	{$mode objfpc}
150	+	interface
151	+
152	+	type
153	+	{ It is not meant to be extended, since it is based on a single article so I sealed the class }
154	+	TMarsagliaPrngs = class sealed
155	+	strict private
156	+	{static variables:}
157	+	class var T: array [0..255] of dword;
158	+	class var a,b,v,w,x,y,z,bro,jcong,jsr:dword;
159	+	class var c:byte;
160	+	class function wnew:dword;inline;static;
161	+	class function znew:dword;inline;static;
162	+	public
163	+	{ You can also call the create explicitly to reset the static
164	+	variables.}
165	+	class constructor create;
166	+	class procedure settable(const i1,i2,i3,i4,i5,i6:dword);static;
167	+	class procedure Randomize;inline;static;
168	+	class function MWC:dword;inline;static;
169	+	class function SHR3:dword;inline;static;
170	+	class function CONG:dword;inline;static;
171	+	class function FIB:dword;inline;static;
172	+	class function KISS:dword;inline;static;
173	+	class function LFIB4:dword;inline;static;
174	+	class function SWB:dword;inline;static;
175	+	// From XSorShift Prng's, 2003
176	+	class function XOS:dword;inline;static;
177	+	end;
178	+
179	+	{ a shorthand alias: }
180	+	Prng = type TMarsagliaPrngs;
181	+
182	+	{ Because many simulations call for uniform
183	+	random variables in 0<x<1 or -1<x<1, here's a utility function.
184	+	Using Signed = false will provide a random single in range (0,1),
185	+	while Signed true will provide one in range (-1,1).
186	+	Simply feed it the output of one of the Prng's}
187	+	function ScaleToFloat(const value:dword; Signed:Boolean=true):single;inline;
188	+
189	+
190	+	implementation
191	+
192	+	class function TMarsagliaPrngs.znew:dword;inline;static;
193	+	begin
194	+	z:=36969*(z and 65535)+(z >> 16);
195	+	Result:= z;
196	+	end;
197	+
198	+	class function TMarsagliaPrngs.wnew:dword;inline;static;
199	+	begin
200	+	w:=18000*(w and 65535)+(w >> 16);
201	+	Result := w;
202	+	end;
203	+
204	+
205	+	class constructor TMarsagliaPrngs.create;
206	+	begin
207	+	a:=224466889; b:=7584631;c:=0;v:=2463534242;w:=521288629;x:=0;y:=0;
208	+	z:=362436069;jsr:=123456789; jcong:=380116160;
209	+	settable(12345,65435,34221,12345,9983651,95746118);
210	+	end;
211	+
212	+	{ Use random seeds to reset z,w,jsr,jcong,a,b, and the table t[256]}
213	+	class procedure TMarsagliaPrngs.settable(const i1,i2,i3,i4,i5,i6:dword);inline;static;
214	+	var
215	+	i:integer;
216	+	begin
217	+	z:=i1; w:=i2; jsr:=i3; jcong:=i4; a:=i5; b:=i6;
218	+	{ Example procedure to set the table using KISS:}
219	+	for i:=0 to 255 do t[i]:=KISS;
220	+	end;
221	+
222	+	class function TMarsagliaPrngs.MWC:dword;inline;static;
223	+	begin
224	+	Result:= (znew << 16)+wnew;
225	+	end;
226	+
227	+	class function TMarsagliaPrngs.SHR3:dword;inline;static;
228	+	begin
229	+	jsr:=jsr xor(jsr << 17); jsr:= jsr xor(jsr >> 13); jsr:=jsr xor(jsr << 5);
230	+	result := jsr;
231	+	end;
232	+
233	+	class function TMarsagliaPrngs.CONG:dword;inline;static;
234	+	begin
235	+	jcong:=69069*jcong+1234567;
236	+	Result := jcong;
237	+	end;
238	+
239	+	class function TMarsagliaPrngs.FIB:dword;inline;static;
240	+	begin
241	+	b:=a+b;a:=b-a;
242	+	Result := a;
243	+	end;
244	+
245	+	class function TMarsagliaPrngs.KISS:dword;inline;static;
246	+	begin
247	+	Result:= (MWC xor CONG)+SHR3;
248	+	end;
249	+
250	+	class function TMarsagliaPrngs.LFIB4:dword;inline;static;
251	+	begin
252	+	inc(c);
253	+	t[c]:=t[c]+t[(c+58) and 255]+t[(c+119) and 255]+t[(c+178) and 255];
254	+	Result:=t[c];
255	+	end;
256	+
257	+	class function TMarsagliaPrngs.SWB:dword;inline;static;
258	+	begin
259	+	inc(c);bro:=ord(x<y);
260	+	x:=t[(c+34) and 255];
261	+	y:=t[(c+19) and 255]+bro;
262	+	t[c] := x-y;
263	+	Result:= t[c];
264	+	end;
265	+
266	+	class function TMarsagliaPrngs.XOS:dword;inline;static;
267	+	var tmp:dword;
268	+	begin
269	+	tmp:=(x xor (x<<15));
270	+	x:=y;
271	+	y:=z;
272	+	z:=w;
273	+	w:=(w xor (w>>21)) xor (tmp xor(tmp>>4));
274	+	Result := w;
275	+	end;
276	+
277	+	class procedure TMarsagliaPrngs.Randomize;inline;static;
278	+	begin
279	+	{$WARNING This is not implemented here yet, so use system.random}
280	+	system.randomize;
281	+	for i:=0 to 255 do t[i]:=system.random(256);
282	+	end;
283	+
284	+	function ScaleToFloat(const value:dword; Signed:Boolean=true):single;inline;
285	+	begin
286	+	case Signed of
287	+	False:Result:= dword(value) * 2.328306e-10;
288	+	True :Result:= longint(value) * 4.656613e-10;
289	+	end;
290	+	end;
291	+
292	+	end.
293	+	```
294	+
295	+	### Example
296	+
297	+	The original C code contained an example, here translated to Object Pascal. It is just a proof that we get the same results as intended.
298	+
299	+	``` pascal
300	+	{$mode objfpc}{$I-}
301	+	{ It is basically Marsaglia's C test program in Pascal }
302	+	uses
303	+	uMarsaglia;
304	+	var
305	+	i:integer;
306	+	k:dword;
307	+	begin
308	+	// Set the table to known values, so we can verify the results.
309	+	// The one million iterations per prng are not a mistake.
310	+	// The whole program will finish in less than .5 seconds even on a Raspberry pi one.
311	+	// Order is important: some seeds are taken from other prng's
312	+	Prng.settable(12345,65435,34221,12345,9983651,95746118);
313	+	for i:=1 to 1000000 do k:=Prng.LFIB4;writeln('LFIB4 ','Pass: ',k-1064612766=0);
314	+	for i:=1 to 1000000 do k:=Prng.SWB; writeln('SWB ','Pass: ',k- 627749721=0);
315	+	for i:=1 to 1000000 do k:=Prng.KISS; writeln('KISS ','Pass: ',k-1372460312=0);
316	+	for i:=1 to 1000000 do k:=Prng.CONG; writeln('CONG ','Pass: ',k-1529210297=0);
317	+	for i:=1 to 1000000 do k:=Prng.SHR3; writeln('SHR3 ','Pass: ',k-2642725982=0);
318	+	for i:=1 to 1000000 do k:=Prng.MWC; writeln('MWC ','Pass: ',k- 904977562=0);
319	+	for i:=1 to 1000000 do k:=Prng.FIB; writeln('FIB ','Pass: ',k-3519793928=0);
320	+	for i:=1 to 1000000 do k:=Prng.XOS; writeln('XOS ','Pass: ',k-1110212780=0);
321	+	// Scale to float demo
322	+	for i := 0 to 9 do writeln(ScaleToFloat(Prng.SHR3):2:10);
323	+	writeln('..........................');
324	+	for i := 0 to 9 do writeln(ScaleToFloat(Prng.SHR3, false):2:10);
325	+	end.
326	+	```
327	+
328	+	### Tests
329	+
330	+	The 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.
331	+
332	+	Note 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.
333	+
334	+	<figure>
335	+	<img src="Selection_018.png" title="Selection_018.png" />
336	+	<figcaption>Selection_018.png</figcaption>
337	+	</figure>
338	+
339	+	### Can we have one of these as the standard random, plz?
340	+
341	+	No you can't.
342	+
343	+	Programming 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..
344	+
345	+	But if you need higher speed or better properties, just use one or a combination of these...
346	+
347	+	### References
348	+
349	+	1. [The original post by George Marsaglia](http://www.ciphersbyritter.com/NEWS4/RANDC.HTM#369B5E30.65A55FD1@stat.fsu.edu)\
350	+	2. [George Marsaglia, XorShift RNG's, Journal of Statistical Software volume 8 issue 14, July 2003](http://www.jstatsoft.org/v08/i14/paper)
351	+	3. [Diehard tests](https://en.wikipedia.org/wiki/Diehard_tests)
352	+	4. [Dieharder tests homepage, diegarder available for most linux distro's](http://www.phy.duke.edu/~rgb/General/dieharder.php)
353	+	5. [FIPS-140-2](https://en.wikipedia.org/wiki/FIPS_140-2)
354	+	6. [rngtest 2, available for most linux distro's](https://wiki.archlinux.org/index.php/Rng-tools)
355	+
356	+	## See also
357	+
358	+	- [[Dev random]]
359	+	- [[Functions for descriptive statistics]]
360	+	- [[Generating Random Numbers]]
361	+
362	+	### TODO
363	+
364	+	Provide X-platform randomize (other than FPC's random) based on timer, Windows RtlGenRandom and Unix /dev/random, /dev/urandom to seed the prng's
365	+
366	+	---
367	+	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).
368	+	Category: Code
369	+	Category: Mathematics
370	+	Category: Modelling and Simulation
371	+	Category: Statistical algorithms

Back to history