#! /usr/bin/env python
#
def r8_erf ( x ):

#*****************************************************************************80
#
## R8_ERF evaluates the error function of an R8 argument.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    23 April 2016
#
#  Author:
#
#    Original FORTRAN77 version by Wayne Fullerton.
#    Python version by John Burkardt.
#
#  Reference:
#
#    Wayne Fullerton,
#    Portable Special Function Routines,
#    in Portability of Numerical Software,
#    edited by Wayne Cowell,
#    Lecture Notes in Computer Science, Volume 57,
#    Springer 1977,
#    ISBN: 978-3-540-08446-4,
#    LC: QA297.W65.
#
#  Parameters:
#
#    Input, real X, the argument.
#
#    Output, real VALUE, the error function of X.
#
  import numpy as np
  from r8_csevl import r8_csevl
  from r8_inits import r8_inits
  from machine import r8_mach

  sqrtpi = 1.77245385090551602729816748334115

  erfcs = np.array ( [ \
      -0.49046121234691808039984544033376E-01, \
      -0.14226120510371364237824741899631E+00, \
      +0.10035582187599795575754676712933E-01, \
      -0.57687646997674847650827025509167E-03, \
      +0.27419931252196061034422160791471E-04, \
      -0.11043175507344507604135381295905E-05, \
      +0.38488755420345036949961311498174E-07, \
      -0.11808582533875466969631751801581E-08, \
      +0.32334215826050909646402930953354E-10, \
      -0.79910159470045487581607374708595E-12, \
      +0.17990725113961455611967245486634E-13, \
      -0.37186354878186926382316828209493E-15, \
      +0.71035990037142529711689908394666E-17, \
      -0.12612455119155225832495424853333E-18, \
      +0.20916406941769294369170500266666E-20, \
      -0.32539731029314072982364160000000E-22, \
      +0.47668672097976748332373333333333E-24, \
      -0.65980120782851343155199999999999E-26, \
      +0.86550114699637626197333333333333E-28, \
      -0.10788925177498064213333333333333E-29, \
      +0.12811883993017002666666666666666E-31 ] )

  nterf = r8_inits ( erfcs, 21, 0.1 * r8_mach ( 3 ) )
  xbig = np.sqrt ( - np.log ( sqrtpi * r8_mach ( 3 ) ) )
  sqeps = np.sqrt ( 2.0 * r8_mach ( 3 ) )

  y = abs ( x )

  if ( y <= sqeps ):
    value = 2.0 * x / sqrtpi
  elif ( y <= 1.0 ):
    value = x * ( 1.0 + r8_csevl ( 2.0 * x * x - 1.0, erfcs, nterf ) )
  elif ( y <= xbig ):
    value = 1.0 - r8_erfc ( y )
    if ( x < 0.0 ):
      value = - value
  else:
    value = 1.0
    if ( x < 0.0 ):
      value = - value

  return value

def r8_erf_test ( ):

#*****************************************************************************80
#
## R8_ERF_TEST tests R8_ERF.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    23 April 2016
#
#  Author:
#
#    John Burkardt
#
  from erf_values import erf_values

  print ( '' )
  print ( 'R8_ERF_TEST:' )
  print ( '  R8_ERF evaluates the error function.' )
  print ( '' )
  print ( '             X          ERF(X)  R8_ERF(X)        Diff' )
  print ( '' )

  n_data = 0

  while ( True ):

    n_data, x, fx1 = erf_values ( n_data )

    if ( n_data == 0 ):
      break

    fx2 = r8_erf ( x )

    print ( '  %14.4g  %14.6g  %14.6g  %14.6g' % ( x, fx1, fx2, abs ( fx1 - fx2 ) ) )
#
#  Terminate.
#
  print ( '' )
  print ( 'R8_ERF_TEST:' )
  print ( '  Normal end of execution.' )
  return

def r8_erfc ( x ):

#*****************************************************************************80
#
## R8_ERFC evaluates the co-error function of an R8 argument.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    23 April 2016
#
#  Author:
#
#    Original FORTRAN77 version by Wayne Fullerton.
#    Python version by John Burkardt.
#
#  Reference:
#
#    Wayne Fullerton,
#    Portable Special Function Routines,
#    in Portability of Numerical Software,
#    edited by Wayne Cowell,
#    Lecture Notes in Computer Science, Volume 57,
#    Springer 1977,
#    ISBN: 978-3-540-08446-4,
#    LC: QA297.W65.
#
#  Parameters:
#
#    Input, real X, the argument.
#
#    Output, real VALUE, the co-error function of X.
#
  import numpy as np
  from r8_csevl import r8_csevl
  from r8_inits import r8_inits
  from machine import r8_mach

  sqrtpi = 1.77245385090551602729816748334115

  erc2cs = np.array ( [ \
    -0.6960134660230950112739150826197E-01, \
    -0.4110133936262089348982212084666E-01, \
    +0.3914495866689626881561143705244E-02, \
    -0.4906395650548979161280935450774E-03, \
    +0.7157479001377036380760894141825E-04, \
    -0.1153071634131232833808232847912E-04, \
    +0.1994670590201997635052314867709E-05, \
    -0.3642666471599222873936118430711E-06, \
    +0.6944372610005012589931277214633E-07, \
    -0.1371220902104366019534605141210E-07, \
    +0.2788389661007137131963860348087E-08, \
    -0.5814164724331161551864791050316E-09, \
    +0.1238920491752753181180168817950E-09, \
    -0.2690639145306743432390424937889E-10, \
    +0.5942614350847910982444709683840E-11, \
    -0.1332386735758119579287754420570E-11, \
    +0.3028046806177132017173697243304E-12, \
    -0.6966648814941032588795867588954E-13, \
    +0.1620854541053922969812893227628E-13, \
    -0.3809934465250491999876913057729E-14, \
    +0.9040487815978831149368971012975E-15, \
    -0.2164006195089607347809812047003E-15, \
    +0.5222102233995854984607980244172E-16, \
    -0.1269729602364555336372415527780E-16, \
    +0.3109145504276197583836227412951E-17, \
    -0.7663762920320385524009566714811E-18, \
    +0.1900819251362745202536929733290E-18, \
    -0.4742207279069039545225655999965E-19, \
    +0.1189649200076528382880683078451E-19, \
    -0.3000035590325780256845271313066E-20, \
    +0.7602993453043246173019385277098E-21, \
    -0.1935909447606872881569811049130E-21, \
    +0.4951399124773337881000042386773E-22, \
    -0.1271807481336371879608621989888E-22, \
    +0.3280049600469513043315841652053E-23, \
    -0.8492320176822896568924792422399E-24, \
    +0.2206917892807560223519879987199E-24, \
    -0.5755617245696528498312819507199E-25, \
    +0.1506191533639234250354144051199E-25, \
    -0.3954502959018796953104285695999E-26, \
    +0.1041529704151500979984645051733E-26, \
    -0.2751487795278765079450178901333E-27, \
    +0.7290058205497557408997703680000E-28, \
    -0.1936939645915947804077501098666E-28, \
    +0.5160357112051487298370054826666E-29, \
    -0.1378419322193094099389644800000E-29, \
    +0.3691326793107069042251093333333E-30, \
    -0.9909389590624365420653226666666E-31, \
    +0.2666491705195388413323946666666E-31 ] )

  erfccs = np.array ( [ \
    +0.715179310202924774503697709496E-01, \
    -0.265324343376067157558893386681E-01, \
    +0.171115397792085588332699194606E-02, \
    -0.163751663458517884163746404749E-03, \
    +0.198712935005520364995974806758E-04, \
    -0.284371241276655508750175183152E-05, \
    +0.460616130896313036969379968464E-06, \
    -0.822775302587920842057766536366E-07, \
    +0.159214187277090112989358340826E-07, \
    -0.329507136225284321486631665072E-08, \
    +0.722343976040055546581261153890E-09, \
    -0.166485581339872959344695966886E-09, \
    +0.401039258823766482077671768814E-10, \
    -0.100481621442573113272170176283E-10, \
    +0.260827591330033380859341009439E-11, \
    -0.699111056040402486557697812476E-12, \
    +0.192949233326170708624205749803E-12, \
    -0.547013118875433106490125085271E-13, \
    +0.158966330976269744839084032762E-13, \
    -0.472689398019755483920369584290E-14, \
    +0.143587337678498478672873997840E-14, \
    -0.444951056181735839417250062829E-15, \
    +0.140481088476823343737305537466E-15, \
    -0.451381838776421089625963281623E-16, \
    +0.147452154104513307787018713262E-16, \
    -0.489262140694577615436841552532E-17, \
    +0.164761214141064673895301522827E-17, \
    -0.562681717632940809299928521323E-18, \
    +0.194744338223207851429197867821E-18, \
    -0.682630564294842072956664144723E-19, \
    +0.242198888729864924018301125438E-19, \
    -0.869341413350307042563800861857E-20, \
    +0.315518034622808557122363401262E-20, \
    -0.115737232404960874261239486742E-20, \
    +0.428894716160565394623737097442E-21, \
    -0.160503074205761685005737770964E-21, \
    +0.606329875745380264495069923027E-22, \
    -0.231140425169795849098840801367E-22, \
    +0.888877854066188552554702955697E-23, \
    -0.344726057665137652230718495566E-23, \
    +0.134786546020696506827582774181E-23, \
    -0.531179407112502173645873201807E-24, \
    +0.210934105861978316828954734537E-24, \
    -0.843836558792378911598133256738E-25, \
    +0.339998252494520890627359576337E-25, \
    -0.137945238807324209002238377110E-25, \
    +0.563449031183325261513392634811E-26, \
    -0.231649043447706544823427752700E-26, \
    +0.958446284460181015263158381226E-27, \
    -0.399072288033010972624224850193E-27, \
    +0.167212922594447736017228709669E-27, \
    -0.704599152276601385638803782587E-28, \
    +0.297976840286420635412357989444E-28, \
    -0.126252246646061929722422632994E-28, \
    +0.539543870454248793985299653154E-29, \
    -0.238099288253145918675346190062E-29, \
    +0.109905283010276157359726683750E-29, \
    -0.486771374164496572732518677435E-30, \
    +0.152587726411035756763200828211E-30 ] )

  erfcs = np.array ( [ \
    -0.49046121234691808039984544033376E-01, \
    -0.14226120510371364237824741899631E+00, \
    +0.10035582187599795575754676712933E-01, \
    -0.57687646997674847650827025509167E-03, \
    +0.27419931252196061034422160791471E-04, \
    -0.11043175507344507604135381295905E-05, \
    +0.38488755420345036949961311498174E-07, \
    -0.11808582533875466969631751801581E-08, \
    +0.32334215826050909646402930953354E-10, \
    -0.79910159470045487581607374708595E-12, \
    +0.17990725113961455611967245486634E-13, \
    -0.37186354878186926382316828209493E-15, \
    +0.71035990037142529711689908394666E-17, \
    -0.12612455119155225832495424853333E-18, \
    +0.20916406941769294369170500266666E-20, \
    -0.32539731029314072982364160000000E-22, \
    +0.47668672097976748332373333333333E-24, \
    -0.65980120782851343155199999999999E-26, \
    +0.86550114699637626197333333333333E-28, \
    -0.10788925177498064213333333333333E-29, \
    +0.12811883993017002666666666666666E-31 ] )

  eta = 0.1 * r8_mach ( 3 )
  nterf = r8_inits ( erfcs, 21, eta )
  nterfc = r8_inits ( erfccs, 59, eta )
  nterc2 = r8_inits ( erc2cs, 49, eta )

  xsml = - np.sqrt ( - np.log ( sqrtpi * r8_mach ( 3 ) ) )
  xmax = np.sqrt (- np.log ( sqrtpi * r8_mach ( 1 ) ) )
  xmax = xmax - 0.5 * np.log ( xmax ) / xmax - 0.01
  sqeps = np.sqrt ( 2.0 * r8_mach ( 3 ) )

  if ( x <= xsml ):

    value = 2.0
    return value

  if ( xmax < x ):
    print ( '' )
    print ( 'R8_ERFC - Warning!' )
    print ( '  X so big that ERFC underflows.' )
    value = 0.0
    return value

  y = abs ( x )

  if ( y < sqeps ):
    value = 1.0 - 2.0 * x / sqrtpi
    return value
  elif ( y <= 1.0 ):
    value = 1.0 - x * ( 1.0 + r8_csevl ( 2.0 * x * x - 1.0, erfcs, nterf ) )
    return value

  y = y * y

  if ( y <= 4.0 ):
    value = np.exp ( - y ) / abs ( x ) * ( 0.5 \
      + r8_csevl ( ( 8.0 / y - 5.0 ) / 3.0, erc2cs, nterc2 ) )
  else:
    value = np.exp ( - y ) / abs ( x ) * ( 0.5 \
      + r8_csevl ( 8.0 / y - 1.0, erfccs, nterfc ) )

  if ( x < 0.0 ):
    value = 2.0 - value

  return value

def r8_erfc_test ( ):

#*****************************************************************************80
#
## R8_ERFC_TEST tests R8_ERFC.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    23 April 2016
#
#  Author:
#
#    John Burkardt
#
  import platform
  from erfc_values import erfc_values

  print ( '' )
  print ( 'R8_ERFC_TEST:' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  R8_ERFC evaluates the complementary error function.' )
  print ( '' )
  print ( '             X         ERFC(X)  R8_ERFC(X)        Diff' )
  print ( '' )

  n_data = 0

  while ( True ):

    n_data, x, fx1 = erfc_values ( n_data )

    if ( n_data == 0 ):
      break

    fx2 = r8_erfc ( x )

    print ( '  %14.4g  %14.6g  %14.6g  %14.6g' % ( x, fx1, fx2, abs ( fx1 - fx2 ) ) )
#
#  Terminate.
#
  print ( '' )
  print ( 'R8_ERFC_TEST:' )
  print ( '  Normal end of execution.' )
  return

if ( __name__ == '__main__' ):
  from timestamp import timestamp
  timestamp ( )
  r8_erf_test ( )
  r8_erfc_test ( )
  timestamp ( )