1 : C CALERF SOURCE CHAT 05/01/12 21:45:49 5004 2 : SUBROUTINE CALERF(ARG,RESULT,JINT) 3 : C------------------------------------------------------------------ 4 : C 5 : C This packet evaluates erf(x), erfc(x), and exp(x*x)*erfc(x) 6 : C for a real argument x. It contains three FUNCTION type 7 : C subprograms: ERF, ERFC, and ERFCX (or DERF, DERFC, and DERFCX), 7 : C and one SUBROUTINE type subprogram, CALERF. The calling 8 : C statements for the primary entries are: 9 : C 10 : C Y=ERF(X) (or Y=DERF(X)), 11 : C 12 : C Y=ERFC(X) (or Y=DERFC(X)), 13 : C and 14 : C Y=ERFCX(X) (or Y=DERFCX(X)). 15 : C 16 : C The routine CALERF is intended for internal packet use only, 17 : C all computations within the packet being concentrated in this 18 : C routine. The function subprograms invoke CALERF with the 19 : C statement 20 : C 20 : C CALL CALERF(ARG,RESULT,JINT) 21 : C 22 : C where the parameter usage is as follows 23 : C 24 : C Function Parameters for CALERF 24 : C call ARG Result JINT 25 : C 26 : C ERF(ARG) ANY REAL ARGUMENT ERF(ARG) 0 27 : C ERFC(ARG) ABS(ARG) .LT. XBIG ERFC(ARG) 1 28 : C ERFCX(ARG) XNEG .LT. ARG .LT. XMAX ERFCX(ARG) 2 29 : C 30 : C The main computation evaluates near-minimax approximations 31 : C from "Rational Chebyshev approximations for the error function" 32 : C by W. J. Cody, Math. Comp., 1969, PP. 631-638. This 32 : C transportable program uses rational functions that theoretically 33 : C approximate erf(x) and erfc(x) to at least 18 significant 34 : C decimal digits. The accuracy achieved depends on the arithmetic 35 : C system, the compiler, the intrinsic functions, and proper 36 : C selection of the machine-dependent constants. 37 : C 38 : C******************************************************************* 39 : C******************************************************************* 40 : C 41 : C Explanation of machine-dependent constants 42 : C 43 : C XMIN = the smallest positive floating-point number. 44 : C XINF = the largest positive finite floating-point number. 45 : C XNEG = the largest negative argument acceptable to ERFCX; 46 : C the negative of the solution to the equation 47 : C 2*exp(x*x) = XINF. 48 : C XSMALL = argument below which erf(x) may be represented by 49 : C 2*x/sqrt(pi) and above which x*x will not underflow. 50 : C A conservative value is the largest machine number X 51 : C such that 1.0 + X = 1.0 to machine precision. 52 : C XBIG = largest argument acceptable to ERFC; solution to 53 : C the equation: W(x) * (1-0.5/x**2) = XMIN, where 54 : C W(x) = exp(-x*x)/[x*sqrt(pi)]. 55 : C XHUGE = argument above which 1.0 - 1/(2*x*x) = 1.0 to 56 : C machine precision. A conservative value is 57 : C 1/[2*sqrt(XSMALL)] 58 : C XMAX = largest acceptable argument to ERFCX; the minimum 59 : C of XINF and 1/[sqrt(pi)*XMIN]. 60 : C 61 : C Approximate values for some important machines are: 62 : C 63 : C XMIN XINF XNEG XSMALL 64 : C 65 : C CDC 7600 (S.P.) 3.13E-294 1.26E+322 -27.220 7.11E-15 66 : C CRAY-1 (S.P.) 4.58E-2467 5.45E+2465 -75.345 7.11E-15 67 : C IEEE (IBM/XT, 68 : C SUN, etc.) (S.P.) 1.18E-38 3.40E+38 -9.382 5.96E-8 69 : C IEEE (IBM/XT, 70 : C SUN, etc.) (D.P.) 2.23D-308 1.79D+308 -26.628 1.11D-16 71 : C IBM 195 (D.P.) 5.40D-79 7.23E+75 -13.190 1.39D-17 72 : C UNIVAC 1108 (D.P.) 2.78D-309 8.98D+307 -26.615 1.73D-18 73 : C VAX D-Format (D.P.) 2.94D-39 1.70D+38 -9.345 1.39D-17 74 : C VAX G-Format (D.P.) 5.56D-309 8.98D+307 -26.615 1.11D-16 75 : C 76 : C 77 : C XBIG XHUGE XMAX 78 : C 79 : C CDC 7600 (S.P.) 25.922 8.39E+6 1.80X+293 80 : C CRAY-1 (S.P.) 75.326 8.39E+6 5.45E+2465 81 : C IEEE (IBM/XT, 82 : C SUN, etc.) (S.P.) 9.194 2.90E+3 4.79E+37 83 : C IEEE (IBM/XT, 84 : C SUN, etc.) (D.P.) 26.543 6.71D+7 2.53D+307 85 : C IBM 195 (D.P.) 13.306 1.90D+8 7.23E+75 86 : C UNIVAC 1108 (D.P.) 26.582 5.37D+8 8.98D+307 87 : C VAX D-Format (D.P.) 9.269 1.90D+8 1.70D+38 88 : C VAX G-Format (D.P.) 26.569 6.71D+7 8.98D+307 89 : C 90 : C******************************************************************* 91 : C******************************************************************* 92 : C 93 : C Error returns 94 : C 95 : C The program returns ERFC = 0 for ARG .GE. XBIG; 96 : C 97 : C ERFCX = XINF for ARG .LT. XNEG; 98 : C and 99 : C ERFCX = 0 for ARG .GE. XMAX. 100 : C 101 : C 102 : C Intrinsic functions required are: 103 : C 104 : C ABS, AINT, EXP 105 : C 106 : C 107 : C Author: W. J. Cody 108 : C Mathematics and Computer Science Division 109 : C Argonne National Laboratory 110 : C Argonne, IL 60439 111 : C 112 : C Latest modification: March 19, 1990 113 : C 114 : C------------------------------------------------------------------ 115 : IMPLICIT INTEGER(I-N) 116 : IMPLICIT REAL*8(A-H,O-Z) 117 : -INC CCREEL 118 : INTEGER I,JINT 119 : REAL*8 120 : 1 A,ARG,B,C,D,DEL,FOUR,HALF,P,ONE,Q,RESULT,SIXTEN,SQRPI, 121 : 2 TWO,THRESH,X,XBIG,XDEN,XHUGE,XINF,XMAX,XNEG,XNUM,XSMALL, 122 : 3 Y,YSQ,ZERO 123 : DIMENSION A(5),B(4),C(9),D(8),P(6),Q(5) 124 : C------------------------------------------------------------------ 125 : C Mathematical constants 126 : C------------------------------------------------------------------ 127 : DATA FOUR,ONE,HALF,TWO,ZERO/4.0D0,1.0D0,0.5D0,2.0D0,0.0D0/, 128 : 1 SQRPI/5.6418958354775628695D-1/,THRESH/0.46875D0/, 129 : 2 SIXTEN/16.0D0/ 130 : C------------------------------------------------------------------ 131 : C Machine-dependent constants 132 : C------------------------------------------------------------------ 133 : C DATA XINF,XNEG,XSMALL/1.79D308,-26.628D0,1.11D-16/, 134 : C 1 XBIG,XHUGE,XMAX/26.543D0,6.71D7,2.53D307/ 135 : XINF=XGRAND 136 : XNEG=-26d0 137 : XSMALL=xzprec 138 : XBIG =26d0 139 : XHUGE=1d0/(2d0*sqrt(XSMALL)) 140 : XMAX=min(xinf,(1d0/(sqrt(Xpi)*Xpetit))) 141 : C------------------------------------------------------------------ 142 : C Coefficients for approximation to erf in first interval 143 : C------------------------------------------------------------------ 144 : DATA A/3.16112374387056560D00,1.13864154151050156D02, 145 : 1 3.77485237685302021D02,3.20937758913846947D03, 146 : 2 1.85777706184603153D-1/ 147 : DATA B/2.36012909523441209D01,2.44024637934444173D02, 148 : 1 1.28261652607737228D03,2.84423683343917062D03/ 149 : C------------------------------------------------------------------ 150 : C Coefficients for approximation to erfc in second interval 151 : C------------------------------------------------------------------ 152 : DATA C/5.64188496988670089D-1,8.88314979438837594D0, 153 : 1 6.61191906371416295D01,2.98635138197400131D02, 154 : 2 8.81952221241769090D02,1.71204761263407058D03, 155 : 3 2.05107837782607147D03,1.23033935479799725D03, 156 : 4 2.15311535474403846D-8/ 157 : DATA D/1.57449261107098347D01,1.17693950891312499D02, 158 : 1 5.37181101862009858D02,1.62138957456669019D03, 159 : 2 3.29079923573345963D03,4.36261909014324716D03, 160 : 3 3.43936767414372164D03,1.23033935480374942D03/ 161 : C------------------------------------------------------------------ 162 : C Coefficients for approximation to erfc in third interval 163 : C------------------------------------------------------------------ 164 : DATA P/3.05326634961232344D-1,3.60344899949804439D-1, 165 : 1 1.25781726111229246D-1,1.60837851487422766D-2, 166 : 2 6.58749161529837803D-4,1.63153871373020978D-2/ 167 : DATA Q/2.56852019228982242D00,1.87295284992346047D00, 168 : 1 5.27905102951428412D-1,6.05183413124413191D-2, 169 : 2 2.33520497626869185D-3/ 170 : C------------------------------------------------------------------ 171 : X = ARG 172 : Y = ABS(X) 173 : IF (Y .LE. THRESH) THEN 174 : C------------------------------------------------------------------ 175 : C Evaluate erf for |X| <= 0.46875 176 : C------------------------------------------------------------------ 177 : YSQ = ZERO 178 : IF (Y .GT. XSMALL) YSQ = Y * Y 179 : XNUM = A(5)*YSQ 180 : XDEN = YSQ 181 : DO 20 I = 1, 3 182 : XNUM = (XNUM + A(I)) * YSQ 183 : XDEN = (XDEN + B(I)) * YSQ 184 : 20 CONTINUE 185 : RESULT = X * (XNUM + A(4)) / (XDEN + B(4)) 186 : IF (JINT .NE. 0) RESULT = ONE - RESULT 187 : IF (JINT .EQ. 2) RESULT = EXP(YSQ) * RESULT 188 : GO TO 800 189 : C------------------------------------------------------------------ 190 : C Evaluate erfc for 0.46875 <= |X| <= 4.0 191 : C------------------------------------------------------------------ 192 : ELSE IF (Y .LE. FOUR) THEN 193 : XNUM = C(9)*Y 194 : XDEN = Y 195 : DO 120 I = 1, 7 196 : XNUM = (XNUM + C(I)) * Y 197 : XDEN = (XDEN + D(I)) * Y 198 : 120 CONTINUE 199 : RESULT = (XNUM + C(8)) / (XDEN + D(8)) 200 : IF (JINT .NE. 2) THEN 201 : YSQ = AINT(Y*SIXTEN)/SIXTEN 202 : DEL = (Y-YSQ)*(Y+YSQ) 203 : RESULT = EXP(-YSQ*YSQ) * EXP(-DEL) * RESULT 204 : END IF 205 : C------------------------------------------------------------------ 206 : C Evaluate erfc for |X| > 4.0 207 : C------------------------------------------------------------------ 208 : ELSE 209 : RESULT = ZERO 210 : IF (Y .GE. XBIG) THEN 211 : IF ((JINT .NE. 2) .OR. (Y .GE. XMAX)) GO TO 300 212 : IF (Y .GE. XHUGE) THEN 213 : RESULT = SQRPI / Y 214 : GO TO 300 215 : END IF 216 : END IF 217 : YSQ = ONE / (Y * Y) 218 : XNUM = P(6)*YSQ 219 : XDEN = YSQ 220 : DO 240 I = 1, 4 221 : XNUM = (XNUM + P(I)) * YSQ 222 : XDEN = (XDEN + Q(I)) * YSQ 223 : 240 CONTINUE 224 : RESULT = YSQ *(XNUM + P(5)) / (XDEN + Q(5)) 225 : RESULT = (SQRPI - RESULT) / Y 226 : IF (JINT .NE. 2) THEN 227 : YSQ = AINT(Y*SIXTEN)/SIXTEN 228 : DEL = (Y-YSQ)*(Y+YSQ) 229 : RESULT = EXP(-YSQ*YSQ) * EXP(-DEL) * RESULT 230 : END IF 231 : END IF 232 : C------------------------------------------------------------------ 233 : C Fix up for negative argument, erf, etc. 234 : C------------------------------------------------------------------ 235 : 300 IF (JINT .EQ. 0) THEN 236 : RESULT = (HALF - RESULT) + HALF 237 : IF (X .LT. ZERO) RESULT = -RESULT 238 : ELSE IF (JINT .EQ. 1) THEN 239 : IF (X .LT. ZERO) RESULT = TWO - RESULT 240 : ELSE 241 : IF (X .LT. ZERO) THEN 242 : IF (X .LT. XNEG) THEN 243 : RESULT = XINF 244 : ELSE 245 : YSQ = AINT(X*SIXTEN)/SIXTEN 246 : DEL = (X-YSQ)*(X+YSQ) 247 : Y = EXP(YSQ*YSQ) * EXP(DEL) 248 : RESULT = (Y+Y) - RESULT 249 : END IF 250 : END IF 251 : END IF 252 : 800 RETURN 253 : C---------- Last card of CALERF ---------- 254 : END 255 : 256 : 257 :
© Cast3M 2003 - All rights reserved.
Disclaimer