{ **********************************************************************
  *                            Unit RANDOM.PAS                         *
  *                              Version 1.3                           *
  *                       (c) J. Debord, April 2000                    *
  **********************************************************************
                           Random number generators
  **********************************************************************
  Reference:
  DIEHARD by G. Marsaglia
  ********************************************************************** }

unit Random;

interface

uses
  FMath, Matrices;

procedure RMarIn(Seed1, Seed2 : Integer);
{ ----------------------------------------------------------------------
  Initializes the random number generator.
  The default initialization corresponds to RMarIn(1802, 9373).
  ---------------------------------------------------------------------- }

function IRanMar : LongInt;
{ ----------------------------------------------------------------------
  Returns a 32 bit random number in [ -2,147,483,648 ; 2,147,483,647 ]
  ---------------------------------------------------------------------- }

function RanMar : Float;
{ ----------------------------------------------------------------------
  Returns a random number in [0, 1[
  ---------------------------------------------------------------------- }

function RanGaussStd : Float;
{ ----------------------------------------------------------------------
  Returns a random number from the standard normal distribution
  (i.e. the Gaussian distribution with zero mean and unit variance)
  ---------------------------------------------------------------------- }

function RanGauss(Mu, Sigma : Float) : Float;
{ ----------------------------------------------------------------------
  Returns a random number from a Gaussian distribution of mean Mu and
  standard deviation Sigma
  ---------------------------------------------------------------------- }

procedure RanMult(M : PVector; L : PMatrix; N : Integer; X : PVector);
{ ----------------------------------------------------------------------
  Samples a vector X from the N-dimensioned multinormal distribution
  with mean vector M. L is the Cholesky factor of the variance-covariance
  matrix.
  ---------------------------------------------------------------------- }

implementation

var
  X1, X2, C1, C2 : LongInt;

  procedure RMarIn(Seed1, Seed2 : Integer);
  begin
    X1 := Seed1;
    X2 := Seed2;
    C1 := 0;
    C2 := 0;
  end;

  function IRanMar : LongInt;
  var
    Y1, Y2 : LongInt;
  begin
    Y1 := 18000 * X1 + C1;
    X1 := Y1 and 65535;
    C1 := Y1 shr 16;
    Y2 := 30903 * X2 + C2;
    X2 := Y2 and 65535;
    C2 := Y2 shr 16;
    IRanMar := (X1 shl 16) + (X2 and 65535);
  end;

  function RanMar : Float;
  begin
    RanMar := (IRanMar + 2147483648.0) / 4294967296.0;
  end;

  function RanGaussStd : Float;
  { Computes 2 random numbers from the standard normal distribution,
    returns one and saves the other for the next call }
  const
    Gauss_Save : Float   = 0.0;    { Saves a random number }
    Gauss_Set  : Boolean = False;  { Flags if a number has been saved }
  var
    X, R, Theta, SinTheta, CosTheta : Float;
  begin
    if not Gauss_Set then
      begin
        X := RanMar;
        if X = 0.0 then X := MACHEP;
        R := Sqrt(- 2.0 * Ln(X));           { Box-Muller algorithm }
        Theta := TWOPI * RanMar;
        SinCos(Theta, SinTheta, CosTheta);
        RanGaussStd := R * CosTheta;        { Return 1st number }
        Gauss_Save := R * SinTheta;         { Save 2nd number }
      end
    else
      RanGaussStd := Gauss_Save;            { Return saved number }
    Gauss_Set := not Gauss_Set;
  end;

  function RanGauss(Mu, Sigma : Float) : Float;
  { Returns a random number from the normal distribution
    with mean Mu and standard deviation Sigma }
  begin
    RanGauss := Mu + Sigma * RanGaussStd;
  end;

  procedure RanMult(M : PVector; L : PMatrix; N : Integer; X : PVector);
  var
    U : PVector;
    I, J : Integer;
  begin
    { Form a vector of N independent standard normal variates }
    DimVector(U, N);
    for I := 1 to N do
      U^[I] := RanGaussStd;

    { Form X = M + L*U, which follows the multinormal distribution }
    for I := 1 to N do
      begin
        X^[I] := M^[I];
        for J := 1 to I do
          X^[I] := X^[I] + L^[I]^[J] * U^[J];
      end;
    DelVector(U, N);
  end;

begin
  RMarIn(1802, 9373);
end.
