program Project1;

uses
  SysUtils,
  spe,
  MathUnit;

// Numerical recipes
function GammaLn(x: double): Double;
{ gibt den Log der vollständigen Gamma-Funktion zurück (x>0).
     Gamma(x) = integral ( t^(x-1) * exp(-t) dt ) (von 0 bis Unendlich)
  (Log, um Floating Point Underflow zu vermeiden). }
const
  stp  = 2.50662827465;
var
  xx,tmp,ser : extended;
begin
  if x<=0 then
    raise Exception.Create('Argument für GammaLn ist negativ.');
  if (x > 1) then begin
    xx  := x  - 1.0;
    tmp := xx + 5.5;
    tmp := (xx+0.5) * ln(tmp) - tmp;
    ser := 1.0 + 76.18009173    /(xx+1.0) - 86.50532033/(xx+2.0)
               + 24.01409822    /(xx+3.0) - 1.231739516/(xx+4.0)
               +  0.120858003E-2/(xx+5.0) - 0.536382E-5/(xx+6.0);
    result := tmp + ln(stp*ser);
  end else
  if (x < 1) then
    result := GammaLn(x+1.0) - ln(x)
  else
  if (x=1) then
    result := 0.0;
end;

const
  xmin = 1e-3;
  xmax = 5.0;

procedure Test;
var
  i: Integer;
  x, y_lazStats, y_numlib, y_numrecip: Double;
begin
  WriteLn('x':20, 'y(lazstats)':20, 'y(numlib)':20, 'y(Num.Recipies)':20);
  x := xmin;
  while (x <= xmax) do begin
    y_numlib := 1.0 - spelga(x);
    y_lazstats := MathUnit.gammaln(x);
    y_numrecip := gammaln(x);

    WriteLn(x:20:6, y_lazstats:20:5, y_numlib:20:5, y_numrecip:20:5);
    x := x*1.25
//    x := x + dx;
  end;
end;

begin
  WriteLn('GammaLn function');
  WriteLn;
  Test;
  ReadLn;
end.