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