unit randomunit;

{$mode objfpc}{$H+}

interface

uses
  Classes, SysUtils;

type
  tMersenneTwister = class
  private
    state: array[0..623] of longword;
    index: longint;
  public
    procedure init(seed: longword);
    function extract_number: longword;
    function random(ma: longword): longword; overload;
    function random: extended; overload;
    procedure generate_numbers;
  end;

implementation

// tMersenneTwister ************************************************************

procedure tMersenneTwister.init(seed: longword);
var
  i: longint;
begin
 index:=0;
 state[0]:=seed;
 for i:=1 to 623 do
   state[i]:=longword($ffffffff and (qword(1812433253) * qword(state[i-1] xor state[i-1] shr 30) + i)); // 0x6c078965
end;

function tMersenneTwister.extract_number: longword;
begin
 if index=0 then
   generate_numbers;

 result:=state[index];
 result:=result xor (result shr 11);
 result:=result xor longword(qword(result shl 7) and 2636928640);  // 0x9d2c5680
 result:=result xor longword(qword(result shl 15) and 4022730752); // 0xefc60000
 result:=result xor (result shr 18);

 index := (index + 1) mod 624;
end;

function tMersenneTwister.random(ma: longword): longword;
var
  i: longword;
begin
 repeat
   i:=extract_number;
 until i<((high(longword)+1) div qword(ma))*ma;
 result:=i mod ma;
end;

function tMersenneTwister.random: extended;
begin
 result:=(extract_number/(high(longword)+1) + extract_number)/(high(longword)+1);
end;

procedure tMersenneTwister.generate_numbers;
var
  i,y: longint;
begin
 for i:=0 to 623 do begin
   y:=longint((state[i] and $80000000) or            // bit 31 (32nd bit) of MT[i]
              (state[(i+1) mod 624] and $7fffffff)); // bits 0-30 (first 31 bits) of MT[...]
   state[i]:=state[(i + 397) mod 624] xor (y shr 1);
   if odd(y) then
     state[i]:=longword(state[i] xor 2567483615); // 0x9908b0df
 end;
end;

end.