ruberon/lib/Math/RandExt.ob07

(* ************************************************************
   Дополнительные алгоритмы генераторов какбыслучайных чисел.
   Вадим Исаев, 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.