C $ f2py -m erfcoremod -h erfcoremod.pyf erfcoremod.f C f2py erfcoremod.pyf erfcoremod.f -c --fcompiler=gnu95 --compiler=mingw32 -lmsvcr71 C $ f2py --fcompiler=gnu95 --compiler=mingw32 -lmsvcr71 -m erfcoremod -c erfcoremod.f C gfortran -fPIC -c erfcoremod.f C f2py -m erfcoremod -c erfcoremod.o erfcoremod_interface.f MODULE ERFCOREMOD C IMPLICIT NONE C INTERFACE CALERF C MODULE PROCEDURE CALERF C END INTERFACE C INTERFACE DERF C MODULE PROCEDURE DERF C END INTERFACE C INTERFACE DERFC C MODULE PROCEDURE DERFC C END INTERFACE C INTERFACE DERFCX C MODULE PROCEDURE DERFCX c END INTERFACE 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 DERFC subprogram computes approximate values for erfc(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, March 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 packet evaluates erf(x), erfc(x), and exp(x*x)*erfc(x) C for a real argument x. It contains three FUNCTION type C subprograms: ERF, ERFC, and ERFCX (or DERF, DERFC, and DERFCX), 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=ERFC(X) (or Y=DERFC(X)), C and C Y=ERFCX(X) (or Y=DERFCX(X)). C C The routine CALERF is intended for internal packet use only, C all computations within the packet being concentrated in this C routine. The function 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 Function Parameters for CALERF C call ARG Result JINT C C ERF(ARG) ANY REAL ARGUMENT ERF(ARG) 0 C ERFC(ARG) ABS(ARG) .LT. XBIG ERFC(ARG) 1 C ERFCX(ARG) XNEG .LT. ARG .LT. XMAX ERFCX(ARG) 2 C C The main computation evaluates near-minimax approximations C from "Rational Chebyshev approximations for the error function" C by W. J. Cody, Math. Comp., 1969, PP. 631-638. This C transportable program uses rational functions that theoretically C approximate erf(x) and erfc(x) to at least 18 significant C decimal digits. The accuracy achieved depends on the arithmetic C system, the compiler, the intrinsic functions, and proper C selection of the machine-dependent constants. C C******************************************************************* C******************************************************************* C C Explanation of machine-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 acceptable to ERFCX; C the negative of the solution to the equation C 2*exp(x*x) = XINF. C XSMALL = argument below which erf(x) may be represented by C 2*x/sqrt(pi) and above which x*x will not underflow. C A conservative value is the largest machine number X C such that 1.0 + X = 1.0 to machine precision. C XBIG = largest argument acceptable to ERFC; 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 which 1.0 - 1/(2*x*x) = 1.0 to C machine precision. A conservative value is C 1/[2*sqrt(XSMALL)] C XMAX = largest acceptable argument to ERFCX; the minimum C of XINF and 1/[sqrt(pi)*XMIN]. C C Approximate values for some important machines 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, etc.) (S.P.) 1.18E-38 3.40E+38 -9.382 5.96E-8 C IEEE (IBM/XT, C SUN, etc.) (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 UNIVAC 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, etc.) (S.P.) 9.194 2.90E+3 4.79E+37 C IEEE (IBM/XT, C SUN, etc.) (D.P.) 26.543 6.71D+7 2.53D+307 C IBM 195 (D.P.) 13.306 1.90D+8 7.23E+75 C UNIVAC 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 ERFC = 0 for ARG .GE. XBIG; C C ERFCX = XINF for ARG .LT. XNEG; C and C ERFCX = 0 for ARG .GE. XMAX. C C C Intrinsic 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 ! Local 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 !--------------------------------------------------------------- ! Coefficents 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 Coefficients 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 Coefficients for approximation to erfc 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 erfc 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