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

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

#DEFINE FOk 0 
' No error

#DEFINE FUnderflow -4
' Underflow range error

#DEFINE FPLoss -5
' Partial loss of precision

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

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

#DEFINE MaxNum 1.797693134862315D+308
' Max. floating point number: 2^1024

#DEFINE MaxLog 709.782712893384#
' Max. argument for Exp

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

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

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

' ------------------------------------------------------------------
' External functions
' ------------------------------------------------------------------

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

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

DECLARE FUNCTION InvNorm(P AS DOUBLE) AS DOUBLE
' Inverse of Normal distribution function

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

FUNCTION InvBeta(A AS DOUBLE, B AS DOUBLE, Y AS DOUBLE) AS DOUBLE

  ErrCode = 0

  IF Y <= 0 THEN InvBeta = 0 : EXIT FUNCTION
  IF Y >= 1 THEN InvBeta = 1 : EXIT FUNCTION

  DIM AS DOUBLE  a1, b1, y0, y1, d, x, x0, x1
  DIM AS DOUBLE  lgm, yp, di, dithresh, yl, yh, xt
  DIM AS INTEGER i, rflg, ndir, nflg

  i = 0
  x0 = 0
  yl = 0
  x1 = 1
  yh = 1
  nflg = 0

  IF A <= 1 OR B <= 1 THEN
    dithresh = 1e-6
    rflg = 0
    a1 = A
    b1 = B
    y0 = Y
    x = a1 / (a1 + b1)
    y1 = IBeta(a1, b1, x)
    GOTO ihalve
  ELSE
    dithresh = 1e-4
  END IF

' approximation to inverse function

  yp = - InvNorm(Y)

  IF Y > 0.5 THEN
    rflg = 1
    a1 = B
    b1 = A
    y0 = 1 - Y
    yp = -yp
  ELSE
    rflg = 0
    a1 = A
    b1 = B
    y0 = Y
  END IF

  lgm = (yp * yp - 3) / 6
  x = 2 / (1 / (2 * a1 - 1) + 1 / (2 * b1 - 1))
  d = yp * SQR(x + lgm) / x _
        - (1 / (2 * b1 - 1) - 1 / (2 * a1 - 1)) _
        * (lgm + 5 / 6 - 2 / (3 * x))
  d = 2 * d

  IF d < MinLog THEN x = 1 : GOTO under

  x = a1 / (a1 + b1 * EXP(d))
  y1 = IBeta(a1, b1, x)
  yp = (y1 - y0) / y0

  IF ABS(yp) < 0.2 THEN GOTO newt

' Resort to interval halving if not close enough.

ihalve:

  ndir = 0
  di = 0.5

  FOR i = 0 TO 99
    IF i <> 0 THEN
      x = x0 + di * (x1 - x0)
      IF x = 1 THEN x = 1 - MachEp
      IF x = 0 THEN
        di = 0.5
        x = x0  +  di * (x1 - x0)
        IF x = 0 THEN GOTO under
      END IF
      y1 = IBeta(a1, b1, x)
      yp = (x1 - x0) / (x1 + x0)
      IF ABS(yp) < dithresh THEN GOTO newt
      yp = (y1 - y0) / y0
      IF ABS(yp) < dithresh THEN GOTO newt
    END IF
    IF y1 < y0 THEN
      x0 = x
      yl = y1
      IF ndir < 0 THEN
        ndir = 0
        di = 0.5
      ELSEIF ndir > 3 THEN
        di = 1 - (1 - di) * (1 - di)
      ELSEIF ndir > 1 THEN
        di = 0.5 * di + 0.5
      ELSE
        di = (y0 - y1) / (yh - yl)
      END IF
      ndir = ndir + 1
      IF x0 > 0.75 THEN
        IF rflg = 1 THEN
          rflg = 0
          a1 = A
          b1 = B
          y0 = Y
        ELSE
          rflg = 1
          a1 = B
          b1 = A
          y0 = 1 - Y
        END IF
        x = 1 - x
        y1 = IBeta(a1, b1, x)
        x0 = 0
        yl = 0
        x1 = 1
        yh = 1
        GOTO ihalve
      END IF
    ELSE
      x1 = x
      IF rflg = 1 AND x1 < MachEp THEN x = 0 : GOTO done
      yh = y1
      IF ndir > 0 THEN
        ndir = 0
        di = 0.5
      ELSEIF ndir < -3 THEN
        di = di * di
      ELSEIF ndir < -1 THEN
        di = 0.5 * di
      ELSE
        di = (y1 - y0) / (yh - yl)
      END IF
      ndir = ndir - 1
    END IF
  NEXT i

  ErrCode = FPLoss

  IF x0 >= 1 THEN x = 1 - MACHEP : GOTO done

  IF x <= 0 THEN
under:
    ErrCode = FUnderflow
    x = 0
    GOTO done
  END IF

newt:

  IF nflg = 1 THEN GOTO done

  nflg = 1
  lgm = LnGamma(a1 + b1) - LnGamma(a1) - LnGamma(b1)

  FOR i = 0 TO 7
    ' Compute the function at this point.
    IF i <> 0 THEN y1 = IBeta(a1, b1, x)
    IF y1 < yl THEN
      x = x0
      y1 = yl
    ELSEIF y1 > yh THEN
      x = x1
      y1 = yh
    ELSEIF y1 < y0 THEN
      x0 = x
      yl = y1
    ELSE
      x1 = x
      yh = y1
    END IF

    IF x = 1 OR x = 0 THEN EXIT FOR

    ' Compute the derivative of the function at this point.
    d = (a1 - 1) * LOG(x) + (b1 - 1) * LOG(1 - x) + lgm
    IF d < MinLog THEN GOTO done
    IF d > MaxLog THEN EXIT FOR
    d = EXP(d)

    ' Compute the step to the next approximation of x.
    d = (y1 - y0) / d
    xt = x - d
    IF xt <= x0 THEN
      y1 = (x - x0) / (x1 - x0)
      xt = x0 + 0.5 * y1 * (x - x0)
      IF xt <= 0 THEN EXIT FOR
    END IF
    IF xt >= x1 THEN
      y1 = (x1 - x) / (x1 - x0)
      xt = x1 - 0.5 * y1 * (x1 - x)
      IF xt >= 1 THEN EXIT FOR
    END IF
    x = xt
    IF ABS(d / x) < 128 * MachEp THEN GOTO done
  NEXT i

  ' Did not converge.
  dithresh = 256 * MachEp
  GOTO ihalve

done:

  IF rflg = 1 THEN
    IF x <= MachEp THEN x = 1 - MachEp ELSE x = 1 - x
  END IF

  InvBeta = x
END FUNCTION

