Difference between revisions of "A simple implementation of the Mersenne twister"

From Lazarus wiki
Jump to navigationJump to search
(Linking to documentation of random function.)
(Adding "See also" links)
Line 87: Line 87:
 
===Reference ===
 
===Reference ===
 
#[https://en.wikipedia.org/wiki/Mersenne_Twister Mersenne twister in WikiPedia]
 
#[https://en.wikipedia.org/wiki/Mersenne_Twister Mersenne twister in WikiPedia]
 +
 +
== See also ==
 +
* [[Dev random]]
 +
* [[Functions for descriptive statistics]]
 +
* [[Generating Random Numbers]]
 +
* [[Marsaglia's pseudo random number generators]]
  
 
[[Category:Statistical algorithms]]
 
[[Category:Statistical algorithms]]
 
[[Category:Modelling and Simulation]]
 
[[Category:Modelling and Simulation]]
 
[[Category: Code]]
 
[[Category: Code]]

Revision as of 01:36, 12 March 2017

When you look at the sourcecode that implements Free Pascal's random it may look overwelmingly complex. The algorithm itself, though, is very understandable. Here is simple 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. Simple call TMersenne.ExtractU32 to obtain a random value.


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;

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.

Reference

  1. Mersenne twister in WikiPedia

See also