module error_function 8,6
! This module provides generic interfaces for functions that evaluate
! erf(x), erfc(x), and exp(x*x)*erfc(x) in either single or double precision.
implicit none
private
save
! Public functions
public :: erf, erfc, erfcx
interface erf
module procedure erf_r4
module procedure derf
end interface
interface erfc
module procedure erfc_r4
module procedure derfc
end interface
interface erfcx
module procedure erfcx_r4
module procedure derfcx
end interface
! Private variables
integer, parameter :: r4 = selected_real_kind(6) ! 4 byte real
integer, parameter :: r8 = selected_real_kind(12) ! 8 byte real
contains
!------------------------------------------------------------------
!
! 6 December 2006 -- B. Eaton
! The following comments are from the original version of CALERF.
! The only changes in implementing this module are that the function
! names previously used for the single precision versions have been
! adopted for the new generic interfaces. To support these interfaces
! there is now both a single precision version (calerf_r4) and a
! double precision version (calerf_r8) of CALERF below. These versions
! are hardcoded to use IEEE arithmetic.
!
!------------------------------------------------------------------
!
! This packet evaluates erf(x), erfc(x), and exp(x*x)*erfc(x)
! for a real argument x. It contains three FUNCTION type
! subprograms: ERF, ERFC, and ERFCX (or DERF, DERFC, and DERFCX),
! and one SUBROUTINE type subprogram, CALERF. The calling
! statements for the primary entries are:
!
! Y=ERF(X) (or Y=DERF(X)),
!
! Y=ERFC(X) (or Y=DERFC(X)),
! and
! Y=ERFCX(X) (or Y=DERFCX(X)).
!
! The routine CALERF is intended for internal packet use only,
! all computations within the packet being concentrated in this
! routine. The function subprograms invoke CALERF with the
! statement
!
! CALL CALERF(ARG,RESULT,JINT)
!
! where the parameter usage is as follows
!
! Function Parameters for CALERF
! call ARG Result JINT
!
! ERF(ARG) ANY REAL ARGUMENT ERF(ARG) 0
! ERFC(ARG) ABS(ARG) .LT. XBIG ERFC(ARG) 1
! ERFCX(ARG) XNEG .LT. ARG .LT. XMAX ERFCX(ARG) 2
!
! The main computation evaluates near-minimax approximations
! from "Rational Chebyshev approximations for the error function"
! by W. J. Cody, Math. Comp., 1969, PP. 631-638. This
! transportable program uses rational functions that theoretically
! approximate erf(x) and erfc(x) to at least 18 significant
! decimal digits. The accuracy achieved depends on the arithmetic
! system, the compiler, the intrinsic functions, and proper
! selection of the machine-dependent constants.
!
!*******************************************************************
!*******************************************************************
!
! Explanation of machine-dependent constants
!
! XMIN = the smallest positive floating-point number.
! XINF = the largest positive finite floating-point number.
! XNEG = the largest negative argument acceptable to ERFCX;
! the negative of the solution to the equation
! 2*exp(x*x) = XINF.
! XSMALL = argument below which erf(x) may be represented by
! 2*x/sqrt(pi) and above which x*x will not underflow.
! A conservative value is the largest machine number X
! such that 1.0 + X = 1.0 to machine precision.
! XBIG = largest argument acceptable to ERFC; solution to
! the equation: W(x) * (1-0.5/x**2) = XMIN, where
! W(x) = exp(-x*x)/[x*sqrt(pi)].
! XHUGE = argument above which 1.0 - 1/(2*x*x) = 1.0 to
! machine precision. A conservative value is
! 1/[2*sqrt(XSMALL)]
! XMAX = largest acceptable argument to ERFCX; the minimum
! of XINF and 1/[sqrt(pi)*XMIN].
!
! Approximate values for some important machines are:
!
! XMIN XINF XNEG XSMALL
!
! CDC 7600 (S.P.) 3.13E-294 1.26E+322 -27.220 7.11E-15
! CRAY-1 (S.P.) 4.58E-2467 5.45E+2465 -75.345 7.11E-15
! IEEE (IBM/XT,
! SUN, etc.) (S.P.) 1.18E-38 3.40E+38 -9.382 5.96E-8
! IEEE (IBM/XT,
! SUN, etc.) (D.P.) 2.23D-308 1.79D+308 -26.628 1.11D-16
! IBM 195 (D.P.) 5.40D-79 7.23E+75 -13.190 1.39D-17
! UNIVAC 1108 (D.P.) 2.78D-309 8.98D+307 -26.615 1.73D-18
! VAX D-Format (D.P.) 2.94D-39 1.70D+38 -9.345 1.39D-17
! VAX G-Format (D.P.) 5.56D-309 8.98D+307 -26.615 1.11D-16
!
!
! XBIG XHUGE XMAX
!
! CDC 7600 (S.P.) 25.922 8.39E+6 1.80X+293
! CRAY-1 (S.P.) 75.326 8.39E+6 5.45E+2465
! IEEE (IBM/XT,
! SUN, etc.) (S.P.) 9.194 2.90E+3 4.79E+37
! IEEE (IBM/XT,
! SUN, etc.) (D.P.) 26.543 6.71D+7 2.53D+307
! IBM 195 (D.P.) 13.306 1.90D+8 7.23E+75
! UNIVAC 1108 (D.P.) 26.582 5.37D+8 8.98D+307
! VAX D-Format (D.P.) 9.269 1.90D+8 1.70D+38
! VAX G-Format (D.P.) 26.569 6.71D+7 8.98D+307
!
!*******************************************************************
!*******************************************************************
!
! Error returns
!
! The program returns ERFC = 0 for ARG .GE. XBIG;
!
! ERFCX = XINF for ARG .LT. XNEG;
! and
! ERFCX = 0 for ARG .GE. XMAX.
!
!
! Intrinsic functions required are:
!
! ABS, AINT, EXP
!
!
! Author: W. J. Cody
! Mathematics and Computer Science Division
! Argonne National Laboratory
! Argonne, IL 60439
!
! Latest modification: March 19, 1990
!
!------------------------------------------------------------------
SUBROUTINE CALERF_r8(ARG, RESULT, JINT) 3
!------------------------------------------------------------------
! This version uses 8-byte reals
!------------------------------------------------------------------
integer, parameter :: rk = r8
! arguments
real(rk), intent(in) :: arg
integer, intent(in) :: jint
real(rk), intent(out) :: result
! local variables
INTEGER :: I
real(rk) :: X, Y, YSQ, XNUM, XDEN, DEL
!------------------------------------------------------------------
! Mathematical constants
!------------------------------------------------------------------
real(rk), parameter :: ZERO = 0.0E0_rk
real(rk), parameter :: FOUR = 4.0E0_rk
real(rk), parameter :: ONE = 1.0E0_rk
real(rk), parameter :: HALF = 0.5E0_rk
real(rk), parameter :: TWO = 2.0E0_rk
real(rk), parameter :: SQRPI = 5.6418958354775628695E-1_rk
real(rk), parameter :: THRESH = 0.46875E0_rk
real(rk), parameter :: SIXTEN = 16.0E0_rk
!------------------------------------------------------------------
! Machine-dependent constants: IEEE single precision values
!------------------------------------------------------------------
!S real, parameter :: XINF = 3.40E+38
!S real, parameter :: XNEG = -9.382E0
!S real, parameter :: XSMALL = 5.96E-8
!S real, parameter :: XBIG = 9.194E0
!S real, parameter :: XHUGE = 2.90E3
!S real, parameter :: XMAX = 4.79E37
!------------------------------------------------------------------
! Machine-dependent constants: IEEE double precision values
!------------------------------------------------------------------
real(rk), parameter :: XINF = 1.79E308_r8
real(rk), parameter :: XNEG = -26.628E0_r8
real(rk), parameter :: XSMALL = 1.11E-16_r8
real(rk), parameter :: XBIG = 26.543E0_r8
real(rk), parameter :: XHUGE = 6.71E7_r8
real(rk), parameter :: XMAX = 2.53E307_r8
!------------------------------------------------------------------
! Coefficients for approximation to erf in first interval
!------------------------------------------------------------------
real(rk), parameter :: A(5) = (/ 3.16112374387056560E00_rk, 1.13864154151050156E02_rk, &
3.77485237685302021E02_rk, 3.20937758913846947E03_rk, &
1.85777706184603153E-1_rk /)
real(rk), parameter :: B(4) = (/ 2.36012909523441209E01_rk, 2.44024637934444173E02_rk, &
1.28261652607737228E03_rk, 2.84423683343917062E03_rk /)
!------------------------------------------------------------------
! Coefficients for approximation to erfc in second interval
!------------------------------------------------------------------
real(rk), parameter :: C(9) = (/ 5.64188496988670089E-1_rk, 8.88314979438837594E00_rk, &
6.61191906371416295E01_rk, 2.98635138197400131E02_rk, &
8.81952221241769090E02_rk, 1.71204761263407058E03_rk, &
2.05107837782607147E03_rk, 1.23033935479799725E03_rk, &
2.15311535474403846E-8_rk /)
real(rk), parameter :: D(8) = (/ 1.57449261107098347E01_rk, 1.17693950891312499E02_rk, &
5.37181101862009858E02_rk, 1.62138957456669019E03_rk, &
3.29079923573345963E03_rk, 4.36261909014324716E03_rk, &
3.43936767414372164E03_rk, 1.23033935480374942E03_rk /)
!------------------------------------------------------------------
! Coefficients for approximation to erfc in third interval
!------------------------------------------------------------------
real(rk), parameter :: P(6) = (/ 3.05326634961232344E-1_rk, 3.60344899949804439E-1_rk, &
1.25781726111229246E-1_rk, 1.60837851487422766E-2_rk, &
6.58749161529837803E-4_rk, 1.63153871373020978E-2_rk /)
real(rk), parameter :: Q(5) = (/ 2.56852019228982242E00_rk, 1.87295284992346047E00_rk, &
5.27905102951428412E-1_rk, 6.05183413124413191E-2_rk, &
2.33520497626869185E-3_rk /)
!------------------------------------------------------------------
X = ARG
Y = ABS(X)
IF (Y .LE. THRESH) THEN
!------------------------------------------------------------------
! Evaluate erf for |X| <= 0.46875
!------------------------------------------------------------------
YSQ = ZERO
IF (Y .GT. XSMALL) YSQ = Y * Y
XNUM = A(5)*YSQ
XDEN = YSQ
DO I = 1, 3
XNUM = (XNUM + A(I)) * YSQ
XDEN = (XDEN + B(I)) * YSQ
end do
RESULT = X * (XNUM + A(4)) / (XDEN + B(4))
IF (JINT .NE. 0) RESULT = ONE - RESULT
IF (JINT .EQ. 2) RESULT = EXP(YSQ) * RESULT
GO TO 80
ELSE IF (Y .LE. FOUR) THEN
!------------------------------------------------------------------
! Evaluate erfc for 0.46875 <= |X| <= 4.0
!------------------------------------------------------------------
XNUM = C(9)*Y
XDEN = Y
DO I = 1, 7
XNUM = (XNUM + C(I)) * Y
XDEN = (XDEN + D(I)) * Y
end do
RESULT = (XNUM + C(8)) / (XDEN + D(8))
IF (JINT .NE. 2) THEN
YSQ = AINT(Y*SIXTEN)/SIXTEN
DEL = (Y-YSQ)*(Y+YSQ)
RESULT = EXP(-YSQ*YSQ) * EXP(-DEL) * RESULT
END IF
ELSE
!------------------------------------------------------------------
! Evaluate erfc for |X| > 4.0
!------------------------------------------------------------------
RESULT = ZERO
IF (Y .GE. XBIG) THEN
IF ((JINT .NE. 2) .OR. (Y .GE. XMAX)) GO TO 30
IF (Y .GE. XHUGE) THEN
RESULT = SQRPI / Y
GO TO 30
END IF
END IF
YSQ = ONE / (Y * Y)
XNUM = P(6)*YSQ
XDEN = YSQ
DO I = 1, 4
XNUM = (XNUM + P(I)) * YSQ
XDEN = (XDEN + Q(I)) * YSQ
end do
RESULT = YSQ *(XNUM + P(5)) / (XDEN + Q(5))
RESULT = (SQRPI - RESULT) / Y
IF (JINT .NE. 2) THEN
YSQ = AINT(Y*SIXTEN)/SIXTEN
DEL = (Y-YSQ)*(Y+YSQ)
RESULT = EXP(-YSQ*YSQ) * EXP(-DEL) * RESULT
END IF
END IF
30 continue
!------------------------------------------------------------------
! Fix up for negative argument, erf, etc.
!------------------------------------------------------------------
IF (JINT .EQ. 0) THEN
RESULT = (HALF - RESULT) + HALF
IF (X .LT. ZERO) RESULT = -RESULT
ELSE IF (JINT .EQ. 1) THEN
IF (X .LT. ZERO) RESULT = TWO - RESULT
ELSE
IF (X .LT. ZERO) THEN
IF (X .LT. XNEG) THEN
RESULT = XINF
ELSE
YSQ = AINT(X*SIXTEN)/SIXTEN
DEL = (X-YSQ)*(X+YSQ)
Y = EXP(YSQ*YSQ) * EXP(DEL)
RESULT = (Y+Y) - RESULT
END IF
END IF
END IF
80 continue
end SUBROUTINE CALERF_r8
!------------------------------------------------------------------------------------------
SUBROUTINE CALERF_r4(ARG, RESULT, JINT) 3
!------------------------------------------------------------------
! This version uses 4-byte reals
!------------------------------------------------------------------
integer, parameter :: rk = r4
! arguments
real(rk), intent(in) :: arg
integer, intent(in) :: jint
real(rk), intent(out) :: result
! local variables
INTEGER :: I
real(rk) :: X, Y, YSQ, XNUM, XDEN, DEL
!------------------------------------------------------------------
! Mathematical constants
!------------------------------------------------------------------
real(rk), parameter :: ZERO = 0.0E0_rk
real(rk), parameter :: FOUR = 4.0E0_rk
real(rk), parameter :: ONE = 1.0E0_rk
real(rk), parameter :: HALF = 0.5E0_rk
real(rk), parameter :: TWO = 2.0E0_rk
real(rk), parameter :: SQRPI = 5.6418958354775628695E-1_rk
real(rk), parameter :: THRESH = 0.46875E0_rk
real(rk), parameter :: SIXTEN = 16.0E0_rk
!------------------------------------------------------------------
! Machine-dependent constants: IEEE single precision values
!------------------------------------------------------------------
real(rk), parameter :: XINF = 3.40E+38_r4
real(rk), parameter :: XNEG = -9.382E0_r4
real(rk), parameter :: XSMALL = 5.96E-8_r4
real(rk), parameter :: XBIG = 9.194E0_r4
real(rk), parameter :: XHUGE = 2.90E3_r4
real(rk), parameter :: XMAX = 4.79E37_r4
!------------------------------------------------------------------
! Coefficients for approximation to erf in first interval
!------------------------------------------------------------------
real(rk), parameter :: A(5) = (/ 3.16112374387056560E00_rk, 1.13864154151050156E02_rk, &
3.77485237685302021E02_rk, 3.20937758913846947E03_rk, &
1.85777706184603153E-1_rk /)
real(rk), parameter :: B(4) = (/ 2.36012909523441209E01_rk, 2.44024637934444173E02_rk, &
1.28261652607737228E03_rk, 2.84423683343917062E03_rk /)
!------------------------------------------------------------------
! Coefficients for approximation to erfc in second interval
!------------------------------------------------------------------
real(rk), parameter :: C(9) = (/ 5.64188496988670089E-1_rk, 8.88314979438837594E00_rk, &
6.61191906371416295E01_rk, 2.98635138197400131E02_rk, &
8.81952221241769090E02_rk, 1.71204761263407058E03_rk, &
2.05107837782607147E03_rk, 1.23033935479799725E03_rk, &
2.15311535474403846E-8_rk /)
real(rk), parameter :: D(8) = (/ 1.57449261107098347E01_rk, 1.17693950891312499E02_rk, &
5.37181101862009858E02_rk, 1.62138957456669019E03_rk, &
3.29079923573345963E03_rk, 4.36261909014324716E03_rk, &
3.43936767414372164E03_rk, 1.23033935480374942E03_rk /)
!------------------------------------------------------------------
! Coefficients for approximation to erfc in third interval
!------------------------------------------------------------------
real(rk), parameter :: P(6) = (/ 3.05326634961232344E-1_rk, 3.60344899949804439E-1_rk, &
1.25781726111229246E-1_rk, 1.60837851487422766E-2_rk, &
6.58749161529837803E-4_rk, 1.63153871373020978E-2_rk /)
real(rk), parameter :: Q(5) = (/ 2.56852019228982242E00_rk, 1.87295284992346047E00_rk, &
5.27905102951428412E-1_rk, 6.05183413124413191E-2_rk, &
2.33520497626869185E-3_rk /)
!------------------------------------------------------------------
X = ARG
Y = ABS(X)
IF (Y .LE. THRESH) THEN
!------------------------------------------------------------------
! Evaluate erf for |X| <= 0.46875
!------------------------------------------------------------------
YSQ = ZERO
IF (Y .GT. XSMALL) YSQ = Y * Y
XNUM = A(5)*YSQ
XDEN = YSQ
DO I = 1, 3
XNUM = (XNUM + A(I)) * YSQ
XDEN = (XDEN + B(I)) * YSQ
end do
RESULT = X * (XNUM + A(4)) / (XDEN + B(4))
IF (JINT .NE. 0) RESULT = ONE - RESULT
IF (JINT .EQ. 2) RESULT = EXP(YSQ) * RESULT
GO TO 80
ELSE IF (Y .LE. FOUR) THEN
!------------------------------------------------------------------
! Evaluate erfc for 0.46875 <= |X| <= 4.0
!------------------------------------------------------------------
XNUM = C(9)*Y
XDEN = Y
DO I = 1, 7
XNUM = (XNUM + C(I)) * Y
XDEN = (XDEN + D(I)) * Y
end do
RESULT = (XNUM + C(8)) / (XDEN + D(8))
IF (JINT .NE. 2) THEN
YSQ = AINT(Y*SIXTEN)/SIXTEN
DEL = (Y-YSQ)*(Y+YSQ)
RESULT = EXP(-YSQ*YSQ) * EXP(-DEL) * RESULT
END IF
ELSE
!------------------------------------------------------------------
! Evaluate erfc for |X| > 4.0
!------------------------------------------------------------------
RESULT = ZERO
IF (Y .GE. XBIG) THEN
IF ((JINT .NE. 2) .OR. (Y .GE. XMAX)) GO TO 30
IF (Y .GE. XHUGE) THEN
RESULT = SQRPI / Y
GO TO 30
END IF
END IF
YSQ = ONE / (Y * Y)
XNUM = P(6)*YSQ
XDEN = YSQ
DO I = 1, 4
XNUM = (XNUM + P(I)) * YSQ
XDEN = (XDEN + Q(I)) * YSQ
end do
RESULT = YSQ *(XNUM + P(5)) / (XDEN + Q(5))
RESULT = (SQRPI - RESULT) / Y
IF (JINT .NE. 2) THEN
YSQ = AINT(Y*SIXTEN)/SIXTEN
DEL = (Y-YSQ)*(Y+YSQ)
RESULT = EXP(-YSQ*YSQ) * EXP(-DEL) * RESULT
END IF
END IF
30 continue
!------------------------------------------------------------------
! Fix up for negative argument, erf, etc.
!------------------------------------------------------------------
IF (JINT .EQ. 0) THEN
RESULT = (HALF - RESULT) + HALF
IF (X .LT. ZERO) RESULT = -RESULT
ELSE IF (JINT .EQ. 1) THEN
IF (X .LT. ZERO) RESULT = TWO - RESULT
ELSE
IF (X .LT. ZERO) THEN
IF (X .LT. XNEG) THEN
RESULT = XINF
ELSE
YSQ = AINT(X*SIXTEN)/SIXTEN
DEL = (X-YSQ)*(X+YSQ)
Y = EXP(YSQ*YSQ) * EXP(DEL)
RESULT = (Y+Y) - RESULT
END IF
END IF
END IF
80 continue
end SUBROUTINE CALERF_r4
!------------------------------------------------------------------------------------------
FUNCTION DERF(X) 1,1
!--------------------------------------------------------------------
!
! This subprogram computes approximate values for erf(x).
! (see comments heading CALERF).
!
! Author/date: W. J. Cody, January 8, 1985
!
!--------------------------------------------------------------------
integer, parameter :: rk = r8 ! 8 byte real
! argument
real(rk), intent(in) :: X
! return value
real(rk) :: DERF
! local variables
INTEGER :: JINT = 0
!------------------------------------------------------------------
CALL CALERF_r8(X, DERF, JINT)
END FUNCTION DERF
!------------------------------------------------------------------------------------------
FUNCTION ERF_r4(X) 1,1
!--------------------------------------------------------------------
!
! This subprogram computes approximate values for erf(x).
! (see comments heading CALERF).
!
! Author/date: W. J. Cody, January 8, 1985
!
!--------------------------------------------------------------------
integer, parameter :: rk = r4 ! 4 byte real
! argument
real(rk), intent(in) :: X
! return value
real(rk) :: ERF_r4
! local variables
INTEGER :: JINT = 0
!------------------------------------------------------------------
CALL CALERF_r4(X, ERF_r4, JINT)
END FUNCTION ERF_r4
!------------------------------------------------------------------------------------------
FUNCTION DERFC(X) 1,1
!--------------------------------------------------------------------
!
! This subprogram computes approximate values for erfc(x).
! (see comments heading CALERF).
!
! Author/date: W. J. Cody, January 8, 1985
!
!--------------------------------------------------------------------
integer, parameter :: rk = r8 ! 8 byte real
! argument
real(rk), intent(in) :: X
! return value
real(rk) :: DERFC
! local variables
INTEGER :: JINT = 1
!------------------------------------------------------------------
CALL CALERF_r8(X, DERFC, JINT)
END FUNCTION DERFC
!------------------------------------------------------------------------------------------
FUNCTION ERFC_r4(X) 1,1
!--------------------------------------------------------------------
!
! This subprogram computes approximate values for erfc(x).
! (see comments heading CALERF).
!
! Author/date: W. J. Cody, January 8, 1985
!
!--------------------------------------------------------------------
integer, parameter :: rk = r4 ! 4 byte real
! argument
real(rk), intent(in) :: X
! return value
real(rk) :: ERFC_r4
! local variables
INTEGER :: JINT = 1
!------------------------------------------------------------------
CALL CALERF_r4(X, ERFC_r4, JINT)
END FUNCTION ERFC_r4
!------------------------------------------------------------------------------------------
FUNCTION DERFCX(X) 1,1
!--------------------------------------------------------------------
!
! This subprogram computes approximate values for exp(x*x) * erfc(x).
! (see comments heading CALERF).
!
! Author/date: W. J. Cody, March 30, 1987
!
!--------------------------------------------------------------------
integer, parameter :: rk = r8 ! 8 byte real
! argument
real(rk), intent(in) :: X
! return value
real(rk) :: DERFCX
! local variables
INTEGER :: JINT = 2
!------------------------------------------------------------------
CALL CALERF_r8(X, DERFCX, JINT)
END FUNCTION DERFCX
!------------------------------------------------------------------------------------------
FUNCTION ERFCX_R4(X) 1,1
!--------------------------------------------------------------------
!
! This subprogram computes approximate values for exp(x*x) * erfc(x).
! (see comments heading CALERF).
!
! Author/date: W. J. Cody, March 30, 1987
!
!--------------------------------------------------------------------
integer, parameter :: rk = r4 ! 8 byte real
! argument
real(rk), intent(in) :: X
! return value
real(rk) :: ERFCX_R4
! local variables
INTEGER :: JINT = 2
!------------------------------------------------------------------
CALL CALERF_r4(X, ERFCX_R4, JINT)
END FUNCTION ERFCX_R4
!------------------------------------------------------------------------------------------
end module error_function