Sign in Join
Wiki copytemporary

A simple implementation of the Mersenne twister: what changed

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

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