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.
|