program erfa
c
c -- gaussian error
c -- figure 11.7
c
      integer out
      real x, e1, e2
c
      out = 6
10    write(out, 101)
      read(*, 102) x
      if (x .lt. 0.0) goto 99
      if (x .ge. 1.5) goto 20
      e1 = erf(x)
      e2 = 1.0 - e1
      goto 30
20    e2 = erfc(x)
      e1 = 1.0 - e2
30    write(out, 104) e1, e2
      goto 10
99    stop
101   format(/' arg? ' )
102   format(e10.0)
104   format(' erf =', 0pf10.5, ', erfc =', 1pe12.4)
      end
      function erf(x)
c
c -- gaussian error function by infinite series
c
      integer i
      real x, x2, sum, sum1, term
      data tol/1.0e-5/, sqrtpi/ 1.772454/
c
      erf = 0.0
      if (x .eq. 0.0) goto 99
      erf = 1.0
      if (x .gt. 4.0) goto 99
      x2 = x * x
      sum = x
      term = x
      i = 0
10    i = i + 1
       sum1 = sum
       term = term * x2 / (i + 0.5)
       sum = term + sum1
       if (term .ge. tol*sum) goto 10
      erf = 2 * sum * exp(-x2) / sqrtpi
99    return
101   format(1x, i7, f9.5)
      end
      function erfc(x)
c
c -- complement of gaussian error function
c -- by asymptotic expansion
c
      integer i, j, terms
      real x, x2, sum, u, v, sqrtpi
      data sqrtpi/ 1.772454/, terms/12/
c
      x2 = x * x
      v = 0.5 / x2
      u = 1.0 + v * (terms + 1)
      do 10 j = 1, terms
       i = terms - j + 1
       sum = 1.0 + i * v / u
       u = sum
10    continue
      erfc = exp(-x2) / (x * sum * sqrtpi)
      return
      end