#! /usr/bin/env python
#
def betain ( x, p, q, beta ):

#*****************************************************************************80
#
## BETAIN computes the incomplete Beta function ratio.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    23 December 2015
#
#  Author:
#
#    Original FORTRAN77 versionby KL Majumder, GP Bhattacharjee.
#    Python version by John Burkardt.
#
#  Reference:
#
#    KL Majumder, GP Bhattacharjee,
#    Algorithm AS 63:
#    The incomplete Beta Integral,
#    Applied Statistics,
#    Volume 22, Number 3, 1973, pages 409-411.
#
#  Parameters:
#
#    Input, real X, the argument, between 0 and 1.
#
#    Input, real P, Q, the parameters, which
#    must be positive.
#
#    Input, real BETA, the logarithm of the complete
#    beta function.
#
#    Output, real BETAIN, the value of the incomplete
#    Beta function ratio.
#
#    Output, integer IFAULT, error flag.
#    0, no error.
#    nonzero, an error occurred.
#
  import numpy as np

  acu = 0.1E-14

  value = x
  ifault = 0
#
#  Check the input arguments.
#
  if ( p <= 0.0 or q <= 0.0 ):
    ifault = 1
    return value, ifault

  if ( x < 0.0 or 1.0 < x ):
    ifault = 2
    return value, ifault
#
#  Special cases.
#
  if ( x == 0.0 or x == 1.0 ):
    return value, ifault
#
#  Change tail if necessary and determine S.
#
  psq = p + q
  cx = 1.0 - x

  if ( p < psq * x ):
    xx = cx
    cx = x
    pp = q
    qq = p
    indx = 1
  else:
    xx = x
    pp = p
    qq = q
    indx = 0

  term = 1.0
  ai = 1.0
  value = 1.0
  ns = np.floor ( qq + cx * psq )
#
#  Use the Soper reduction formula.
#
  rx = xx / cx
  temp = qq - ai
  if ( ns == 0 ):
    rx = xx

  while ( True ):

    term = term * temp * rx / ( pp + ai )
    value = value + term
    temp = abs ( term )

    if ( temp <= acu and temp <= acu * value ):

      value = value * np.exp ( pp * np.log ( xx ) \
      + ( qq - 1.0 ) * np.log ( cx ) - beta ) / pp

      if ( indx ):
        value = 1.0 - value

      break

    ai = ai + 1.0
    ns = ns - 1

    if ( 0 <= ns ):
      temp = qq - ai
      if ( ns == 0 ):
        rx = xx
    else:
      temp = psq
      psq = psq + 1.0

  return value, ifault

def betain_test ( ):

#*****************************************************************************80
#
## ASA063_TEST01 demonstrates the use of BETAIN.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    23 December 2015
#
#  Author:
#
#    John Burkardt
#
  import platform
  from alogam import alogam
  from beta_inc_values import beta_inc_values

  print ( '' )
  print ( 'BETAIN_TEST' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  BETAIN computes the incomplete beta function.' )
  print ( '  Compare to tabulated values.' )
  print ( '' )
  print ( '      A       B       X     FX                       FX2' )
  print ( '                            (Tabulated)              (BETAIN)                DIFF' )
  print ( '' )

  n_data = 0

  while ( True ):

    n_data, a, b, x, fx = beta_inc_values ( n_data )

    if ( n_data == 0 ):
      break

    beta_log = alogam ( a ) \
             + alogam ( b ) \
             - alogam ( a + b )

    fx2, ifault = betain ( x, a, b, beta_log )

    print ( '  %6.2f  %6.2f  %6.2f  %24.16e  %24.16e  %10.4e' \
      % ( a, b, x, fx, fx2, abs ( fx - fx2 ) ) )
#
#  Terminate.
#
  print ( '' )
  print ( 'BETAIN_TEST:' )
  print ( '  Normal end of execution.' )

  return

if ( __name__ == '__main__' ):
  from timestamp import timestamp
  timestamp ( )
  betain_test ( )
  timestamp ( )