' ******************************************************************
' Incomplete Beta function
' Translated from C code in Cephes library (http://www.moshier.net)
' ******************************************************************

' ------------------------------------------------------------------
' Error codes
' ------------------------------------------------------------------

#DEFINE FOk 0  
' No error

#DEFINE FDomain -1  
' Argument domain error

' ------------------------------------------------------------------
' Machine-dependent constants
' ------------------------------------------------------------------

#DEFINE MachEp 2.220446049250313D-16  
' Floating point precision: 2^(-52)

#DEFINE MinLog -708.3964185322641#       
' Min. argument for Exp

' ------------------------------------------------------------------
' Other constants
' ------------------------------------------------------------------

#DEFINE Big    1.44115188075856D+17
#DEFINE BigInv 6.93889390390722D-18

' ------------------------------------------------------------------
' Global variable
' ------------------------------------------------------------------

COMMON SHARED ErrCode AS INTEGER
' Error code from the latest function evaluation

' ------------------------------------------------------------------
' External function
' ------------------------------------------------------------------

DECLARE FUNCTION LnGamma (X AS DOUBLE) AS DOUBLE
' Ln of Gamma function

' ******************************************************************

FUNCTION IbCf1 (A AS DOUBLE, B AS DOUBLE, X AS DOUBLE) AS DOUBLE
' Incomplete beta function - Continuous fraction expansion #1

DIM AS DOUBLE Xk, R, T, Ans, Pk, Pkm1, Pkm2, Qk, Qkm1, Qkm2
DIM AS DOUBLE K1, K2, K3, K4, K5, K6, K7, K8

DIM AS INTEGER N

K1 = A
K2 = A + B
K3 = A
K4 = A + 1
K5 = 1
K6 = B - 1
K7 = K4
K8 = A + 2

Pkm2 = 0
Qkm2 = 1
Pkm1 = 1
Qkm1 = 1
Ans = 1
N = 0

DO
    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 THEN R = Pk / Qk
    IF R <> 0 THEN
        T = ABS((Ans - R) / R)
        Ans = R
    ELSE
        T = 1
    END IF

    IF T < MachEp THEN GOTO Done1

    K1 = K1 + 1
    K2 = K2 + 1
    K3 = K3 + 2
    K4 = K4 + 2
    K5 = K5 + 1
    K6 = K6 - 1
    K7 = K7 + 2
    K8 = K8 + 2

    IF ABS(Qk) + ABS(Pk) > Big THEN
        Pkm2 = Pkm2 * BigInv
        Pkm1 = Pkm1 * BigInv
        Qkm2 = Qkm2 * BigInv
        Qkm1 = Qkm1 * BigInv
    END IF
  
    IF ABS(Qk) < MachEp OR ABS(Pk) < MachEp THEN
        Pkm2 = Pkm2 * Big
        Pkm1 = Pkm1 * Big
        Qkm2 = Qkm2 * Big
        Qkm1 = Qkm1 * Big
    END IF

    N = N + 1
LOOP WHILE (N < 100)

Done1:
IbCf1 = Ans
END FUNCTION

FUNCTION IbCf2 (A AS DOUBLE, B AS DOUBLE, X AS DOUBLE) AS DOUBLE
' Incomplete beta function - Continuous fraction expansion #2

DIM AS DOUBLE Xk, R, T, Ans, Pk, Pkm1, Pkm2, Qk, Qkm1, Qkm2
DIM AS DOUBLE K1, K2, K3, K4, K5, K6, K7, K8, Z

DIM AS INTEGER N

K1 = A
K2 = B - 1
K3 = A
K4 = A + 1
K5 = 1
K6 = A + B
K7 = A + 1
K8 = A + 2

Pkm2 = 0
Qkm2 = 1
Pkm1 = 1
Qkm1 = 1
Z = X / (1 - X)
Ans = 1
N = 0

DO
    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 THEN R = Pk / Qk
    IF R <> 0 THEN
        T = ABS((Ans - R) / R)
        Ans = R
    ELSE
        T = 1
    END IF

    IF T < MachEp THEN GOTO Done2

    K1 = K1 + 1
    K2 = K2 - 1
    K3 = K3 + 2
    K4 = K4 + 2
    K5 = K5 + 1
    K6 = K6 + 1
    K7 = K7 + 2
    K8 = K8 + 2

    IF ABS(Qk) + ABS(Pk) > Big THEN
        Pkm2 = Pkm2 * BigInv
        Pkm1 = Pkm1 * BigInv
        Qkm2 = Qkm2 * BigInv
        Qkm1 = Qkm1 * BigInv
    END IF
    
    IF ABS(Qk) < MachEp OR ABS(Pk) < MachEp THEN
        Pkm2 = Pkm2 * Big
        Pkm1 = Pkm1 * Big
        Qkm2 = Qkm2 * Big
        Qkm1 = Qkm1 * Big
    END IF
    
    N = N + 1
LOOP WHILE (N < 100)

Done2:
IbCf2 = Ans
END FUNCTION

FUNCTION IBeta (A AS DOUBLE, B AS DOUBLE, X AS DOUBLE) AS DOUBLE
' Incomplete beta function

DIM AS DOUBLE Ans, A1, B1, T, X1

DIM AS INTEGER Flag 

IF A <= 0 OR B <= 0 OR X < 0 THEN
    ErrCode = FDomain
    IBeta = 0
    EXIT FUNCTION
END IF

IF X > 1 THEN
    ErrCode = FDomain
    IBeta = 1
    EXIT FUNCTION
END IF

IF X = 0 OR X = 1 THEN
    ErrCode = FOk
    IBeta = X
    EXIT FUNCTION
END IF

ErrCode = FOk
IF A <= 1 THEN
    Ans = IBeta(A + 1, B, X)
    T = A * LOG(X) + B * LOG(1 - X)
    T = T + LnGamma(A + B) - LnGamma(A + 1) - LnGamma(B)
    IF T > MinLog THEN Ans = Ans + EXP(T)
    IBeta = Ans
ELSE
    IF X > A / (A + B) THEN
        Flag = 1
        A1 = B
        B1 = A
        T = X
        X1 = 1 - X
    ELSE
        Flag = 0
        A1 = A
        B1 = B
        T = 1 - X
        X1 = X
    END IF
    Ans = X1 * (A1 + B1 - 2) / (A1 - 1)
    IF Ans < 1 THEN
        Ans = IbCf1(A1, B1, X1)
        T = B1 * LOG(T)
    ELSE
        Ans = IbCf2(A1, B1, X1)
        T = (B1 - 1) * LOG(T)
    END IF
    T = T + A1 * LOG(X1) + LnGamma(A1 + B1) - LnGamma(A1) - LnGamma(B1) + LOG(Ans / A1)
    IF T < MinLog THEN
        IF Flag = 0 THEN
            IBeta = 0
        ELSE
            IBeta = 1
        END IF
    END IF
    T = EXP(T)
    IF Flag = 1 THEN T = 1 - T
    IBeta = T
END IF
END FUNCTION
