{
#
   use Maurer's test to measure randomness of random number generators
   NB: a single test on its own is not enough to form a reliable judgement,
       but maurer's test is generally regarded as a good one.
   This version gives each RNG a score in the range 0 .. 100, where
   a score of 50 or more indicates a good RNG.
}

program rtest(output);

#include <clib.inc.pas> // include rand()

var seed: integer;


#if defined __unix__
{
/dev/urandom random numbers collected from hardware noise
to test if your environment generates hw random nrs, try one of these comamnds:
               xxd /dev/urandom
               hexdump /dev/urandom
}
var   rf : file of integer;
function drandom:integer;
var
   n  : integer;
begin
   read(rf, n);
   drandom := n;
end; { drandom }
#else
{
generate an initial seed from random contents of uninitialised memory
note that global memory is initialised to zero,
so we must use variables that are inside a function or procedure
}
function iseed: integer;
const max =  1024;
var
   a : array[1..max] of integer;
   i : integer;
   s : integer;
begin
   for i := 1 to max do begin
//      writeln('s is ', s, ', a[', i, '] is ', a[i]);
      s := bitxor(s, a[i]);
   end;
   iseed := s;
end; { iseed }
#endif


{
  random - return random real numbers uniformly distributed over (0,1)
  based on minimum standard random number generator from
  "Random Number Generators: Good Ones Are Hard to Find",
  S.K. Park and K.W. Miller, Communications of the ACM 31:10 (Oct 1988)
  Assumes integers are at least 32 bits.
}
function minstd : integer;
const
   { 16807, 48271 & 69621 are all good multipliers}
   a = 16807;         { multiplier = 7**5 }
   m = 2147483647;    { 2**31 - 1 }
   q = 127773;        { m div a }
   r = 2836;          { m mod a }
var
   lo, hi : integer; { need to be at least 32 bits }
begin
   hi := seed div q;
   lo := seed mod q;

   if a * lo > r * hi then
      seed := a * lo - r * hi
   else
      seed := (m - r * hi) + a * lo;

   minstd := seed;
end; { minstd }


{
 Wichmann–Hill random number generator
}
var
   s1, s2, s3 : integer; {initialise these to random numbers from 1 to 30,000}

function whrand: integer;
var
   r :  real;
/*
[r, s1, s2, s3] = function(s1, s2, s3)
    % s1, s2, s3 should be random from 1 to 30,000. Use clock if available
    s1 = mod(171 × s1, 30269)
    s2 = mod(172 × s2, 30307)
    s3 = mod(170 × s3, 30323)

    r = mod(s1/30269.0 + s2/30307.0 + s3/30323.0, 1)
*/
begin
   s1 := (171 * s1) mod 30269;
   s2 := (172 * s2) mod 30307;
   s3 := (170 * s3) mod 30323;

   r := s1/30269.0 + s2/30307.0 + s3/30323.0;
   r := r - trunc(r);

   whrand := round(maxint*r); {convert from 0..1 to 0..maxint}

end; { whrand() }


{
 ansi random nr generator
}
function arand: integer;
const  ARAND_MAX  =  32768;
begin
{$d-  -- need to allow overflow}
   seed := (seed*1103515245 + 12345) mod ARAND_MAX;
{$d+}
   arand := seed;
end;


function prand:integer;
begin
   seed := (seed*125) mod 2796203;
   prand := seed;
end; { prand }


function randhp35: integer;
{ random number 0-maxint }
{ adapted from HP-35 applications programs }
{seed should be initialised to 4}
const
 pi	= 3.1415926535897932384626433832795;

var
 x	: real;

begin
 x:=seed/maxint+pi;
 x:=exp(5.0*ln(x));
 x:=x-trunc(x);
 seed:=trunc(x*maxint);
 randhp35:=seed;
end; { randhp35 }


{ ******************************************************************
    Mersenne Twister Random Number Generator for pascal
  ******************************************************************
  From http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/VERSIONS/PASCAL/Alex.pas
}

const
   N	      = 624;
   N_1	      = 623;
   M	      = 397;
   MATRIX_A   = $9908b0df;  { constant vector a }
   UPPER_MASK = $80000000;  { most significant w-r bits }
   LOWER_MASK = $7fffffff;  { least significant r bits }

type
   Int32 = integer;
   Int16 = 0..65535;
   Cardinal = Int32;
   arrNLong = array [0..N] of Int32;

var
   mt	 : array[0..N_1] of Cardinal{Int32};  { the array for the state vector }
   mti	 : Int16;                        { mti == N+1 means mt[N] is not initialized }
   mag01 : array[0..1] of Cardinal{Int32};
   init	 : arrNLong;
   i	 : integer;

procedure InitMT(Seed : Int32);
var
  i : Int16;
begin
  mt[0] := bitand(Seed, $ffffffff);
  for i := 1 to N_1 do
    begin
      {$d-  --allow multiply overflow}
      mt[i] := (1812433253 * bitxor(mt[i-1], rshiftu(mt[i-1], 30)) + i);
      {$d+}
        { 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] := bitand(mt[i], $ffffffff);
        { For >32 Bit machines }
    end;
  mti := N;
end; { InitMT }


procedure InitMTbyArray(InitKey : arrNLong ; KeyLength : Int16);
var
  i, j, k, k1 : Int16;
begin
  InitMT(seed {19650218});

  i := 1;
  j := 0;

  if N > KeyLength then k1 := N else k1 := KeyLength;

  for k := k1 downto 1 do
  begin
      {$d-   --allow multiply overflow}
      mt[i] := bitxor(mt[i], (bitxor(mt[i-1], rshiftu(mt[i-1], 30)) * 1664525)) + InitKey[j] + j; { non linear }
//      mt[i] := mt[i] and $ffffffff; { for WORDSIZE > 32 machines }
      {$d+}
      i := i + 1;
      j := j + 1;
      if i >= N then
        begin
          mt[0] := mt[N-1];
          i := 1;
        end;
      if j >= KeyLength then j := 0;
    end;

  for k := N-1 downto 1 do
    begin
      {$d-   --allow multiply overflow}
      mt[i] := bitxor(mt[i], (bitxor(mt[i-1], rshiftu(mt[i-1], 30)) * 1566083941)) - i; { non linear }
      mt[i] := bitand(mt[i], $ffffffff); { for WORDSIZE > 32 machines }
      {$d+}
      i := i + 1;
      if i >= N then
        begin
          mt[0] := mt[N-1];
          i := 1;
        end;
    end;

    mt[0] := $80000000; { MSB is 1; assuring non-zero initial array }
end;

function IRanMT : Int32;
var
  y : Int32;
  k : Int16;
begin
  if mti >= N then  { generate N words at one Time }
    begin
      { If IRanMT() has not been called, a default initial seed is used }
      if mti = N + 1 then InitMT(5489);

      for k := 0 to (N-M)-1 do
        begin
          y := bitor(bitand(mt[k], UPPER_MASK), bitand(mt[k+1], LOWER_MASK));
          mt[k] := bitxor(mt[k+M], bitxor(rshiftu(y, 1), mag01[bitand(y, $1)]));
        end;

      for k := (N-M) to (N-2) do
        begin
          y := bitor(bitand(mt[k], UPPER_MASK), bitand(mt[k+1], LOWER_MASK));
          mt[k] := bitxor(mt[k - (N - M)], bitxor(rshiftu(y, 1), mag01[bitand(y, $1)]));
        end;

      y := bitor(bitand(mt[N-1], UPPER_MASK), bitand(mt[0], LOWER_MASK));
      mt[N-1] := bitxor(mt[M-1], bitxor(rshiftu(y, 1), mag01[bitand(y, $1)]));

      mti := 0;
    end;

  y := mt[mti];
  mti := mti + 1;

  { Tempering }
  y := bitxor(y, rshiftu(y, 11));
  y := bitxor(y, bitand(rshiftu(y, -7), $9d2c5680));
  y := bitxor(y, bitand(rshiftu(y, -15), $efc60000));
  y := bitxor(y, rshiftu(y, 18));

  IRanMT := y
end;


{ ******************************************************************
    Mersenne Twister Random Number Generator end
  ******************************************************************}



{ ------------ Maurer's random number test ---------- }
var
   tab : array [0..255] of integer;

function maurer(function gen : integer ): real;
label 99;
const
   Q         = 10000;
   K         = 100*Q;
   MEAN      = 7.1836656;
   DEVIATION = 0.00271744365167 {1.5*sqrt(3.282/K)};

var
   sum, ftu : real ;
   i        : integer;
   n        : integer;

begin
   //writeln('initialising, K is ', K );
   for i:=0 to 255 do tab[i] := 0;

   for n:=1 to 11 do i := gen {writeln( 'r is ', gen)};
   for n:=1 to Q do
      tab[gen mod 256] := n;
   for i:=0 to 255 do
      if tab[i] = 0 then begin
         writeln( 'not enough randomness to start up' );
         maurer := 0;//break, return 0;
         goto 99;
      end; {if}


   //writeln( 'start accumulating data ...' );
   sum := 0.0;
   for n:=Q+1 to Q+K do begin
      i := gen;
      sum := sum + ln(n - tab[i mod 256]);
      tab[i mod 256] := n;
      if n mod 512 = 0 then begin
         ftu := sum/(n-Q)/ln(2.0) - MEAN;
         //writeln( chr(13), n:1, ': ft = ', ftu, ',  merit = ', 100*sqr(DEVIATION)/(sqr(DEVIATION)+sqr(ftu)) );
      end; {if}
   end; {for}
   ftu := sum/K/ln(2.0) - MEAN;
   //writeln;
   //writeln( 'ftu = ', ftu, ', merit = ', 100*sqr(DEVIATION)/(sqr(DEVIATION) + sqr(ftu)) );

    {ftu <= DEVIATION and ftu >= -DEVIATION;}
   maurer := 100*sqr(DEVIATION)/(sqr(DEVIATION) + sqr(ftu)) ;
99:
end; { maurer }

begin

#if defined __unix__
   assign(rf, '/dev/urandom');
   reset(rf);
   writeln('/dev/urandom score is ', maurer(drandom):1:1);
   seed := drandom;
#else
   seed := iseed;
#endif

   while seed <=0 do seed := seed + maxint;
   writeln('using seed = ', seed:1);

   srand(seed);
   writeln('c library rand score is ', maurer(rand):1:1);

   writeln('minstd score is ', maurer(minstd):1:1);
{   if maurer(minstd) < 50 then
      writeln('there are better choices of random number generator')
   else
      writeln('this seems to be a good random bit generator');
}

   mag01[0] := 0;
   mag01[1] := MATRIX_A;
   init[0] := $123;
   init[1] := $234;
   init[2] := $345;
   init[3] := $456;
   InitMTbyArray(init, 4);

   writeln('mersenne twister score is ', maurer(IRanMT):1:1);

   seed := seed mod 2796203;
   writeln('prand score is ', maurer(prand):1:1);

   s1 := rand mod 30000; s2 := rand mod 30000; s3 := rand mod 30000;
   writeln('whrand score is ', maurer(whrand):1:1);

   writeln('arand score is ', maurer(arand):1:1);

   seed := 4;
   writeln('hp35 rand score is ', maurer(randhp35):1:1);

end.


{%%%%%%%%%%%%%%%%%%%%%%%%%%%%% end of rtest.pas %%%%%%%%%%%%%%%%%%%%%%%%%%%%%}
