Sign in Join
Wiki copytemporary

A simple implementation of the Mersenne twister: what changed

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

OldNewLine
11### TMersenne32
22
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.
3+This is an implementation of the Mersenne Twister that is used in FreePascal versions prior to 3.3.1
4+It is verified against several mt19937 implementations in different languages, including Freepascal and C/C++.
5+The implementation is fully static and threadsafe in the sense that given the same seed the sequence is the same for all threads.
6+Simply call Rnd32.URand, URandf, Srand of Srandf to obtain a random value.
177
18-``` pascal
19-{ new version 2017-02-14 }
8+If you have data that relies on mt19937 and you have a FPC version > 3.2.2 you can use this instead of the new PRNG
9+```
2010unit 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.
11+ {$mode objfpc}{$A8}{$R-}{$Q-}
12+ interface
13+ type
14+ TMersenne32 = class sealed
15+ strict private
16+ const
17+ // Define MT19937 constants (32-bit RNG)
18+ N = 624;M = 397;R = 31;A = $9908B0DF;F = 1812433253;
19+ U = 11;S = 7;B = $9D2C5680;T = 15;C = $EFC60000;L = 18;
20+ MASK_LOWER = (QWORD(1) << R) - 1;
21+ MASK_UPPER = QWORD(1) << R;
22+ class var mt:array[0..N-1] of dword;
23+ class var index:word;
24+ class procedure twist;inline;static;
25+ public
26+ class constructor create;
27+ class procedure initialize(const seed:dword);inline;static;
28+ { 32 bit unsigned integer }
29+ class function URand:dword;inline;static;
30+ { 32 bit float in range 0..1 }
31+ class function URandf:single;inline;static;
32+ { 32 bit signed integer }
33+ class function SRand:integer;inline;static;
34+ {32 bit float in the range -1..1 }
35+ class function SRandf:single;inline;static;
36+ end;
37+
38+ Rnd32 = type TMersenne32;
39+
40+ implementation
41+
42+ class constructor TMersenne32.Create;
43+ begin
44+ initialize(5489); // default seed for mt19937
45+ end;
46+
47+ class procedure TMersenne32.Initialize(const seed:dword);inline;static;
48+ var
49+ i:dword;
50+ begin
51+ mt[0] := seed;
52+ for i := 1 to pred(N) do
53+ mt[i] := F * (mt[i - 1] xor (mt[i - 1] >> 30)) + i;
54+ index := N;
55+ end;
56+
57+ class procedure TMersenne32.Twist;inline;static;
58+ var
59+ i:integer;
60+ begin
61+ for i:=0 to N-M-1 do
62+ mt[i]:=mt[i+M] xor {twist} (((mt[i] and MASK_UPPER) or
63+ (mt[i+1] and MASK_LOWER)) shr 1)xor(dword(-(mt[i+1] and 1)) and A);
64+ for i:=N-M to N-2 do
65+ mt[i]:=mt[i+(M-N)]xor{twist}(((mt[i] and MASK_UPPER) or
66+ (mt[i+1] and MASK_LOWER)) shr 1)xor(dword(-(mt[i+1] and 1)) and A);
67+ mt[N-1]:=mt[M-1] xor {twist} (((mt[n-1] and MASK_UPPER) or (mt[0] and
68+ MASK_LOWER)) shr 1)xor(dword(-(mt[0] and 1)) and A);
69+ index:=0;
70+ end;
71+
72+
73+ class function TMersenne32.URand:dword;inline;static;
74+ var
75+ i:integer;
76+ begin
77+ i := index;
78+ if index >= N then
79+ begin
80+ Twist;
81+ i := index;
82+ end;
83+ Result := mt[i];index := i + 1;
84+ Result := Result xor (mt[i] >> U);
85+ Result := Result xor (Result << S) and B;
86+ Result := Result xor (Result << T) and C;
87+ Result := Result xor(Result >> L);
88+ end;
89+
90+ class function TMersenne32.URandf:single;inline;static;
91+ begin
92+ Result := URand * 2.32830643653869628906e-10;
93+ end;
94+
95+ class function TMersenne32.SRandf:single;inline;static;
96+ begin
97+ Result :=URand * 4.6566128730773926e-010;
98+ end;
99+
100+ class function TMersenne32.SRand:integer;inline;static;
101+ begin
102+ {$ifopt D+}Result := 0;{$endif}
103+ Result := URand;
104+ end;
105+ end.
115106```
116107
117108### C++ verification test:
118109
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
110+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:
111+Compile this with g++ mt_cpp_float32.cpp -std=gnu++11 -o mt_cpp_float32
112+```
123113#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-}
114+ #include <random>
115+ using namespace std;
116+
117+ int main()
118+ {
119+ mt19937 mt_rand(5489ul);
120+ for(int i=0;i<5;i++){
121+ cout << (float) mt_rand() * 2.32830643653869628906e-10 << endl;};
122+
123+ return 0;
124+ }
135125```
136126
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
127+Compile this with g++ mt_cpp_uint32.cpp -std=gnu++11 -o mt_cpp_uint32
128+Run ./mersplusplus and compare the output with TMersenne32. This code is for 32 bit unsigned. It should be (is) equal.
129+```
141130#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-}
131+ #include <random>
132+ using namespace std;
133+ int main()
134+ {
135+ mt19937 mt_rand(5489u);
136+ for(int i=0;i<5;i++){
137+ cout << mt_rand() << endl;};
138+
139+ return 0;
140+ }
152141```
153142
154-Left FPC, right C++\
155-![](/wiki-images/6/68/Selection_019.png)
143+Left FPC, right C++
144+[![Selection 019.png](/wiki-images/6/68/Selection_019.png)](</File:Selection_019.png>)
156145
157-### Verification against FPC's Random
146+### Verification against FPC's Random prior to 3.3.1
158147
159-``` pascal
148+This test is only valid for fpc versions prior to 3.3.1
149+```
160150program 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.
151+ {$mode objfpc}
152+ uses sysutils,dateutils,mersenne;
153+ var
154+ i:integer;
155+ begin
156+ Randseed :=5489; // default seed for mt19937
157+ for i := 0 to 9999999 do Rnd32.urandf;write(Rnd32.urandf :2:10,' ');writeln('Mersenne32');
158+ for i := 0 to 9999999 do Random ;write(single(Random) :2:10,' ');writeln('FPC Random');
159+ end.
170160```
171161
172162## Reference
173163
174-1. [Mersenne twister in WikiPedia](https://en.wikipedia.org/wiki/Mersenne_Twister)
164+ 1. [Mersenne twister in WikiPedia](<https://en.wikipedia.org/wiki/Mersenne_Twister>)
175165
176166## See also
177167
178-- [[Dev random]]
179-- [[Functions for descriptive statistics]]
180-- [[Generating Random Numbers]]
181-- [[Marsaglia&|Marsaglia's pseudo random number generators]]
168+ * [[Dev random]]
169+ * [[Functions for descriptive statistics]]
170+ * [[Generating Random Numbers]]
171+ * [[Marsaglia's pseudo random number generators]]
182172
183173---
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
174+*From the Free Pascal wiki, page "A simple implementation of the Mersenne twister", as the Internet Archive kept it on 2026-02-16, last changed on 2024-12-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