Another retro (16-bit) Pseudo Random Number Generators written in Pascal before. Do anything with it.

PRNGS.PAS

unit prngs;
{ Pseudo Random Number Generators written in Pascal }
{ http://wiki.freepascal.org/Generating_Random_Numbers }

interface

type
  algorithm = (rtl, lcg, pcg, xors, mt);
  dist = (uniform, normal, exponential, gamma, erlang, poisson);

{ Default uniform distribution random integer range in RTL }
{ Should be Linear congruential generator (LCG) algorithm }
function rnd(const min, max: word): word;

{ LCG algorithm used by Borland Pascal/Delphi and other old languages }
{ Should be consistent with RTL one }
function rndLcg(const min, max: word): word;
procedure randomizeLcg(const seed: longint);

{ Tiny Mersenne Twister algorithm, MT used by Free Pascal/Python/PHP7 }
function rndMt(const min, max: word): word;
procedure randomizeMt(const seed: longint);

{ Xorshift algorithm 16-bit variant, its improvements used by CUDA/V8JS }
function rndXors(const min, max: word): word;
procedure randomizeXors(const seed: longint);

{ Permuted congruential generator (PCG) algorithm 16-bit, improved LCG }
function rndPcg(const min, max: word): word;
procedure randomizePcg(const seed: longint);

{ Scale the random integer to float value ranged in [0, 1) }
function scaleToFloat(const value: word): real;

{ Normal(Gaussian) distribution random real according the Box-Muller approach }
function rndNorm(const mean, stddev: real): real;

{ Exponential distribution random real }
function rndExp(const a, rate: real): real;

{ Gamma distribution random real }
{function rndGamma(const a, b, c: real): real;}

{ Erlang distribution random real }
function rndErlang(const mean: real; const k: integer): real;

{ Poisson distribution random integer }
function rndPoisson(const mean: integer): integer;


{ Return a 32-bit system clock timestamp }
{ But not real timestamp since MSDOS service $2C return no millisecond }
function getDosTime: longint;

{ Return the timestamp when this unit was initialized }
function getInitTime: longint;


implementation

uses dos;

const
  WordMax = 65535;

{ States records for corresponding algorithm }
{ see http://wiki.freepascal.org/Marsaglia%27s_pseudo_random_number_generators }type
  pcgStates = record
    state, inc: longint;
  end;
  tinyMtStates = record
    status: array[0..3] of word;
    mat1, mat2, tmat: word;
  end;

var
  initTime: longint;
  lcgSeed: longint;
  xors16: word;
  pcgs: pcgStates;
  tmts: tinyMtStates;

function rnd;
begin
  rnd := min + random(max - min);
end;

function getDosTime;
var
  regs: registers;
begin
  regs.ah := $2C;
  msdos(regs);
  { CH: hours, CL: minutes; DH: seconds, DL: hundredthsecs }
  { Some system cannot offer even accurate hudredthsecs (0 returned) }
  getDosTime := (longint(regs.cx) shl 16) + regs.dx;
end;

{ Swap timestamp high 16-bit with low 16-bit for some reason }
function timeToRandSeed(const dosTime: longint): longint;
begin
  timeToRandSeed := (dosTime shl 16) + (dosTime shr 16);
end;

function getInitTime;
begin
  getInitTime := initTime;
end;

function scaleToFloat;
begin
  scaleToFloat := value / WordMax;
end;

procedure randomizeLcg;
begin
  lcgSeed := seed;
end;

function rndLcgNext: longint;
begin
  { Multiply factor in RTLSYS\RAND.ASM: $08088405 }
  { Parameters used by other implementations:
      https://en.wikipedia.org/wiki/Linear_congruential_generator }
  lcgSeed := lcgSeed * $8088405 + 1;
  rndLcgNext := lcgSeed;
end;

function rndLcg;
var
  next: longint;
  tmp1, tmp2, range: word;
begin
  { To keep high 16-bit output for DWORD nextrand * range shr 16, }
  { To get overflowed bits without 64-bit QWORD type support }
  next := rndLcgNext;
  range := max - min;
  tmp1 := next shr 16;
  next := longint(next and $FFFF) * range;
  tmp2 := next shr 16;
  next := longint(tmp1) * range + tmp2;
  rndLcg := min + word(next shr 16);
end;

procedure randomizeXors;
begin
  { Initial state must not be 0 }
  xors16 := seed mod WordMax;
  repeat
    inc(xors16);
  until xors16 > 0;
end;

{ Orginal paper assumed 32-bit output from XSorShift Prng's, 2003 }
{ This simple 16-bit one from: http://www.retroprogramming.com/2017/07/xorshift-pseudorandom-numbers-in-z80.html }
function rndXorsNext: word;
begin
  xors16 := xors16 xor (xors16 shl 7);
  xors16 := xors16 xor (xors16 shr 9);
  xors16 := xors16 xor (xors16 shl 8);
  rndXorsNext := xors16;
end;

function rndXors;
var
  next, range, threshold: word;
begin
  { To get limited output from scaled float, but this quite slow }
{
  rndXors := min + trunc(scaleToFloat(rndXorsNext) * (max - min));
}
  { Borrowed from PCG bound output method }
  range := max - min;
  threshold := word(-range) mod range;
  repeat
    next := rndXorsNext;
  until next >= threshold;
  rndXors := min + (next mod range);
end;

{ Simple oneseq PCG next step, so inc unused }
{ See https://github.com/imneme/pcg-c }
{ See http://www.pcg-random.org }
procedure pcgOneseqStep;
const
  multi32 = $2C9277B5;
  inc32 = $AC564B05;
begin
  { pcg_oneseq_32_step_r }
  pcgs.state := pcgs.state * multi32 + inc32;
end;

procedure randomizePcg;
begin
  with pcgs do
  begin
    state := 0;
    inc := 1;
    pcgOneseqStep;
    state := state + seed;
    pcgOneseqStep;
  end;
end;

function rndPcgNext: word;
var
  oldstate: longint;
  xorshifted, rot: word;
begin
  oldstate := pcgs.state;
  pcgOneseqStep;
  { pcg_output_xsh_rr_32_16 }
  xorshifted := ((oldstate shr 10) xor oldstate) shr 12;
  { pcg_rotr_16 }
  rot := oldstate shr 28;
  rndPcgNext := (xorshifted shr rot) or (xorshifted shl (word(-rot) and 15));
end;

function rndPcg;
var
  bound, threshold, next: word;
begin
  bound := max - min;
  threshold := word(-bound) mod bound;
  repeat
    next := rndPcgNext;
  until next >= threshold;
  rndPcg := min + (next mod bound);
end;

{ TinyMT 32-bit see http://github.com/MersenneTwister-Lab/TinyMT }
procedure randomizeMt;
var
  i: word;
begin
  with tmts do
  begin
    { 32-bit initial values from check32.out.txt }
    mat1 := word($8F7011EE mod WordMax);
    mat2 := word($FC78FF1F mod WordMax);
    tmat := word($3793FDFF mod WordMax);
    status[0] := seed mod WordMax;
    status[1] := mat1;
    status[2] := mat2;
    status[3] := tmat;
    for i := 1 to 8 do
      status[i and 3] := status[i and 3] xor
        (i + ($6C078965 mod WordMax) *
          (status[(i - 1) and 3] xor
          (status[(i - 1) and 3] shr 30)));
  end;
end;

function rndMtNext: word;
var
  x, y: word;
  t0, t1: word;
begin
  with tmts do
  begin
    { next state, still using 32-bit shift values }
    y := status[3];
    x := (status[0] and $7FFF) xor status[1] xor status[2];
    x := x xor (x shl 1);
    y := y xor ((y shr 1) xor x);
    status[0] := status[1];
    status[1] := status[2];
    status[2] := x xor (y shl 10);
    status[3] := y;
    status[1] := status[1] xor (word(-(y and 1)) and mat1);
    status[2] := status[2] xor (word(-(y and 1)) and mat2);
    { temper, internal states are already 16-bit }
    t0 := status[3];
    t1 := status[0] + (status[2] shr 8);
    t0 := t0 xor t1;
    t0 := t0 xor (word(-(t1 and 1)) and tmat);
    rndMtNext := t0;
  end;
end;

function rndMt;
var
  next, range, threshold: word;
begin
  range := max - min;
  threshold := word(-range) mod range;
  repeat
    next := rndMtNext;
  until next >= threshold;
  rndMt := min + (next mod range);
end;


{ Simulating other distributions with RTL random goes here }

function rndNorm;
var
  u1, u2: real;
begin
  u1 := random;
  u2 := random;
  rndNorm := mean * abs(1 + sqrt(-2 * (ln(u1))) * cos(2 * pi * u2) * stddev);
end;

function rndExp;
const
  reso = 1000;
var
  unif: real;
begin
  if rate = 0 then
    rndExp := -1
  else begin
    repeat
      unif := random(reso) / reso;
    until unif <> 0;
    rndExp := a - rate * ln(unif);
  end;
end;

{
function rndGamma;
const
  reso = 1000;
  t = 4.5;
var
  d: real;
  unif: real;
  a2, b2, c2, q, p, y: real;
  p1, p2, v, w, z: real;
  found: boolean;
begin
  d := 1 + ln(t);
  a2 := 1 / sqrt(2 * c - 1);
  b2 := c - ln(4);
  q := c + 1 / a2;
  c2 := 1 + c / exp(1);
  found := false;
  if c < 1 then
  begin
    repeat
      repeat
        unif := random(reso) / reso;
      until unif > 0;
      p := c2 * unif;
      if p > 1 then
      begin
        repeat
          unif := random(reso) / reso;
        until unif > 0;
        y := -ln((c2 - p) / c);
        if unif <= power(y, c - 1) then
        begin
          rndGamma := a + b * y;
          found := true;
        end;
      end else
      begin
        y := power(p, 1 / c);
        if unif <= exp(-y) then
        begin
          rndGamma := a + b * y;
          found := true;
        end;
      end;
    until found;
  end else if c = 1 then
  begin
    rndGamma := rndExp(a, b);
  end else
  begin
    repeat
      repeat
        p1 := random(reso) / reso;
      until p1 > 0;
      repeat
        p2 := random(reso) / reso;
      until p2 > 0;
      v := a2 * ln(p1 / (1 - p1));
      y := c * exp(v);
      z := p1 * p1 * p2;
      w := b2 + q * v - y;
      if (w + d - t * z >= 0) or (w >= ln(z)) then
      begin
        rndGamma := a + b * y;
        found := true;
      end;
    until found;
  end;
end;
}

function rndErlang;
const
  reso = 1000;
var
  i: integer;
  unif, prod: real;
begin
  if (mean <= 0) or (k < 1) then
    rndErlang := -1
  else
  begin
    prod := 1;
    for i := 1 to k do
    begin
      repeat
        unif := random(reso) / reso;
      until unif <> 0;
      prod := prod * unif;
    end;
    rndErlang := -mean * ln(prod);
  end;
end;

function rndPoisson;
const
  reso = 1000;
var
  k: integer;
  b, l: real;
begin
  if mean <= 0 then
    rndPoisson := -1
  else
  begin
    k := 0;
    b := 1;
    l := exp(-mean);
    while b > l do
    begin
      inc(k);
      b := b * random(reso) / reso;
    end;
    rndPoisson := k - 1;
  end;
end;

begin
  initTime := getDosTime;
  { Use a swapped system timestamp as entropy }
  RandSeed := timeToRandSeed(initTime);
  randomizeLcg(timeToRandSeed(initTime));
  randomizeXors(timeToRandSeed(initTime));
  randomizePcg(timeToRandSeed(initTime));
  randomizeMt(timeToRandSeed(initTime));
end.