Marsaglia's pseudo random number generators: what changed

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

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

Back to history