function erfc (x)
c july 1977 edition.  w. fullerton, c3, los alamos scientific lab.
      dimension erfcs(13), erfccs(24), erc2cs(23)
      external alog, csevl, exp, inits, r1mach, sqrt
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
      data erf cs( 1) /   -.0490461212 34691808e0 /
      data erf cs( 2) /   -.1422612051 0371364e0 /
      data erf cs( 3) /    .0100355821 87599796e0 /
      data erf cs( 4) /   -.0005768764 69976748e0 /
      data erf cs( 5) /    .0000274199 31252196e0 /
      data erf cs( 6) /   -.0000011043 17550734e0 /
      data erf cs( 7) /    .0000000384 88755420e0 /
      data erf cs( 8) /   -.0000000011 80858253e0 /
      data erf cs( 9) /    .0000000000 32334215e0 /
      data erf cs(10) /   -.0000000000 00799101e0 /
      data erf cs(11) /    .0000000000 00017990e0 /
      data erf cs(12) /   -.0000000000 00000371e0 /
      data erf cs(13) /    .0000000000 00000007e0 /
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
      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 /
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                                    decimal places required  17.01
c
      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 /
c
      data sqrtpi /1.772453850 9055160e0/
      data nterf, nterfc, nterc2, xsml, xmax, sqeps /3*0, 3*0./
c
      if (nterf.ne.0) go to 10
      eta = 0.1*r1mach(3)
      nterf = inits (erfcs, 13, eta)
      nterfc = inits (erfccs, 24, eta)
      nterc2 = inits (erc2cs, 23, eta)
c
      xsml = -sqrt (-alog(sqrtpi*r1mach(3)))
      xmax = sqrt (-alog(sqrtpi*r1mach(1)))
      xmax = xmax - 0.5*alog(xmax)/xmax - 0.01
      sqeps = sqrt (2.0*r1mach(3))
c
 10   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   call seteru (32herfc    x so big erfc underflows, 32, 1, 0)
      erfc = 0.
      return
c
      end