{ Delphi Translation developed by ESB Consultancy
 http://www.esbconsult.com/
 mailto:support@esbconsult.com

 Also a part of ESBPCS: http://www.esbconsult.com/esbpcs2

 v3.0.4 28-Jan-2004
 - Synchronised versions with ESBPCS
 - Fixed bug in Random_Gamma thanks to Mykhalevych Kyryl.
 - Other minor improvements.

 v1.1.3 17-May-2002
 - Fixed bug Demo which showed wrong paramaters.

 v1.1.2 9-Jan-2002
 - Borland C++ Builder Support added

 v1.1.1 19-Sep-2001
 - Fixed bug in Multivariate Normal Distribution - thanks to Haris Baltzakis

 v1.1 19-Sep-2001
 - now tested with Delphi 6
 - Tester now built with Delphi 6 and ESBPCS v2

 v1.00 9-Jan-2001
 - Fixed a bug that Alan Miller reported in Random_T
 - Some Improvements to the TestRandom Project
 - Now also incorporated into our ESBPCS with the release of v1.5

 Public Beta 3 Release 23-Oct-2000
 - Fixed problem where Poisson would get divide by zero - incorrect initialisation
  on my part.

 Public Beta 2 Release 12-Oct-2000
 - Exit condition for Random_Normal modified slightly
 - Changed routines to be overloaded to allow you to pass in the Uniform
  Random Number Generator you want to use. If omitted then Delphi's
  Random Number Generator will be used.
 - Added routines to call the Delphi Random and the MRNG Random Generators
  in the form of TRandomGenFunction.
 - Added AMRandomRS Unit to store strings as Resource Strings. Replace this
  unit with other language translations. Feel free to send to us so
  we can share with all users.
 - Added a Test Utility - source is included but ESBPCS (Shareware) and
  EFD Systems' MRNG.PAS (Freeware) is required.
  ESBPCS: http://www.esbconsult.com.au/esbpcs.html
  EFD: http://efd.home.mindspring.com/tools.htm
 - Test Utility still doesn't test Multi-variate.

 Public Beta 1: 30-Sep-2000
}

{ A unit for random number generation from the following distributions:

 Distribution                    Function/subroutine name

 Normal (Gaussian)               Random_Normal
 Gamma                           Random_Gamma
 Chi-squared                     Random_ChiSq
 Exponential                     Random_Exponential
 Weibull                         Random_Weibull
 Beta                            Random_Beta
 t                               Random_T
 Multivariate normal             Random_MVNorm
 Generalized inverse Gaussian    Random_Inv_Gauss
 Poisson                         Random_Poisson
 Binomial                        Random_Binomial1   *
         Random_Binomial2   *
 Negative binomial               Random_neg_Binomial
 von Mises                       Random_von_Mises
 Cauchy                          Random_Cauchy

 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.

  N.B. 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

 Some of the code is adapted from Dagpunar's book:
 Dagpunar, J. 'Principles of random variate generation'
 Clarendon Press, Oxford, 1988.   ISBN 0-19-852202-9

 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.

 First Delphi Version by Glenn Crouch of ESB Consultancy
     info@esbconsult.com.au
     http://www.esbconsult.com.au/

 FORTRAN Version History
 Version 1.12, 2 July 2000
 Changes from version 1.01
 1. The random_order, random_Poisson & random_binomial routines have been
    replaced with more efficient routines.
 2. A routine, seed_random_number, has been added to seed the uniform random
    number generator.   This requires input of the required number of seeds
    for the particular compiler from a specified I/O unit such as a keyboard.
 3. Made compatible with Lahey's ELF90.
 4. Marsaglia & Tsang algorithm used for random_gamma when shape parameter > 1.

 Author: Alan Miller
     CSIRO Division of Mathematical & Information Sciences
     Private Bag 10, Clayton South MDC
     Clayton 3169, Victoria, Australia
 Phone: (+61) 3 9545-8036      Fax: (+61) 3 9545-8080
 e-mail: Alan.Miller @ vic.cmis.csiro.au

}
unit AMRandom;

interface

{.$Define UseMRNG}// Remove this if you don't have EFD's free MRNG.PAS

type
     TAMFloatVector = array of Extended;

type
     TRandomGenFunction = function: Extended;

     {: Used to call Delphi's Inbuilt Random Number Generator. Remember to use
      Randomize if you don't want repeatable sequences. }
function DelphiRandom: Extended;

{$IFDEF UseMRNG}
{: Used to call EFD System's MRNG Random Number Generator. }
function MRNGRandom: Extended;
{$ENDIF}

{ 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 function random_normal() 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: Extended; overload;
function Random_Normal (RandomGenerator: TRandomGenFunction): Extended; overload;

{ Adapted from Fortran 77 code from the book:
 Dagpunar, J. 'Principles of random variate generation'
 Clarendon Press, Oxford, 1988.   ISBN 0-19-852202-9

 FUNCTION 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;
function Random_Gamma (const Shape: Extended;
     RandomGenerator: TRandomGenFunction): Extended; overload;

{  	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;

{ Adapted from Fortran 77 code from the book:
  Dagpunar, J. 'Principles of random variate generation'
  Clarendon Press, Oxford, 1988.   ISBN 0-19-852202-9

  FUNCTION 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;

{ Adapted from Fortran 77 code from the book:
  Dagpunar, J. 'Principles of random variate generation'
  Clarendon Press, Oxford, 1988.   ISBN 0-19-852202-9

  FUNCTION 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;

{  Adapted from Fortran 77 code from the book:
  Dagpunar, J. 'Principles of random variate generation'
  Clarendon Press, Oxford, 1988.   ISBN 0-19-852202-9

  FUNCTION 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;

{  Adapted from Fortran 77 code from the book:
 Dagpunar, J. 'Principles of random variate generation'
  Clarendon Press, Oxford, 1988.   ISBN 0-19-852202-9

  N.B. An extra argument, ier, has been added to Dagpunar's routine

  SUBROUTINE 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;

{ Adapted from Fortran 77 code from the book:
  Dagpunar, J. 'Principles of random variate generation'
  Clarendon Press, Oxford, 1988.   ISBN 0-19-852202-9

  FUNCTION 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;

{  FUNCTION 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;

{**********************************************************************
  Translated to Fortran 90 by Alan Miller from:
       RANLIB

  Library of Fortran Routines for Random Number Generation

       Compiled and Written by:

        Barry W. Brown
         James Lovato

    Department of Biomathematics, Box 237
    The University of Texas, M.D. Anderson Cancer Center
    1515 Holcombe Boulevard
    Houston, TX      77030

  This work was supported by grant CA-16672 from the National Cancer Institute.

     GENerate BINomial random deviate

       Function

  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;
function Random_Binomial2 (const n: Int64; const pp: Extended;
     RandomGenerator: TRandomGenFunction): Int64; overload;

{ Adapted from Fortran 77 code from the book:
  Dagpunar, J. 'Principles of random variate generation'
  Clarendon Press, Oxford, 1988.   ISBN 0-19-852202-9

  FUNCTION 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 from:
  Dagpunar, J.S. (1990) `Sampling from the von Mises distribution via a
  comparison of random numbers', J. of Appl. Statist., 17, 165-168.

  Fortran 90 code by Alan Miller
  CSIRO Division of Mathematical & Information Sciences

  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;

{**********************************************************************
  Translated to Fortran 90 by Alan Miller from:
        RANLIB

  Library of Fortran Routines for Random Number Generation

     Compiled and Written by:

      Barry W. Brown
       James Lovato

      Department of Biomathematics, Box 237
      The University of Texas, M.D. Anderson Cancer Center
      1515 Holcombe Boulevard
      Houston, TX      77030

  This work was supported by grant CA-16672 from the National Cancer Institute.

     GENerate POIsson random deviate

         Function

  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

uses
     {$IFDEF UseMRNG}
     MRNG,
     {$ENDIF}
     Math, SysUtils,
     AMRandomRS;

const
     VSmall = MinDouble;
     VLarge = MaxDouble;

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

{$IFDEF UseMRNG}

function MRNGRandom: Extended;
begin
     Result := MRandom;
end;
{$ENDIF}

{ 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 function random_normal() 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;
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 (rsDivideByZero);

     Result := v / u;
end;

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

{ 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.

  Generates a random gamma deviate for shape parameter Shape >= 1.
}

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 (rsInvalidShape);

     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;

{ Adapted from Fortran 77 code from the book:
  Dagpunar, J. 'Principles of random variate generation'
  Clarendon Press, Oxford, 1988.   ISBN 0-19-852202-9

  FUNCTION 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 (rsInvalidShape);
     if (Shape < VSmall) then
          raise EMathError.Create (rsShapeTooSmall);

     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;

{ Adapted from Fortran 77 code from the book:
 Dagpunar, J. 'Principles of random variate generation'
 Clarendon Press, Oxford, 1988.   ISBN 0-19-852202-9

 FUNCTION 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;
begin
     if (Shape <= 0.0) then
          raise EMathError.Create (rsInvalidShape);

     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, DelphiRandom);
end;

{  	Generates a random variate from the chi-squared distribution with
 dff degrees of freedom.
}

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, DelphiRandom);
end;

{ Adapted from Fortran 77 code from the book:
  Dagpunar, J. 'Principles of random variate generation'
  Clarendon Press, Oxford, 1988.   ISBN 0-19-852202-9

  FUNCTION GENERATES A RANDOM VARIATE IN [0,INFINITY) FROM
  A NEGATIVE EXPONENTIAL DlSTRIBUTION WlTH DENSITY PROPORTIONAL
  TO EXP(-random_exponential), USING INVERSION.
}

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 (DelphiRandom);
end;

{     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;
     RandomGenerator: TRandomGenFunction): Extended;
begin
     if Abs (a) <= VSmall then
          raise EMathError.Create (rsInvalidValue);

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

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

{ Adapted from Fortran 77 code from the book:
  Dagpunar, J. 'Principles of random variate generation'
  Clarendon Press, Oxford, 1988.   ISBN 0-19-852202-9

  FUNCTION 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;
     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 (rsInvalidValue);

     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, DelphiRandom);
end;

{  Adapted from Fortran 77 code from the book:
  Dagpunar, J. 'Principles of random variate generation'
  Clarendon Press, Oxford, 1988.   ISBN 0-19-852202-9

  FUNCTION 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;
     RandomGenerator: TRandomGenFunction): Extended;
var
     s, c, a, f, g: Extended;
     r, x, v: Extended;
     Done: Boolean;
begin
     if (DF < 1) then
          raise EMathError.Create (rsInvalidValue);

     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 >= 3) 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, DelphiRandom);
end;

{  Adapted from Fortran 77 code from the book:
 Dagpunar, J. 'Principles of random variate generation'
  Clarendon Press, Oxford, 1988.   ISBN 0-19-852202-9

  N.B. An extra argument, ier, has been added to Dagpunar's routine

  SUBROUTINE 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;
     RandomGenerator: TRandomGenFunction);
var
     i, j, m, n, n2, fn: Integer;
     y, v: Extended;
begin
     n := Length (h);
     if (n = 0) then
          raise EMathError.Create (rsVectorEmpty);

     fn := n * (n + 1) div 2;
     if (Length (D) < fn) then
          raise EMathError.Create (rsVectorTooSmall);

     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, DelphiRandom);
end;

{ Adapted from Fortran 77 code from the book:
  Dagpunar, J. 'Principles of random variate generation'
  Clarendon Press, Oxford, 1988.   ISBN 0-19-852202-9

  FUNCTION 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;
     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 (rsInvalidValue);

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

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

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

     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 (rsInvalidValue);

     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, DelphiRandom);
end;

{ 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;
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 (rsInvalidValue);

     // 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 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;
     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, DelphiRandom);
end;
{**********************************************************************
  Translated to Fortran 90 by Alan Miller from:
       RANLIB

  Library of Fortran Routines for Random Number Generation

       Compiled and Written by:

        Barry W. Brown
         James Lovato

    Department of Biomathematics, Box 237
    The University of Texas, M.D. Anderson Cancer Center
    1515 Holcombe Boulevard
    Houston, TX      77030

  This work was supported by grant CA-16672 from the National Cancer Institute.

     GENerate BINomial random deviate

       Function

  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

!     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;
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);
end;

{ Adapted from Fortran 77 code from the book:
  Dagpunar, J. 'Principles of random variate generation'
  Clarendon Press, Oxford, 1988.   ISBN 0-19-852202-9

  FUNCTION 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;
     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 (rsInvalidValue);

     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 (rsInvalidValue);

     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, DelphiRandom);
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;

{     Algorithm VMD from:
  Dagpunar, J.S. (1990) `Sampling from the von Mises distribution via a
  comparison of random numbers', J. of Appl. Statist., 17, 165-168.

  Fortran 90 code by Alan Miller
  CSIRO Division of Mathematical & Information Sciences

  Arguments:
  k (real)        parameter of the von Mises distribution.
}

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 (rsInvalidValue);

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

     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, DelphiRandom);
end;

{      Generate a random deviate from the standard Cauchy distribution
}

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 (DelphiRandom);
end;

{**********************************************************************
  Translated to Fortran 90 by Alan Miller from:
        RANLIB

  Library of Fortran Routines for Random Number Generation

     Compiled and Written by:

      Barry W. Brown
       James Lovato

      Department of Biomathematics, Box 237
      The University of Texas, M.D. Anderson Cancer Center
      1515 Holcombe Boulevard
      Houston, TX      77030

  This work was supported by grant CA-16672 from the National Cancer Institute.

     GENerate POIsson random deviate

         Function

  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;
     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, DelphiRandom);
end;

end.

