{$N+}{$E-}
UNIT Gamma;

INTERFACE

FUNCTION gammln(xx: real): real;
  {Vypocet ln (Gamma(xx))}
PROCEDURE gser(a,x: real;
    VAR gamser,gln: real);
   {Vypocet Gammp(a,x) = 1/Gamma(a)* Int_0^x exp(-t) t^(a-1) dt
                               radou - vhodne pro x < a + 1}
PROCEDURE gcf(a,x: real;
   VAR gammcf,gln: real);
   {Vypocet Gammq(a,x) = 1/Gamma(a)* Int_x^inf
             nekonecnymi zlomky - vhodne pro x > a + 1}
FUNCTION gammq(a,x: real): real;
    {Vypocet Gammq(a,x) pro vsechna x s vyuzitim GSER a GCF}
FUNCTION gammp(a,x: real): real;
    {Vypocet Gammp(a,x) pro vsechna x s vyuzitim GSER a GCF}
FUNCTION betacf(a,b,x: real): real;
    {Pomocny vypocet pro BETAI pomoci nekonecnych zlomku}
FUNCTION betai(a,b,x: real): real;
    {Vypocet neuplne BETA funkce I_x(a,b) =
               1/B(a,b)*Int_0^x t^(a-1) * (1 - t)^(b-1)  dt}


IMPLEMENTATION

FUNCTION gammln(xx: real): real;
CONST
   stp  =  2.50662827465;
VAR
   x,tmp,ser: double;
BEGIN
   x := xx-1.0;
   tmp := x+5.5;
   tmp := (x+0.5)*ln(tmp)-tmp;
   ser := 1.0+76.18009173/(x+1.0)-86.50532033/(x+2.0)+24.01409822/(x+3.0)
            -1.231739516/(x+4.0)+0.120858003e-2/(x+5.0)-0.536382e-5/(x+6.0);
   gammln := tmp+ln(stp*ser)
END;
{end GAMMLN}
PROCEDURE gser(a,x: real;
    VAR gamser,gln: real);
LABEL 99;
CONST
   itmax = 100;
   eps = 3.0e-7;
VAR
   n: integer;
   sum,del,ap: real;
BEGIN
   gln := gammln(a);
   IF x <= 0.0 THEN BEGIN
      IF x < 0.0 THEN BEGIN
         writeln('pause in GSER - x less than 0');
         readln
      END;
      gamser := 0.0
   END
   ELSE BEGIN
      ap := a;
      sum := 1.0/a;
      del := sum;
      FOR n := 1 TO itmax DO BEGIN
         ap := ap+1.0;
         del := del*x/ap;
         sum := sum+del;
         IF abs(del) < abs(sum)*eps THEN GOTO 99
      END;
      writeln('pause in GSER - a too large, itmax too small');
      readln;
99:   gamser := sum*exp(-x+a*ln(x)-gln)
   END
END;
{end GSER}
PROCEDURE gcf(a,x: real;
   VAR gammcf,gln: real);
LABEL 99;
CONST
   itmax = 100;
   eps = 3.0e-7;
VAR
   n: integer;
   gold,g,fac,b1,b0,anf,ana,an,a1,a0: real;
BEGIN
   gln := gammln(a);
   gold := 0.0;
   a0 := 1.0;
   a1 := x;
   b0 := 0.0;
   b1 := 1.0;
   fac := 1.0;
   FOR n := 1 TO itmax DO BEGIN
      an := 1.0*n;
      ana := an-a;
      a0 := (a1+a0*ana)*fac;
      b0 := (b1+b0*ana)*fac;
      anf := an*fac;
      a1 := x*a0+anf*a1;
      b1 := x*b0+anf*b1;
      IF a1 <> 0.0 THEN BEGIN
         fac := 1.0/a1;
         g := b1*fac;
         IF abs((g-gold)/g) < eps THEN GOTO 99;
         gold := g
      END
   END;
   writeln('pause in GCF - a too large, itmax too small');
   readln;
99:
   gammcf := exp(-x+a*ln(x)-gln)*g
END;
{end GCF}
FUNCTION gammp(a,x: real): real;
VAR
   gamser,gammcf,gln: real;
BEGIN
   IF (x < 0.0) OR (a <= 0.0) THEN BEGIN
      writeln('pause in GAMMP - invalid arguments');
      readln
   END;
   IF x < a+1.0 THEN BEGIN
      gser(a,x,gamser,gln);
      gammp := gamser
   END
   ELSE BEGIN
      gcf(a,x,gammcf,gln);
      gammp := 1.0-gammcf
   END
END;
{end GAMMP}
FUNCTION gammq(a,x: real): real;
VAR
   gamser,gammcf,gln: real;
BEGIN
   IF (x < 0.0) OR (a <= 0.0) THEN BEGIN
      writeln('pause in GAMMQ - invalid arguments');
      readln
   END;
   IF x < a+1.0 THEN BEGIN
      gser(a,x,gamser,gln);
      gammq := 1.0-gamser
   END
   ELSE BEGIN
      gcf(a,x,gammcf,gln);
      gammq := gammcf
   END
END;
{end GAMMQ}
FUNCTION betacf(a,b,x: real): real;
LABEL 99;
CONST
   itmax = 100;
   eps = 3.0e-7;
VAR
   tem,qap,qam,qab,em,d: real;
   bz,bpp,bp,bm,az,app: real;
   am,aold,ap: real;
   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 GOTO 99
   END;
   writeln('pause in BETACF');
   writeln('a or b too big, or itmax too small');
   readln;
99:
   betacf := az
END;
{END betacf}
FUNCTION betai(a,b,x: real): real;
VAR
   bt: real;
BEGIN
   IF (x < 0.0) OR (x > 1.0) THEN BEGIN
      writeln('pause in routine BETAI');
      readln
   END;
   IF (x = 0.0) OR (x = 1.0) THEN bt := 0.0
   ELSE bt := exp(gammln(a+b)-gammln(a)-gammln(b)
      +a*ln(x)+b*ln(1.0-x));
   IF x < (a+1.0)/(a+b+2.0) THEN
      betai := bt*betacf(a,b,x)/a
   ELSE
      betai := 1.0-bt*betacf(b,a,1.0-x)/b
END;

END.