LnGamma 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 |
X | Value to process. |
Category
Arithmetic Routines for FloatsImplementation
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; |
|