Difference between revisions of "A simple implementation of the Mersenne twister"
m |
|||
Line 144: | Line 144: | ||
begin | begin | ||
Randseed :=5489; // default seed for mt199937 | Randseed :=5489; // default seed for mt199937 | ||
− | start := Now;for i := 0 to 9999999 do | + | start := Now;for i := 0 to 9999999 do Rnd32.urandf;write(Rnd32.urandf :2:10,' ');writeln('Mersenne32');m:=MilliSecondsBetween(Now,start); |
for i := 0 to 9999999 do Random ;write(single(Random) :2:10,' ');writeln('FPC Random');f:=MilliSecondsBetween(Now,start); | for i := 0 to 9999999 do Random ;write(single(Random) :2:10,' ');writeln('FPC Random');f:=MilliSecondsBetween(Now,start); | ||
writeln('Mersenne32: ',m,' FPC Random ',f, ', speed difference (f *100)/m: ', 'Mersenne32 is ',f*100/m -100 :4:0 ,'% Faster than FPC's Random'); | writeln('Mersenne32: ',m,' FPC Random ',f, ', speed difference (f *100)/m: ', 'Mersenne32 is ',f*100/m -100 :4:0 ,'% Faster than FPC's Random'); |
Revision as of 08:32, 13 March 2017
TMersenne32
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 a simple 32 bit implementation according to MT19937. It is compatible with FPC's Random for 32 bit values. It is also about 45% faster than FPC's random.
Enjoy,
Thaddy.
It is verified against several mt19937 implementations in different languages, including Freepascal and C/C++.
Of course I added this one to wikipedia too. You don't have to create this class. it is fully static.
Simply call Rnd32.URand, URandf, Srand of Srandf to obtain a random value.
{ new version 2017-02-12 }
unit mersenne;
{$mode objfpc}
interface
type
TMersenne32 = 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;
{ 32 bit unsigned integer }
class function URand:dword;inline;static;
{ 32 bit float in range 0..1 }
class function URandf:single;inline;static;
{ 32 bit signed integer }
class function SRand:integer;inline;static;
{32 bit float in the range -1..1 }
class function SRandf:single;inline;static;
end;
Rnd32 = type TMersenne32;
implementation
class constructor TMersenne32.Create;
begin
initialize(5489); // default seed for mt199937
end;
class procedure TMersenne32.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 TMersenne32.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 TMersenne32.URand: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;
class function TMersenne32.URandf:single;inline;static;
begin
Result := URand * 2.32830643653869628906e-10;
end;
class function TMersenne32.SRandf:single;inline;static;
begin
Result :=URand * 4.6566128730773926e-010;
end;
class function TMersenne32.SRand:integer;inline;static;
begin
{$ifopt D+}Result := 0;{$endif}
Result := URand;
end;
end.
C++ verification test:
Compile this with g++ mersplusplus.cpp -std=gnu++11 -o mersplusplus
Run ./mersplusplus and compare the output with TMersenne32. It should be (is) equal.
#include <iostream>
#include <random>
using namespace std;
int main()
{
mt19937 mt_rand(5489u);
for(int i=0;i<5;i++){
cout << mt_rand() << endl;};
return 0;
}
Verification against FPC's Random
program mersennetest2;
{$mode objfpc}
uses sysutils,dateutils,mersenne;
var
i:integer;
start:TdateTime;
m,f:int64;
begin
Randseed :=5489; // default seed for mt199937
start := Now;for i := 0 to 9999999 do Rnd32.urandf;write(Rnd32.urandf :2:10,' ');writeln('Mersenne32');m:=MilliSecondsBetween(Now,start);
for i := 0 to 9999999 do Random ;write(single(Random) :2:10,' ');writeln('FPC Random');f:=MilliSecondsBetween(Now,start);
writeln('Mersenne32: ',m,' FPC Random ',f, ', speed difference (f *100)/m: ', 'Mersenne32 is ',f*100/m -100 :4:0 ,'% Faster than FPC's Random');
end.