erfc.in.f

Go to the documentation of this file.

*DECK ERFC
      FUNCTION erfc (X)
C***BEGIN PROLOGUE  ERFC
C***PURPOSE  Compute the complementary error function.
C***LIBRARY   SLATEC (FNLIB)
C***CATEGORY  C8A, L5A1E
C***TYPE      SINGLE PRECISION (ERFC-S, DERFC-D)
C***KEYWORDS  COMPLEMENTARY ERROR FUNCTION, ERFC, FNLIB,
C             SPECIAL FUNCTIONS
C***AUTHOR  Fullerton, W., (LANL)
C***DESCRIPTION
C
C ERFC(X) calculates the single precision complementary error
C function for single precision argument X.
C
C Series for ERF        on the interval  0.          to  1.00000D+00
C                                        with weighted error   7.10E-18
C                                         log weighted error  17.15
C                               significant figures required  16.31
C                                    decimal places required  17.71
C
C Series for ERFC       on the interval  0.          to  2.50000D-01
C                                        with weighted error   4.81E-17
C                                         log weighted error  16.32
C                        approx significant figures required  15.0
C
C
C Series for ERC2       on the interval  2.50000D-01 to  1.00000D+00
C                                        with weighted error   5.22E-17
C                                         log weighted error  16.28
C                        approx significant figures required  15.0
C                                    decimal places required  16.96
C
C***REFERENCES  (NONE)
C***ROUTINES CALLED  CSEVL, INITS, R1MACH, XERMSG
C***REVISION HISTORY  (YYMMDD)
C   770701  DATE WRITTEN
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   890531  REVISION DATE from Version 3.2
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
C   920618  Removed space from variable names.  (RWC, WRB)
C***END PROLOGUE  ERFC
      dimension erfcs(13), erfccs(24), erc2cs(23)
      LOGICAL FIRST
      SAVE erfcs, erc2cs, erfccs, sqrtpi, nterf, nterfc,
     1 nterc2, xsml, xmax, sqeps, first
      DATA erfcs( 1) /   -.0490461212 34691808e0 /
      DATA erfcs( 2) /   -.1422612051 0371364e0 /
      DATA erfcs( 3) /    .0100355821 87599796e0 /
      DATA erfcs( 4) /   -.0005768764 69976748e0 /
      DATA erfcs( 5) /    .0000274199 31252196e0 /
      DATA erfcs( 6) /   -.0000011043 17550734e0 /
      DATA erfcs( 7) /    .0000000384 88755420e0 /
      DATA erfcs( 8) /   -.0000000011 80858253e0 /
      DATA erfcs( 9) /    .0000000000 32334215e0 /
      DATA erfcs(10) /   -.0000000000 00799101e0 /
      DATA erfcs(11) /    .0000000000 00017990e0 /
      DATA erfcs(12) /   -.0000000000 00000371e0 /
      DATA erfcs(13) /    .0000000000 00000007e0 /
      DATA erc2cs( 1) /   -.0696013466 02309501e0 /
      DATA erc2cs( 2) /   -.0411013393 62620893e0 /
      DATA erc2cs( 3) /    .0039144958 66689626e0 /
      DATA erc2cs( 4) /   -.0004906395 65054897e0 /
      DATA erc2cs( 5) /    .0000715747 90013770e0 /
      DATA erc2cs( 6) /   -.0000115307 16341312e0 /
      DATA erc2cs( 7) /    .0000019946 70590201e0 /
      DATA erc2cs( 8) /   -.0000003642 66647159e0 /
      DATA erc2cs( 9) /    .0000000694 43726100e0 /
      DATA erc2cs(10) /   -.0000000137 12209021e0 /
      DATA erc2cs(11) /    .0000000027 88389661e0 /
      DATA erc2cs(12) /   -.0000000005 81416472e0 /
      DATA erc2cs(13) /    .0000000001 23892049e0 /
      DATA erc2cs(14) /   -.0000000000 26906391e0 /
      DATA erc2cs(15) /    .0000000000 05942614e0 /
      DATA erc2cs(16) /   -.0000000000 01332386e0 /
      DATA erc2cs(17) /    .0000000000 00302804e0 /
      DATA erc2cs(18) /   -.0000000000 00069666e0 /
      DATA erc2cs(19) /    .0000000000 00016208e0 /
      DATA erc2cs(20) /   -.0000000000 00003809e0 /
      DATA erc2cs(21) /    .0000000000 00000904e0 /
      DATA erc2cs(22) /   -.0000000000 00000216e0 /
      DATA erc2cs(23) /    .0000000000 00000052e0 /
      DATA erfccs( 1) /   0.0715179310 202925e0 /
      DATA erfccs( 2) /   -.0265324343 37606719e0 /
      DATA erfccs( 3) /    .0017111539 77920853e0 /
      DATA erfccs( 4) /   -.0001637516 63458512e0 /
      DATA erfccs( 5) /    .0000198712 93500549e0 /
      DATA erfccs( 6) /   -.0000028437 12412769e0 /
      DATA erfccs( 7) /    .0000004606 16130901e0 /
      DATA erfccs( 8) /   -.0000000822 77530261e0 /
      DATA erfccs( 9) /    .0000000159 21418724e0 /
      DATA erfccs(10) /   -.0000000032 95071356e0 /
      DATA erfccs(11) /    .0000000007 22343973e0 /
      DATA erfccs(12) /   -.0000000001 66485584e0 /
      DATA erfccs(13) /    .0000000000 40103931e0 /
      DATA erfccs(14) /   -.0000000000 10048164e0 /
      DATA erfccs(15) /    .0000000000 02608272e0 /
      DATA erfccs(16) /   -.0000000000 00699105e0 /
      DATA erfccs(17) /    .0000000000 00192946e0 /
      DATA erfccs(18) /   -.0000000000 00054704e0 /
      DATA erfccs(19) /    .0000000000 00015901e0 /
      DATA erfccs(20) /   -.0000000000 00004729e0 /
      DATA erfccs(21) /    .0000000000 00001432e0 /
      DATA erfccs(22) /   -.0000000000 00000439e0 /
      DATA erfccs(23) /    .0000000000 00000138e0 /
      DATA erfccs(24) /   -.0000000000 00000048e0 /
      DATA sqrtpi /1.772453850 9055160e0/
      DATA first /.true./
C***FIRST EXECUTABLE STATEMENT  ERFC
      IF (first) THEN
         eta = 0.1*r1mach(3)
         nterf = inits(erfcs, 13, eta)
         nterfc = inits(erfccs, 24, eta)
         nterc2 = inits(erc2cs, 23, eta)
C
         xsml = -sqrt(-log(sqrtpi*r1mach(3)))
         txmax = sqrt(-log(sqrtpi*r1mach(1)))
         xmax = txmax - 0.5*log(txmax)/txmax - 0.01
         sqeps = sqrt(2.0*r1mach(3))
      ENDIF
      first = .false.
C
      IF (isnan(x)) THEN
         erfc = x
         RETURN
      ENDIF
C
      IF (x.GT.xsml) go to 20
C
C ERFC(X) = 1.0 - ERF(X) FOR X .LT. XSML
C
      erfc = 2.
      RETURN
C
 20   IF (x.GT.xmax) go to 40
      y = abs(x)
      IF (y.GT.1.0) go to 30
C
C ERFC(X) = 1.0 - ERF(X) FOR -1. .LE. X .LE. 1.
C
      IF (y.LT.sqeps) erfc = 1.0 - 2.0*x/sqrtpi
      IF (y.GE.sqeps) erfc = 1.0 -
     1  x*(1.0 + csevl(2.*x*x-1., erfcs, nterf) )
      RETURN
C
C ERFC(X) = 1.0 - ERF(X) FOR 1. .LT. ABS(X) .LE. XMAX
C
 30   y = y*y
      IF (y.LE.4.) erfc = exp(-y)/abs(x) * (0.5 + csevl((8./y-5.)/3.,
     1  erc2cs, nterc2) )
      IF (y.GT.4.) erfc = exp(-y)/abs(x) * (0.5 + csevl(8./y-1.,
     1  erfccs, nterfc) )
      IF (x.LT.0.) erfc = 2.0 - erfc
      RETURN
C
 40   erfc = 0.
      RETURN
C
      END