!*******************************************************************************
! Obtain erfc(x) (and erf(x)) to the maximum machine accuracy possible.
!
! REFERENCE: Irene A. Stegun and Ruth Zucker, J. Res. NBS  74B, 211-224 (1970).
!            ``Automatic Computing Methods for Special Functions'',
!            Part 1. Error, Probability, and Related Functions
!
!            CARDS GIVEN AS DOUBLE PRECISION CONSTANTS, YIELDS DOUBLE
!            PRECISION RESULTS ON THE UNIVAC 1108 SYSTEM.
!
!            SINGLE PRECISION RESULTS ARE OBTAINED BY
!            (1) DELETING THE TYPE STATEMENT
!            (2) SETTING NBC=8, NBM=27 ON THE FIRST DATA CARD
!            (3) CHANGING THE D'S TO E'S ON THE SECOND AND THIRD DATA
!            CARDS.
!            (4) CHANGING THE FUNCTION TYPE - DABS TO ABS, DEXP TO EXP
!
!            FOR OTHER VALUES OF NBM, THE CONSTANT ON THE THIRD DATA
!            CARD SHOULD BE GIVEN TO THE CORRESPONDING PRECISION.
!
! NBC:       Number of binary digits in the characteristic of a floating 
!            point number.
!            Single precision: NBC = 8, Double precision, NBC = 11

! NBM:       Accuracy desired or maximum number of binary digits in the 
!            mantissa of a floating point number.
!            Single precision: NBM = 11, Double precision, NBM = 60
!
! METHOD:    Power series for ABS(X) <= 1 (= ULPS, upper limit for series).
!            Continued fraction for  ULPS < ABS(X) <= ULCF 
!            (ULCF = upper limit of continued fraction)
!
! RANGE:     ABS(X) <= ULCF ; ULCF = 0.83*(2.0**((NBC-1)/2))
!            ABS(ULCF) approximately    9.3 for NBC = 8
!                                      26.5 for NBC = 11
!            Beyond this range the limiting values of the functions are
!            supplied: ERF(INF) = 1.0,   ERFC(INF) = 0.0
!
! ACCURACY:  Routine will return (NBM-I-3) significant binary digits
!            where I is the number of binary digits representing the
!            integer part of X**2. (This is essentially the accuracy of
!            the exponential routine.)
!
! PRECISION: Variable - by setting desired NBM
!*******************************************************************************

      function erfc05(x) result(erfc)

      use precision, only: wp

      implicit none
      real(wp), intent(in)    :: x
      real(wp)                :: erf,erfc
      real(wp)                :: an,bn,c1,dn,f,fn,  &
                                 fnm1,fnm2,gn,gnm1,gnm2,  &
                                 prev,pwr,rnbc,scf,sum,tn,toler, &
                                 ulcf,wn,y,ysq
      integer, parameter      :: nbc = 11
      integer, parameter      :: nbm = 60
      real(wp), parameter     :: zero = 0.0_wp
      real(wp), parameter     :: one = 1.0_wp
      real(wp), parameter     :: two = 2.0_wp
      real(wp), parameter     :: four = 4.0_wp
      real(wp), parameter     :: ulps = 1.0_wp
      real(wp), parameter     :: cons = 0.83_wp
      real(wp), parameter     :: trrtpi = 1.128379167095512574_wp


      rnbc = nbc
      toler = two**(-nbm)

!     Test on zero
      if (x == zero) then
          erf = zero
          erfc = one
          return
      endif

      y = abs(x)
      ysq = y**2

      if (y <= ulps) then
!         Method: power series
          sum = zero
          dn = one
          tn = one
          pwr = two*ysq
          do
              dn = dn + two
              tn = pwr*tn/dn
              sum = tn + sum
!             Tolerance check
              if (tn < toler) exit
          enddo
          erf = (sum + one)*trrtpi*y*exp(-ysq)
          erfc = one - erf
      else
!         Maximum argument
          c1 = two**((rnbc-one)/two)
          ulcf = cons*c1

!         Scale factor
          scf = two**(c1**2 - rnbc)

!         Limiting value
          if (y > ulcf) then
              erf = one
              erfc = zero
          else
!             Method: continued fraction
              fnm2 = zero
              gnm2 = one
              fnm1 = two*y
              gnm1 = two*ysq + one
              prev = fnm1/gnm1
              wn = one
              bn = gnm1 + four
              do
                  an = -wn*(wn + one)
                  fn = bn*fnm1 + an*fnm2
                  gn = bn*gnm1 + an*gnm2
                  f = fn/gn

!                 Tolerance check
                  if (abs(one - (f/prev)) <= toler) exit
                  if (prev > f) then
                      f = prev
                      exit
                  endif

!                 Both fn and gn must be tested if abs(x) < 0.61
                  if (gn >= scf) then
!                     Scaling
                      fn = fn/scf
                      gn = gn/scf
                      fnm1 = fnm1/scf
                      gnm1 = gnm1/scf
                  endif

                  fnm2 = fnm1
                  gnm2 = gnm1
                  fnm1 = fn
                  gnm1 = gn

                  wn = wn + two
                  bn = bn + four
                  prev = f
              enddo
              erfc = f*exp(-ysq)*trrtpi/two
              erf = one - erfc
          endif
      endif

!     Negative argument
      if (x < zero) erfc = two - erfc
!     if (x < zero) erf = -erf

      end function erfc05