(*-------------------------------------------------------------------------*)
(*               ALGama  -- Logarithm of Gamma Distribution                *)
(*-------------------------------------------------------------------------*)

FUNCTION ALGama( Arg : REAL ) : REAL;

(*-------------------------------------------------------------------------*)
(*                                                                         *)
(*       Function:  ALGama                                                 *)
(*                                                                         *)
(*       Purpose:  Calculates Log (base E) of Gamma function               *)
(*                                                                         *)
(*       Calling Sequence:                                                 *)
(*                                                                         *)
(*            Val := ALGama( Arg )                                         *)
(*                                                                         *)
(*                 Arg   --- Gamma distribution parameter (Input)          *)
(*                 Val   --- output Log Gamma value                        *)
(*                                                                         *)
(*       Calls:  None                                                      *)
(*                                                                         *)
(*       Called By:  Many (CDBeta, etc.)                                   *)
(*                                                                         *)
(*       Remarks:                                                          *)
(*                                                                         *)
(*            Minimax polynomial approximations are used over the          *)
(*            intervals [-inf,0], [0,.5], [.5,1.5], [1.5,4.0],             *)
(*            [4.0,12.0], [12.0,+inf].                                     *)
(*                                                                         *)
(*            See Hart et al, "Computer Approximations",                   *)
(*            Wiley(1968), p. 130F, and also                               *)
(*                                                                         *)
(*            Cody and Hillstrom, "Chebyshev approximations for            *)
(*            the natural logarithm of the Gamma function",                *)
(*            Mathematics of Computation, 21, April, 1967, P. 198F.        *)
(*                                                                         *)
(*                                                                         *)
(*            There are some machine-dependent constants used --           *)
(*                                                                         *)
(*               Rmax   --- Largest value for which ALGama                 *)
(*                          can be safely computed.                        *)
(*               Rinf   --- Largest floating-point number.                 *)
(*               Zeta   --- Smallest floating-point number                 *)
(*                          such that (1 + Zeta) = 1.                      *)
(*                                                                         *)
(*-------------------------------------------------------------------------*)

VAR
   Rarg        : REAL;
   Alinc       : REAL;
   Scale       : REAL;
   Top         : REAL;
   Bot         : REAL;
   Frac        : REAL;
   Algval      : REAL;

   I           : INTEGER;
   Iapprox     : INTEGER;
   Iof         : INTEGER;
   Ilo         : INTEGER;
   Ihi         : INTEGER;

   Qminus      : BOOLEAN;
   Qdoit       : BOOLEAN;


(* Structured *) CONST

   P           : ARRAY [ 1 .. 29 ] OF REAL =
                 ( 4.12084318584770E+00 ,
                   8.56898206283132E+01 ,
                   2.43175243524421E+02 ,
                  -2.61721858385614E+02 ,
                  -9.22261372880152E+02 ,
                  -5.17638349802321E+02 ,
                  -7.74106407133295E+01 ,
                  -2.20884399721618E+00 ,

                   5.15505761764082E+00 ,
                   3.77510679797217E+02 ,
                   5.26898325591498E+03 ,
                   1.95536055406304E+04 ,
                   1.20431738098716E+04 ,
                  -2.06482942053253E+04 ,
                  -1.50863022876672E+04 ,
                  -1.51383183411507E+03 ,

                  -1.03770165173298E+04 ,
                  -9.82710228142049E+05 ,
                  -1.97183011586092E+07 ,
                  -8.73167543823839E+07 ,
                   1.11938535429986E+08 ,
                   4.81807710277363E+08 ,
                  -2.44832176903288E+08 ,
                  -2.40798698017337E+08 ,

                   8.06588089900001E-04 ,
                  -5.94997310888900E-04 ,
                   7.93650067542790E-04 ,
                  -2.77777777688189E-03 ,
                   8.33333333333330E-02   );

   Q           : ARRAY [ 1 .. 24 ] OF REAL =
                 ( 1.00000000000000E+00 ,
                   4.56467718758591E+01 ,
                   3.77837248482394E+02 ,
                   9.51323597679706E+02 ,
                   8.46075536202078E+02 ,
                   2.62308347026946E+02 ,
                   2.44351966250631E+01 ,
                   4.09779292109262E-01 ,

                   1.00000000000000E+00 ,
                   1.28909318901296E+02 ,
                   3.03990304143943E+03 ,
                   2.20295621441566E+04 ,
                   5.71202553960250E+04 ,
                   5.26228638384119E+04 ,
                   1.44020903717009E+04 ,
                   6.98327414057351E+02 ,

                   1.00000000000000E+00 ,
                  -2.01527519550048E+03 ,
                  -3.11406284734067E+05 ,
                  -1.04857758304994E+07 ,
                  -1.11925411626332E+08 ,
                  -4.04435928291436E+08 ,
                  -4.35370714804374E+08 ,
                  -7.90261111418763E+07   );

BEGIN (* ALGama *)


                                   (* Initialize *)

      Algval := Rinf;
      Scale  := 1.0;
      Alinc  := 0.0;
      Frac   := 0.0;
      Rarg   := Arg;
      Iof    := 1;
      Qminus := FALSE;
      Qdoit  := TRUE;
                                   (* Adjust for negative argument *)

      IF( Rarg < 0.0 ) THEN
         BEGIN

            Qminus := TRUE;
            Rarg   := -Rarg;
            Top    := Int( Rarg );
            Bot    := 1.0;

            IF( ( INT( Top / 2.0 ) * 2.0 ) = 0.0 ) THEN Bot := -1.0;

            Top    := Rarg - Top;

            IF( Top = 0.0 ) THEN
               Qdoit := FALSE
            ELSE
               BEGIN
                  Frac   := Bot * PI / SIN( Top * PI );
                  Rarg   := Rarg + 1.0;
                  Frac   := LN( ABS( Frac ) );
               END;

         END;
                                   (* Choose approximation interval *)
                                   (* based upon argument range     *)

      IF( Rarg = 0.0 ) THEN
         Qdoit := FALSE

      ELSE IF( Rarg <= 0.5 ) THEN
         BEGIN

            Alinc  := -LN( Rarg );
            Scale  := Rarg;
            Rarg   := Rarg + 1.0;

            IF( Scale < Zeta ) THEN
               BEGIN
                  Algval := Alinc;
                  Qdoit  := FALSE;
               END;

         END

      ELSE IF ( Rarg <= 1.5 ) THEN
         Scale := Rarg - 1.0

      ELSE IF( Rarg <= 4.0 ) THEN
         BEGIN
            Scale  := Rarg - 2.0;
            Iof    := 9;
         END

      ELSE IF( Rarg <= 12.0 ) THEN
         Iof := 17

      ELSE IF( Rarg <= RMAX   ) THEN
         BEGIN

            Alinc  := ( Rarg - 0.5 ) * LN( Rarg ) - Rarg + Xln2sp;
            Scale  := 1.0 / Rarg;
            Rarg   := Scale * Scale;

            Top    := P[ 25 ];

            FOR I := 26 TO 29 DO
               Top    := Top * Rarg + P[ I ];

            Algval := Scale * Top + Alinc;

            Qdoit  := FALSE;

         END;

                                   (* Common evaluation code for Arg <= 12. *)
                                   (* Horner's method is used, which seems  *)
                                   (* to give better accuracy than          *)
                                   (* continued fractions.                  *)

      IF Qdoit THEN
         BEGIN

            Ilo    := Iof + 1;
            Ihi    := Iof + 7;

            Top    := P[ Iof ];
            Bot    := Q[ Iof ];

            FOR I := Ilo TO Ihi DO
               BEGIN
                  Top    := Top * Rarg + P[ I ];
                  Bot    := Bot * Rarg + Q[ I ];
               END;

            Algval := Scale * ( Top / Bot ) + Alinc;

         END;

      IF( Qminus ) THEN Algval := Frac - Algval;

      ALGama := Algval;

END   (* ALGama *);