module QUADPREC_ERRORFUNC
use QUADRUPLE_PRECISION
implicit none

private
public :: Q_ERF

contains

subroutine Q_ERF(X, ERF, ERFC)

!   Calculate the error function & its complement in quadruple precision.
!   0 <= x <= 1 uses series 7.1.5 from Abramowitz & Stegun
!   1 < x < 4 uses Chebyshev series
!   4 <= x uses a continued fraction
!
!   If x < 0 then erfc(x) = 1 + erf(|x|) and erf(x) = - erf(|x|)
!
!   WARNING: x must NOT occupy the same locations as any of
!            the other arguments if x < 0.
!
!   Programmer: Alan Miller    Alan.Miller @ vic.cmis.csiro.au
!               www.ozemail.com.au/~milleraj   FAX: (+61) 3-9545-8080
!   Latest revision - (Fortran 77 version) 11 May 1988
!   Latest revision - 4 February 1998

type (QUAD), intent(in)            :: X
type (QUAD), intent(out), optional :: ERF, ERFC

!     Chebyshev coefficients
type (QUAD), dimension(40), parameter :: COEFF =  &
(/ QUAD(0.4888515056984716e+00_DP,-0.1253349490474187e-15_DP), QUAD(-0.1358986289159460e+00_DP, 0.1828677302762691e-16_DP), &
   QUAD(0.3563511414853497e-01_DP,-0.1026978994630672e-17_DP), QUAD(-0.8884098618222492e-02_DP, 0.1869288268188562e-17_DP), &
   QUAD(0.2118423997553177e-02_DP, 0.1709987556862464e-18_DP), QUAD(-0.4853987640800370e-03_DP,-0.6960364726525852e-19_DP), &
   QUAD(0.1072737545377177e-03_DP, 0.1534550435363124e-19_DP), QUAD(-0.2293663448190487e-04_DP, 0.6928883131080608e-20_DP), &
   QUAD(0.4756874691575556e-05_DP,-0.1627977086918653e-20_DP), QUAD(-0.9589944832662770e-06_DP,-0.2249363452839194e-22_DP), &
   QUAD(0.1882900688295897e-06_DP,-0.1968876883111176e-22_DP), QUAD(-0.3606320662143092e-07_DP, 0.7483138202418590e-23_DP), &
   QUAD(0.6747602051937456e-08_DP,-0.1005646109123554e-23_DP), QUAD(-0.1234907042333986e-08_DP,-0.8553919151789978e-25_DP), &
   QUAD(0.2213150119337267e-09_DP,-0.8378636504284192e-27_DP), QUAD(-0.3887955444797815e-10_DP, 0.4599438293044550e-27_DP), &
   QUAD(0.6701388253274512e-11_DP,-0.4533239266031100e-27_DP), QUAD(-0.1134237426103359e-11_DP,-0.2592968654684961e-27_DP), &
   QUAD(0.1886555920176928e-12_DP,-0.9535048043067878e-29_DP), QUAD(-0.3085780489354876e-13_DP, 0.7801557018723816e-29_DP), &
   QUAD(0.4966710292058390e-14_DP,-0.1626370800837245e-29_DP), QUAD(-0.7871154824211008e-15_DP, 0.8389350837750800e-31_DP), &
   QUAD(0.1228887884514987e-15_DP,-0.5265492467954702e-31_DP), QUAD(-0.1891087633190591e-16_DP,-0.2388153131040173e-32_DP), &
   QUAD(0.2869742042877464e-17_DP, 0.7703719777548944e-34_DP), QUAD(-0.4296338289538412e-18_DP, 0.9629649721936182e-35_DP), &
   QUAD(0.6348319360536490e-19_DP,-0.7222237291452136e-35_DP), QUAD(-0.9261746469308854e-20_DP,-0.1805559322863034e-35_DP), &
   QUAD(0.1334626827242897e-20_DP, 0.2633107345841924e-36_DP), QUAD(-0.1900244302571750e-21_DP, 0.3761581922631320e-37_DP), &
   QUAD(0.2674133030089899e-22_DP,-0.4701977403289150e-38_DP), QUAD(-0.3720610299437181e-23_DP,-0.2938735877055719e-39_DP), &
   QUAD(0.5119522849758822e-24_DP,-0.1578156674465586e-56_DP), QUAD(-0.6968654298053806e-25_DP,-0.2295887403949780e-41_DP), &
   QUAD(0.9386117844916958e-26_DP,-0.2008901478456058e-41_DP), QUAD(-0.1251273661158375e-26_DP, 0.5022253696140144e-42_DP), &
   QUAD(0.1651286941713772e-27_DP,-0.3587324068671532e-43_DP), QUAD(-0.2157623168556909e-28_DP,-0.5605193857299268e-45_DP), &
   QUAD(0.2771990968956549e-29_DP,-0.3503246160812043e-45_DP), QUAD(-0.3468792422835851e-30_DP, 0.7006492321624086e-46_DP) /)

!       Local variables

type (QUAD)            :: D, DD, TERM, ARG, Y2, SV, XX, Y, TEMP, ERFF, ERFCC
real (kind=DP)         :: A
integer                :: M, J
logical                :: SMALL
type (QUAD), parameter :: QONE = QUAD(1.0_DP, 0.0_DP),  &
                          QZERO = QUAD(0.0_DP, 0.0_DP)

if (X%HI < 0.0_DP) then
  ARG = -X
else
  ARG = X
end if

if (ARG%HI <= 1.0_DP) then
!----------------------------------------------------------------------
!
!                  2x          x^2   x^4    x^6
!       erf(x) = --------.(1 - --- + ---- - ---- + .. )
!                sqrt(pi)      1.3   2!.5   3!.7
!
!
  SMALL = .false.
  TERM = TWO_ON_RTPI * ARG
  ERFF = TERM
  XX = ARG * ARG
  M = 1
  do
    if (SMALL) then
      TERM%HI = - TERM%HI * XX%HI / M
      TEMP%HI = TERM%HI / (2*M + 1)
      TEMP%LO = 0.0_DP
    else
      TERM = - TERM * XX
      select case (M)
        case (1)
          TERM = TERM
        case (2)
          TERM = scale(TERM, -1)
        case (4)
          TERM = scale(TERM, -2)
        case (8)
          TERM = scale(TERM, -3)
        case default
          TERM = TERM / real(M, kind=DP)
      end select
      TEMP = TERM / real(2*M+1, kind=DP)
    end if
    ERFF = ERFF + TEMP
    M = M + 1
    if (.not. SMALL .and. abs(TEMP%HI) < ERFF%HI * 1.0e-16_DP) then
      SMALL = .true.
    end if
    if (abs(TEMP%HI) < ERFF%HI * 1.0e-30_DP) then
      exit
    end if
  end do

  ERFCC = QONE - ERFF

!----------------------------------------------------------------------
!
!    Use Chebyshev series for 1 < x < 4.
!    Adapted from the routine CHEBEV from 'Numerical Recipes' by Press,
!    Flannery, Teukolsky & Vetterling, Cambridge Uni Press, 1986.
!
else if (ARG%HI > 1.0_DP .and. ARG%HI < 4.0_DP) then
  D = QZERO
  DD = QZERO
  Y = (ARG - 2.5_DP) / 1.5_DP
  Y2 = scale(Y, 1)
  do J = 40, 2, -1
    SV = D
    D = Y2 * D - DD + COEFF(J)
    DD = SV
  end do
  TERM = Y * D - DD + scale(COEFF(1), -1)

!     Finally multiply by exp(-x^2)

  ERFCC = TERM * exp(-X*X)
  ERFF = QONE - ERFCC

!----------------------------------------------------------------------
!
!       Use a continued fraction for x >= 4.0

else
  A = 0.5_DP * int(230.0_DP/ARG%HI) + 2.0_DP
  ERFCC = QZERO
  do
    ERFCC = ERFCC + ARG
    ERFCC = QUAD(A, 0.0_DP) / ERFCC
    A = A - 0.5_DP
    if (A < 0.1_DP) then
      exit
    end if
  end do
  ERFCC = ERFCC + ARG
  ERFCC = QONE / ERFCC

!       Multiply result by exp(-x**2) / sqrt(pi)

  ERFCC = ERFCC * exp(-X*X) / SQRTPI
  ERFF = QONE - ERFCC
end if

!----------------------------------------------------------------------
!       Replace erf & erfc if argument was negative.

if (X%HI < 0.0_DP) then
  ERFCC = QONE + ERFF
  ERFF = -ERFF
end if

if (present(ERF)) then
  ERF = ERFF
end if
if (present(ERFC)) then
  ERFC = ERFCC
end if

return
end subroutine Q_ERF

end module QUADPREC_ERRORFUNC


program TEST_Q_ERF
! Test quadruple-precision erf function by approximating its first derivative
! using divided differences.
! f'(x) = [f(x-2h) - 8.f(x-h) + 8.f(x+h) - f(x-2h)] / (12.h)  approx.

use QUADRUPLE_PRECISION
use QUADPREC_ERRORFUNC
implicit none

type (QUAD)            :: X, H, ARG, ERF_2, ERF_1, ERF1, ERF2, ERFC, AV_ERFC, &
                          D_EST, DERIV, ERROR
type (QUAD), parameter :: QONE = QUAD(1.0_DP, 0.0_DP)

do
  write(unit=*, fmt="(a)", advance="NO") " Enter x: "
  read(unit=*, fmt=*) X%HI
  X%LO = 0.0_DP
  H = scale(QONE, -24)
  if (abs(X%HI) < 0.477_DP) then       ! Use erf if erf(x) > 0.5
    ARG = X - 2.0_DP * H
    call Q_ERF(ARG, ERF_2, AV_ERFC)
    ARG = X - H
    call Q_ERF(ARG, ERF_1, ERFC)
    AV_ERFC = AV_ERFC + ERFC
    ARG = X + H
    call Q_ERF(ARG, ERF1, ERFC)
    AV_ERFC = AV_ERFC + ERFC
    ARG = X + 2.0_DP * H
    call Q_ERF(ARG, ERF2, ERFC)
    AV_ERFC = (AV_ERFC + ERFC) / 4.0_DP
    D_EST = (ERF_2 - ERF2 + scale( ERF1 - ERF_1, 3)) / (12.0_DP * H)

  else                                 ! Else use erfc(x)
    ARG = X - 2.0_DP * H
    call Q_ERF(ARG, ERFC=AV_ERFC)
    ERF_2 = AV_ERFC
    ARG = X - H
    call Q_ERF(ARG, ERFC=ERF_1)
    AV_ERFC = AV_ERFC + ERF_1
    ARG = X + H
    call Q_ERF(ARG, ERFC=ERF1)
    AV_ERFC = AV_ERFC + ERF1
    ARG = X + 2.0_DP * H
    call Q_ERF(ARG, ERFC=ERF2)
    AV_ERFC = (AV_ERFC + ERF2) / 4.0_DP
    D_EST = (ERF2 - ERF_2 + scale( ERF_1 - ERF1, 3)) / (12.0_DP * H)
  end if

  DERIV = TWO_ON_RTPI * exp(-X*X)
  ERROR = D_EST - DERIV
  write(unit=*, fmt="(a)") "      Estimate        Derivative         Error          erfc"
  write(unit=*, fmt="(2f20.16, g12.4, f20.16)")  &
                    D_EST%HI, DERIV%HI, ERROR%HI, AV_ERFC%HI
end do

stop

end program TEST_Q_ERF