LnGamma Function
Logarithm to base e of the gamma function.

Unit
QESBPCSMath

Declaration
Function LnGamma(const X: Extended): Extended;

Description
Defined for all positive values of X.

Accurate to about 14 digits.

Programmer: Alan Miller - developed for Fortan 77, converted with permission.

Parameters
Value to process.

Category
Arithmetic Routines for Floats

Implementation

function LnGamma (const X: Extended): Extended;
const
     A1 = -4.166666666554424E-02;
     A2 = 2.430554511376954E-03;
     A3 = -7.685928044064347E-04;
     A4 = 5.660478426014386E-04;
var
     Temp, Arg, Product: Extended;
     Reflect: Boolean;
begin
     //  lngamma is not defined if x = 0 or a negative integer.
     if FloatIsZero (X) or (FloatIsNegative (X) and SameFloat (X, Int (X))) then
          raise EMathError.Create (rsNotDefinedForValue);

     // 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 (ESBPi * X);
          Result := Ln (ESBPi / Temp) - Result;
     end;
End;


HTML generated by Time2HELP
http://www.time2help.com