Sign in Join
Wiki copytemporary

A simple implementation of the Mersenne twister: what changed

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

Title: A simple implementation of the Mersenne twister

OldNewLine
1+### TMersenne32
2+
3+When you look at the source code that implements Free Pascal's [[doc:rtl/system/random.html|random]] it may look overwelmingly complex. The algorithm itself, though, is very understandable.\
4+Here is a simple 32 bit implementation according to MT19937. It is compatible with FPC's Random for 32 bit values.\
5+You may actually prefer this one over the one in system. That is because system is usually compiled with conservative optimization settings.\
6+Although the core implementation is about the same, you can compile this with more agressive compiler optimization settings.\
7+E.g. on armhf I got a 50% better performance using fpc -CX -XXs -OpARMV7a -CpARMV7A -CfVFPV4 -Sv -OoFASTMATH -O4 -B mersenne.pas\
8+E.g. on AMD64 I got a 30% better performance using fpc -CX -XXs -OpATHLON64 -CpATHLON64 -CfSSE41 -Sv -OoFASTMATH -B mersenne.pas\
9+The performance gain is solely because of the compiler settings.
10+
11+Enjoy,\
12+Thaddy.
13+
14+It is verified against several mt19937 implementations in different languages, including Freepascal and C/C++.\
15+Of course I added this one to wikipedia too. You don't have to create this class. it is fully static.\
16+Simply call Rnd32.URand, URandf, Srand of Srandf to obtain a random value.
17+
18+``` pascal
19+{ new version 2017-02-14 }
20+unit mersenne;
21+{$mode objfpc}{$A8}{$R-}{$Q-}
22+interface
23+type
24+ TMersenne32 = class sealed
25+ strict private
26+ const
27+ // Define MT19937 constants (32-bit RNG)
28+ N = 624;M = 397;R = 31;A = $9908B0DF;F = 1812433253;
29+ U = 11;S = 7;B = $9D2C5680;T = 15;C = $EFC60000;L = 18;
30+ MASK_LOWER = (QWORD(1) << R) - 1;
31+ MASK_UPPER = QWORD(1) << R;
32+ class var mt:array[0..N-1] of dword;
33+ class var index:word;
34+ class procedure twist;inline;static;
35+ public
36+ class constructor create;
37+ class procedure initialize(const seed:dword);inline;static;
38+ { 32 bit unsigned integer }
39+ class function URand:dword;inline;static;
40+ { 32 bit float in range 0..1 }
41+ class function URandf:single;inline;static;
42+ { 32 bit signed integer }
43+ class function SRand:integer;inline;static;
44+ {32 bit float in the range -1..1 }
45+ class function SRandf:single;inline;static;
46+ end;
47+
48+ Rnd32 = type TMersenne32;
49+
50+implementation
51+
52+class constructor TMersenne32.Create;
53+begin
54+ initialize(5489); // default seed for mt199937
55+end;
56+
57+class procedure TMersenne32.Initialize(const seed:dword);inline;static;
58+var
59+ i:dword;
60+begin
61+ mt[0] := seed;
62+ for i := 1 to pred(N) do
63+ mt[i] := F * (mt[i - 1] xor (mt[i - 1] >> 30)) + i;
64+ index := N;
65+end;
66+
67+class procedure TMersenne32.Twist;inline;static;
68+var
69+ i:integer;
70+begin
71+ for i:=0 to N-M-1 do
72+ mt[i]:=mt[i+M] xor {twist} (((mt[i] and MASK_UPPER) or
73+ (mt[i+1] and MASK_LOWER)) shr 1)xor(dword(-(mt[i+1] and 1)) and A);
74+ for i:=N-M to N-2 do
75+ mt[i]:=mt[i+(M-N)]xor{twist}(((mt[i] and MASK_UPPER) or
76+ (mt[i+1] and MASK_LOWER)) shr 1)xor(dword(-(mt[i+1] and 1)) and A);
77+ mt[N-1]:=mt[M-1] xor {twist} (((mt[n-1] and MASK_UPPER) or (mt[0] and
78+ MASK_LOWER)) shr 1)xor(dword(-(mt[0] and 1)) and A);
79+ index:=0;
80+end;
81+
82+class function TMersenne32.URand:dword;inline;static;
83+var
84+ i:integer;
85+begin
86+ i := index;
87+ if index >= N then
88+ begin
89+ Twist;
90+ i := index;
91+ end;
92+ Result := mt[i];index := i + 1;
93+ Result := Result xor (mt[i] >> U);
94+ Result := Result xor (Result << S) and B;
95+ Result := Result xor (Result << T) and C;
96+ Result := Result xor(Result >> L);
97+end;
98+
99+class function TMersenne32.URandf:single;inline;static;
100+begin
101+ Result := URand * 2.32830643653869628906e-10;
102+end;
103+
104+class function TMersenne32.SRandf:single;inline;static;
105+begin
106+ Result :=URand * 4.6566128730773926e-010;
107+end;
108+
109+class function TMersenne32.SRand:integer;inline;static;
110+begin
111+ {$ifopt D+}Result := 0;{$endif}
112+ Result := URand;
113+end;
114+end.
115+```
116+
117+### C++ verification test:
118+
119+This code is for a 32 bit float (single) and can be used to check against FPC's Random as well, provided you cast the output to single:\
120+Compile this with g++ mt_cpp_float32.cpp -std=gnu++11 -o mt_cpp_float32
121+
122+``` c
123+#include <iostream>
124+#include <random>
125+using namespace std;
126+
127+int main()
128+{
129+ mt19937 mt_rand(5489ul);
130+ for(int i=0;i<5;i++){
131+ cout << (float) mt_rand() * 2.32830643653869628906e-10 << endl;};
132+
133+ return 0;
134+}
135+```
136+
137+Compile this with g++ mt_cpp_uint32.cpp -std=gnu++11 -o mt_cpp_uint32\
138+Run ./mersplusplus and compare the output with TMersenne32. This code is for 32 bit unsigned. It should be (is) equal.
139+
140+``` c
141+#include <iostream>
142+#include <random>
143+using namespace std;
144+ int main()
145+{
146+ mt19937 mt_rand(5489u);
147+ for(int i=0;i<5;i++){
148+ cout << mt_rand() << endl;};
149+
150+ return 0;
151+}
152+```
153+
154+Left FPC, right C++\
155+![](/wiki-images/6/68/Selection_019.png)
156+
157+### Verification against FPC's Random
158+
159+``` pascal
160+program mersennetest2;
161+{$mode objfpc}
162+uses sysutils,dateutils,mersenne;
163+var
164+ i:integer;
165+begin
166+ Randseed :=5489; // default seed for mt199937
167+ for i := 0 to 9999999 do Rnd32.urandf;write(Rnd32.urandf :2:10,' ');writeln('Mersenne32');
168+ for i := 0 to 9999999 do Random ;write(single(Random) :2:10,' ');writeln('FPC Random');
169+end.
170+```
171+
172+## Reference
173+
174+1. [Mersenne twister in WikiPedia](https://en.wikipedia.org/wiki/Mersenne_Twister)
175+
176+## See also
177+
178+- [[Dev random]]
179+- [[Functions for descriptive statistics]]
180+- [[Generating Random Numbers]]
181+- [[Marsaglia&|Marsaglia's pseudo random number generators]]
182+
183+---
184+*From the Free Pascal wiki, page "A simple implementation of the Mersenne twister", as it stood on 2020-02-01 after a change by GAN, 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).*
185+Category: Code
186+Category: Modelling and Simulation
187+Category: Statistical algorithms

Back to history