label
    Done;
  begin
    MathError := FN_OK;

    if (A <= 0.0) or (B <= 0.0) or (X < 0.0) then
      begin
        IBeta := DefaultVal(FN_DOMAIN, 0.0);
        Exit;
      end;

    if X > 1.0 then
      begin
        IBeta := DefaultVal(FN_DOMAIN, 1.0);
        Exit;
      end;

    if (X = 0.0) or (X = 1.0) then
      begin
        IBeta := X;
        Exit;
      end;

    Flag := False;
    if (B * X <= 1.0) and (X <= 0.95) then
      begin
        T := PSeries(A, B, X);
        goto Done;
      end;

    W := 1.0 - X;

    { Reverse a and b if x is greater than the mean. }
    if X > A / (A + B) then
      begin
        Flag := True;
        A1 := B;
        B1 := A;
        Xc := X;
        X1 := W;
      end
    else
      begin
        A1 := A;
        B1 := B;
        Xc := W;
        X1 := X;
      end;

    if Flag and (B1 * X1 <= 1.0) and (X1 <= 0.95) then
      begin
        T := PSeries(A1, B1, X1);
        goto Done;
      end;

    { Choose expansion for optimal convergence }
    Y := X1 * (A1 + B1 - 2.0) - (A1 - 1.0);
    if Y < 0.0 then
      W := CFrac1(A1, B1, X1)
    else
      W := CFrac2(A1, B1, X1) / Xc;

    { Multiply w by the factor
     a      b   _             _     _
    x  (1-x)   | (a+b) / ( a | (a) | (b) )    }

    Y := A1 * Ln(X1);
    T := B1 * Ln(Xc);
    if (A1 + B1 < MAXGAM) and (Abs(Y) < MAXLOG) and (Abs(T) < MAXLOG) then
      begin
        T := Power(Xc, B1) ;
        T := T * Power(X1, A1);
        T := T / A1;
        T := T * W;
        T := T * Gamma(A1 + B1) / (Gamma(A1) * Gamma(B1));
      end
    else
      begin
        { Resort to logarithms }
        Y := Y + T + LnGamma(A1 + B1) - LnGamma(A1) - LnGamma(B1) + Ln(W / A1);
        if Y < MINLOG then
          T := 0.0
        else
          T := Exp(Y);
      end;

Done:
    if Flag then
      if T <= MACHEP then
        T := 1.0 - MACHEP
      else
        T := 1.0 - T;

    IBeta := T;
  end;

  function Erf(X : Float) : Float;
  begin
    if X < 0.0 then
      Erf := - IGamma(0.5, Sqr(X))
    else
      Erf := IGamma(0.5, Sqr(X));
  end;

  function Erfc(X : Float) : Float;
  begin
    if X < 0.0 then
      Erfc := 1.0 + IGamma(0.5, Sqr(X))
    else
      Erfc := JGamma(0.5, Sqr(X));
  end;

{ ----------------------------------------------------------------------
  Probability functions
  ---------------------------------------------------------------------- }

  function PBinom(N : Integer; P : Float; K : Integer) : Float;
  begin
    MathError := FN_OK;
    if (P < 0.0) or (P > 1.0) or (N <= 0) or (N < K) then
      PBinom := DefaultVal(FN_DOMAIN, 0.0)
    else if K = 0 then
      PBinom := Power(1.0 - P, N)
    else if K = N then
      PBinom := Power(P, N)
    else
      PBinom := Binomial(N, K) * Power(P, K) * Power(1.0 - P, N - K);
  end;

  function FBinom(N : Integer; P : Float; K : Integer) : Float;
  begin
    MathError := FN_OK;
    if (P < 0.0) or (P > 1.0) or (N <= 0) or (N < K) then
      FBinom := DefaultVal(FN_DOMAIN, 0.0)
    else if K = 0 then
      FBinom := Power(1.0 - P, N)
    else if K = N then
      FBinom := 1.0
    else
      FBinom := 1.0 - IBeta(K + 1, N - K, P);
  end;

  function PPoisson(Mu : Float; K : Integer) : Float;
  var
    P : Float;
    I : Integer;
  begin
    MathError := FN_OK;
    if (Mu <= 0.0) or (K < 0) then
      PPoisson := DefaultVal(FN_DOMAIN, 0.0)
    else if K = 0 then
      PPoisson := Expo(- Mu)
    else
      begin
        P := Mu;
        for I := 2 to K do
          P := P * Mu / I;
        PPoisson := Expo(- Mu) * P;
      end;
  end;

  function FPoisson(Mu : Float; K : Integer) : Float;
  begin
    MathError := FN_OK;
    if (Mu <= 0.0) or (K < 0) then
      FPoisson := DefaultVal(FN_DOMAIN, 0.0)
    else if K = 0 then
      FPoisson := Expo(- Mu)
    else
      FPoisson := JGamma(K + 1, Mu);
  end;

  function DNorm(X : Float) : Float;
  begin
    DNorm := INVSQRT2PI * Expo(- 0.5 * Sqr(X));
  end;

  function FNorm(X : Float) : Float;
  begin
    FNorm := 0.5 * (1.0 + Erf(X * SQRT2DIV2));
  end;

  function InvNorm(P : Float) : Float;
{ ----------------------------------------------------------------------
  Inverse of Normal distribution function

  Returns the argument, X, for which the area under the Gaussian
  probability density function (integrated from minus infinity to X)
  is equal to P.
  ---------------------------------------------------------------------- }
  const
    P0 : TabCoef = (
        8.779679420055069160496E-3,
      - 7.649544967784380691785E-1,
        2.971493676711545292135E0,
      - 4.144980036933753828858E0,
        2.765359913000830285937E0,
      - 9.570456817794268907847E-1,
        1.659219375097958322098E-1,
      - 1.140013969885358273307E-2,
        0, 0);

    Q0 : TabCoef = (
      - 5.303846964603721860329E0,
        9.908875375256718220854E0,
      - 9.031318655459381388888E0,
        4.496118508523213950686E0,
      - 1.250016921424819972516E0,
        1.823840725000038842075E-1,
      - 1.088633151006419263153E-2,
        0, 0, 0);

    P1 : TabCoef = (
      4.302849750435552180717E0,
      4.360209451837096682600E1,
      9.454613328844768318162E1,
      9.336735653151873871756E1,
      5.305046472191852391737E1,
      1.775851836288460008093E1,
      3.640308340137013109859E0,
      3.691354900171224122390E-1,
      1.403530274998072987187E-2,
      1.377145111380960566197E-4);

    Q1 : TabCoef = (
      2.001425109170530136741E1,
      7.079893963891488254284E1,
      8.033277265194672063478E1,
      5.034715121553662712917E1,
      1.779820137342627204153E1,
      3.845554944954699547539E0,
      3.993627390181238962857E-1,
      1.526870689522191191380E-2,
      1.498700676286675466900E-4,
      0);

    P2 : TabCoef = (
      3.244525725312906932464E0,
      6.856256488128415760904E0,
      3.765479340423144482796E0,
      1.240893301734538935324E0,
      1.740282292791367834724E-1,
      9.082834200993107441750E-3,
      1.617870121822776093899E-4,
      7.377405643054504178605E-7,
      0, 0);

    Q2 : TabCoef = (
      6.021509481727510630722E0,
      3.528463857156936773982E0,
      1.289185315656302878699E0,
      1.874290142615703609510E-1,
      9.867655920899636109122E-3,
      1.760452434084258930442E-4,
      8.028288500688538331773E-7,
      0, 0, 0);

    P3 : TabCoef = (
        2.020331091302772535752E0,
        2.133020661587413053144E0,
        2.114822217898707063183E-1,
      - 6.500909615246067985872E-3,
      - 7.279315200737344309241E-4,
      - 1.275404675610280787619E-5,
      - 6.433966387613344714022E-8,
      - 7.772828380948163386917E-11,
        0, 0);

    Q3 : TabCoef = (
        2.278210997153449199574E0,
        2.345321838870438196534E-1,
      - 6.916708899719964982855E-3,
      - 7.908542088737858288849E-4,
      - 1.387652389480217178984E-5,
      - 7.001476867559193780666E-8,
      - 8.458494263787680376729E-11,
        0, 0, 0);

  var
    X, Y, Z, Y2, X0, X1 : Float;
    Code : Integer;
  begin
    if (P <= 0.0) or (P >= 1.0) then
      begin
        InvNorm := DefaultVal(FN_DOMAIN, 0.0);
        Exit;
      end;

    Code := 1;
    Y := P;
    if Y > (1.0 - 0.13533528323661269189) then { 0.135... = exp(-2) }
      begin
        Y := 1.0 - Y;
        Code := 0;
      end;
    if Y > 0.13533528323661269189 then
      begin
        Y := Y - 0.5;
        Y2 := Y * Y;
        X := Y + Y * (Y2 * PolEvl(Y2, P0, 7) / P1Evl(Y2, Q0, 7));
        X := X * SQRT2PI;
        InvNorm := X;
        Exit;
      end;

    X := Sqrt(- 2.0 * Ln(Y));
    X0 := X - Ln(X) / X;
    Z := 1.0 / X;
    if X < 8.0 then
      X1 := Z * PolEvl(Z, P1, 9) / P1Evl(Z, Q1, 9)
    else if X < 32.0 then
      X1 := Z * PolEvl(Z, P2, 7) / P1Evl(Z, Q2, 7)
    else
      X1 := Z * PolEvl(Z, P3, 7) / P1Evl(Z, Q3, 7);
    X := X0 - X1;
    if Code <> 0 then
      X := - X;
    InvNorm := X;
  end;

  function PNorm(X : Float) : Float;
  var
    A : Float;
  begin
    A := Abs(X);
    MathError := FN_OK;
    if A = 0.0 then
      PNorm := 1.0
    else if A < 1.0 then
      PNorm := 1.0 - Erf(A * SQRT2DIV2)
    else
      PNorm := Erfc(A * SQRT2DIV2);
  end;

  function DStudent(Nu : Integer; X : Float) : Float;
  var
    L, P, Q : Float;
  begin
    MathError := FN_OK;
    if Nu < 1 then
      DStudent := DefaultVal(FN_DOMAIN, 0.0)
    else
      begin
        P := 0.5 * (Nu + 1);
    	Q := 0.5 * Nu;
    	L := LnGamma(P) - LnGamma(Q) - 0.5 * Ln(Nu * PI) - P * Ln(1.0 + Sqr(X) / Nu);
    	DStudent := Expo(L);
      end;
  end;

  function FStudent(Nu : Integer; X : Float) : Float;
  var
    F : Float;
  begin
    MathError := FN_OK;
    if Nu < 1 then
      FStudent := DefaultVal(FN_DOMAIN, 0.0)
    else if X = 0.0 then
      FStudent := 0.5
    else
      begin
        F := 0.5 * IBeta(0.5 * Nu, 0.5, Nu / (Nu + Sqr(X)));
        if X < 0.0 then FStudent := F else FStudent := 1.0 - F;
      end;
  end;

  function PStudent(Nu : Integer; X : Float) : Float;
  begin
    MathError := FN_OK;
    if Nu < 1 then
      PStudent := DefaultVal(FN_DOMAIN, 0.0)
    else
      PStudent := IBeta(0.5 * Nu, 0.5, Nu / (Nu + Sqr(X)));
  end;

  function DKhi2(Nu : Integer; X : Float) : Float;
  begin
    MathError := FN_OK;
    DKhi2 := DGamma(0.5 * Nu, 0.5, X);
  end;

  function FKhi2(Nu : Integer; X : Float) : Float;
  begin
    MathError := FN_OK;
    if (Nu < 1) or (X <= 0.0) then
      FKhi2 := DefaultVal(FN_DOMAIN, 0.0)
    else
      FKhi2 := IGamma(0.5 * Nu, 0.5 * X);
  end;

  function PKhi2(Nu : Integer; X : Float) : Float;
  begin
    MathError := FN_OK;
    if (Nu < 1) or (X <= 0.0) then
      PKhi2 := DefaultVal(FN_DOMAIN, 0.0)
    else
      PKhi2 := JGamma(0.5 * Nu, 0.5 * X);
  end;

  function DSnedecor(Nu1, Nu2 : Integer; X : Float) : Float;
  var
    P1, P2, R, S, L : Float;
  begin
    MathError := FN_OK;
    if (Nu1 < 1) or (Nu2 < 1) or (X <= 0.0) then
      DSnedecor := DefaultVal(FN_DOMAIN, 0.0)
    else
      begin
        R := Int(Nu1) / Int(Nu2);
        P1 := 0.5 * Nu1;
        P2 := 0.5 * Nu2;
        S := P1 + P2;
        L := LnGamma(S) - LnGamma(P1) - LnGamma(P2) + P1 * Ln(R) +
               (P1 - 1.0) * Ln(X) - S * Ln(1.0 + R * X);
        DSnedecor := Expo(L);
      end;
  end;

  function FSnedecor(Nu1, Nu2 : Integer; X : Float) : Float;
  begin
    MathError := FN_OK;
    if (Nu1 < 1) or (Nu2 < 1) or (X <= 0.0) then
      FSnedecor := DefaultVal(FN_DOMAIN, 0.0)
    else
      FSnedecor := 1.0 - IBeta(0.5 * Nu2, 0.5 * Nu1, Nu2 / (Nu2 + Nu1 * X));
  end;

  function PSnedecor(Nu1, Nu2 : Integer; X : Float) : Float;
  begin
    MathError := FN_OK;
    if (Nu1 < 1) or (Nu2 < 1) or (X <= 0.0) then
      PSnedecor := DefaultVal(FN_DOMAIN, 0.0)
    else
      PSnedecor := IBeta(0.5 * Nu2, 0.5 * Nu1, Nu2 / (Nu2 + Nu1 * X));
  end;

  function DExpo(A, X : Float) : Float;
  begin
    if (A <= 0.0) or (X < 0.0) then
      DExpo := DefaultVal(FN_DOMAIN, 0.0)
    else
      DExpo := A * Expo(- A * X);
  end;

  function FExpo(A, X : Float) : Float;
  begin
    if (A <= 0.0) or (X < 0.0) then
      FExpo := DefaultVal(FN_DOMAIN, 0.0)
    else
      FExpo := 1.0 - Expo(- A * X);
  end;

  function DBeta(A, B, X : Float) : Float;
  var
    L : Float;
  begin
    MathError := FN_OK;
    if (A <= 0.0) or (B <= 0.0) or (X < 0.0) or (X > 1.0) then
      DBeta := DefaultVal(FN_DOMAIN, 0.0)
    else if X = 0.0 then
      if A < 1.0 then DBeta := DefaultVal(FN_SING, MAXNUM) else DBeta := 0.0
    else if X = 1.0 then
      if B < 1.0 then DBeta := DefaultVal(FN_SING, MAXNUM) else DBeta := 0.0
    else
      begin
        L := LnGamma(A + B) - LnGamma(A) - LnGamma(B) +
               (A - 1.0) * Ln(X) + (B - 1.0) * Ln(1.0 - X);
        DBeta := Expo(L);
      end;
  end;

  function FBeta(A, B, X : Float) : Float;
  begin
    FBeta := IBeta(A, B, X);
  end;

  function DGamma(A, B, X : Float) : Float;
  var
    L : Float;
  begin
    MathError := FN_OK;
    if (A <= 0.0) or (B <= 0.0) or (X < 0.0) then
      DGamma := DefaultVal(FN_DOMAIN, 0.0)
    else if X = 0.0 then
      if A < 1.0 then
        DGamma := DefaultVal(FN_SING, MAXNUM)
      else if A = 1.0 then
        DGamma := B
      else
        DGamma := 0.0
    else
      begin
        L := A * Ln(B) - LnGamma(A) + (A - 1.0) * Ln(X) - B * X;
        DGamma := Expo(L);
      end;
  end;

  function FGamma(A, B, X : Float) : Float;
  begin
    FGamma := IGamma(A, B * X);
  end;

{ ----------------------------------------------------------------------
  Initialization code
  ---------------------------------------------------------------------- }

begin
  MathError := FN_OK;
end.