' ******************************************************************
'               Inverse of incomplete Gamma 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 JGamma (A AS DOUBLE, X AS DOUBLE) AS DOUBLE
' Complement of incomplete Gamma function

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

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

FUNCTION InvGamma(A AS DOUBLE, P AS DOUBLE) AS DOUBLE
' ------------------------------------------------------------------
' Given P, the function finds X such that
'
' IGamma(A, X) = P
'
' It is best valid in the right-hand tail of the distribution, P > 0.5.
' ------------------------------------------------------------------

  DIM AS DOUBLE  x0, x1, x, Y, yl, yh, y1, d, lgm, dithresh
  DIM AS INTEGER i, ndir

  IF P > 0.5 THEN ErrCode = FPLoss ELSE ErrCode = 0

  Y = 1 - P

  ' bound the solution
  x0 = MaxNum
  yl = 0
  x1 = 0               
  yh = 1
  dithresh = 5 * MachEp 

  ' approximation to inverse function
  d = 1 / (9 * a)
  y1 = 1 - d - InvNorm(Y) * SQR(d)
  x = a * y1 * y1 * y1

  lgm = LnGamma(a)

  FOR i = 0 TO 9
    IF x > x0 OR x < x1 THEN GOTO ihalve
    y1 = JGamma(a, x)
    IF y1 < yl OR y1 > yh THEN GOTO ihalve
    IF y1 < Y THEN
      x0 = x
      yl = y1
    ELSE
      x1 = x
      yh = y1
    END IF

    ' compute the derivative of the function at this point
    d = (a - 1) * LOG(x) - x - lgm
    IF d < MinLog THEN GOTO ihalve
    d = -EXP(d)

    ' compute the step to the next approximation of x
    d = (y1 - Y) / d
    IF ABS(d / x) < MachEp THEN GOTO done
    x = x - d 
  NEXT i 

  ' Resort to interval halving if Newton iteration did not converge.
ihalve:

  d = 0.0625
  IF x0 = MaxNum THEN
    IF x <= 0 THEN x = 1
    DO WHILE x0 = MaxNum
      x = (1 + d) * x
      y1 = JGamma(a, x)
      IF y1 < Y THEN
        x0 = x
        yl = y1
        EXIT DO
      END IF
      d = d + d
    LOOP
  END IF
  d = 0.5
  ndir = 0

  FOR i = 0 TO 399
    x = x1 + d * (x0 - x1)
    y1 = JGamma(a, x)
    lgm = (x0 - x1) / (x1 + x0)
    IF ABS(lgm) < dithresh THEN EXIT FOR
    lgm = (y1 - Y) / Y
    IF ABS(lgm) < dithresh THEN EXIT FOR
    IF x <= 0 THEN EXIT FOR
    IF y1 >= Y THEN
      x1 = x
      yh = y1
      IF ndir < 0 THEN 
        ndir = 0
        d = 0.5
      ELSEIF ndir > 1 THEN
        d = 0.5 * d + 0.5
      ELSE
        d = (Y - yl) / (yh - yl)
      END IF
      ndir = ndir + 1
    ELSE
      x0 = x
      yl = y1
      IF ndir > 0 THEN
        ndir = 0
        d = 0.5
      ELSEIF ndir < -1 THEN
        d = 0.5 * d
      ELSE
        d = (Y - yl) / (yh - yl)
      END IF
      ndir = ndir - 1
    END IF
  NEXT i

  IF x = 0 THEN ErrCode = FUnderflow

done:
  InvGamma = x
END FUNCTION
