      implicit real*8 (a-h,o-z)
c
c This program computes P values for tau and z tests using
c approximations from the paper:
c
c James G. MacKinnon, ``Approximate asymptotic distribution functions
c for unit root and cointegration tests,'' JBES, 12, April 1994.
c
  100 write(6,*) 'Enter 1 for tau test or 2 for z test: '
      read(5,*) itest
      if(itest.lt.1.or.itest.gt.2) then
        write(6,*) 'Invalid response. Try again.'
        go to 100
      end if
  101 write(6,*) 'Enter 1, 2, 3 or 4 for nc, c, ct, or ctt test: '
      read(5,*) it
      if(it.lt.1.or.it.gt.4) then
        write(6,*) 'Invalid response. Try again.'
        go to 101
      end if
  102 write(6,*) 'Enter number of variables (1 through 6): '
      read(5,*) k
      if(k.lt.1.or.k.gt.6) then
        write(6,*) 'Invalid response. Try again.'
        go to 102
      end if
      write(6,*) 'Enter value of test statistic: '
      read(5,*) tstat
c
      if (itest.eq.1) call pvtau(it,k,tstat,pval)
      if (itest.eq.2) call pvz(it,k,tstat,pval)
c
      if (itest.eq.1.and.it.eq.1) write(6,90) k, tstat, pval
   90 format(' tau_nc(',i1,') test statistic = ', f9.3,   ' Approximate
     1asymptotic P value = ', f7.4)
      if (itest.eq.2.and.it.eq.1) write(6,91) k, tstat, pval
   91 format(' z_nc(',i1,') test statistic = ',   f9.3,   ' Approximate
     1asymptotic P value = ', f7.4)
      if (itest.eq.1.and.it.eq.2) write(6,92) k, tstat, pval
   92 format(' tau_c(',i1,') test statistic = ', f9.3,    ' Approximate
     1asymptotic P value = ', f7.4)
      if (itest.eq.2.and.it.eq.2) write(6,93) k, tstat, pval
   93 format(' z_c(',i1,') test statistic = ',  f9.3,     ' Approximate
     1asymptotic P value = ', f7.4)
      if (itest.eq.1.and.it.eq.3) write(6,94) k, tstat, pval
   94 format(' tau_ct(',i1,') test statistic = ', f9.3,   ' Approximate
     1asymptotic P value = ', f7.4)
      if (itest.eq.2.and.it.eq.3) write(6,95) k, tstat, pval
   95 format(' z_ct(',i1,') test statistic = ',  f9.3,    ' Approximate
     1asymptotic P value = ', f7.4)
      if (itest.eq.1.and.it.eq.4) write(6,96) k, tstat, pval
   96 format(' tau_ctt(',i1,') test statistic = ', f9.3,  ' Approximate
     1asymptotic P value = ', f7.4)
      if (itest.eq.2.and.it.eq.4) write(6,97) k, tstat, pval
   97 format(' z_ctt(',i1,') test statistic = ',  f9.3,   ' Approximate
     1asymptotic P value = ', f7.4)
      stop
      end
      subroutine pvtau(it,k,tstat,pval)
c
c This subroutine computes P values for tau tests using approximations
c from the paper:
c
c James G. MacKinnon, ``Approximate asymptotic distribution functions
c for unit root and cointegration tests,'' JBES, 1993 or 1994.
c
c it = 1 ==> \tau_nc test
c it = 2 ==> \tau_c test
c it = 3 ==> \tau_ct test
c it = 4 ==> \tau_ctt test
c
c k = 1 for unit root tests
c k = 2, 3, 4, 5, or 6 for cointegration tests
c
      implicit real*8 (a-h,o-z)
      real*8 gamnc(6,4), gamc(6,4), gamct(6,4), gamctt(6,4)
      real*8 gmsnc(6,3), gmsc(6,3), gmsct(6,3), gmsctt(6,3)
      real*8 tmxnc(6), tmxc(6), tmxct(6), tmxctt(6)
      real*8 tmnnc(6), tmnc(6), tmnct(6), tmnctt(6)
      real*8 tsnc(6), tsc(6), tsct(6), tsctt(6)
      data gamnc /0.4797d0, 1.5578d0, 2.2268d0, 2.7654d0,
     1 3.2684d0, 3.7268d0,
     2 0.93557d0, 0.85580d0, 0.68093d0, 0.64502d0, 0.68051d0, 0.71670d0,
     3 -0.06999d0, -0.20830d0, -0.32362d0, -0.30811d0,
     4 -0.26778d0, -0.23648d0,
     5 0.033066d0, -0.033549d0, -0.054448d0, -0.044946d0,
     6 -0.034972d0, -0.028288d0 /
      data gamc /1.7339d0, 2.1945d0, 2.5893d0, 3.0387d0,
     1 3.5049d0, 3.9489d0,
     2 0.93202d0, 0.64695d0, 0.45168d0, 0.45452d0, 0.52098d0, 0.58933d0,
     3 -0.12745d0, -0.29198d0, -0.36529d0, -0.33666d0,
     4 -0.29158d0, -0.25359d0,
     5 -0.010368d0, -0.042377d0, -0.050074d0, -0.041921d0,
     6 -0.033468d0, -0.027210d0 /
      data gamct /2.5261d0, 2.8500d0, 3.2210d0, 3.6520d0,
     1 4.0712d0, 4.4735d0,
     2 0.61654d0, 0.52720d0, 0.52550d0, 0.59758d0, 0.66428d0, 0.71757d0,
     3 -0.37956d0, -0.36622d0, -0.32685d0, -0.27483d0,
     4 -0.23464d0, -0.20681d0,
     5 -0.060285d0, -0.051695d0, -0.041501d0, -0.032081d0,
     6 -0.025460d0, -0.021196d0 /
      data gamctt /3.0778d0, 3.4713d0, 3.8637d0, 4.2736d0,
     1 4.6679d0, 5.0009d0,
     2 0.49529d0, 0.59670d0, 0.67852d0, 0.76199d0, 0.82618d0, 0.83735d0,
     3 -0.41477d0, -0.32507d0, -0.26286d0, -0.21534d0,
     4 -0.18220d0, -0.16994d0,
     5 -0.059359d0, -0.042286d0, -0.031381d0, -0.024026d0,
     6 -0.019147d0, -0.016928d0 /
      data gmsnc /0.6344d0, 1.9129d0, 2.7648d0, 3.4336d0, 4.0999d0,
     1 4.5388d0,
     2 1.2378d0, 1.3857d0, 1.4502d0, 1.4835d0, 1.5533d0, 1.5344d0,
     3 0.032496d0, 0.035322d0, 0.034186d0, 0.031900d0,
     4 0.035900d0, 0.029807d0 /
      data gmsc /2.1659d0, 2.9200d0, 3.4699d0, 3.9673d0,
     1 4.5509d0, 5.1399d0,
     2 1.4412d0, 1.5012d0, 1.4856d0, 1.4777d0, 1.5338d0, 1.6036d0,
     3 0.038269d0, 0.039796d0, 0.031640d0, 0.026315d0,
     4 0.029545d0, 0.034445d0 /
      data gmsct /3.2512d0, 3.6646d0, 4.0983d0, 4.5844d0,
     1 5.0722d0, 5.5300d0,
     2 1.6047d0, 1.5419d0, 1.5173d0, 1.5338d0, 1.5634d0, 1.5914d0,
     3 0.049588d0, 0.036448d0, 0.029898d0, 0.028796d0,
     4 0.029472d0, 0.030392d0 /
      data gmsctt /4.0003d0, 4.3534d0, 4.7343d0, 5.2140d0,
     1 5.6481d0, 5.9296d0,
     2 1.6580d0, 1.6016d0, 1.5768d0, 1.6077d0, 1.6274d0, 1.5929d0,
     3 0.048288d0, 0.037947d0, 0.032396d0, 0.033449d0,
     4 0.033455d0, 0.028223d0 /
      data tmxnc /99.0d0, 1.51d0, 0.86d0, 0.88d0, 1.05d0, 1.24d0 /
      data tmxc /2.74d0, 0.92d0, 0.55d0, 0.61d0, 0.79d0, 1.00d0 /
      data tmxct /0.70d0, 0.63d0, 0.71d0, 0.93d0, 1.19d0, 1.42d0 /
      data tmxctt /0.54d0, 0.79d0, 1.08d0, 1.43d0, 3.49d0, 1.92d0 /
      data tmnnc /-19.04d0, -19.62d0, -21.21d0, -23.25d0,
     1 -21.63d0, -25.74d0 /
      data tmnc /-18.83d0, -18.86d0, -23.48d0, -28.07d0,
     1 -25.96d0, -23.27d0 /
      data tmnct /-16.18d0, -21.15d0, -25.37d0, -26.63d0,
     1 -26.53d0, -26.18d0 /
      data tmnctt /-17.17d0, -21.10d0, -24.33d0, -24.03d0,
     1 -24.33d0, -28.22d0 /
      data tsnc /-1.04d0, -1.53d0, -2.68d0, -3.09d0, -3.07d0, -3.77d0 /
      data tsc /-1.61d0, -2.62d0, -3.13d0, -3.47d0, -3.78d0, -3.93d0 /
      data tsct /-2.89d0, -3.19d0, -3.50d0, -3.65d0, -3.80d0, -4.36d0 /
      data tsctt /-3.21d0, -3.51d0, -3.81d0, -3.83d0, -4.12d0, -4.63d0 /
c
c check that number of cointegrating variables is correct
c
      if (k.lt.1.or.k.gt.6) then
        write(6,*) 'Invalid argument. k must be between 1 and 6.'
        return
      end if
c
c find out which type of test is being used
c
      if (it.eq.1) goto 10
      if (it.eq.2) goto 20
      if (it.eq.3) goto 30
      if (it.eq.4) goto 40
      write(6,*) 'Invalid argument. It must be 1, 2, 3, or 4.'
      return
c
   10 continue
c \tau_nc test
      arg = tstat
      if (tstat.gt.tmxnc(k)) arg = tmxnc(k)
      if (tstat.lt.tmnnc(k)) arg = tmnnc(k)
      if (arg.lt.tsnc(k)) then
         h = gmsnc(k,1) + gmsnc(k,2)*arg + gmsnc(k,3)*arg**2
      else
         h = gamnc(k,1) + gamnc(k,2)*arg + gamnc(k,3)*arg**2
     1   + gamnc(k,4)*arg**3
      end if
      go to 50
   20 continue
c \tau_c test
      arg = tstat
      if (tstat.gt.tmxc(k)) arg = tmxc(k)
      if (tstat.lt.tmnc(k)) arg = tmnc(k)
      if (arg.lt.tsc(k)) then
         h = gmsc(k,1) + gmsc(k,2)*arg + gmsc(k,3)*arg**2
      else
         h = gamc(k,1) + gamc(k,2)*arg + gamc(k,3)*arg**2
     1   + gamc(k,4)*arg**3
      end if
      go to 50
   30 continue
c \tau_ct test
      arg = tstat
      if (tstat.gt.tmxct(k)) arg = tmxct(k)
      if (tstat.lt.tmnct(k)) arg = tmnct(k)
      if (arg.lt.tsct(k)) then
         h = gmsct(k,1) + gmsct(k,2)*arg + gmsct(k,3)*arg**2
      else
         h = gamct(k,1) + gamct(k,2)*arg + gamct(k,3)*arg**2
     1   + gamct(k,4)*arg**3
      end if
      go to 50
   40 continue
c \tau_ctt test
      arg = tstat
      if (tstat.gt.tmxctt(k)) arg = tmxctt(k)
      if (tstat.lt.tmnctt(k)) arg = tmnctt(k)
      if (arg.lt.tsctt(k)) then
         h = gmsctt(k,1) + gmsctt(k,2)*arg + gmsctt(k,3)*arg**2
      else
         h = gamctt(k,1) + gamctt(k,2)*arg + gamctt(k,3)*arg**2
     1   + gamctt(k,4)*arg**3
      end if
   50 continue
c
c now call routine to compute cumulative normal
c
      call ddnor(h,gauss)
      pval = gauss
c     write(6,*) tstat, h, pval
      return
      end
      subroutine pvz(it,k,tstat,pval)
c
c This subroutine computes P values for z tests using approximations
c from the paper:
c
c James G. MacKinnon, ``Approximate asymptotic distribution functions
c for unit root and cointegration tests,'' JBES, 1993 or 1994.
c
c it = 1 ==> z_nc test
c it = 2 ==> z_c test
c it = 3 ==> z_ct test
c it = 4 ==> z_ctt test
c
c k = 1 for unit root tests
c k = 2, 3, 4, 5, or 6 for cointegration tests
c
      implicit real*8 (a-h,o-z)
      real*8 gamnc(6,5), gamc(6,5), gamct(6,5), gamctt(6,5)
      real*8 gmsnc(6,4), gmsc(6,4), gmsct(6,4), gmsctt(6,4)
      real*8 zsnc(6), zsc(6), zsct(6), zsctt(6)
      data gamnc /0.4927d0, 1.5167d0, 2.2347d0, 2.8239d0,
     1 3.3174d0, 3.7290d0,
     2 0.69060d0, 0.46859d0, 0.39465d0, 0.36265d0, 0.33492d0, 0.30611d0,
     3 .132331d0, 4.2401d-2, 2.2406d-2, 1.6738d-2, 1.2792d-2, 0.9579d-2,
     4 12.099d-3, 2.7939d-3, 0.8746d-3, 0.5408d-3, 0.3416d-3, 0.2087d-3,
     5 0.d0, 7.9601d-5, 1.4239d-5, 0.7449d-5, 0.3894d-5, 0.1943d-5 /
      data gamc /1.7170d0, 2.2394d0, 2.7430d0, 3.2280d0,
     1 3.6583d0, 4.0379d0,
     2 0.55243d0, 0.42377d0, 0.36260d0, 0.33399d0, 0.30934d0, 0.28735d0,
     3 4.3463d-2, 2.4320d-2, 1.5703d-2, 1.2319d-2, 0.9681d-2, 0.7694d-2,
     4 1.6671d-3, 0.9241d-3, 0.4612d-3, 0.3162d-3, 0.2111d-3, 0.1433d-3,
     5 0.d0, 1.4364d-5, 0.5670d-5, 0.3482d-5, 0.1979d-5, 0.1146d-5 /
      data gamct /2.7117d0, 3.0972d0, 3.4594d0, 3.8060d0,
     1 4.1402d0, 4.4497d0,
     2 0.45731d0, 0.40873d0, 0.36326d0, 0.32634d0, 0.29867d0, 0.27534d0,
     3 2.2868d-2, 1.8982d-2, 1.4284d-2, 1.0689d-2, 0.8323d-2, 0.6582d-2,
     4 0.6362d-3, 0.5796d-3, 0.3813d-3, 0.2402d-3, 0.1600d-3, 0.1089d-3,
     5 0.500d-5, 0.7384d-5, 0.4325d-5, 0.2304d-5, 0.1315d-5, 0.0773d-5 /
      data gamctt /3.4671d0, 3.7827d0, 4.0520d0, 4.3311d0,
     1 4.5940d0, 4.8479d0,
     2 0.43476d0, 0.39421d0, 0.34947d0, 0.31625d0, 0.28739d0, 0.26447d0,
     3 1.9231d-2, 1.5699d-2, 1.1772d-2, 0.9126d-2, 0.70695d-2,0.5647d-2,
     4 0.5381d-3, 0.4093d-3, 0.2642d-3, 0.1775d-3, 0.11812d-3,0.0827d-3,
     5 0.6216d-5, 0.4485d-5, 0.2502d-5, 0.1462d-5, 0.0838d-5,0.0518d-5 /
c .70695d-2 above appears in the table as .7070d-2
c and .11812d-3 appears as .1181d-3
c max error declines from .00067 to .00028 when this is done.
      data gmsnc /0.0342d0, 1.3426d0, 3.8607d0, 6.1072d0, 7.7800d0,
     1 4.0253d0,
     2 -0.6376d0, -0.7680d0, -2.4159d0, -3.7250d0, -4.4579d0, -0.8815d0,
     3 0.d0, 0.d0, 0.51293d0, 0.85887d0, 1.00056d0, 0.d0,
     4 -0.03872d0, -0.04104d0, -0.09835d0, -0.13102d0, -0.14014d0,
     5 -0.04887d0 /
      data gmsc /2.2142d0, 1.1662d0, 6.6584d0, 3.3249d0,
     1 4.0356d0, 13.9959d0,
     2 -1.7863d0, 0.1814d0, -4.3486d0, -0.8456d0, -0.9306d0, -8.4314d0,
     3 0.32828d0, -0.36707d0, 1.04705d0, 0.d0, 0.d0, 1.97411d0,
     4 -0.07727d0, 0.d0, -0.15011d0, -0.04818d0, -0.04776d0,-0.22234d0 /
      data gmsct /4.6476d0, 7.2453d0, 3.4893d0, 1.6604d0,
     1 2.0060d0, 11.1626d0,
     2 -2.8932d0, -4.7021d0, -0.8914d0, 1.0375d0, 1.1197d0, -5.6858d0,
     3 0.58320d0, 1.12700d0, 0.d0, -0.53377d0, -0.55315d0, 1.21479d0,
     4 -0.09990d0, -0.15665d0, -0.04755d0, 0.d0, 0.d0, -0.15428d0 /
      data gmsctt /3.6739d0, 3.9783d0, 2.0062d0, 4.9218d0,
     1 5.1433d0, 23.6812d0,
     2 -1.1549d0, -1.0619d0, 0.8907d0, -1.0663d0, -0.9877d0, -14.6485d0,
     3 0.d0, 0.d0, -0.51708d0, 0.d0, 0.d0, 3.42909d0,
     4 -0.03947d0, -0.04394d0, 0.d0, -0.04691d0, -0.04993d0,-0.33794d0 /
      data zsnc /-2.9d0, -8.7d0, -14.8d0, -20.9d0, -25.7d0, -30.5d0 /
      data zsc /-8.9d0, -14.3d0, -19.5d0, -25.1d0, -29.6d0, -34.4d0 /
      data zsct /-15.0d0, -19.6d0, -25.3d0, -29.6d0, -31.8d0, -38.4d0 /
      data zsctt /-20.7d0, -25.3d0, -29.9d0, -34.4d0, -38.5d0, -44.2d0 /
c
c check that number of cointegrating variables is correct
c
      if (k.lt.1.or.k.gt.6) then
        write(6,*) 'Invalid argument. k must be between 1 and 6.'
        return
      end if
c
c find out which type of test is being used
c
      if (it.eq.1) goto 10
      if (it.eq.2) goto 20
      if (it.eq.3) goto 30
      if (it.eq.4) goto 40
      write(6,*) 'Invalid argument. It must be 1, 2, 3, or 4.'
      return
c
   10 continue
c z_nc test
      arg = tstat
      if (arg.lt.zsnc(k)) then
         arg = log(-arg)
         h = gmsnc(k,1) + gmsnc(k,2)*arg + gmsnc(k,3)*arg**2
     1   + gmsnc(k,4)*arg**3
      else
         h = gamnc(k,1) + gamnc(k,2)*arg + gamnc(k,3)*arg**2
     1   + gamnc(k,4)*arg**3 + gamnc(k,5)*arg**4
      end if
      go to 50
   20 continue
c z_c test
      arg = tstat
      if (arg.lt.zsc(k)) then
         arg = log(-arg)
         h = gmsc(k,1) + gmsc(k,2)*arg + gmsc(k,3)*arg**2
     1   + gmsc(k,4)*arg**3
      else
         h = gamc(k,1) + gamc(k,2)*arg + gamc(k,3)*arg**2
     1   + gamc(k,4)*arg**3 + gamc(k,5)*arg**4
      end if
      go to 50
   30 continue
c z_ct test
      arg = tstat
      if (arg.lt.zsct(k)) then
         arg = log(-arg)
         h = gmsct(k,1) + gmsct(k,2)*arg + gmsct(k,3)*arg**2
     1   + gmsct(k,4)*arg**3
      else
         h = gamct(k,1) + gamct(k,2)*arg + gamct(k,3)*arg**2
     1   + gamct(k,4)*arg**3 + gamct(k,5)*arg**4
      end if
      go to 50
   40 continue
c z_ctt test
      arg = tstat
      if (arg.lt.zsctt(k)) then
         arg = log(-arg)
         h = gmsctt(k,1) + gmsctt(k,2)*arg + gmsctt(k,3)*arg**2
     1   + gmsctt(k,4)*arg**3
      else
         h = gamctt(k,1) + gamctt(k,2)*arg + gamctt(k,3)*arg**2
     1   + gamctt(k,4)*arg**3 + gamctt(k,5)*arg**4
      end if
   50 continue
c
c now call routine to compute cumulative normal
c
      call ddnor(h,gauss)
      pval = gauss
c     write(6,*) tstat, h, pval
      return
      end
      subroutine ddnor(ystar,gauss)
c
c written by James G. MacKinnon, 1977.
c
c This subroutine uses Cody's method to evaluate the cumulative
c normal distribution. it is probably accurate to 19 or 20
c significant digits on a machine with that much accuracy. It is
c based on the Cody article referred to in the documentation
c for IMSL subroutine mdnor.
c
      implicit real*8(a-h,o-z)
      real*8 p(6), q(5), a(9), b(8), c(5), d(4)
      data p(1)/-6.58749161529837803157d-04/,
     1     p(2)/-1.60837851487422766278d-02/,
     2     p(3)/-1.25781726111229246204d-01/,
     3     p(4)/-3.60344899949804439429d-01/,
     4     p(5)/-3.05326634961232344035d-01/,
     5     p(6)/-1.63153871373020978498d-02/
      data q(1)/2.33520497626869185443d-03/,
     1     q(2)/6.05183413124413191178d-02/,
     2     q(3)/5.27905102951428412248d-01/,
     3     q(4)/1.87295284992346047209d00/,
     4     q(5)/2.56852019228982242072d00/
      data a(1)/1.23033935479799725272d03/,
     1     a(2)/2.05107837782607146532d03/,
     2     a(3)/1.71204761263407058314d03/,
     3     a(4)/8.81952221241769090411d02/,
     4     a(5)/2.98635138197400131132d02/,
     5     a(6)/6.61191906371416294775d01/,
     6     a(7)/8.88314979438837594118d00/,
     7     a(8)/5.64188496988670089180d-01/,
     8     a(9)/2.15311535474403846343d-08/
      data b(1)/1.23033935480374942043d03/,
     1     b(2)/3.43936767414372163696d03/,
     2     b(3)/4.36261909014324715820d03/,
     3     b(4)/3.29079923573345962678d03/,
     4     b(5)/1.62138957456669018874d03/,
     5     b(6)/5.37181101862009857509d02/,
     6     b(7)/1.17693950891312499305d02/,
     7     b(8)/1.57449261107098347253d01/
      data c(1)/3.209377589138469472562d03/,
     1     c(2)/3.774852376853020208137d02/,
     2     c(3)/1.138641541510501556495d02/,
     3     c(4)/3.161123743870565596947d00/,
     4     c(5)/1.857777061846031526730d-01/
      data d(1)/2.844236833439170622273d03/,
     1     d(2)/1.282616526077372275645d03/,
     2     d(3)/2.440246379344441733056d02/,
     3     d(4)/2.360129095234412093499d01/
      data orpi/.5641895835477562869483d0/,
     1   root2/.70710678118654752440083d0/
c
      isw = 1
      y = ystar
      if (ystar.lt.-16.d0) y = -16.d0
      if (ystar.gt.16.d0) y = 16.d0
      x = -y*root2
      if(x.gt.0.d0) go to 1
      if(x.lt.0.d0) go to 2
      gauss = .5d0
      return
    2 continue
      x = - x
      isw = -1
    1 continue
      if(x.lt..477d0) go to 10
      if(x.le.4.d0) go to 20
c  evaluate erfc for x.gt.4.0
      x2 = x*x
      xm2 = 1.d0/x2
      xm4 = xm2*xm2
      xm6 = xm4*xm2
      xm8 = xm4*xm4
      xm10 = xm6*xm4
      top = p(1) + p(2)*xm2 + p(3)*xm4 + p(4)*xm6 + p(5)*xm8 + p(6)*xm10
      bot = q(1) + q(2)*xm2 + q(3)*xm4 + q(4)*xm6 + q(5)*xm8 + xm10
      crap = orpi + top/(bot*x2)
      erfc = dexp(-x2)*crap/x
c
      if(isw.eq.-1) erfc = 2.d0 - erfc
      gauss = erfc*.5d0
      return
   20 continue
c  evaluate erfc for .477.lt.x.le.4.0
      x2 = x*x
      x3 = x2*x
      x4 = x2*x2
      x5 = x3*x2
      x6 = x3*x3
      x7 = x3*x4
      x8 = x4*x4
      top = a(1) + a(2)*x + a(3)*x2 + a(4)*x3 + a(5)*x4 + a(6)*x5 +
     1 a(7)*x6 + a(8)*x7 + a(9)*x8
      bot = b(1) + b(2)*x + b(3)*x2 + b(4)*x3 + b(5)*x4 + b(6)*x5 +
     1 b(7)*x6 + b(8)*x7 + x8
      erfc = dexp(-x2)*top/bot
c
      if(isw.eq.-1) erfc = 2.d0 - erfc
      gauss = erfc*.5d0
      return
   10 continue
c  evaluate erf for x.lt..477
      x2 = x*x
      x4 = x2*x2
      x6 = x4*x2
      x8 = x4*x4
      top = c(1) + c(2)*x2 + c(3)*x4 + c(4)*x6 + c(5)*x8
      bot = d(1) + d(2)*x2 + d(3)*x4 + d(4)*x6 + x8
      erf = x*top/bot
c
      erf = erf*isw
      erfc = 1.d0 - erf
      gauss = erfc*.5d0
      return
      end
