{  (C) Copr. 1986-92 Numerical Recipes Software .5-28.}

{ Гамма-функция и связанные с ней функции }
unit uGamma;

interface

{ Натуральный логарифм Гамма-функции }
function GammLn(x:Double):Double;

{ Неполная Гамма-функция }
function GammP(a,x:Double):Double;

implementation

const
  cEps=3.0e-14;    { 3.0e-7;  - for Real }
  cFPMin=1.0e-300; { 1.0e-30; - for Real }
  cItMax=1000;     { 100;     - for Real }

{ внутренняя процедура }
procedure GCF(var gammcf:Double; a,x:Double; var gln:Double);
  var i:word; an,b,c,d,del,h:Extended;
begin
  gln:=gammln(a);
  b:=x+(1.0-a);
  c:=1.0/cFPMin;
  d:=1.0/b;
  h:=d;
  i:=1;
  repeat
    an:=-i*(i-a);
    b:=b+2.0;
    d:=an*d+b;
    if abs(d)<cFPMin then d:=cFPMin;
    c:=b+an/c;
    if abs(c)<cFPMin then c:=cFPMin;
    d:=1.0/d;
    del:=d*c;
    h:=h*del;
    Inc(i);
{$IfOpt R+}
    if (i>CItMAx) then RunError(201);
{$EndIf R+}
  until (abs(del-1) < cEps);
  gammcf:=exp(-x+a*ln(x)-gln)*h;
end;

{ внутренняя процедура }
procedure GSer(var gamser:Double; a,x:Double; var gln:Double);
  var ap,del,sum: Extended;
{$IfOpt R+}
  n:word;
{$EndIf R+}
begin
  gln:=GammLn(a);
  if x<=0.0 then begin
{$IfOpt R+}
    if x<0.0 then RunError(201);
{$EndIf R+}
    gamser:=0.0;
    Exit;
  end;
  ap:=a;
  sum:=1.0/a;
  del:=sum;
{$IfOpt R+}
  n:=cItMAx;
{$EndIf R+}
  repeat
    ap:=ap+1.0;
    del:=del*x/ap;
    sum:=sum+del;
{$IfOpt R+}
    Dec(n);
    if (n<=0) then RunError(201);
{$EndIf R+}
  until abs(del) < (abs(sum)*cEps);
  gamser:=sum*exp(-x+a*ln(x)-gln);
end;

function GammLn(x:Double):Double;
  const
    cof:array[1..6] of Extended=(76.18009172947146,
				-86.50532032941677,
				 24.01409824083091,
				- 1.231739572450155,
				  0.1208650973866179e-2,
				- 0.5395239384953e-5);
   stp=2.5066282746310005;

  var j:byte;  ser,tmp,y:Extended;
begin
  y:=x;
  tmp:=x+5.5;
  tmp:=(x+0.5)*ln(tmp)-tmp;
  ser:=1.000000000190015;
  for j:=1 to 6 do begin
    y:=y+1.0;
    ser:=ser+cof[j]/y;
  end;
  GammLn:=tmp+ln(stp*ser/x);
end;

function GammP(a,x: Double):Double;
  var gamser,gammcf,gln: Double;
begin
{$IfOpt R+}
  if (x<0) or (a<=0) then RunError(201);
{$EndIf R+}
  if x < (a+1) then begin
    GSer(gamser,a,x,gln);
    GammP:=gamser;
  end else begin
    GCF(gammcf,a,x,gln);
    GammP:=1.0-gammcf;
  end;
end;

end.