unit UnitZufall;
(*******************************************************************************)
(** **)
(** TRNGenerator: **)
(** Generiert Zufallszahlen nach dem Mersenne Twister. **)
(** **)
(*******************************************************************************)
{$R-} {range checking off}
{$Q-} {overflow checking off}
interface
const rtsRNGN = 624;
type
TRNGStateArray =
array [0..rtsRNGN-1]
of Integer;
// the array for the state vector
type TRNGenerator =
class
private
mt : TRNGStateArray;
mti: Integer;
// Generiert eine Zufallszahl von Low(Integer) bis High(Integer)
function GetFullRangeRandom: Integer;
public
constructor Create;
procedure SetSeed(Seed: Integer);
overload;
procedure SetSeed(
const SeedArray: TRNGStateArray);
overload;
procedure Randomize;
// Liefert eine Zufallszahl im Bereich von 0 bis +Range zurück:
function GetRandom(Range: Integer): Integer;
overload;
// Liefert eine Zufallszahl im Bereich von 0 bis 1 zurück:
function GetRandom: Double;
overload;
// Liefert eine Zufallszahl im Bereich von rFrom bis rTo zurück:
function GetRandomRange(rFrom, rTo: Integer): Integer;
// Liefert eine Gauss-Verteilte Zufallszahl zurück:
function Gauss(Mean, Variance: Extended): Extended;
end;
implementation
const
{ Period parameters }
rtsRNGM=397;
rtsRNGMATRIX_A =$9908b0df;
rtsRNGUPPER_MASK=$80000000;
rtsRNGLOWER_MASK=$7fffffff;
{ Tempering parameters }
TEMPERING_MASK_B=$9d2c5680;
TEMPERING_MASK_C=$efc60000;
constructor TRNGenerator.Create;
begin
mti := rtsRNGN+1;
inherited;
end;
procedure TRNGenerator.SetSeed(Seed: Integer);
var i: Integer;
begin
for i := 0
to rtsRNGN-1
do begin
mt[i]:= seed
and $ffff0000;
seed := 69069 * seed + 1;
mt[i]:= mt[i]
or ((seed
and $ffff0000)
shr 16);
seed := 69069 * seed + 1;
end;
mti := rtsRNGN;
end;
procedure TRNGenerator.SetSeed(
const SeedArray: TRNGStateArray);
var i: Integer;
begin
for i := 0
to rtsRNGN-1
do
mt[i] := SeedArray[i];
mti := rtsRNGN;
end;
function TRNGenerator.GetFullRangeRandom: Integer;
const mag01 :
array [0..1]
of Integer =(0, rtsRNGMATRIX_A);
var y, kk: Integer;
begin
if mti >= rtsRNGN
then begin
if mti = (rtsRNGN+1)
then SetSeed(4357);
for kk := 0
to rtsRNGN-rtsRNGM-1
do begin
y := (mt[kk]
and rtsRNGUPPER_MASK)
or (mt[kk+1]
and rtsRNGLOWER_MASK);
mt[kk] := mt[kk+rtsRNGM]
xor (y
shr 1)
xor mag01[y
and $00000001];
end;
for kk:= rtsRNGN-rtsRNGM
to rtsRNGN-2
do begin
y := (mt[kk]
and rtsRNGUPPER_MASK)
or (mt[kk+1]
and rtsRNGLOWER_MASK);
mt[kk] := mt[kk+(rtsRNGM-rtsRNGN)]
xor (y
shr 1)
xor mag01[y
and $00000001];
end;
y := (mt[rtsRNGN-1]
and rtsRNGUPPER_MASK)
or (mt[0]
and rtsRNGLOWER_MASK);
mt[rtsRNGN-1] := mt[rtsRNGM-1]
xor (y
shr 1)
xor mag01[y
and $00000001];
mti := 0;
end;
y := mt[mti]; inc(mti);
y := y
xor (y
shr 11);
y := y
xor (y
shl 7)
and TEMPERING_MASK_B;
y := y
xor (y
shl 15)
and TEMPERING_MASK_C;
y := y
xor (y
shr 18);
result := y;
end;
procedure TRNGenerator.Randomize;
var alt: Integer;
begin
alt := System.Randseed;
System.Randomize;
SetSeed(System.RandSeed);
System.RandSeed := alt;
end;
function TRNGenerator.GetRandom(Range: Integer): Integer;
begin
// 0 <= result < Range
result := Trunc(Range * GetRandom);
end;
function TRNGenerator.GetRandom: Double;
begin
result := (Abs(GetFullRangeRandom / High(Integer)))
end;
function TRNGenerator.GetRandomRange(rFrom, rTo: Integer): Integer;
begin
// rFrom <= result <= rTo
if rFrom > rTo
then result := GetRandom(rFrom - rTo + 1) + rTo
else result := GetRandom(rTo - rFrom + 1) + rFrom;
end;
function TRNGenerator.Gauss(Mean, Variance: Extended): Extended;
{ Marsaglia-Bray algorithm }
var a, b: Extended;
begin
repeat
a := 2 * GetRandom - 1;
b := Sqr(a) + Sqr(2*GetRandom-1);
until b < 1;
result := Sqrt(-2*Ln(b)/b) * a * Variance + Mean;
end;
end.