(* ************************************************************
Дополнительные алгоритмы генераторов какбыслучайных чисел.
Вадим Исаев, 2020
Additional generators of pseudorandom numbers.
Vadim Isaev, 2020
************************************************************ *)
MODULE RandExt;
IMPORT HOST, MathRound, MathBits;
CONST
(* Для алгоритма Мерсена-Твистера *)
N = 624;
M = 397;
MATRIX_A = 9908B0DFH; (* constant vector a *)
UPPER_MASK = 80000000H; (* most significant w-r bits *)
LOWER_MASK = 7FFFFFFFH; (* least significant r bits *)
INT_MAX = 4294967295;
TYPE
(* структура служебных данных, для алгоритма mrg32k3a *)
random_t = RECORD
mrg32k3a_seed : REAL;
mrg32k3a_x : ARRAY 3 OF REAL;
mrg32k3a_y : ARRAY 3 OF REAL
END;
(* Для алгоритма Мерсена-Твистера *)
MTKeyArray = ARRAY N OF INTEGER;
VAR
(* Для алгоритма mrg32k3a *)
prndl: random_t;
(* Для алгоритма Мерсена-Твистера *)
mt : MTKeyArray; (* the array for the state vector *)
mti : INTEGER; (* mti == N+1 means mt[N] is not initialized *)
(* ---------------------------------------------------------------------------
Генератор какбыслучайных чисел в диапазоне [a,b].
Алгоритм 133б из книги "Агеев и др. - Бибилотека алгоритмов 101б-150б",
стр. 53.
Переделка из Algol на Oberon и доработка, Вадим Исаев, 2020
Generator pseudorandom numbers, algorithm 133b from
Comm ACM 5,10 (Oct 1962) 553.
Convert from Algol to Oberon Vadim Isaev, 2020.
Входные параметры:
a - начальное вычисляемое значение, тип REAL;
b - конечное вычисляемое значение, тип REAL;
seed - начальное значение для генерации случайного числа.
Должно быть в диапазоне от 10 000 000 000 до 34 359 738 368 (2^35),
нечётное.
--------------------------------------------------------------------------- *)
PROCEDURE alg133b* (a, b: REAL; VAR seed: INTEGER): REAL;
CONST
m35 = 34359738368;
m36 = 68719476736;
m37 = 137438953472;
VAR
x: INTEGER;
BEGIN
IF seed # 0 THEN
IF (seed MOD 2 = 0) THEN
seed := seed + 1
END;
x:=seed;
seed:=0;
END;
x:=5*x;
IF x>=m37 THEN
x:=x-m37
END;
IF x>=m36 THEN
x:=x-m36
END;
IF x>=m35 THEN
x:=x-m35
END;
RETURN FLT(x) / FLT(m35) * (b - a) + a
END alg133b;
(* ----------------------------------------------------------
Генератор почти равномерно распределённых
какбыслучайных чисел mrg32k3a
(Combined Multiple Recursive Generator) от 0 до 1.
Период повторения последовательности = 2^127
Generator pseudorandom numbers,
algorithm mrg32k3a.
Переделка из FreePascal на Oberon, Вадим Исаев, 2020
Convert from FreePascal to Oberon, Vadim Isaev, 2020
---------------------------------------------------------- *)
(* Инициализация генератора.
Входные параметры:
seed - значение для инициализации. Любое. Если передать
ноль, то вместо ноля будет подставлено кол-во
процессорных тиков. *)
PROCEDURE mrg32k3a_init* (seed: REAL);
BEGIN
prndl.mrg32k3a_x[0] := 1.0;
prndl.mrg32k3a_x[1] := 1.0;
prndl.mrg32k3a_y[0] := 1.0;
prndl.mrg32k3a_y[1] := 1.0;
prndl.mrg32k3a_y[2] := 1.0;
IF seed # 0.0 THEN
prndl.mrg32k3a_x[2] := seed;
ELSE
prndl.mrg32k3a_x[2] := FLT(HOST.GetTickCount());
END;
END mrg32k3a_init;
(* Генератор какбыслучайных чисел от 0.0 до 1.0. *)
PROCEDURE mrg32k3a* (): REAL;
CONST
(* random MRG32K3A algorithm constants *)
MRG32K3A_NORM = 2.328306549295728E-10;
MRG32K3A_M1 = 4294967087.0;
MRG32K3A_M2 = 4294944443.0;
MRG32K3A_A12 = 1403580.0;
MRG32K3A_A13 = 810728.0;
MRG32K3A_A21 = 527612.0;
MRG32K3A_A23 = 1370589.0;
RAND_BUFSIZE = 512;
VAR
xn, yn, result: REAL;
BEGIN
(* Часть 1 *)
xn := MRG32K3A_A12 * prndl.mrg32k3a_x[1] - MRG32K3A_A13 * prndl.mrg32k3a_x[2];
xn := xn - MathRound.trunc(xn / MRG32K3A_M1) * MRG32K3A_M1;
IF xn < 0.0 THEN
xn := xn + MRG32K3A_M1;
END;
prndl.mrg32k3a_x[2] := prndl.mrg32k3a_x[1];
prndl.mrg32k3a_x[1] := prndl.mrg32k3a_x[0];
prndl.mrg32k3a_x[0] := xn;
(* Часть 2 *)
yn := MRG32K3A_A21 * prndl.mrg32k3a_y[0] - MRG32K3A_A23 * prndl.mrg32k3a_y[2];
yn := yn - MathRound.trunc(yn / MRG32K3A_M2) * MRG32K3A_M2;
IF yn < 0.0 THEN
yn := yn + MRG32K3A_M2;
END;
prndl.mrg32k3a_y[2] := prndl.mrg32k3a_y[1];
prndl.mrg32k3a_y[1] := prndl.mrg32k3a_y[0];
prndl.mrg32k3a_y[0] := yn;
(* Смешение частей *)
IF xn <= yn THEN
result := ((xn - yn + MRG32K3A_M1) * MRG32K3A_NORM)
ELSE
result := (xn - yn) * MRG32K3A_NORM;
END;
RETURN result
END mrg32k3a;
(* -------------------------------------------------------------------
Генератор какбыслучайных чисел, алгоритм Мерсена-Твистера (MT19937).
Переделка из Delphi в Oberon Вадим Исаев, 2020.
Mersenne Twister Random Number Generator.
A C-program for MT19937, with initialization improved 2002/1/26.
Coded by Takuji Nishimura and Makoto Matsumoto.
Adapted for DMath by Jean Debord - Feb. 2007
Adapted for Oberon-07 by Vadim Isaev - May 2020
------------------------------------------------------------ *)
(* Initializes MT generator with a seed *)
PROCEDURE InitMT(Seed : INTEGER);
VAR
i : INTEGER;
BEGIN
mt[0] := MathBits.iand(Seed, INT_MAX);
FOR i := 1 TO N-1 DO
mt[i] := (1812433253 * MathBits.ixor(mt[i-1], LSR(mt[i-1], 30)) + i);
(* See Knuth TAOCP Vol2. 3rd Ed. P.106 For multiplier.
In the previous versions, MSBs of the seed affect
only MSBs of the array mt[].
2002/01/09 modified by Makoto Matsumoto *)
mt[i] := MathBits.iand(mt[i], INT_MAX);
(* For >32 Bit machines *)
END;
mti := N;
END InitMT;
(* Initialize MT generator with an array InitKey[0..(KeyLength - 1)] *)
PROCEDURE InitMTbyArray(InitKey : MTKeyArray; KeyLength : INTEGER);
VAR
i, j, k, k1 : INTEGER;
BEGIN
InitMT(19650218);
i := 1;
j := 0;
IF N > KeyLength THEN
k1 := N
ELSE
k1 := KeyLength;
END;
FOR k := k1 TO 1 BY -1 DO
(* non linear *)
mt[i] := MathBits.ixor(mt[i], (MathBits.ixor(mt[i-1], LSR(mt[i-1], 30)) * 1664525)) + InitKey[j] + j;
mt[i] := MathBits.iand(mt[i], INT_MAX); (* for WORDSIZE > 32 machines *)
INC(i);
INC(j);
IF i >= N THEN
mt[0] := mt[N-1];
i := 1;
END;
IF j >= KeyLength THEN
j := 0;
END;
END;
FOR k := N-1 TO 1 BY -1 DO
(* non linear *)
mt[i] := MathBits.ixor(mt[i], (MathBits.ixor(mt[i-1], LSR(mt[i-1], 30)) * 1566083941)) - i;
mt[i] := MathBits.iand(mt[i], INT_MAX); (* for WORDSIZE > 32 machines *)
INC(i);
IF i >= N THEN
mt[0] := mt[N-1];
i := 1;
END;
END;
mt[0] := UPPER_MASK; (* MSB is 1; assuring non-zero initial array *)
END InitMTbyArray;
(* Generates a integer Random number on [-2^31 .. 2^31 - 1] interval *)
PROCEDURE IRanMT(): INTEGER;
VAR
mag01 : ARRAY 2 OF INTEGER;
y,k : INTEGER;
BEGIN
IF mti >= N THEN (* generate N words at one Time *)
(* If IRanMT() has not been called, a default initial seed is used *)
IF mti = N + 1 THEN
InitMT(5489);
END;
FOR k := 0 TO (N-M)-1 DO
y := MathBits.ior(MathBits.iand(mt[k], UPPER_MASK), MathBits.iand(mt[k+1], LOWER_MASK));
mt[k] := MathBits.ixor(MathBits.ixor(mt[k+M], LSR(y, 1)), mag01[MathBits.iand(y, 1H)]);
END;
FOR k := (N-M) TO (N-2) DO
y := MathBits.ior(MathBits.iand(mt[k], UPPER_MASK), MathBits.iand(mt[k+1], LOWER_MASK));
mt[k] := MathBits.ixor(mt[k - (N - M)], MathBits.ixor(LSR(y, 1), mag01[MathBits.iand(y, 1H)]));
END;
y := MathBits.ior(MathBits.iand(mt[N-1], UPPER_MASK), MathBits.iand(mt[0], LOWER_MASK));
mt[N-1] := MathBits.ixor(mt[M-1], MathBits.ixor(LSR(y, 1), mag01[MathBits.iand(y, 1H)]));
mti := 0;
END;
y := mt[mti];
INC(mti);
(* Tempering *)
y := MathBits.ixor(y, LSR(y, 11));
y := MathBits.ixor(y, MathBits.iand(LSL(y, 7), 9D2C5680H));
y := MathBits.ixor(y, MathBits.iand(LSL(y, 15), 4022730752));
y := MathBits.ixor(y, LSR(y, 18));
RETURN y
END IRanMT;
(* Generates a real Random number on [0..1] interval *)
PROCEDURE RRanMT(): REAL;
BEGIN
RETURN FLT(IRanMT())/FLT(INT_MAX)
END RRanMT;
END RandExt.