program Project1;

uses
  SysUtils,
  spe,
  MathUnit;

// Numerical recipes

function BetaCf(a,b,x: Double) : Double;
{ Kettenbruch-Entwicklung der unvollständigen Beta-Funktion }
const
  ItMax = 100;
  EPS = 1E-9;
var
  tem, em, d            : Double;
  qap, qam, qab         : Double;
  bz, bpp, BP, bm       : Double;
  az, app, am, aold, ap : Double;
  m                     : integer;
begin
  am  := 1.0;
  bm  := 1.0;
  az  := 1.0;
  qab := a+b;
  qap := a+1.0;
  qam := a-1.0;
  bz  := 1.0 - qab*x/qap;
  for m:=1 to ItMax do begin
    em  := m;
    tem := em+em;
    d   := em*(b-m)*x/((qam+tem)*(a+tem));
    ap  := az + d*am;
    BP  := bz + d*bm;
    d   := -(a+em)*(qab+em)*x/((a+tem)*(qap+tem));
    app := ap + d*az;
    bpp := BP + d*bz;
    aOld:= az;
    am  := ap/bpp;
    bm  := BP/bpp;
    az  := app/bpp;
    bz  := 1.0;
    if abs(az-aOld) < EPS*abs(az) then begin
      result := az;
      Exit;
    end;
  end;
  raise Exception.Create('Kettenbruchentwicklung der unvollst.Beta-Funktion divergiert.');
  // a oder b zu groß, oder ItMax zu klein
end;

function BetaI(a,b,x: Double): Double;
{ berechnet die Unvollständige Beta-Funktion Ix(a,b) = Bx(a,b)/B(a,b),
                  /x  a-1         b-1
  wobei Bx(a,b) = |  t     * (1-t)    dt
                  /0
}
var
  bt : Double;
begin
  if (x<0.0) then
    raise EMathError.Create('Argument der unvollst.Beta-Funktion ist <0.');
  if (x>1.0) then
    raise EMathError.Create('Argument der unvollst.Beta-Funktion ist >1.');
  if (x=0.0) or (x=1.0) then
    bt := 0.0
  else
    bt := exp(GammaLn(a+b) - GammaLn(a) - GammaLn(b) + a*ln(x) + b*ln(1.0-x));
  if x < (a+1.0)/(a+b+2.0)
    then result := bt * BetaCf(a,b,x)/a
    else result := 1.0 - bt*BetaCf(b,a,1.0-x)/b;
end;


const
  xmin = 0.0;
  xmax = 1.0;
  dx = 0.2;

procedure Test(a, b: Double);
var
  i: Integer;
  x, y_lazStats, y_numlib, y_numrecip: Double;
begin
  WriteLn('a = ', a:0:3, ' b = ', b:0:3);
  WriteLn;
  WriteLn('x':20, 'y(lazstats)':20, 'y(numlib)':20, 'y(Num.Recip)':20);
  x := xmin;
  while (x <= xmax) do begin
    y_numlib := spe.betai(a, b, x);
    y_lazstats := mathunit.betai(a, b, x);
    y_numrecip := betai(a, b, x);

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

begin
  WriteLn('incomplete beta function');
  WriteLn;
  Test(0.5, 0.5);
  Test(1.0, 2.0);
  Test(2.0, 1.0);
  ReadLn;
end.