GLScene/Src/Stage.RandomLib.pas

//
// The graphics engine GLScene
//
unit Stage.RandomLib;
(*
  Some of the code is translated from:
  RANLIB - Library of Fortran Routines for Random Number Generation
  OR
  Adapted from Fortran 77 code from the book:
  Dagpunar, J. 'Principles of random variate generation' Clarendon Press,
  Oxford, 1988. ISBN 0-19-852202-9.
  First Delphi Version prepared by Glenn Crouch of ESB Consultancy
  In most of Dagpunar's routines, there is a test to see
  whether the value of one or two floating-point parameters
  has changed since the last call. These tests have been
  replaced by using a logical variable FIRST. This should be
  set to .TRUE. on the first call using new values of the
  parameters, and .FALSE. if the parameter values are the same
  as for the previous call.

  Generate a random ordering of the integers 1 .. N
  Random_Order Initialize (seed) the uniform random number
  generator for ANY compiler Seed_Rrandom_Number

  * Two functions are provided for the binomial distribution.
  If the parameter values remain constant, it is recommended
  that the first function is used (random_binomial1). If one or
  both of the parameters change, use the second function
  (random_binomial2).

  Delphi's own random number generator, Random, is used to
  provide a source of uniformly distributed random numbers.
  At this stage, only one random number is generated at
  each call to one of the functions above.

  The unit uses the following functions which are included
  here: bin_prob to calculate a single binomial probability
  lngamma to calculate the logarithm to base e of the gamma
  function
*)

interface

uses
  System.Math,
  System.SysUtils;

type
  // Used as a Dynamic Vector
  TAMFloatVector = array of Extended;

type
  // Used for passing Random Number Generators
  TRandomGenFunction = function: Extended;

(*
  Used to call Random Number Generator.
  Remember to use Randomize if you don't want repeatable sequences.
*)
function GenRandom: Extended;
(*
Common random distributions  TEST :
                           Unif.      Expo.    Norm.
Total Numbers             100000     100000   100000
Given Mean                    50         25       25
Computed Mean               49.9       24.9     24.9
Expect.Strd_Dev               29         25        1
*)
// Uniform Random single numbers in RLow..RHigh range
function RandUniform(RLow, RHigh: Single): Single; overload;
function RandUniform(RLow, RHigh: Integer): Integer; overload;
function RandExponent(Mean: Single): single;
function RandNorm(Mean, StDev: single): single;
function RandLogNorm(Mean, StDev: single): single;
//  Marsaglia-Bray algorithm
function RandGauss(Mean, StdDev: single): single;
(*
  Random Number for Empiricle Distribution
  X[I]  - middle of I-interval
  Fx[I] - cummulative empiricle distribution with values in [0..1]
*)
function RandEmpiric(X: array of single; Fx: array of single; NClass: integer): single;

(*
  Adapted from the following Fortran 77 code
  ALGORITHM 712, COLLECTED ALGORITHMS FROM ACM.
  THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE,
  VOL. 18, NO. 4, DECEMBER, 1992, PP. 434-435.
  The algorithm uses the ratio of uniforms method of A.J. Kinderman
  and J.F. Monahan augmented with quadratic bounding curves
  Returns a normally distributed pseudo-random number with zero mean and unit variance
*)
function Random_Normal: Extended; overload;
(*
  Adapted from COLLECTED ALGORITHMS FROM ACM.
  THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE,
  VOL. 18, NO. 4, DECEMBER, 1992, PP. 434-435.
  Returns a normally distributed pseudo-random number with zero mean and unit variance.
  The algorithm uses the ratio of uniforms method of A.J. Kinderman
  and J.F. Monahan augmented with quadratic bounding curves
*)
function Random_Normal(RandomGenerator: TRandomGenFunction): Extended; overload;
function Random_Normal(const Mean, StdDev: Extended): Extended; overload;
function Random_Normal(const Mean, StdDev: Extended;
  RandomGenerator: TRandomGenFunction): Extended; overload;

// Adapted from approach suggested by Greg Hood, used with Permission.
function Random_LogNormal(const Mean, StdDev: Extended): Extended; overload;
// Adapted from approach suggested by Greg Hood, used with Permission.
function Random_LogNormal(const Mean, StdDev: Extended;
  RandomGenerator: TRandomGenFunction): Extended; overload;

(*
  GENERATES A RANDOM GAMMA VARIATE.
  CALLS EITHER random_gamma1 (Shape > 1.0) OR random_exponential (Shape = 1.0)
  OR random_gamma2 (Shape < 1.0).
  Shape = Shape PARAMETER OF DISTRIBUTION (0 < REAL).
*)
function Random_Gamma(const Shape: Extended): Extended; overload;
(*
  GENERATES A RANDOM GAMMA VARIATE.
  CALLS EITHER random_gamma1 (Shape > 1.0)
  OR random_exponential (Shape = 1.0)
  OR random_gamma2 (Shape < 1.0).
  Shape = Shape PARAMETER OF DISTRIBUTION (0 < REAL)
*)
function Random_Gamma(const Shape: Extended;
  RandomGenerator: TRandomGenFunction): Extended; overload;
(*
  Uses the algorithm in Marsaglia, G. and Tsang, W.W. (2000)
  'A simple method for generating gamma variables', Trans. om Math.
  Software (TOMS), vol.26(3), pp.363-372.
  Generates a random gamma deviate for shape parameter Shape >= 1.
*)
function Random_Gamma1(const Shape: Extended;
  RandomGenerator: TRandomGenFunction): Extended;

(*
  Generates a random variate from the chi-squared distribution
  with dff degrees of freedom.
*)
function Random_ChiSq(const DF: Integer): Extended; overload;
function Random_ChiSq(const DF: Integer; RandomGenerator: TRandomGenFunction)
  : Extended; overload;

(*
  Generates a random variate in [0,INFINITY) FROM
  A NEGATIVE EXPONENTIAL DlSTRIBUTION WlTH DENSITY PROPORTIONAL
  TO EXP(-random_exponential), USING INVERSION.
*)
function Random_Exponential: Extended; overload;
function Random_Exponential(RandomGenerator: TRandomGenFunction)
  : Extended; overload;

(*
  Generates a random variate from the Weibull distribution with
  probability density:
  a
  a-1  -x
  f(x) = a.x    e
*)
function Random_Weibull(const a: Extended): Extended; overload;
function Random_Weibull(const a: Extended; RandomGenerator: TRandomGenFunction)
  : Extended; overload;

(*
  Generates a random variate in [0,1]
  FROM A BETA DISTRIBUTION WITH DENSITY
  PROPORTIONAL TO BETA**(AA-1) * (1-BETA)**(BB-1).
  USING CHENG'S LOG LOGISTIC METHOD.
  AA = Shape PARAMETER FROM DISTRIBUTION (0 < REAL)
  BB = Shape PARAMETER FROM DISTRIBUTION (0 < REAL)
*)
function Random_Beta(const aa, bb: Extended): Extended; overload;
function Random_Beta(const aa, bb: Extended;
  RandomGenerator: TRandomGenFunction): Extended; overload;

(*
  Generates a random variate FROM A
  T DISTRIBUTION USING KINDERMAN AND MONAHAN'S RATIO METHOD.
  DF = DEGREES OF FREEDOM OF DISTRIBUTION (DF >= 1)
*)
function Random_T(const DF: Integer): Extended; overload;
function Random_T(const DF: Integer; RandomGenerator: TRandomGenFunction)
  : Extended; overload;

(*
  GENERATES AN N VARIATE RANDOM NORMAL
  Vector USING A CHOLESKY DECOMPOSITION.

  ARGUMENTS:
  N = NUMBER OF VARIATES IN Vector (INPUT,INTEGER >= 1)
  H(J) = J'TH ELEMENT OF Vector OF MEANS (INPUT,REAL)
  X(J) = J'TH ELEMENT OF DELIVERED Vector (OUTPUT,REAL)

  D(J*(J-1)/2+I) = (I,J)'TH ELEMENT OF VARIANCE MATRIX (J> = I) (INPUT,REAL)
  F((J-1)*(2*N-J)/2+I) = (I,J)'TH ELEMENT OF LOWER TRIANGULAR
  DECOMPOSITION OF VARIANCE MATRIX (J <= I) (OUTPUT,REAL)
  ier = 1 if the input covariance matrix is not +ve definite = 0 otherwise
*)
procedure Random_MVNorm(const h, d: TAMFloatVector; var f, x: TAMFloatVector;
  var ier: Integer); overload;
procedure Random_MVNorm(const h, d: TAMFloatVector; var f, x: TAMFloatVector;
  var ier: Integer; RandomGenerator: TRandomGenFunction); overload;

(*
  Generates a random variate in [0,INFINITY] FROM
  A REPARAMETERISED GENERALISED INVERSE GAUSSIAN (GIG) DISTRIBUTION
  WITH DENSITY PROPORTIONAL TO  GIG**(H-1) * EXP(-0.5*B*(GIG+1/GIG))
  USING A RATIO METHOD.

  h = PARAMETER OF DISTRIBUTION (0 <= REAL)
  b = PARAMETER OF DISTRIBUTION (0 < REAL)
*)
function Random_Inv_Gauss(const h, b: Extended): Extended; overload;
function Random_Inv_Gauss(const h, b: Extended;
  RandomGenerator: TRandomGenFunction): Extended; overload;

(*
  GENERATES A RANDOM BINOMIAL VARIATE USING C.D.Kemp's method.
  This algorithm is suitable when many random variates are required
  with the SAME parameter values for n & p.
  P = Bernoulli SUCCESS PROBABILITY
  (0 <= REAL <= 1)
  N = Number of Bernoulli trials
  (1 <= INTEGER)
  Reference: Kemp, C.D. (1986). `A modal method for generating binomial
  variables', Commun. Statist. - Theor. Meth. 15(3), 805-813.
*)
function Random_Binomial1(const n: Integer; const p: Extended): Int64; overload;
function Random_Binomial1(const n: Integer; const p: Extended;
  RandomGenerator: TRandomGenFunction): Int64; overload;

(*
  Generates a single random deviate from a binomial
  distribution whose number of trials is N and whose
  probability of an event in each trial is P.

  Arguments
  N  --> The number of trials in the binomial distribution
  from which a random deviate is to be generated.
  INTEGER N
  P  --> The probability of an event in each trial of the
  binomial distribution from which a random deviate
  is to be generated.
  REAL P
*)
function Random_Binomial2(const n: Int64; const pp: Extended): Int64; overload;
(*
  FIRST --> Set FIRST = .TRUE. for the first call to perform initialization
  the set FIRST = .FALSE. for further calls using the same pair
  of parameter values (N, P).
  LOGICAL FIRST
  random_binomial2 <-- A random deviate yielding the number of events
  from N independent trials, each of which has a probability of event P.
  INTEGER --> random_binomial
  Method:
  This is algorithm BTPE from:
  Kachitvichyanukul, V. and Schmeiser, B. W.
  Binomial Random Variate Generation.
  Communications of the ACM, 31, 2 (February, 1988) 216.
*)
function Random_Binomial2(const n: Int64; const pp: Extended;
  RandomGenerator: TRandomGenFunction): Int64; overload;
(*
  GENERATES A RANDOM NEGATIVE BINOMIAL VARIATE USING UNSTORED
  INVERSION AND/OR THE REPRODUCTIVE PROPERTY.
  SK = NUMBER OF FAILURES REQUIRED (Dagpunar's words!)
  = the `power' parameter of the negative binomial
  (0 < REAL)
  P = BERNOULLI SUCCESS PROBABILITY
  (0 < REAL < 1)
  THE PARAMETER H IS SET SO THAT UNSTORED INVERSION ONLY IS USED WHEN P <= H,
  OTHERWISE A COMBINATION OF UNSTORED INVERSION AND
  THE REPRODUCTIVE PROPERTY IS USED.
*)
function Random_Neg_Binomial(const sk, p: Extended): Int64; overload;
function Random_Neg_Binomial(const sk, p: Extended;
  RandomGenerator: TRandomGenFunction): Int64; overload;
(*
  Algorithm VMD
  Arguments:  k (real) -  parameter of the von Mises distribution.
*)
function Random_von_Mises(const k: Extended): Extended; overload;
function Random_von_Mises(const k: Extended;
  RandomGenerator: TRandomGenFunction): Extended; overload;
// Generate a random deviate from the standard Cauchy distribution
function Random_Cauchy: Extended; overload;
function Random_Cauchy(RandomGenerator: TRandomGenFunction): Extended; overload;
(*
  Generates a single random deviate from a Poisson distribution with mean mu.
  Arguments:
  mu --> The mean of the Poisson distribution from which
  a random deviate is to be generated.
  REAL mu
  Method: For details see Ahrens, J.H. and Dieter, U.
  Computer Generation of Poisson Deviates
  From Modified Normal Distributions. ACM Trans. Math. Software, 8, 2
  (June 1982),163-179

  TABLES: COEFFICIENTS A0-A7 FOR STEP F. FACTORIALS FACT
  COEFFICIENTS A(K) - FOR PX = FK*V*V*SUM(A(K)*V**K)-DEL
*)
function Random_Poisson(var mu: Extended): Int64; overload;
function Random_Poisson(var mu: Extended; RandomGenerator: TRandomGenFunction)
  : Int64; overload;

(*
  Logarithm to base e of the gamma function.
  Accurate to about 1.e-14.
  Programmer: Alan Miller
  Latest revision of Fortran 77 version - 28 February 1988
*)
function LnGamma(const x: Extended): Extended;

// =====================================================================
implementation
// =====================================================================

const
  VSmall = MinDouble;
  VLarge = MaxDouble;

function GenRandom: Extended;
begin
  Result := Random;
end;

//---------------------------------------------------------------

function RandUniform(RLow, RHigh: Single): Single; overload;
begin
  Result := RLow + (RHigh - RLow) * Random();
end;


function RandUniform(RLow, RHigh: Integer): Integer; overload;
begin
  Result := RLow + Round((RHigh - RLow) * Random());
end;

//------------------------------------------------------
//                  RANDOM EXPONENTIAL
//------------------------------------------------------
function RandExponent(Mean: single): single;
begin
  RandExponent := -Mean * Ln(Random());
end;

//------------------------------------------------------
//                   RANDOM  NORMAL
//------------------------------------------------------
function RandNorm(Mean, StDev: single): single;
var
  S, Y, U1, U2: single;
begin
  repeat
    U1 := 2 * Random() - 1;
    U2 := 2 * Random() - 1;
    S  := U1 * U1 + U2 * U2;
  until (S < 1);
  Y := U2 * Sqrt(-2 * Ln(S) / S);
  RandNorm := Mean + StDev * Y;
end; // of Normal

//-------------------------------------------------------
//              RANDOM GAUSSIAN
//-------------------------------------------------------

function RandGauss(Mean, StdDev: single): single;
var
  U1, S2: single;
begin
  repeat
    U1 := 2 * Random - 1;
    S2 := Sqr(U1) + Sqr(2 * Random - 1);
  until S2 < 1;
  RandGauss := Sqrt(-2 * Ln(S2) / S2) * U1 * StdDev + Mean;
end;

//------------------------------------------------------
//                RANDOM  LOGNORMAL
//------------------------------------------------------
function RandLogNorm(Mean, StDev: single): single;
var
  S, Y, U1, U2: single;
begin
  S := 0.0;
  repeat
    U1 := 2 * Random - 1;
    U2 := 2 * Random - 1;
    S  := U1 * U1 + U2 * U2;
  until (S < 1);
  Y := Sqrt(-2 * Ln(S) / S) * U2;
  RandLogNorm := Exp(Y) + Mean;
end; // of LogNorm

function RandEmpiric(X: array of single; Fx: array of single; NClass: integer): single;
var
  I: integer;
  R: single;
begin
  RandEmpiric := 0;
  R := Random;
  for I := 1 to NClass do
  begin
    if (R <= Fx[I]) then
      RandEmpiric := X[I];
  end;
end;

//---------------------------------------------------------------

function Random_Normal(RandomGenerator: TRandomGenFunction): Extended;
const
  s = 0.449871;
  t = -0.386595;
  a = 0.19600;
  b = 0.25472;
  r1 = 0.27597;
  r2 = 0.27846;
var
  u, v, x, y, q: Extended;
  Done: Boolean;

begin
  // Generate P = (u,v) uniform in rectangle enclosing acceptance region
  Done := False;
  repeat
    u := RandomGenerator;
    v := RandomGenerator;
    v := 1.7156 * (v - 0.5);

    // Evaluate the quadratic form
    x := u - s;
    y := abs(v) - t;
    q := Sqr(x) + y * (a * y - b * x);

    // Accept P if inside inner ellipse
    if (q < r1) then
      Done := True
    else if (q <= r2) and (Sqr(v) < -4.0 * Ln(u) * Sqr(u)) then
      Done := True;
  until Done;
  // Return ratio of P's coordinates as the normal deviate
  if u < VSmall then
    raise EMathError.Create('Divided By Zero');

  Result := v / u;
end;

//---------------------------------------------------------------

function Random_Normal: Extended;
begin
  Result := Random_Normal(GenRandom);
end;

function Random_Normal(const Mean, StdDev: Extended;
  RandomGenerator: TRandomGenFunction): Extended;
begin
  Result := Random_Normal(RandomGenerator) * StdDev + Mean;
end;

function Random_Normal(const Mean, StdDev: Extended): Extended;
begin
  Result := Random_Normal(Mean, StdDev, GenRandom);
end;

//---------------------------------------------------------------

function Random_LogNormal(const Mean, StdDev: Extended;
  RandomGenerator: TRandomGenFunction): Extended;
var
  C, M, s: Extended;
begin
  if Mean <= 0 then
    raise EMathError.Create('Invalid Mean');
  if StdDev <= 0 then
    raise EMathError.Create('Invalid StdDev');

  C := StdDev / Mean;
  M := Ln(Mean) - 0.5 * Ln(Sqr(C) + 1);
  s := Sqrt(Ln(Sqr(C) + 1));
  Result := Exp(Random_Normal(M, s, RandomGenerator));
end;

function Random_LogNormal(const Mean, StdDev: Extended): Extended;
begin
  Result := Random_LogNormal(Mean, StdDev, GenRandom);
end;

//---------------------------------------------------------------

function Random_Gamma1(const Shape: Extended;
  RandomGenerator: TRandomGenFunction): Extended;
var
  C, d: Extended;
  u, v, x: Extended;
  Done: Boolean;
begin
  if (Shape <= 1.0) then
    raise EMathError.Create('Invalid Shape');
  d := Shape - 1 / 3;
  C := 1 / Sqrt(9.0 * d);
  Result := 0.0;
  Done := False;
  repeat
    // Generate v = (1+cx)^3 where x is random normal; repeat if v <= 0.
    repeat
      x := Random_Normal(RandomGenerator);
      v := Power(1 + C * x, 3);
    until v > 0;
    // Generate uniform variable U
    u := RandomGenerator;
    if (u < 1.0 - 0.0331 * Sqr(Sqr(x))) then
    begin
      Result := d * v;
      Done := True;
    end
    else if Ln(u) < 0.5 * Sqr(x) + d * (1.0 - v + Ln(v)) then
    begin
      Result := d * v;
      Done := True;
    end
    until Done;
  end;

  (*
    Generates a RANDOM VARIATE in [0,INFINITY) from
    A GAMMA DISTRIBUTION with DENSITY proportional to
    GAMMA2**(S-1) * EXP(-GAMMA2),
    using a switching method.
    S = Shape PARAMETER of distribution
    (real < 1.0)
   *)
  function Random_Gamma2(const Shape: Extended;
    RandomGenerator: TRandomGenFunction): Extended;
  var
    r, x, w: Extended;
    a, p, C, uf, vr, d: Extended;
    Done: Boolean;
  begin
    if (Shape <= 0.0) and (Shape >= 1.0) then
      raise EMathError.Create('Invalid Shape');
    if (Shape < VSmall) then
      raise EMathError.Create('Shape Too Small');

    a := 1.0 - Shape;
    p := a / (a + Shape * Exp(-a));
    C := 1.0 / Shape;
    uf := p * Power(VSmall / a, Shape);
    vr := 1.0 - VSmall;
    d := a * Ln(a);
    w := 0;
    x := 0;
    Done := False;
    repeat
      r := RandomGenerator;
      if (r >= vr) then
        Continue;
      if (r > p) then
      begin
        x := a - Ln((1.0 - r) / (1.0 - p));
        w := a * Ln(x) - d;
      end
      else if (r > uf) then
      begin
        x := a * Power(r / p, C);
        w := x;
      end
      else
      begin
        Result := 0.0;
        Exit;
      end;
      r := RandomGenerator;
      if (1.0 - r <= w) and (r > 0.0) then
      begin
        if (r * (w + 1.0) >= 1.0) then
          Continue;
        if (-Ln(r) <= w) then
          Continue;
      end;
      Done := True;
    until Done;
    Result := x;
  end;

//---------------------------------------------------------------

  function Random_Gamma(const Shape: Extended;
    RandomGenerator: TRandomGenFunction): Extended;
  begin
    if (Shape <= 0.0) then
      raise EMathError.Create('Invalid Shape');
    if (Shape > 1.0) then
      Result := Random_Gamma1(Shape, RandomGenerator)
    else if (Shape < 1.0) then
      Result := Random_Gamma2(Shape, RandomGenerator)
    else
      Result := Random_Exponential(RandomGenerator);
  end;

  function Random_Gamma(const Shape: Extended): Extended;
  begin
    Result := Random_Gamma(Shape, GenRandom);
  end;

  function Random_ChiSq(const DF: Integer; RandomGenerator: TRandomGenFunction)
    : Extended;
  begin
    Result := 2.0 * Random_Gamma(0.5 * DF, RandomGenerator)
  end;

  function Random_ChiSq(const DF: Integer): Extended;
  begin
    Result := Random_ChiSq(DF, GenRandom);
  end;

//---------------------------------------------------------------

  function Random_Exponential(RandomGenerator: TRandomGenFunction): Extended;
  var
    r: Extended;
    Done: Boolean;
  begin
    Done := False;
    repeat
      r := RandomGenerator;
      if (r > VSmall) then
        Done := True;
    until Done;
    Result := -Ln(r)
  end;

  function Random_Exponential: Extended;
  begin
    Result := Random_Exponential(GenRandom);
  end;

  function Random_Weibull(const a: Extended;
    RandomGenerator: TRandomGenFunction): Extended;
  begin
    if abs(a) <= VSmall then
      raise EMathError.Create('Invalid Value');

    Result := Power(Random_Exponential(RandomGenerator), 1.0 / a);
  end;

 //---------------------------------------------------------------

  function Random_Weibull(const a: Extended): Extended;
  begin
    Result := Random_Weibull(a, GenRandom);
  end;

  function Random_Beta(const aa, bb: Extended;
    RandomGenerator: TRandomGenFunction): Extended;
  const
    aln4 = 1.3862943611198906;
  var
    a, b, g, r, s, x, y, z: Extended;
    d, f, h, t, C: Extended;
    Swap: Boolean;
    Done: Boolean;
  begin
    if (aa <= 0.0) or (bb <= 0.0) then
      raise EMathError.Create('Invalid Value');
    a := aa;
    b := bb;
    Swap := b > a;
    if Swap then
    begin
      g := b;
      b := a;
      a := g;
    end;
    d := a / b;
    f := a + b;
    if (b > 1.0) then
    begin
      h := Sqrt((2.0 * a * b - f) / (f - 2.0));
      t := 1.0;
    end
    else
    begin
      h := b;
      t := 1.0 / (1.0 + Power(a / (VLarge * b), b));
    end;
    C := a + h;

    Done := False;
    Result := 0.0;
    repeat
      r := RandomGenerator;
      x := RandomGenerator;
      s := Sqr(r) * x;
      if (r < VSmall) or (s <= 0.0) then
        Continue;

      if (r < t) then
      begin
        x := Ln(r / (1.0 - r)) / h;
        y := d * Exp(x);
        z := C * x + f * Ln((1.0 + d) / (1.0 + y)) - aln4;
        if (s - 1.0 > z) then
        begin
          if (s - s * z > 1.0) then
            Continue;
          if (Ln(s) > z) then
            Continue;
        end;
        Result := y / (1.0 + y);
      end
      else
      begin
        if 4.0 * s > Power(1.0 + 1.0 / d, f) then
          Continue;
        Result := 1.0
      end;
      Done := True
    until Done;

    if Swap then
      Result := 1 - Result;
  end;

  //---------------------------------------------------------------

  function Random_Beta(const aa, bb: Extended): Extended;
  begin
    Result := Random_Beta(aa, bb, GenRandom);
  end;

  function Random_T(const DF: Integer; RandomGenerator: TRandomGenFunction)
    : Extended;
  var
    s, C, a, f, g: Extended;
    r, x, v: Extended;
    Done: Boolean;
  begin
    if (DF < 1) then
      raise EMathError.Create('Invalid Value');
    s := DF;
    C := -0.25 * (s + 1.0);
    a := 4.0 / Power(1.0 + 1.0 / s, C);
    f := 16.0 / a;
    if (DF > 1) then
    begin
      g := s - 1.0;
      g := Power((s + 1.0) / g, C) * Sqrt((s + s) / g);
    end
    else
      g := 1.0;
    Done := False;
    x := 0;
    repeat
      r := RandomGenerator;
      if (r <= 0.0) then
        Continue;
      v := RandomGenerator;
      x := (2.0 * v - 1.0) * g / r;
      v := Sqr(x);
      if (v > 5.0 - a * r) then
      begin
        if (DF >= 1) and (r * (v + 3.0) > f) then
          Continue;
        if r > Power(1.0 + v / s, C) then
          Continue;
      end;
      Done := True;
    until Done;
    Result := x;
  end;

//---------------------------------------------------------------

  function Random_T(const DF: Integer): Extended;
  begin
    Result := Random_T(DF, GenRandom);
  end;

  procedure Random_MVNorm(const h, d: TAMFloatVector; var f, x: TAMFloatVector;
    var ier: Integer; RandomGenerator: TRandomGenFunction);
  var
    i, j, M, n, n2, fn: Integer;
    y, v: Extended;
  begin
    n := Length(h);
    if (n = 0) then
      raise EMathError.Create('Vector Empty');
    fn := n * (n + 1) div 2;
    if (Length(d) < fn) then
      raise EMathError.Create('Vector Too Small');

    SetLength(f, fn);
    ier := 0;
    n2 := 2 * n;
    if (d[0] <= 0.0) then
    begin
      ier := 1;
      Exit;
    end;
    f[0] := Sqrt(d[0]);
    y := 1.0 / f[0];
    for j := 2 to n do
      f[j - 1] := d[j * (j - 1) div 2] * y;
    for i := 2 to n do
    begin
      v := d[i * (i - 1) div 2 + i - 1];
      for M := 1 to i - 1 do
        v := v - Sqr(f[(M - 1) * (n2 - M) div 2 + i - 1]);

      if (v <= 0.0) then
      begin
        ier := 1;
        Exit;
      end;
      v := Sqrt(v);
      y := 1.0 / v;
      f[(i - 1) * (n2 - i) div 2 + i - 1] := v;

      for j := i + 1 to n do
      begin
        v := d[j * (j - 1) div 2 + i - 1];
        for M := 1 to i - 1 do
        begin
          v := v - f[(M - 1) * (n2 - M) div 2 + i - 1] *
            f[(M - 1) * (n2 - M) div 2 + j - 1]
        end;
        f[(i - 1) * (n2 - i) div 2 + j - 1] := v * y;
      end;
    end;
    x := Copy(h, 0, n);
    for j := 1 to n do
    begin
      y := Random_Normal(RandomGenerator);
      for i := j to n do
        x[i - 1] := x[i - 1] + f[(j - 1) * (n2 - j) div 2 + i - 1] * y;
    end;
  end;

  procedure Random_MVNorm(const h, d: TAMFloatVector; var f, x: TAMFloatVector;
    var ier: Integer);
  begin
    Random_MVNorm(h, d, f, x, ier, GenRandom);
  end;

 //---------------------------------------------------------------

  function Random_Inv_Gauss(const h, b: Extended;
    RandomGenerator: TRandomGenFunction): Extended;
  var
    ym, xm, r, w, r1, r2, x: Extended;
    a, C, d, e: Extended;
    Done: Boolean;
  begin
    if (h < 0.0) or (b <= 0.0) then
      raise EMathError.Create('Invalid Value');

    if (h > 0.25 * b * Sqrt(VLarge)) then
      raise EMathError.Create('Invalid Value');

    e := Sqr(b);
    d := h + 1.0;
    ym := (-d + Sqrt(Sqr(d) + e)) / b;
    if (ym < VSmall) then
      raise EMathError.Create('Invalid Value');

    d := h - 1.0;
    xm := (d + Sqrt(Sqr(d) + e)) / b;
    if (xm < VSmall) then
      raise EMathError.Create('Invalid Value');

    d := 0.5 * d;
    e := -0.25 * b;
    r := xm + 1.0 / xm;
    w := xm * ym;
    a := Power(w, -0.5 * h) * Sqrt(xm / ym) * Exp(-e * (r - ym - 1.0 / ym));
    if (a < VSmall) then
      raise EMathError.Create('Invalid Value');
    C := -d * Ln(xm) - e * r;
    Done := False;
    x := 0.0;
    repeat
      r1 := RandomGenerator;
      if (r1 > 0.0) then
      begin
        r2 := RandomGenerator;
        x := a * r2 / r1;
        if (x > 0.0) then
        begin
          if (Ln(r1) < d * Ln(x) + e * (x + 1.0 / x) + C) then
            Done := True
        end;
      end
      until Done;
      Result := x;
    end;

    function Random_Inv_Gauss(const h, b: Extended): Extended;
    begin
      Result := Random_Inv_Gauss(h, b, GenRandom);
    end;

//---------------------------------------------------------------

    function LnGamma(const x: Extended): Extended;
    const
      a1 = -4.166666666554424E-02;
      a2 = 2.430554511376954E-03;
      a3 = -7.685928044064347E-04;
      a4 = 5.660478426014386E-04;
      lnrt2pi = 9.189385332046727E-1;
      pi = 3.141592653589793E0;
    var
      temp, arg, product: Extended;
      reflect: Boolean;
    begin
      // lngamma is not defined if x = 0 or a negative integer.
      if (x = 0.0) or ((x < 0.0) and (abs(x - Int(x)) < VSmall)) then
        raise EMathError.Create('Invalid Value');

      // If x < 0, use the reflection formula:
      // gamma(x) * gamma(1-x) = pi * cosec(pi.x)

      reflect := x < 0.0;
      if reflect then
        arg := 1.0 - x
      else
        arg := x;

      // Increase the argument, if necessary, to make it > 10.
      product := 1.0;
      while (arg <= 10.0) do
      begin
        product := product * arg;
        arg := arg + 1.0;
      end;

      // Use a polynomial approximation to Stirling's formula.
      // N.B. The real Stirling's formula is used here, not the simpler, but less
      // accurate formula given by De Moivre in a letter to Stirling, which
      // is the one usually quoted.

      arg := arg - 0.5;
      temp := 1.0 / Sqr(arg);
      Result := lnrt2pi + arg *
        (Ln(arg) - 1.0 + (((a4 * temp + a3) * temp + a2) * temp + a1) * temp) -
        Ln(product);

      if reflect then
      begin
        temp := Sin(pi * x);
        Result := Ln(pi / temp) - Result;
      end;
    end;

    // Calculate a binomial probability

    function Bin_Prob(const n: Int64; const p: Extended; const r: Int64)
      : Extended;
    begin
      Result := Exp(LnGamma(n + 1.0) - LnGamma(r + 1.0) - LnGamma(n - r + 1.0) +
        r * Ln(p) + (n - r) * Ln(1.0 - p));
    end;

   //---------------------------------------------------------------

    function Random_Binomial1(const n: Integer; const p: Extended;
      RandomGenerator: TRandomGenFunction): Int64;
    var
      ru, rd: Int64;
      r0: Int64;
      u, pd, pu: Real;
      odds_ratio, p_r: Real;
    begin
      r0 := Trunc((n + 1) * p);
      p_r := Bin_Prob(n, p, r0);
      if p < 1 then
        odds_ratio := p / (1.0 - p)
      else
        odds_ratio := VLarge;

      u := RandomGenerator;
      u := u - p_r;
      if (u < 0.0) then
      begin
        Result := r0;
        Exit;
      end;

      pu := p_r;
      ru := r0;
      pd := p_r;
      rd := r0;

      repeat
        Dec(rd);
        if (rd >= 0) then
        begin
          pd := pd * (rd + 1.0) / (odds_ratio * (n - rd));
          u := u - pd;
          if (u < 0.0) then
          begin
            Result := rd;
            Exit;
          end;
        end;

        Inc(ru);
        if (ru <= n) then
        begin
          pu := pu * (n - ru + 1.0) * odds_ratio / ru;
          u := u - pu;
          if (u < 0.0) then
          begin
            Result := ru;
            Exit;
          end;
        end;
      until False;
    end;

   //---------------------------------------------------------------

    function Random_Binomial1(const n: Integer; const p: Extended): Int64;
    begin
      Result := Random_Binomial1(n, p, GenRandom);
    end;

    function Random_Binomial2(const n: Int64; const pp: Extended;
      RandomGenerator: TRandomGenFunction): Int64;
    var
      alv, amaxp, f, f1, f2, u, v, w, w2, x, x1, x2, ynorm, z, z2: Extended;
      i: Integer;
      ix, ix1, k, mp: Int64;
      M: Int64;
      p, q, xnp, ffm, fm, xnpq, p1, xm, xl, xr, C, al, xll, xlr, p2, p3, p4, qn,
        r, g: Extended;
      Done: Boolean;
    begin
      // SETUP
      p := Min(pp, 1.0 - pp);
      q := 1.0 - p;
      xnp := n * p;

      if (xnp > 30.0) then
      begin
        ffm := xnp + p;
        M := Trunc(ffm);
        fm := M;
        xnpq := xnp * q;
        p1 := Int(2.195 * Sqrt(xnpq) - 4.6 * q) + 0.5;
        xm := fm + 0.5;
        xl := xm - p1;
        xr := xm + p1;
        C := 0.134 + 20.5 / (15.3 + fm);
        al := (ffm - xl) / (ffm - xl * p);
        xll := al * (1.0 + 0.5 * al);
        al := (xr - ffm) / (xr * q);
        xlr := al * (1.0 + 0.5 * al);
        p2 := p1 * (1.0 + C + C);
        p3 := p2 + C / xll;
        p4 := p3 + C / xlr;

        // GENERATE VARIATE, Binomial mean at least 30.
        ix := 0;
        repeat
          ;
          u := RandomGenerator; // 20
          u := u * p4;
          v := RandomGenerator;

          // TRIANGULAR REGION
          if (u <= p1) then
          begin
            ix := Trunc(xm - p1 * v + u);
            Break;
          end;

          // PARALLELOGRAM REGION
          if (u <= p2) then
          begin
            x := xl + (u - p1) / C;
            v := v * C + 1.0 - abs(xm - x) / p1;
            if (v > 1.0) or (v <= 0) then
              Continue;
            ix := Trunc(x);
          end
          else
          begin
            // LEFT TAIL
            if (u <= p3) then
            begin
              ix := Trunc(xl + Ln(v) / xll);
              if (ix < 0) then
                Continue;
              v := v * (u - p2) * xll
            end
            // RIGHT TAIL
            else
            begin
              ix := Trunc(xr - Ln(v) / xlr);
              if (ix > n) then
                Continue;
              v := v * (u - p3) * xlr;
            end;
          end;

          // DETERMINE APPROPRIATE WAY TO PERFORM ACCEPT/REJECT TEST
          k := abs(ix - M);
          if (k <= 20) or (k >= xnpq / 2 - 1) then
          begin
            // EXPLICIT EVALUATION
            f := 1.0;
            r := p / q;
            g := (n + 1) * r;
            if (M < ix) then
            begin
              mp := M + 1;
              for i := mp to ix do
                f := f * (g / i - r);
            end
            else if (M > ix) then
            begin
              ix1 := ix + 1;
              for i := ix1 to M do
                f := f / (g / i - r);
            end;

            if (v > f) then
              Continue
            else
              Break
          end;

          // SQUEEZING USING UPPER AND LOWER BOUNDS ON LOG(F(X))
          amaxp := (k / xnpq) * ((k * (k / 3.0 + 0.625) + 0.1666666666666) /
            xnpq + 0.5);
          ynorm := -Sqr(k) / (2.0 * xnpq);
          alv := Ln(v);
          if (alv < ynorm - amaxp) then
            Break;
          if (alv > ynorm + amaxp) then
            Continue;
          // STIRLING'S (actually de Moivre's) FORMULA TO MACHINE ACCURACY FOR
          // THE FINAL ACCEPTANCE/REJECTION TEST

          x1 := ix + 1;
          f1 := fm + 1.0;
          z := n + 1 - fm;
          w := n - ix + 1.0;
          z2 := Sqr(z);
          x2 := Sqr(x1);
          f2 := Sqr(f1);
          w2 := Sqr(w);
          if (alv - (xm * Ln(f1 / x1) + (n - M + 0.5) * Ln(z / w) + (ix - M) *
            Ln(w * p / (x1 * q)) + (13860.0 -
            (462.0 - (132.0 - (99.0 - 140.0 / f2) / f2) / f2) / f2) / f1 /
            166320.0 + (13860.0 - (462.0 - (132.0 - (99.0 - 140.0 / z2) / z2) /
            z2) / z2) / z / 166320.0 +
            (13860.0 - (462.0 - (132.0 - (99.0 - 140.0 / x2) / x2) / x2) / x2) /
            x1 / 166320.0 + (13860.0 - (462.0 - (132.0 - (99.0 - 140.0 / w2) /
            w2) / w2) / w2) / w / 166320.0) > 0.0) then
          begin
            Continue;
          end
          else
            Break;
        until False;
      end
      else
      begin
        // INVERSE CDF LOGIC FOR MEAN LESS THAN 30
        qn := Power(q, n);
        r := p / q;
        g := r * (n + 1);

        Done := False;
        repeat
          ix := 0;
          f := qn;
          u := RandomGenerator;
          while (u >= f) do
          begin
            if (ix > 110) then
              Done := False
            else
            begin
              Done := True;
              u := u - f;
              ix := ix + 1;
              f := f * (g / ix - r);
            end;
          end;
        until Done;
      end;
      if (pp > 0.5) then
        ix := n - ix;
      Result := ix;
    end;

  //---------------------------------------------------------------

    function Random_Binomial2(const n: Int64; const pp: Extended): Int64;
    begin
      Result := Random_Binomial2(n, pp, GenRandom);
    end;

    function Random_Neg_Binomial(const sk, p: Extended;
      RandomGenerator: TRandomGenFunction): Int64;
    const
      h = 0.7;
    var
      q, x, st, uln, v, r, s, y, g: Extended;
      k, n: Int64;
      i: Integer;
    begin
      if (sk <= 0.0) or (p <= 0.0) or (p >= 1.0) then
        raise EMathError.Create('Invalid Value');

      q := 1.0 - p;
      x := 0.0;
      st := sk;
      if (p > h) then
      begin
        v := 1.0 / Ln(p);
        k := Trunc(st);
        for i := 1 to k do
        begin
          repeat
            r := RandomGenerator;
          until r > 0.0;
          n := Trunc(v * Ln(r));
          x := x + n;
        end;
        st := st - k;
      end;

      s := 0.0;
      uln := -Ln(VSmall);
      if (st > -uln / Ln(q)) then
        raise EMathError.Create('Invalid Value');

      y := Power(q, st);
      g := st;
      r := RandomGenerator;
      while (y <= r) do
      begin
        r := r - y;
        s := s + 1.0;
        y := y * p * g / s;
        g := g + 1.0;
      end;
      Result := Trunc(x + s + 0.5)
    end;

  //---------------------------------------------------------------

    function Random_Neg_Binomial(const sk, p: Extended): Int64;
    begin
      Result := Random_Neg_Binomial(sk, p, GenRandom);
    end;

    // Gaussian integration of exp(k.cosx) from a to b.
    procedure Integral(const a, b: Extended; var res: Extended;
      const dk: Extended);
    var
      xmid, range, x1, x2: Extended;
      x, w: array [1 .. 3] of Extended;
      i: Integer;
    begin
      x[1] := 0.238619186083197;
      x[2] := 0.661209386466265;
      x[3] := 0.932469514203152;
      w[1] := 0.467913934572691;
      w[2] := 0.360761573048139;
      w[3] := 0.171324492379170;

      xmid := (a + b) / 2.0;
      range := (b - a) / 2.0;

      res := 0.0;
      for i := 1 to 3 do
      begin
        x1 := xmid + x[i] * range;
        x2 := xmid - x[i] * range;
        res := res + w[i] * (Exp(dk * Cos(x1)) + Exp(dk * Cos(x2)))
      end;

      res := res * range;
    end;

    function Random_von_Mises(const k: Extended;
      RandomGenerator: TRandomGenFunction): Extended;
    var
      j, n, jj: Integer;
      nk: Integer;
      p: array [1 .. 20] of Extended;
      theta: array [0 .. 20] of Extended;
      xx, sump, r, th, lambda, rlast: Extended;
      dk: Extended;
    begin
      if (k < 0.0) then
        raise EMathError.Create('Invalid Value');

      nk := Trunc(k + k + 1.0);
      if (nk > 20) then
        raise EMathError.Create('Invalid Value');

      dk := k;
      theta[0] := 0.0;
      if (k > 0.5) then
      begin
        // Set up array p of probabilities.
        sump := 0.0;
        for j := 1 to nk do
        begin
          if (j < nk) then
            theta[j] := ArcCos(1.0 - j / k)
          else
            theta[nk] := pi;

          // Numerical integration of e^[k.cos(x)] from theta(j-1) to theta(j)
          Integral(theta[j - 1], theta[j], p[j], dk);
          sump := sump + p[j];
        end;
        for j := 1 to nk do
          p[j] := p[j] / sump;
      end
      else
      begin
        p[1] := 1.0;
        theta[1] := pi;
      end;

      r := RandomGenerator;
      jj := nk;
      for j := 1 to nk do
      begin
        r := r - p[j];
        if (r < 0.0) then
        begin
          jj := j;
          Break;
        end;
      end;
      r := -r / p[jj];

      repeat
        th := theta[jj - 1] + r * (theta[jj] - theta[j - 1]);
        lambda := k - jj + 1.0 - k * Cos(th);
        n := 1;
        rlast := lambda;

        repeat
          r := RandomGenerator;
          if r <= rlast then
          begin
            Inc(n);
            rlast := r;
          end;
        until (r > rlast);

        if not Odd(n) then
          r := RandomGenerator until Odd(n);

        th := abs(th);
        xx := (r - rlast) / (1.0 - rlast) - 0.5;
        if xx < 0 then
          Result := -1 * th
        else
          Result := th;
      end;

//---------------------------------------------------------------

      function Random_von_Mises(const k: Extended): Extended;
      begin
        Result := Random_von_Mises(k, GenRandom);
      end;

      function Random_Cauchy(RandomGenerator: TRandomGenFunction): Extended;
      var
        V1, V2: Extended;
      begin
        repeat
          V1 := 2.0 * (RandomGenerator - 0.5);
          V2 := 2.0 * (RandomGenerator - 0.5);
        until (abs(V2) > VSmall) and (Sqr(V1) + Sqr(V2) < 1.0);
        Result := V1 / V2;
      end;

      function Random_Cauchy: Extended;
      begin
        Result := Random_Cauchy(GenRandom);
      end;

//---------------------------------------------------------------

      function Random_Poisson(var mu: Extended;
        RandomGenerator: TRandomGenFunction): Int64;
      const
        a0 = -0.5;
        a1 = 0.3333333;
        a2 = -0.2500068;
        a3 = 0.2000118;
        a4 = -0.1661269;
        a5 = 0.1421878;
        a6 = -0.1384794;
        a7 = 0.1250060;

      var
        b1, b2, C, c0, c1, c2, c3, del, difmuk, e, fk, fx, fy, g, omega, px, py,
          t, u, v, x, xx: Extended;
        s, d, p, q, p0: Extended;
        j, k, kflag: Integer;
        l, M: Integer;
        pp: array [1 .. 35] of Extended;
        fact: array [1 .. 10] of Extended;
        FirstK0: Boolean;
      begin
        u := 0;
        e := 0;
        fk := 0;
        difmuk := 0;

        fact[1] := 1; // Factorial 0
        fact[2] := 1;
        fact[3] := 2;
        fact[4] := 6;
        fact[5] := 24;
        fact[6] := 120;
        fact[7] := 720;
        fact[8] := 5040;
        fact[9] := 40320;
        fact[10] := 362880; // Factorial 9

        Result := 0;
        if (mu > 10.0) then
        begin
          // C A S E  A. (RECALCULATION OF S, D, L IF MU HAS CHANGED)
          FirstK0 := False;

          s := Sqrt(mu);
          d := 6.0 * Sqr(mu);

          // THE POISSON PROBABILITIES PK EXCEED THE DISCRETE NORMAL
          // PROBABILITIES FK WHENEVER K >= M(MU). L=IFIX(MU-1.1484)
          // IS AN UPPER BOUND TO M(MU) FOR ALL MU >= 10 .

          l := Trunc(mu - 1.1484);

          // STEP N. NORMAL SAMPLE - random_normal() FOR STANDARD NORMAL DEVIATE

          g := mu + s * Random_Normal(RandomGenerator);
          if (g > 0.0) then
          begin
            Result := Trunc(g);
            // STEP I. IMMEDIATE ACCEPTANCE IF ival IS LARGE ENOUGH
            if (Result >= l) then
              Exit;
            // STEP S. SQUEEZE ACCEPTANCE - SAMPLE U
            fk := Result;
            difmuk := mu - fk;
            u := RandomGenerator;
            if (d * u >= Sqr(difmuk) * difmuk) then
              Exit;
          end;

          // STEP P. PREPARATIONS FOR STEPS Q AND H.
          // (RECALCULATIONS OF PARAMETERS IF NECESSARY)
          // .3989423=(2*Pi)**(-.5)  .416667E-1=1./24.  .1428571=1./7.
          // THE QUANTITIES B1, B2, C3, C2, C1, C0 ARE FOR THE HERMITE
          // APPROXIMATIONS TO THE DISCRETE NORMAL PROBABILITIES FK.
          // C=.1069/MU GUARANTEES MAJORIZATION BY THE 'HAT'-FUNCTION.

          omega := 0.3989423 / s;
          b1 := 0.4166667E-1 / mu;
          b2 := 0.3 * Sqr(b1);
          c3 := 0.1428571 * b1 * b2;
          c2 := b2 - 15.0 * c3;
          c1 := b1 - 6.0 * b2 + 45.0 * c3;
          c0 := 1.0 - b1 + 3.0 * b2 - 15.0 * c3;
          C := 0.1069 / mu;

          kflag := -1;
          if (g >= 0.0) then
          begin
            // 'SUBROUTINE' F IS CALLED (KFLAG=0 FOR CORRECT RETURN)
            kflag := 0;
            FirstK0 := True;
          end;
          // STEP E. EXPONENTIAL SAMPLE - random_exponential() FOR STANDARD EXPONENTIAL
          // DEVIATE E AND SAMPLE T FROM THE LAPLACE 'HAT'
          // (IF T <= -.6744 THEN PK < FK FOR ALL MU >= 10.)
          repeat // 50
            if not FirstK0 then
            begin
              repeat // 50
                e := Random_Exponential(RandomGenerator);
                u := RandomGenerator;
                u := u + u - 1.0;
                if u < 0 then
                  t := 1.8 - abs(e)
                else
                  t := 1.8 + abs(e);
              until t > -0.6744;

              Result := Trunc(mu + s * t);
              fk := Result;
              difmuk := mu - fk;

              // 'SUBROUTINE' F IS CALLED (KFLAG=1 FOR CORRECT RETURN)

              kflag := 1
            end
            else
              FirstK0 := False;

            // STEP F. 'SUBROUTINE' F. CALCULATION OF PX, PY, FX, FY.
            // CASE ival < 10 USES FACTORIALS FROM TABLE FACT

            if (Result < 10) then // 70
            begin
              px := -mu;
              py := Power(mu, Result) / fact[Result + 1];
            end
            else
            begin
              // CASE ival >= 10 USES POLYNOMIAL APPROXIMATION
              // A0-A7 FOR ACCURACY WHEN ADVISABLE
              // .8333333E-1=1./12.  .3989423=(2*Pi)**(-.5)

              del := 0.8333333E-1 / fk; // 80
              del := del - 4.8 * Sqr(del) * del;
              v := difmuk / fk;
              if abs(v) > 0.25 then
                px := fk * Ln(1.0 + v) - difmuk - del
              else
                px := fk * Sqr(v) *
                  (((((((a7 * v + a6) * v + a5) * v + a4) * v + a3) * v + a2) *
                  v + a1) * v + a0) - del;
              py := 0.3989423 / Sqrt(fk);
            end;
            x := (0.5 - difmuk) / s; // 110
            xx := Sqr(x);
            fx := -0.5 * xx;
            fy := omega * (((c3 * xx + c2) * xx + c1) * xx + c0);
            if (kflag <= 0) then
            begin
              // STEP Q. QUOTIENT ACCEPTANCE (RARE CASE)

              if (fy - u * fy <= py * Exp(px - fx)) then // 40
                Exit;
            end
            else
            begin
              // STEP H. HAT ACCEPTANCE (E IS REPEATED ON REJECTION)

              if (C * abs(u) <= py * Exp(px + e) - fy * Exp(fx + e)) then // 60
                Exit;
            end;
          until False;
        end
        else
        begin
          // C A S E  B.    mu < 10
          // START NEW TABLE AND CALCULATE P0 IF NECESSARY

          M := MAX(1, Trunc(mu));
          l := 0;
          p := Exp(-mu);
          q := p;
          p0 := p;

          // STEP U. UNIFORM SAMPLE FOR INVERSION METHOD
          repeat
            u := RandomGenerator;
            Result := 0;
            if (u <= p0) then
              Exit;
            // STEP T. TABLE COMPARISON UNTIL THE END PP(L) OF THE
            // PP-TABLE OF CUMULATIVE POISSON PROBABILITIES
            // (0.458=PP(9) FOR MU=10)
            if l <> 0 then
            begin
              j := 1;
              if (u > 0.458) then
                j := Min(l, M);
              for k := j to l do
                if (u <= pp[k]) then
                begin
                  Result := k;
                  Exit;
                end;

              if l = 35 then
                Continue;
            end;

            // STEP C. CREATION OF NEW POISSON PROBABILITIES P
            // AND THEIR CUMULATIVES Q = PP(K)

            l := l + 1; // 150
            for k := l to 35 do
            begin
              p := p * mu / k;
              q := q + p;
              pp[k] := q;
              if (u <= q) then
              begin
                Result := k;
                Exit;
              end;
            end;
            l := 35;
          until False end;
        end;

//---------------------------------------------------------------

        function Random_Poisson(var mu: Extended): Int64;
        begin
          Result := Random_Poisson(mu, GenRandom);
        end;

end.