Sign in Join
Wiki copytemporary

A simple implementation of the Mersenne twister: what changed

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

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.
1+## Contents
132
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.
3+ * 1 TMersenne32
4+ * 2 C++ verification test:
5+ * 3 Verification against FPC's Random
6+ * 4 Reference
7+ * 5 See also
178
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
9+### TMersenne32
5110
52-class constructor TMersenne32.Create;
53-begin
54- initialize(5489); // default seed for mt199937
55-end;
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.
5618
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;
19+Enjoy,
20+Thaddy.
8121
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;
22+It is verified against several mt19937 implementations in different languages, including Freepascal and C/C++.
9823
99-class function TMersenne32.URandf:single;inline;static;
100-begin
101- Result := URand * 2.32830643653869628906e-10;
102-end;
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.
10326
104-class function TMersenne32.SRandf:single;inline;static;
105-begin
106- Result :=URand * 4.6566128730773926e-010;
107-end;
27+
10828
109-class function TMersenne32.SRand:integer;inline;static;
110-begin
111- {$ifopt D+}Result := 0;{$endif}
112- Result := URand;
113-end;
114-end.
11529```
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.
127+```
116128
117129### C++ verification test:
118130
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
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+```
123134#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+ #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+ }
135146```
136147
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
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+```
141151#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+ #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+ }
152162```
153163
154-Left FPC, right C++\
155-![](/wiki-images/6/68/Selection_019.png)
164+Left FPC, right C++
165+[![Selection 019.png](/wiki-images/6/68/Selection_019.png)](</File:Selection_019.png>)
156166
157167### Verification against FPC's Random
158-
159-``` pascal
168+```
160169program 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+ {$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.
170179```
171180
172181## Reference
173182
174-1. [Mersenne twister in WikiPedia](https://en.wikipedia.org/wiki/Mersenne_Twister)
183+ 1. [Mersenne twister in WikiPedia](<https://en.wikipedia.org/wiki/Mersenne_Twister>)
175184
176185## See also
177186
178-- [[Dev random]]
179-- [[Functions for descriptive statistics]]
180-- [[Generating Random Numbers]]
181-- [[Marsaglia&|Marsaglia's pseudo random number generators]]
187+ * [[Dev random]]
188+ * [[Functions for descriptive statistics]]
189+ * [[Generating Random Numbers]]
190+ * [[Marsaglia's pseudo random number generators]]
182191
183192---
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
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).*

Back to history