function exp (x)
c may 1980 edition.  w. fullerton c3, los alamos scientific lab.
      dimension expcs(8), twon16(17)
      external aint, alog, csevl, inits, r1mach, r9pak
c
c series for exp        on the interval -1.00000d+00 to  1.00000d+00
c                                        with weighted error   4.81e-18
c                                         log weighted error  17.32
c                               significant figures required  15.95
c                                    decimal places required  17.77
c
      data exp cs( 1) /    .0866569493 31498571e0 /
      data exp cs( 2) /    .0009384948 69299839e0 /
      data exp cs( 3) /    .0000067760 39709981e0 /
      data exp cs( 4) /    .0000000366 93120039e0 /
      data exp cs( 5) /    .0000000001 58959053e0 /
      data exp cs( 6) /    .0000000000 00573859e0 /
      data exp cs( 7) /    .0000000000 00001775e0 /
      data exp cs( 8) /    .0000000000 00000004e0 /
c
      data twon16( 1) /  0.e0                  /
      data twon16( 2) /   .44273782427413840e-1 /
      data twon16( 3) /   .90507732665257659e-1 /
      data twon16( 4) /   .13878863475669165e0 /
      data twon16( 5) /   .18920711500272107e0 /
      data twon16( 6) /   .24185781207348405e0 /
      data twon16( 7) /   .29683955465100967e0 /
      data twon16( 8) /   .35425554693689273e0 /
      data twon16( 9) /   .41421356237309505e0 /
      data twon16(10) /   .47682614593949931e0 /
      data twon16(11) /   .54221082540794082e0 /
      data twon16(12) /   .61049033194925431e0 /
      data twon16(13) /   .68179283050742909e0 /
      data twon16(14) /   .75625216037329948e0 /
      data twon16(15) /   .83400808640934246e0 /
      data twon16(16) /   .91520656139714729e0 /
      data twon16(17) /  1.0e0                 /
c
c aln216 is 16.0/alog(2.) - 23.0
      data aln216 /  .083120654223414518e0 /
      data nterms, xmin, xmax /0, 0., 0./
c
      if (nterms.ne.0) go to 10
      nterms = inits (expcs, 8, 0.1*r1mach(3))
      xmin = alog (r1mach(1)) + 0.01
      xmax = alog(r1mach(2)) - 0.001
c
 10   if (x.lt.xmin) go to 20
      if (x.gt.xmax) call seteru (
     1  30hexp     x so big exp overflows, 30, 2, 2)
c
      xint = aint(x)
      y = x - xint
c
      y = 23.0*y + x*aln216
      n = y
      f = y - float(n)
      n = 23.0*xint + float(n)
      n16 = n/16
      if (n.lt.0) n16 = n16 - 1
      ndx = n - 16*n16 + 1
c
      exp = 1.0 + (twon16(ndx) + f*(1.0+twon16(ndx)) *
     1  csevl (f, expcs, nterms) )
c
      exp = r9pak (exp, n16)
      return
c
 20   call seteru (33hexp     x so small exp underflows, 33, 1, 0)
      exp = 0.0
      return
c
      end