Difference between revisions of "A simple implementation of the Mersenne twister"
From Lazarus wiki
Jump to navigationJump to searchJwdietrich (talk | contribs) (Typo) |
m |
||
| Line 1: | Line 1: | ||
When you look at the source code that implements Free Pascal's <code>[[doc:rtl/system/random.html|random]]</code> it may look overwelmingly complex. | When you look at the source code that implements Free Pascal's <code>[[doc:rtl/system/random.html|random]]</code> it may look overwelmingly complex. | ||
| − | The algorithm itself, though, is very understandable. Here is simple implementation according to MT19937. | + | The algorithm itself, though, is very understandable. Here is simple 32 bit implementation according to MT19937. |
Enjoy,<br/> | Enjoy,<br/> | ||
| Line 6: | Line 6: | ||
| − | It is basically taken from the C version in wikipedia. Of course I added this one there too. You don't have to create this class. it is fully static. | + | It is basically taken from the C version in wikipedia. Of course I added this one there too. You don't have to create this class. it is fully static.<br> |
| + | Simply call Mer32.ExtractU32 to obtain a random value. Cast to integer for signed output. | ||
| Line 30: | Line 31: | ||
class function Extractu32:dword;inline;static; | class function Extractu32:dword;inline;static; | ||
end; | end; | ||
| + | |||
| + | mer32 = type TMersenne; | ||
implementation | implementation | ||
Revision as of 13:11, 12 March 2017
When you look at the source code that implements Free Pascal's random it may look overwelmingly complex.
The algorithm itself, though, is very understandable. Here is simple 32 bit implementation according to MT19937.
Enjoy,
Thaddy.
It is basically taken from the C version in wikipedia. Of course I added this one there too. You don't have to create this class. it is fully static.
Simply call Mer32.ExtractU32 to obtain a random value. Cast to integer for signed output.
unit mersenne;
{$mode objfpc}
interface
type
TMersenne = class sealed
strict private
const
// Define MT19937 constants (32-bit RNG)
N = 624;M = 397;R = 31;A = $9908B0DF;F = 1812433253;
U = 11;S = 7;B = $9D2C5680;T = 15;C = $EFC60000;L = 18;
MASK_LOWER = (QWORD(1) << R) - 1;
MASK_UPPER = QWORD(1) << R;
class var mt:array[0..N-1] of dword;
class var index:word;
class procedure twist;inline;static;
public
class constructor create;
class procedure initialize(const seed:dword);inline;static;
class function Extractu32:dword;inline;static;
end;
mer32 = type TMersenne;
implementation
class constructor TMersenne.Create;
begin
initialize(12345678);
end;
class procedure TMersenne.Initialize(const seed:dword);inline;static;
var
i:dword;
begin
mt[0] := seed;
for i := 1 to pred(N) do
mt[i] := F * (mt[i - 1] xor (mt[i - 1] >> 30)) + i;
index := N;
end;
class procedure TMersenne.Twist;inline;static;
var
i,x,xA:dword;
begin
for i := 0 to pred(N) do
begin
x := (mt[i] and MASK_UPPER) + (mt[(i + 1) mod N] and MASK_LOWER);
xA := x >> 1;
if ( x and $1 ) <> 0 then
xA := xA xor A;
mt[i] := mt[(i + M) mod N] xor xA;
end;
index := 0;
end;
class function TMersenne.ExtractU32:dword;inline;static;
var
y:dword;
i:integer;
begin
i := index;
if index >= N then
begin
Twist;
i := index;
end;
y := mt[i];
index := i + 1;
y := y xor (mt[i] >> U);
y := y xor (y << S) and B;
y := y xor (y << T) and C;
y := y xor(y >> L);
Result :=y;
end;
end.