Z := A * Sin(PI * Z);
            if Z = 0.0 then
              begin
                LnGamma := DefaultVal(FN_OVERFLOW, SgnGam * MAXNUM);
                Exit;
              end;
            Z := LNPI - Ln(Z) - StirfL(A);
          end
        else
          Z := StirfL(X);
        LnGamma := Z;
      end
    else if X < 13.0 then
      begin
        Z := 1.0;
        X1 := X;
        while X1 >= 3 do
          begin
            X1 := X1 - 1.0;
            Z := Z * X1;
          end;
        while X1 < 2.0 do
          begin
            if Abs(X1) <= 0.03125 then
              begin
                LnGamma := Ln(Abs(GamSmall(X1, Z)));
                Exit;
              end;
            Z := Z / X1;
            X1 := X1 + 1.0;
          end;
        if Z < 0.0 then Z := - Z;
        if X1 = 2.0 then
          LnGamma := Ln(Z)
        else
          begin
            X1 := X1 - 2.0;
            LnGamma := X1 * PolEvl(X1, P, 6) / P1Evl(X1, Q, 7) + Ln(Z);
          end;
      end
    else
      LnGamma := StirfL(X);
  end;

  function ApproxLnGamma(Z : Complex) : Complex;
  { This is the approximation used in the National Bureau of
    Standards "Table of the Gamma Function for Complex Arguments,"
    Applied Mathematics Series 34, 1954. The NBS table was created
    using this approximation over the area 9 < Re(z) < 10 and
    0 < Im(z) < 10. Other table values were computed using the
    relationship:
        _                   _
    ln | (z+1) = ln z + ln | (z) }

  const
    C : array[1..8] of Float =
    (8.33333333333333E-02, - 2.77777777777778E-03,
     7.93650793650794E-04, - 5.95238095238095E-04,
     8.41750841750842E-04, - 1.91752691752692E-03,
     6.41025641025641E-03, - 2.95506535947712E-02);
  var
    I : Integer;
    Powers : array[1..8] of Complex;
    Temp1, Temp2, Sum : Complex;
  begin
    Temp1 := Log(Z);                       { Ln(Z) }
    Temp2 := Cmplx(Z.X - 0.5, Z.Y);        { Z - 0.5 }
    Sum := CMult(Temp1, Temp2);            { (Z - 0.5)*Ln(Z) }
    Sum := CSub(Sum, Z);                   { (Z - 0.5)*ln(Z) - Z }
    Sum.X := Sum.X + LN2PIDIV2;
    Temp1 := C_one;
    Powers[1] := CDiv(Temp1, Z);           { Z^(-1) }
    Temp2 := CMult(Powers[1], Powers[1]);  { Z^(-2) }
    for I := 2 to 8 do
      Powers[I] := CMult(Powers[I - 1], Temp2);
    for I := 8 downto 1 do
      Sum := CAdd(Sum, Cmplx(C[I] * Powers[I].X, C[I] * Powers[I].Y));
    ApproxLnGamma := Sum;
  end;

  function LnGamma(Z : Complex) : Complex; overload;
  var
    A, LnZ, Temp : Complex;
  begin
    if (Z.X <= 0.0) and (Z.Y = 0.0) then
      if (Int(Z.X - 1E-8) - Z.X) = 0.0 then  { Negative integer? }
        begin
          MathError := FN_SING;
          LnGamma := C_infinity;
          Exit
        end;
    if Z.Y < 0.0 then            { 3rd or 4th quadrant? }
      begin
        A := LnGamma(CConj(Z));  { Try again in 1st or 2nd quadrant }
        LnGamma := CConj(A);     { Left this out! 1/3/91 }
      end
    else
      begin
        if Z.X < 9.0 then  { "left" of NBS table range }
          begin
            LnZ := Log(Z);
            Z := Cmplx(Z.X + 1.0, Z.Y);
            Temp := LnGamma(Z);
            LnGamma := CSub(Temp, LnZ)
          end
        else
          LnGamma := ApproxLnGamma(Z)  { NBS table range: 9 < Re(z) < 10 }
      end
  end;

  function IGamma(A, X : Float) : Float;
  var
    Ans, Ax, C, R : Float;
  begin
    MathError := FN_OK;

    if (X <= 0.0) or (A <= 0.0) then
      begin
        IGamma := 0.0;
        Exit;
      end;

    if (X > 1.0) and (X > A) then
      begin
        IGamma := 1.0 - JGamma(A, X);
        Exit;
      end;

    Ax := A * Ln(X) - X - LnGamma(A);
    if Ax < MINLOG then
      begin
        IGamma := DefaultVal(FN_UNDERFLOW, 0.0);
        Exit;
      end;

    Ax := Exp(Ax);

    { Power series }
    R := A;
    C := 1.0;
    Ans := 1.0;

    repeat
      R := R + 1.0;
      C := C * X / R;
      Ans := Ans + C;
    until C / Ans <= MACHEP;

    IGamma := Ans * Ax / A;
  end;

  function JGamma(A, X : Float) : Float;
  var
    Ans, C, Yc, Ax, Y, Z, R, T,
    Pk, Pkm1, Pkm2, Qk, Qkm1, Qkm2 : Float;
  begin
    MathError := FN_OK;

    if (X <= 0.0) or (A <= 0.0) then
      begin
        JGamma := 1.0;
        Exit;
      end;

    if (X < 1.0) or (X < A) then
      begin
        JGamma := 1.0 - IGamma(A, X);
        Exit;
      end;

    Ax := A * Ln(X) - X - LnGamma(A);

    if Ax < MINLOG then
      begin
        JGamma := DefaultVal(FN_UNDERFLOW, 0.0);
        Exit;
      end;

    Ax := Exp(Ax);

    { Continued fraction }
    Y := 1.0 - A;
    Z := X + Y + 1.0;
    C := 0.0;
    Pkm2 := 1.0;
    Qkm2 := X;
    Pkm1 := X + 1.0;
    Qkm1 := Z * X;
    Ans := Pkm1 / Qkm1;

    repeat
      C := C + 1.0;
      Y := Y + 1.0;
      Z := Z + 2.0;
      Yc := Y * C;
      Pk := Pkm1 * Z - Pkm2 * Yc;
      Qk := Qkm1 * Z - Qkm2 * Yc;
      if Qk <> 0.0 then
        begin
          R := Pk / Qk;
          T := Abs((Ans - R) / R);
          Ans := R;
        end
      else
        T := 1.0;
      Pkm2 := Pkm1;
      Pkm1 := Pk;
      Qkm2 := Qkm1;
      Qkm1 := Qk;
      if Abs(Pk) > BIG then
        begin
          Pkm2 := Pkm2 / BIG;
          Pkm1 := Pkm1 / BIG;
          Qkm2 := Qkm2 / BIG;
          Qkm1 := Qkm1 / BIG;
        end;
    until T <= MACHEP;

    JGamma := Ans * Ax;
  end;

  function Fact(N : Integer) : Float;
  begin
    MathError := FN_OK;
    if N < 0 then
      Fact := DefaultVal(FN_DOMAIN, 1.0)
    else if N > MAXFAC then
      Fact := DefaultVal(FN_OVERFLOW, MAXNUM)
    else if N <= NFACT then
      Fact := FactArray[N]
    else
      Fact := Gamma(N + 1);
  end;

  function Binomial(N, K : Integer) : Float;
  var
    I, N1 : Integer;
    Prod : Float;
  begin
    MathError := FN_OK;
    if K < 0 then
      Binomial := 0.0
    else if (K = 0) or (K = N) then
      Binomial := 1.0
    else if (K = 1) or (K = N - 1) then
      Binomial := N
    else
      begin
        if K > N - K then K := N - K;
        N1 := Succ(N);
        Prod := N;
        for I := 2 to K do
          Prod := Prod * ((N1 - I) / I);
        Binomial := Int(0.5 + Prod);
      end;
  end;

  function Beta(X, Y : Float) : Float;
  { Computes Beta(X, Y) = Gamma(X) * Gamma(Y) / Gamma(X + Y) }
  var
    Lx, Ly, Lxy : Float;
    SgnBeta : Integer;
  begin
    MathError := FN_OK;
    SgnBeta := SgnGamma(X) * SgnGamma(Y) * SgnGamma(X + Y);
    Lxy := LnGamma(X + Y);
    if MathError <> FN_OK then
      begin
        Beta := 0.0;
        Exit;
      end;
    Lx := LnGamma(X);
    if MathError <> FN_OK then
      begin
        Beta := SgnBeta * MAXNUM;
        Exit;
      end;
    Ly := LnGamma(Y);
    if MathError <> FN_OK then
      begin
        Beta := SgnBeta * MAXNUM;
        Exit;
      end;
    Beta := SgnBeta * Exp(Lx + Ly - Lxy);
  end;

  function PSeries(A, B, X : Float) : Float;
  { Power series for incomplete beta integral. Use when B*X is small }
  var
    S, T, U, V, T1, Z, Ai : Float;
    N : Integer;
  begin
    Ai := 1.0 / A;
    U := (1.0 - B) * X;
    V := U / (A + 1.0);
    T1 := V;
    T := U;
    N := 2;
    S := 0.0;
    Z := MACHEP * Ai;
    while Abs(V) > Z do
      begin
        U := (N - B) * X / N;
        T := T * U;
        V := T / (A + N);
        S := S + V;
        N := N + 1;
      end;
    S := S + T1;
    S := S + Ai;

    U := A * Ln(X);
    if (A + B < MAXGAM) and (Abs(U) < MAXLOG) then
      begin
        T := Gamma(A + B) / (Gamma(A) * Gamma(B));
        S := S * T * Power(X, A);
      end
    else
      begin
        T := LnGamma(A + B) - LnGamma(A) - LnGamma(B) + U + Ln(S);
        if T < MINLOG then
          S := 0.0
        else
          S := Exp(T);
      end;
    PSeries := S;
  end;

  function CFrac1(A, B, X : Float) : Float;
  { Continued fraction expansion #1 for incomplete beta integral }
  var
    Xk, Pk, Pkm1, Pkm2, Qk, Qkm1, Qkm2,
    K1, K2, K3, K4, K5, K6, K7, K8,
    R, T, Ans, Thresh : Float;
    N : Integer;
  label
    CDone;
  begin
    K1 := A;
    K2 := A + B;
    K3 := A;
    K4 := A + 1.0;
    K5 := 1.0;
    K6 := B - 1.0;
    K7 := K4;
    K8 := A + 2.0;

    Pkm2 := 0.0;
    Qkm2 := 1.0;
    Pkm1 := 1.0;
    Qkm1 := 1.0;
    Ans := 1.0;
    R := 1.0;
    N := 0;
    Thresh := 3.0 * MACHEP;

    repeat
      Xk := - (X * K1 * K2) / (K3 * K4);
      Pk := Pkm1 + Pkm2 * Xk;
      Qk := Qkm1 + Qkm2 * Xk;
      Pkm2 := Pkm1;
      Pkm1 := Pk;
      Qkm2 := Qkm1;
      Qkm1 := Qk;

      Xk := (X * K5 * K6) / (K7 * K8);
      Pk := Pkm1 + Pkm2 * Xk;
      Qk := Qkm1 + Qkm2 * Xk;
      Pkm2 := Pkm1;
      Pkm1 := Pk;
      Qkm2 := Qkm1;
      Qkm1 := Qk;

      if Qk <> 0.0 then R := Pk / Qk;

      if R <> 0.0 then
        begin
          T := Abs((Ans - R) / R);
          Ans := R;
        end
      else
        T := 1.0;

      if T < Thresh then goto CDone;

      K1 := K1 + 1.0;
      K2 := K2 + 1.0;
      K3 := K3 + 2.0;
      K4 := K4 + 2.0;
      K5 := K5 + 1.0;
      K6 := K6 - 1.0;
      K7 := K7 + 2.0;
      K8 := K8 + 2.0;

      if Abs(Qk) + Abs(Pk) > BIG then
        begin
          Pkm2 := Pkm2 * BIGINV;
          Pkm1 := Pkm1 * BIGINV;
          Qkm2 := Qkm2 * BIGINV;
          Qkm1 := Qkm1 * BIGINV;
        end;

      if (Abs(Qk) < BIGINV) or (Abs(Pk) < BIGINV) then
        begin
          Pkm2 := Pkm2 * BIG;
          Pkm1 := Pkm1 * BIG;
          Qkm2 := Qkm2 * BIG;
          Qkm1 := Qkm1 * BIG;
        end;
      N := N + 1;
    until N > 400;
    CFrac1 := DefaultVal(FN_PLOSS, 0.0);

CDone:
    CFrac1 := Ans;
  end;

  function CFrac2(A, B, X : Float) : Float;
  { Continued fraction expansion #2 for incomplete beta integral }
  var
    Xk, Pk, Pkm1, Pkm2, Qk, Qkm1, Qkm2,
    K1, K2, K3, K4, K5, K6, K7, K8,
    R, T, Z, Ans, Thresh : Float;
    N : Integer;
  label
    CDone;
  begin
    K1 := A;
    K2 := B - 1.0;
    K3 := A;
    K4 := A + 1.0;
    K5 := 1.0;
    K6 := A + B;
    K7 := A + 1.0;
    K8 := A + 2.0;

    Pkm2 := 0.0;
    Qkm2 := 1.0;
    Pkm1 := 1.0;
    Qkm1 := 1.0;
    Z := X / (1.0 - X);
    Ans := 1.0;
    R := 1.0;
    N := 0;
    Thresh := 3.0 * MACHEP;

    repeat
      Xk := - (Z * K1 * K2) / (K3 * K4);
      Pk := Pkm1 + Pkm2 * Xk;
      Qk := Qkm1 + Qkm2 * Xk;
      Pkm2 := Pkm1;
      Pkm1 := Pk;
      Qkm2 := Qkm1;
      Qkm1 := Qk;

      Xk := (Z * K5 * K6) / (K7 * K8);
      Pk := Pkm1 + Pkm2 * Xk;
      Qk := Qkm1 + Qkm2 * Xk;
      Pkm2 := Pkm1;
      Pkm1 := Pk;
      Qkm2 := Qkm1;
      Qkm1 := Qk;

      if Qk <> 0.0 then R := Pk / Qk;

      if R <> 0.0 then
        begin
          T := Abs((Ans - R) / R);
          Ans := R;
        end
      else
        T := 1.0;

      if T < Thresh then goto CDone;

      K1 := K1 + 1.0;
      K2 := K2 - 1.0;
      K3 := K3 + 2.0;
      K4 := K4 + 2.0;
      K5 := K5 + 1.0;
      K6 := K6 + 1.0;
      K7 := K7 + 2.0;
      K8 := K8 + 2.0;

      if Abs(Qk) + Abs(Pk) > BIG then
        begin
          Pkm2 := Pkm2 * BIGINV;
          Pkm1 := Pkm1 * BIGINV;
          Qkm2 := Qkm2 * BIGINV;
          Qkm1 := Qkm1 * BIGINV;
        end;

      if (Abs(Qk) < BIGINV) or (Abs(Pk) < BIGINV) then
        begin
          Pkm2 := Pkm2 * BIG;
          Pkm1 := Pkm1 * BIG;
          Qkm2 := Qkm2 * BIG;
          Qkm1 := Qkm1 * BIG;
        end;
      N := N + 1;
    until N > 400;
    MathError := FN_PLOSS;

CDone:
    CFrac2 := Ans;
  end;

  function IBeta(A, B, X : Float) : Float;
  var
    A1, B1, X1, T, W, Xc, Y : Float;
    Flag : Boolean;