forked from KolibriOS/kolibrios
299 lines
8.9 KiB
Plaintext
299 lines
8.9 KiB
Plaintext
|
(* ************************************************************
|
|||
|
|
Дополнительные алгоритмы генераторов какбыслучайных чисел.
|
|||
|
|
Вадим Исаев, 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.
|