!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