1 real function SERF(X) 2c Copyright (c) 1996 California Institute of Technology, Pasadena, CA. 3c ALL RIGHTS RESERVED. 4c Based on Government Sponsored Research NAS7-03001. 5c>> 1999-12-27 SERF Krogh Added SAVE XMAX in SERFCE. 6c>> 1997-06-05 SERF Snyder Added SAVE EPS3 in SERFC 7c>> 1996-04-27 SERF Krogh Changes to use .C. and C%%. 8c>> 1996-03-30 SERF Krogh Added external statements. 9C>> 1995-11-28 SERF Krogh Converted SFTRAN to Fortran 10C>> 1995-10-24 SERF Krogh Removed blanks in numbers for C conversion. 11C>> 1994-10-19 SERF Krogh Changes to use M77CON 12C>> 1991-01-14 SERF CLL Changed to use generic names: abs and sign. 13C>> 1990-11-29 CLL 14c>> 1989-08-16 SERF CLL Changes for Cray 15C>> 1989-06-30 SERF CLL More digits in sqrt(PI), S.P. & D.P. more alike. 16C>> 1986-03-18 SERF Lawson Initial code. 17c--S replaces "?": ?ERF, ?ERFC, ?ERFCE, ?INITS, ?CSEVL, ?ERM1, ?ERF1, 18c--& ?ERFC1, ?ERFE1, ?ERFE2, ?ERFSP 19C 20C JULY 1977 EDITION. W.FULLERTON, C3, 21C LOS ALAMOS SCIENTIFIC LAB. 22C 23C REORGANIZATION OF FULLERTON'S SERF & SERFC 24C C.L.LAWSON & S.CHAN, JUNE 2 ,1983, JPL. 25c When SERFC is called with 0.5 .lt. X .le. 1.0 the original code 26c computed erfc from erfce, computing erfce from Cody's rational 27c formula whose coeffs (PS() and QS()) are only given to 18 decimal 28c digits. For Cray D.P. of 29D precision we can do better by 29c computing SERFC as 1 - erf in this interval, since the coeffs for 30c SERF are good up to 30D. 31c With this change SERFC will also be good up to 30 or 29D for all 32c X up to XMAX, where it would underflow. 33c We compute SERFCE only for X .ge. 0. SERFCE is computed from 34c SERF in [0, 1] and from SERFC in [1, XMAX] and so is good to 35c about 30D there. For X > XMAX we use the Cody rational 36c approximations for SERFCE which have coeffs good only to 18D. 37c Added variable EPS3. 38C ------------------------------------------------------------------ 39C-- Begin mask code changes 40C MACHINE DEPENDENT VALUES ARE SET ON THE FIRST ENTRY 41C TO THIS CODE. EXAMPLES OF THESE VALUES FOLLOW: 42C 43C SYSTEM NTERF NTERFC NTERC2 SQEPS XBIG XMAX 44C ------ ----- ------ ------ ----- ---- ---- 45c IEEE S.P. 7 9 10 .34E-3 4.01 9.18 46c IEEE D.P. 12 25 24 .15E-7 6.01 26.53 47C UNIVAC S.P. 8 11 11 .12E-3 4.26 9.29 48C UNIVAC D.P. 14 30 27 .13E-8 6.40 26.58 49c Cray S.P. 11 21 21 .12E-6 5.66 75.30 50c Cray D.P. 19 55 45 .10E-13 8.04 74.89 51C 52c The above values are derived from the following values: 53C 54C SYSTEM R1MACH(1) R1MACH(3) D1MACH(1) D1MACH(3) 55C ------ --------- --------- --------- --------- 56c IEEE 2**(-126) 2**(-24) 2**(-1022) 2**(-53) 57c VAX 2**(-128) 2**(-27) 2**(-128) 2**(-56) 58c UNIVAC 2**(-129) 2**(-27) 2**(-1025) 2**(-60) 59C Cray 2**(-8190) 2**(-47) 2**(-8100) 2**(-94) 60C 61C-- End mask code changes 62C ------------------------------------------------------------------ 63 external SERF1, SERFC1 64 real X 65 real FAC, U, SQEPS, XBIG, Y 66 real SERF1, SERFC1 67 logical FIRST 68C 69 save FIRST, SQEPS, XBIG 70C FAC = 2/sqrt(PI) 71 parameter(FAC = 1.12837916709551257389615890312154516E0) 72 data FIRST / .TRUE. / 73c 74C ------------------------------------------------------------------ 75C 76 if (FIRST) call SERFSP(FIRST, U, SQEPS, XBIG, Y) 77 Y = abs(X) 78 IF (Y .LE. SQEPS) THEN 79 U = FAC * X 80 ELSE IF (Y .LE. 1.0E0) THEN 81 U = SERF1(X) 82 ELSE IF (Y .LT. XBIG) THEN 83 U = sign( 0.5E0+(0.5E0 - SERFC1(Y) * exp(-Y*Y) / Y), X ) 84 ELSE 85 U = sign ( 1.0E0,X ) 86 END IF 87 SERF = U 88 RETURN 89 END 90C 91C 92c ------------------------------------------------------------------ 93 real function SERFC(X) 94C 95 external SERF1, SERFE1, SERFC1 96 real X 97 real EPS3, FAC, U, SQEPS, XBIG, XMAX, Y 98 real SERF1, SERFC1, SERFE1 99 logical FIRST 100C 101 save EPS3, FIRST, SQEPS, XBIG, XMAX 102C FAC = 2/sqrt(PI) 103 parameter(FAC = 1.12837916709551257389615890312154516E0) 104 data FIRST / .TRUE. / 105c 106 if (FIRST) call SERFSP(FIRST, EPS3, SQEPS, XBIG, XMAX) 107 IF (X .LE. -XBIG) THEN 108 U = 2.0E0 109 ELSE 110 Y = abs(X) 111 IF (Y .LE. SQEPS) THEN 112 U = 0.5E0 + (0.5E0 - FAC*X) 113 ELSE IF (Y .LE. 1.0E0) THEN 114 IF (X .LT. 0.5E0) THEN 115C 116C *** Here X is in [-1,-SQEPS] or [SQEPS,0.5] 117C 118 U = 0.5E0 + (0.5E0-SERF1(X)) 119 ELSE 120C 121C . *** Here X is in [0.5,1.0] 122c . Following test on 1.0E-19 added 8/89 to make better use 123c . of 28 D precision available on a Cray. The stored 124c . approximation for erfce is only good to about 18 D. 125c . On machines having precision greater than 18 D we get 126c . a more accurate value for erfc by computing U = erf 127c . and then 0.5E0 - (0.5E0 - U). 128c 129 if(EPS3 .gt. 1.0E-19) then 130 U = EXP(-X**2) * SERFE1(Y) 131 else 132 U = 0.5E0 + (0.5E0-SERF1(X)) 133 end if 134 END IF 135 ELSE IF (Y .LE. XMAX) THEN 136 U = SERFC1(Y) * exp(-Y*Y) / Y 137 IF (X .LT. 0.0E0) U = 2.0E0 - U 138 ELSE 139 U = 0.0E0 140 CALL SERM1('SERFC',1,0,'X SO BIG SERFC UNDERFLOWS', 141 * 'X',X,'.') 142 END IF 143 END IF 144 SERFC = U 145 RETURN 146 END 147C 148c ------------------------------------------------------------------ 149 real function SERFCE(X) 150C 151 external SERF1, SERFC1, SERFE1, SERFE2 152 real X 153 real CUT2, CUT3, PII, U, XMAX, Y 154 real SERF1, SERFC1, SERFE1, SERFE2 155 logical FIRST 156 save FIRST, XMAX 157C PII = 1/sqrt(PI) 158 parameter(PII = 0.56418958354775628694807945156077258E0) 159C 160c CUT2 = 4.0 161C CUT3 = 10.**16 WILL BE SATISFACTORY FOR ALL 162C MACHINES HAVING RELATIVE PRECISION LESS 163C THAN 32 DECIMAL PLACES AND OVERFLOW LIMIT 164C GREATER THAN 10.**32 165 parameter(CUT2 = 4.0E0, CUT3 = 1.0E16) 166 data FIRST / .TRUE. / 167c 168 if (FIRST) call SERFSP(FIRST, U, Y, Y, XMAX) 169 Y = X 170 IF (X .LT.0.0E0) THEN 171 CALL SERM1('SERFCE',1,0,'X MUST BE .GE. ZERO', 172 * 'X',X,'.') 173 U = 0.0E0 174 ELSE IF (X .LT. 1.0E0) THEN 175 U = EXP(Y*Y)*(0.5E0 + (0.5E0-SERF1(X))) 176 ELSE IF (X .LT. XMAX) THEN 177 U = SERFC1(Y) / X 178 ELSE IF ( X .lt. CUT2) then 179 U = SERFE1(Y) 180 ELSE IF (X .LT. CUT3) THEN 181 U = SERFE2(Y) 182 ELSE 183 U = PII / X 184 END IF 185 SERFCE = U 186 RETURN 187 END 188C 189C 190 subroutine SERFSP(FIRST, EPS3, SQEPS, XBIG, XMAX) 191c PROCEDURE (SET PARAMETERS) 192 external R1MACH 193 real EPS3, SQEPS, XBIG, XMAX 194 logical FIRST 195 real EPS3I,R1MACH,SQEPSI,SQRTPI,XBIGI,XMAXI 196 logical FIRSTI 197 parameter(SQRTPI = 1.77245385090551602729816748334114518E0) 198 save EPS3I, FIRSTI, SQEPSI, XBIGI, XMAXI 199 data FIRSTI / .true. / 200c 201 if (FIRSTI) then 202 FIRSTI = .false. 203C 204 EPS3I = R1MACH(3) 205 SQEPSI = sqrt (2.0E0*EPS3I) 206 XBIGI = sqrt (-log(SQRTPI*EPS3I)) 207 XMAXI = sqrt (-log(SQRTPI*R1MACH(1)) ) 208 XMAXI = XMAXI - 0.5E0*log(XMAXI)/XMAXI - 0.01E0 209 end if 210 EPS3 = EPS3I 211 SQEPS = SQEPSI 212 XBIG = XBIGI 213 XMAX = XMAXI 214c PRINT*,'SQEPS = ',SQEPS, SQEPSI 215c PRINT*,'XBIG = ',XBIG, XBIGI 216c PRINT*,'XMAX = ',XMAX, XMAXI 217 FIRST = .false. 218 return 219 end 220C 221C 222 real function SERF1(X) 223c PROCEDURE (SERF(X) FOR ABS(X) .LE. 1) 224 external R1MACH, SCSEVL 225 real X 226 real R1MACH, SCSEVL, ERFCS(21) 227 integer NTERF 228 save NTERF 229C 230C SERIES FOR ERF ON THE INTERVAL 0. TO 1.00000E+00 231C WITH WEIGHTED ERROR 1.28E-32 232C LOG WEIGHTED ERROR 31.89 233C SIGNIFICANT FIGURES REQUIRED 31.05 234C DECIMAL PLACES REQUIRED 32.55 235C 236c++ Save data by elements if ~.C. 237 DATA ERFCS(1) / -.49046121234691808039984544033376E-1 / 238 DATA ERFCS(2) / -.14226120510371364237824741899631E+0 / 239 DATA ERFCS(3) / +.10035582187599795575754676712933E-1 / 240 DATA ERFCS(4) / -.57687646997674847650827025509167E-3 / 241 DATA ERFCS(5) / +.27419931252196061034422160791471E-4 / 242 DATA ERFCS(6) / -.11043175507344507604135381295905E-5 / 243 DATA ERFCS(7) / +.38488755420345036949961311498174E-7 / 244 DATA ERFCS(8) / -.11808582533875466969631751801581E-8 / 245 DATA ERFCS(9) / +.32334215826050909646402930953354E-10 / 246 DATA ERFCS(10) / -.79910159470045487581607374708595E-12 / 247 DATA ERFCS(11) / +.17990725113961455611967245486634E-13 / 248 DATA ERFCS(12) / -.37186354878186926382316828209493E-15 / 249 DATA ERFCS(13) / +.71035990037142529711689908394666E-17 / 250 DATA ERFCS(14) / -.12612455119155225832495424853333E-18 / 251 DATA ERFCS(15) / +.20916406941769294369170500266666E-20 / 252 DATA ERFCS(16) / -.32539731029314072982364160000000E-22 / 253 DATA ERFCS(17) / +.47668672097976748332373333333333E-24 / 254 DATA ERFCS(18) / -.65980120782851343155199999999999E-26 / 255 DATA ERFCS(19) / +.86550114699637626197333333333333E-28 / 256 DATA ERFCS(20) / -.10788925177498064213333333333333E-29 / 257 DATA ERFCS(21) / +.12811883993017002666666666666666E-31 / 258c 259 DATA NTERF / 0 / 260c 261 if (NTERF .eq. 0) call SINITS(ERFCS, 21, 0.1E0*R1MACH(3), NTERF) 262c 263 SERF1 = X + X * SCSEVL ( 2.0E0 * X * X - 1.0E0, ERFCS, NTERF ) 264 RETURN 265 END 266C 267C 268 real function SERFC1(Y) 269C *(exp(-Y*Y) / Y) = PROCEDURE (SERFC(Y) FOR 1 .LT. Y .LT. XMAX) 270C *(1 / X) = PROCEDURE (SERFCE(X) FOR 1 .LT. X .LT. XMAX) 271 external R1MACH, SCSEVL 272 real Y 273 integer NTERC2, NTERFC 274 real R1MACH, SCSEVL, ERFCCS(59), ERC2CS(49), ETA, YSQ 275 save NTERC2, NTERFC 276C 277C 278C SERIES FOR ERC2 ON THE INTERVAL 2.50000E-01 TO 1.00000E+00 279C WITH WEIGHTED ERROR 2.67E-32 280C LOG WEIGHTED ERROR 31.57 281C SIGNIFICANT FIGURES REQUIRED 30.31 282C DECIMAL PLACES REQUIRED 32.42 283C 284c++ Save data by elements if ~.C. 285 DATA ERC2CS(1) / -.6960134660230950112739150826197E-1 / 286 DATA ERC2CS(2) / -.4110133936262089348982212084666E-1 / 287 DATA ERC2CS(3) / +.3914495866689626881561143705244E-2 / 288 DATA ERC2CS(4) / -.4906395650548979161280935450774E-3 / 289 DATA ERC2CS(5) / +.7157479001377036380760894141825E-4 / 290 DATA ERC2CS(6) / -.1153071634131232833808232847912E-4 / 291 DATA ERC2CS(7) / +.1994670590201997635052314867709E-5 / 292 DATA ERC2CS(8) / -.3642666471599222873936118430711E-6 / 293 DATA ERC2CS(9) / +.6944372610005012589931277214633E-7 / 294 DATA ERC2CS(10) / -.1371220902104366019534605141210E-7 / 295 DATA ERC2CS(11) / +.2788389661007137131963860348087E-8 / 296 DATA ERC2CS(12) / -.5814164724331161551864791050316E-9 / 297 DATA ERC2CS(13) / +.1238920491752753181180168817950E-9 / 298 DATA ERC2CS(14) / -.2690639145306743432390424937889E-10 / 299 DATA ERC2CS(15) / +.5942614350847910982444709683840E-11 / 300 DATA ERC2CS(16) / -.1332386735758119579287754420570E-11 / 301 DATA ERC2CS(17) / +.3028046806177132017173697243304E-12 / 302 DATA ERC2CS(18) / -.6966648814941032588795867588954E-13 / 303 DATA ERC2CS(19) / +.1620854541053922969812893227628E-13 / 304 DATA ERC2CS(20) / -.3809934465250491999876913057729E-14 / 305 DATA ERC2CS(21) / +.9040487815978831149368971012975E-15 / 306 DATA ERC2CS(22) / -.2164006195089607347809812047003E-15 / 307 DATA ERC2CS(23) / +.5222102233995854984607980244172E-16 / 308 DATA ERC2CS(24) / -.1269729602364555336372415527780E-16 / 309 DATA ERC2CS(25) / +.3109145504276197583836227412951E-17 / 310 DATA ERC2CS(26) / -.7663762920320385524009566714811E-18 / 311 DATA ERC2CS(27) / +.1900819251362745202536929733290E-18 / 312 DATA ERC2CS(28) / -.4742207279069039545225655999965E-19 / 313 DATA ERC2CS(29) / +.1189649200076528382880683078451E-19 / 314 DATA ERC2CS(30) / -.3000035590325780256845271313066E-20 / 315 DATA ERC2CS(31) / +.7602993453043246173019385277098E-21 / 316 DATA ERC2CS(32) / -.1935909447606872881569811049130E-21 / 317 DATA ERC2CS(33) / +.4951399124773337881000042386773E-22 / 318 DATA ERC2CS(34) / -.1271807481336371879608621989888E-22 / 319 DATA ERC2CS(35) / +.3280049600469513043315841652053E-23 / 320 DATA ERC2CS(36) / -.8492320176822896568924792422399E-24 / 321 DATA ERC2CS(37) / +.2206917892807560223519879987199E-24 / 322 DATA ERC2CS(38) / -.5755617245696528498312819507199E-25 / 323 DATA ERC2CS(39) / +.1506191533639234250354144051199E-25 / 324 DATA ERC2CS(40) / -.3954502959018796953104285695999E-26 / 325 DATA ERC2CS(41) / +.1041529704151500979984645051733E-26 / 326 DATA ERC2CS(42) / -.2751487795278765079450178901333E-27 / 327 DATA ERC2CS(43) / +.7290058205497557408997703680000E-28 / 328 DATA ERC2CS(44) / -.1936939645915947804077501098666E-28 / 329 DATA ERC2CS(45) / +.5160357112051487298370054826666E-29 / 330 DATA ERC2CS(46) / -.1378419322193094099389644800000E-29 / 331 DATA ERC2CS(47) / +.3691326793107069042251093333333E-30 / 332 DATA ERC2CS(48) / -.9909389590624365420653226666666E-31 / 333 DATA ERC2CS(49) / +.2666491705195388413323946666666E-31 / 334C 335C SERIES FOR ERFC ON THE INTERVAL 0. TO 2.50000E-01 336C WITH WEIGHTED ERROR 1.53E-31 337C LOG WEIGHTED ERROR 30.82 338C SIGNIFICANT FIGURES REQUIRED 29.47 339C DECIMAL PLACES REQUIRED 31.70 340C 341c++ Save data by elements if ~.C. 342 DATA ERFCCS(1) / +.715179310202924774503697709496E-1 / 343 DATA ERFCCS(2) / -.265324343376067157558893386681E-1 / 344 DATA ERFCCS(3) / +.171115397792085588332699194606E-2 / 345 DATA ERFCCS(4) / -.163751663458517884163746404749E-3 / 346 DATA ERFCCS(5) / +.198712935005520364995974806758E-4 / 347 DATA ERFCCS(6) / -.284371241276655508750175183152E-5 / 348 DATA ERFCCS(7) / +.460616130896313036969379968464E-6 / 349 DATA ERFCCS(8) / -.822775302587920842057766536366E-7 / 350 DATA ERFCCS(9) / +.159214187277090112989358340826E-7 / 351 DATA ERFCCS(10) / -.329507136225284321486631665072E-8 / 352 DATA ERFCCS(11) / +.722343976040055546581261153890E-9 / 353 DATA ERFCCS(12) / -.166485581339872959344695966886E-9 / 354 DATA ERFCCS(13) / +.401039258823766482077671768814E-10 / 355 DATA ERFCCS(14) / -.100481621442573113272170176283E-10 / 356 DATA ERFCCS(15) / +.260827591330033380859341009439E-11 / 357 DATA ERFCCS(16) / -.699111056040402486557697812476E-12 / 358 DATA ERFCCS(17) / +.192949233326170708624205749803E-12 / 359 DATA ERFCCS(18) / -.547013118875433106490125085271E-13 / 360 DATA ERFCCS(19) / +.158966330976269744839084032762E-13 / 361 DATA ERFCCS(20) / -.472689398019755483920369584290E-14 / 362 DATA ERFCCS(21) / +.143587337678498478672873997840E-14 / 363 DATA ERFCCS(22) / -.444951056181735839417250062829E-15 / 364 DATA ERFCCS(23) / +.140481088476823343737305537466E-15 / 365 DATA ERFCCS(24) / -.451381838776421089625963281623E-16 / 366 DATA ERFCCS(25) / +.147452154104513307787018713262E-16 / 367 DATA ERFCCS(26) / -.489262140694577615436841552532E-17 / 368 DATA ERFCCS(27) / +.164761214141064673895301522827E-17 / 369 DATA ERFCCS(28) / -.562681717632940809299928521323E-18 / 370 DATA ERFCCS(29) / +.194744338223207851429197867821E-18 / 371 DATA ERFCCS(30) / -.682630564294842072956664144723E-19 / 372 DATA ERFCCS(31) / +.242198888729864924018301125438E-19 / 373 DATA ERFCCS(32) / -.869341413350307042563800861857E-20 / 374 DATA ERFCCS(33) / +.315518034622808557122363401262E-20 / 375 DATA ERFCCS(34) / -.115737232404960874261239486742E-20 / 376 DATA ERFCCS(35) / +.428894716160565394623737097442E-21 / 377 DATA ERFCCS(36) / -.160503074205761685005737770964E-21 / 378 DATA ERFCCS(37) / +.606329875745380264495069923027E-22 / 379 DATA ERFCCS(38) / -.231140425169795849098840801367E-22 / 380 DATA ERFCCS(39) / +.888877854066188552554702955697E-23 / 381 DATA ERFCCS(40) / -.344726057665137652230718495566E-23 / 382 DATA ERFCCS(41) / +.134786546020696506827582774181E-23 / 383 DATA ERFCCS(42) / -.531179407112502173645873201807E-24 / 384 DATA ERFCCS(43) / +.210934105861978316828954734537E-24 / 385 DATA ERFCCS(44) / -.843836558792378911598133256738E-25 / 386 DATA ERFCCS(45) / +.339998252494520890627359576337E-25 / 387 DATA ERFCCS(46) / -.137945238807324209002238377110E-25 / 388 DATA ERFCCS(47) / +.563449031183325261513392634811E-26 / 389 DATA ERFCCS(48) / -.231649043447706544823427752700E-26 / 390 DATA ERFCCS(49) / +.958446284460181015263158381226E-27 / 391 DATA ERFCCS(50) / -.399072288033010972624224850193E-27 / 392 DATA ERFCCS(51) / +.167212922594447736017228709669E-27 / 393 DATA ERFCCS(52) / -.704599152276601385638803782587E-28 / 394 DATA ERFCCS(53) / +.297976840286420635412357989444E-28 / 395 DATA ERFCCS(54) / -.126252246646061929722422632994E-28 / 396 DATA ERFCCS(55) / +.539543870454248793985299653154E-29 / 397 DATA ERFCCS(56) / -.238099288253145918675346190062E-29 / 398 DATA ERFCCS(57) / +.109905283010276157359726683750E-29 / 399 DATA ERFCCS(58) / -.486771374164496572732518677435E-30 / 400 DATA ERFCCS(59) / +.152587726411035756763200828211E-30 / 401c 402 DATA NTERC2 / 0 / 403c 404 if (NTERC2 .eq. 0) then 405 ETA = 0.1E0 * R1MACH(3) 406 call SINITS (ERC2CS, 49, ETA, NTERC2) 407 call SINITS (ERFCCS, 59, ETA, NTERFC) 408 end if 409 YSQ = Y * Y 410 IF (YSQ .LE. 4.0E0) THEN 411 SERFC1 = 0.5E0+SCSEVL((8.0E0/YSQ-5.0E0)/3.0E0,ERC2CS,NTERC2) 412 ELSE 413 SERFC1 = 0.5E0+SCSEVL(8.0E0/YSQ-1.0E0,ERFCCS,NTERFC) 414 END IF 415 RETURN 416 END 417C 418C 419 real function SERFE1(Y) 420c PROCEDURE (SERFCE(Y) FOR 15/32 .LE. Y .LE. 4) 421 real Y 422 real PS(9), PSUM, QS(10), QSUM 423 DATA PS/ 424 * 1.00000000000036828E+0, 425 * 1.87051017604560834E+0, 426 * 1.74642369370058320E+0, 427 * 1.02438464807598001E+0, 428 * 4.07413180167223764E-1, 429 * 1.11870870991098165E-1, 430 * 2.07045775788719818E-2, 431 * 2.37133372752999036E-3, 432 * 1.29992515945788642E-4/ 433C 434 DATA QS/ 435 * 1.00000000000000000E+0, 436 * 2.99888934314798253E+0, 437 * 4.13030795287321183E+0, 438 * 3.43830153103630866E+0, 439 * 1.91273588328781533E+0, 440 * 7.40352738163508723E-1, 441 * 2.00387662412610424E-1, 442 * 3.68131014202168126E-2, 443 * 4.20307996290648223E-3, 444 * 2.30405728794132537E-4/ 445C 446 PSUM=PS(1)+Y*(PS(2)+Y*(PS(3)+Y*(PS(4)+Y*(PS(5)+Y*(PS(6) 447 * +Y*(PS(7)+Y*(PS(8)+Y*PS(9)))))))) 448 QSUM=QS(1)+Y*(QS(2)+Y*(QS(3)+Y*(QS(4)+Y*(QS(5)+Y*(QS(6) 449 * +Y*(QS(7)+Y*(QS(8)+Y*(QS(9)+Y*QS(10))))))))) 450 SERFE1=PSUM/QSUM 451 RETURN 452 END 453C 454C 455 real function SERFE2(Y) 456C PROCEDURE (SERFCE(Y) FOR 4 .LE. Y .LE. 10**16) 457 real Y 458 real PL(6), PII, PSUM, QL(6), QSUM, YSQ 459C PII = 1/sqrt(PI) 460 parameter(PII = 0.56418958354775628694807945156077258E0) 461 DATA PL/ 462 * -2.82094791773876911E-1, 463 * -6.88752606824337911E+0, 464 * -5.38632485753818598E+1, 465 * -1.54309751654255924E+2, 466 * -1.30749393763753716E+2, 467 * -6.98670450907125305E+0/ 468C 469 DATA QL/ 470 * 1.00000000000000000E+0, 471 * 2.59156442057350244E+1, 472 * 2.26063711030173368E+2, 473 * 8.02050727433854173E+2, 474 * 1.09991209268230208E+3, 475 * 4.28227932949102556E+2/ 476C 477 YSQ=1.E0/Y/Y 478 PSUM=PL(1)+YSQ*(PL(2)+YSQ*(PL(3)+YSQ*(PL(4)+YSQ*(PL(5) 479 * +YSQ*PL(6))))) 480 QSUM=QL(1)+YSQ*(QL(2)+YSQ*(QL(3)+YSQ*(QL(4)+YSQ*(QL(5) 481 * +YSQ*QL(6))))) 482 SERFE2=(PII+YSQ*PSUM/QSUM)/Y 483 RETURN 484 END 485