!C $ f2py -m erf!Coremod -h erf!Coremod.pyf erf!Coremod.f
!C f2py erf!Coremod.pyf erf!Coremod.f -!C --f!Compiler=gnu95 --!Compiler=mingw32 -lmsv!Cr71
!C $ f2py --f!Compiler=gnu95 --!Compiler=mingw32 -lmsv!Cr71 -m erf!Coremod -!C erf!Coremod.f
!C gfortran -fPI!C -!C erf!Coremod.f
!C f2py -m erf!Coremod -DUPPER!CASE_FORTRAN -!C erf!Coremod.o erf!Coremod_interfa!Ce.f
MODULE ERFCOREMOD
!C IMPLI!CIT NONE
!C INTERFA!CE !CALERF
!C MODULE PRO!CEDURE !CALERF
!C END INTERFA!CE
!C INTERFA!CE DERF
!C MODULE PRO!CEDURE DERF
!C END INTERFA!CE
!C INTERFA!CE DERF!C
!C MODULE PRO!CEDURE DERF!C
!C END INTERFA!CE
!C INTERFA!CE DERF!CX
!C MODULE PRO!CEDURE DERF!CX
!C END INTERFA!CE
CONTAINS
!C--------------------------------------------------------------------
!C
!C DERF subprogram !Computes approximate values for erf(x).
!C (see !Comments heading !CALERF).
!C
!C Author/date: W. J. !Cody, January 8, 1985
!C
!C--------------------------------------------------------------------
FUNCTION DERF( X ) RESULT (VALUE)
IMPLICIT NONE
DOUBLE PRECISION, INTENT(IN) :: X
DOUBLE PRECISION :: VALUE
INTEGER, PARAMETER :: JINT = 0
CALL CALERF(X,VALUE,JINT)
RETURN
END FUNCTION DERF
!C--------------------------------------------------------------------
!C
!C DERF!C subprogram !Computes approximate values for erf!C(x).
!C (see !Comments heading !CALERF).
!C
!C Author/date: W. J. !Cody, January 8, 1985
!C
!C--------------------------------------------------------------------
FUNCTION DERFC( X ) RESULT (VALUE)
IMPLICIT NONE
DOUBLE PRECISION, INTENT(IN) :: X
DOUBLE PRECISION :: VALUE
INTEGER, PARAMETER :: JINT = 1
CALL CALERF(X,VALUE,JINT)
RETURN
END FUNCTION DERFC
!C------------------------------------------------------------------
!C
!C DERFCX subprogram Computes approximate values for exp(x*x) * erfC(x).
!C (see !Comments heading !CALERF).
!C
!C Author/date: W. J. !Cody, Mar!Ch 30, 1987
!C
!C------------------------------------------------------------------
FUNCTION DERFCX( X ) RESULT (VALUE)
IMPLICIT NONE
DOUBLE PRECISION, INTENT(IN) :: X
DOUBLE PRECISION :: VALUE
INTEGER, PARAMETER :: JINT = 2
CALL CALERF(X,VALUE,JINT)
RETURN
END FUNCTION DERFCX
SUBROUTINE CALERF(ARG,RESULT,JINT)
IMPLICIT NONE
!C------------------------------------------------------------------
!C
!C !CALERF pa!Cket evaluates erf(x), erf!C(x), and exp(x*x)*erf!C(x)
!C for a real argument x. It !Contains three FUN!CTION type
!C subprograms: ERF, ERF!C, and ERF!CX (or DERF, DERF!C, and DERF!CX),
!C and one SUBROUTINE type subprogram, !CALERF. The !Calling
!C statements for the primary entries are:
!C
!C Y=ERF(X) (or Y=DERF(X)),
!C
!C Y=ERF!C(X) (or Y=DERF!C(X)),
!C and
!C Y=ERF!CX(X) (or Y=DERF!CX(X)).
!C
!C The routine !CALERF is intended for internal pa!Cket use only,
!C all !Computations within the pa!Cket being !Con!Centrated in this
!C routine. The fun!Ction subprograms invoke !CALERF with the
!C statement
!C
!C !CALL !CALERF(ARG,RESULT,JINT)
!C
!C where the parameter usage is as follows
!C
!C Fun!Ction Parameters for !CALERF
!C !Call ARG Result JINT
!C
!C ERF(ARG) ANY REAL ARGUMENT ERF(ARG) 0
!C ERF!C(ARG) ABS(ARG) .LT. XBIG ERF!C(ARG) 1
!C ERF!CX(ARG) XNEG .LT. ARG .LT. XMAX ERF!CX(ARG) 2
!C
!C The main !Computation evaluates near-minimax approximations
!C from "Rational !Chebyshev approximations for the error fun!Ction"
!C by W. J. !Cody, Math. !Comp., 1969, PP. 631-638. This
!C transportable program uses rational fun!Ctions that theoreti!Cally
!C approximate erf(x) and erf!C(x) to at least 18 signifi!Cant
!C de!Cimal digits. The a!C!Cura!Cy a!Chieved depends on the arithmeti!C
!C system, the !Compiler, the intrinsi!C fun!Ctions, and proper
!C sele!Ction of the ma!Chine-dependent !Constants.
!C
!C*******************************************************************
!C*******************************************************************
!C
!C Explanation of ma!Chine-dependent !Constants
!C
!C XMIN = the smallest positive floating-point number.
!C XINF = the largest positive finite floating-point number.
!C XNEG = the largest negative argument a!C!Ceptable to ERF!CX;
!C the negative of the solution to the equation
!C 2*exp(x*x) = XINF.
!C XSMALL = argument below whi!Ch erf(x) may be represented by
!C 2*x/sqrt(pi) and above whi!Ch x*x will not underflow.
!C A !Conservative value is the largest ma!Chine number X
!C su!Ch that 1.0 + X = 1.0 to ma!Chine pre!Cision.
!C XBIG = largest argument a!C!Ceptable to ERF!C; solution to
!C the equation: W(x) * (1-0.5/x**2) = XMIN, where
!C W(x) = exp(-x*x)/[x*sqrt(pi)].
!C XHUGE = argument above whi!Ch 1.0 - 1/(2*x*x) = 1.0 to
!C ma!Chine pre!Cision. A !Conservative value is
!C 1/[2*sqrt(XSMALL)]
!C XMAX = largest a!C!Ceptable argument to ERF!CX; the minimum
!C of XINF and 1/[sqrt(pi)*XMIN].
!C
!C Approximate values for some important ma!Chines are:
!C
!C XMIN XINF XNEG XSMALL
!C
!C !C 7600 (S.P.) 3.13E-294 1.26E+322 -27.220 7.11E-15
!C !CRAY-1 (S.P.) 4.58E-2467 5.45E+2465 -75.345 7.11E-15
!C IEEE (IBM/XT,
!C SUN, et!C.) (S.P.) 1.18E-38 3.40E+38 -9.382 5.96E-8
!C IEEE (IBM/XT,
!C SUN, et!C.) (D.P.) 2.23D-308 1.79D+308 -26.628 1.11D-16
!C IBM 195 (D.P.) 5.40D-79 7.23E+75 -13.190 1.39D-17
!C UNIVA!C 1108 (D.P.) 2.78D-309 8.98D+307 -26.615 1.73D-18
!C VAX D-Format (D.P.) 2.94D-39 1.70D+38 -9.345 1.39D-17
!C VAX G-Format (D.P.) 5.56D-309 8.98D+307 -26.615 1.11D-16
!C
!C
!C XBIG XHUGE XMAX
!C
!C !C 7600 (S.P.) 25.922 8.39E+6 1.80X+293
!C !CRAY-1 (S.P.) 75.326 8.39E+6 5.45E+2465
!C IEEE (IBM/XT,
!C SUN, et!C.) (S.P.) 9.194 2.90E+3 4.79E+37
!C IEEE (IBM/XT,
!C SUN, et!C.) (D.P.) 26.543 6.71D+7 2.53D+307
!C IBM 195 (D.P.) 13.306 1.90D+8 7.23E+75
!C UNIVA!C 1108 (D.P.) 26.582 5.37D+8 8.98D+307
!C VAX D-Format (D.P.) 9.269 1.90D+8 1.70D+38
!C VAX G-Format (D.P.) 26.569 6.71D+7 8.98D+307
!C
!C*******************************************************************
!C*******************************************************************
!C
!C Error returns
!C
!C The program returns ERF!C = 0 for ARG .GE. XBIG;
!C
!C ERF!CX = XINF for ARG .LT. XNEG;
!C and
!C ERF!CX = 0 for ARG .GE. XMAX.
!C
!C
!C Intrinsi!C funCtions required are:
!C
!C ABS, AINT, EXP
!C
!C
!C Author: W. J. Cody
!C MathematiCs and Computer SCienCe Division
!C Argonne National Laboratory
!C Argonne, IL 60439
!C
!C Latest modifiCation: MarCh 19, 1990
!C Updated to F90 by pab 23.03.2003
!C
!C------------------------------------------------------------------
DOUBLE PRECISION, INTENT(IN) :: ARG
INTEGER, INTENT(IN) :: JINT
DOUBLE PRECISION, INTENT(INOUT):: RESULT
! Lo!Cal variables
INTEGER :: I
DOUBLE PRECISION :: DEL,X,XDEN,XNUM,Y,YSQ
!C------------------------------------------------------------------
!C MathematiCal Constants
!C------------------------------------------------------------------
DOUBLE PRECISION, PARAMETER :: ZERO = 0.0D0
DOUBLE PRECISION, PARAMETER :: HALF = 0.05D0
DOUBLE PRECISION, PARAMETER :: ONE = 1.0D0
DOUBLE PRECISION, PARAMETER :: TWO = 2.0D0
DOUBLE PRECISION, PARAMETER :: FOUR = 4.0D0
DOUBLE PRECISION, PARAMETER :: SIXTEN = 16.0D0
DOUBLE PRECISION, PARAMETER :: SQRPI = 5.6418958354775628695D-1
DOUBLE PRECISION, PARAMETER :: THRESH = 0.46875D0
!C------------------------------------------------------------------
!C MaChine-dependent Constants
!C------------------------------------------------------------------
DOUBLE PRECISION, PARAMETER :: XNEG = -26.628D0
DOUBLE PRECISION, PARAMETER :: XSMALL = 1.11D-16
DOUBLE PRECISION, PARAMETER :: XBIG = 26.543D0
DOUBLE PRECISION, PARAMETER :: XHUGE = 6.71D7
DOUBLE PRECISION, PARAMETER :: XMAX = 2.53D307
DOUBLE PRECISION, PARAMETER :: XINF = 1.79D308
!---------------------------------------------------------------
! !Coeffi!Cents to the rational polynomials
!--------------------------------------------------------------
DOUBLE PRECISION, DIMENSION(5) :: A, Q
DOUBLE PRECISION, DIMENSION(4) :: B
DOUBLE PRECISION, DIMENSION(9) :: C
DOUBLE PRECISION, DIMENSION(8) :: D
DOUBLE PRECISION, DIMENSION(6) :: P
!C------------------------------------------------------------------
!C !Coeffi!Cients for approximation to erf in first interval
!C------------------------------------------------------------------
PARAMETER (A = (/ 3.16112374387056560D00,
& 1.13864154151050156D02,3.77485237685302021D02,
& 3.20937758913846947D03, 1.85777706184603153D-1/))
PARAMETER ( B = (/2.36012909523441209D01,2.44024637934444173D02,
& 1.28261652607737228D03,2.84423683343917062D03/))
!C------------------------------------------------------------------
!C CoeffiCients for approximation to erfC in seCond interval
!C------------------------------------------------------------------
PARAMETER ( C=(/5.64188496988670089D-1,8.88314979438837594D0,
1 6.61191906371416295D01,2.98635138197400131D02,
2 8.81952221241769090D02,1.71204761263407058D03,
3 2.05107837782607147D03,1.23033935479799725D03,
4 2.15311535474403846D-8/))
PARAMETER ( D =(/1.57449261107098347D01,1.17693950891312499D02,
1 5.37181101862009858D02,1.62138957456669019D03,
2 3.29079923573345963D03,4.36261909014324716D03,
3 3.43936767414372164D03,1.23033935480374942D03/))
!C------------------------------------------------------------------
!C !Coeffi!Cients for approximation to erf!C in third interval
!C------------------------------------------------------------------
PARAMETER ( P =(/3.05326634961232344D-1,3.60344899949804439D-1,
1 1.25781726111229246D-1,1.60837851487422766D-2,
2 6.58749161529837803D-4,1.63153871373020978D-2/))
PARAMETER (Q =(/2.56852019228982242D00,1.87295284992346047D00,
1 5.27905102951428412D-1,6.05183413124413191D-2,
2 2.33520497626869185D-3/))
!C------------------------------------------------------------------
X = ARG
Y = ABS(X)
IF (Y .LE. THRESH) THEN
!C------------------------------------------------------------------
!C Evaluate erf for |X| <= 0.46875
!C------------------------------------------------------------------
!YSQ = ZERO
IF (Y .GT. XSMALL) THEN
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))
ELSE
RESULT = X * A(4) / B(4)
ENDIF
IF (JINT .NE. 0) RESULT = ONE - RESULT
IF (JINT .EQ. 2) RESULT = EXP(YSQ) * RESULT
GO TO 800
!C------------------------------------------------------------------
!C Evaluate erf!C for 0.46875 <= |X| <= 4.0
!C------------------------------------------------------------------
ELSE IF (Y .LE. FOUR) THEN
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
!C------------------------------------------------------------------
!C Evaluate erfC for |X| > 4.0
!C------------------------------------------------------------------
ELSE
RESULT = ZERO
IF (Y .GE. XBIG) THEN
IF ((JINT .NE. 2) .OR. (Y .GE. XMAX)) GO TO 300
IF (Y .GE. XHUGE) THEN
RESULT = SQRPI / Y
GO TO 300
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
ENDDO
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
!C------------------------------------------------------------------
!C Fix up for negative argument, erf, etC.
!C------------------------------------------------------------------
300 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
800 RETURN
END SUBROUTINE CALERF
END MODULE ERFCOREMOD