function cos (x)
c august 1980 edition.  w. fullerton, los alamos scientific lab.
c
c this routine is based on the algorithm of cody and waite
c in argonne tm-321, software manual working note number 1.
c
      dimension sincs(10)
      external aint, csevl, inits, r1mach, sqrt
c
c series for sin    on the interval  0.00000e+00 to  2.46740e+00
c                                        with weighted error   2.06e-19
c                                         log weighted error  18.69
c                               significant figures required  18.10
c                                    decimal places required  19.19
c
      data sin cs(  1) /    -0.3749911549 5587317584 0e0/
      data sin cs(  2) /    -0.1816031552 3725020186 4e0/
      data sin cs(  3) /     0.0058047092 7459863355 9e0/
      data sin cs(  4) /    -0.0000869543 1177934075 7e0/
      data sin cs(  5) /     0.0000007543 7014808885 1e0/
      data sin cs(  6) /    -0.0000000042 6712966505 6e0/
      data sin cs(  7) /     0.0000000000 1698042294 5e0/
      data sin cs(  8) /    -0.0000000000 0005012057 9e0/
      data sin cs(  9) /     0.0000000000 0000011410 1e0/
      data sin cs( 10) /    -0.0000000000 0000000020 6e0/
c
c pihi + pilo = pi.  pihi is exactly representable on all machines
c with at least 8 bits of precision.  whether it is exactly
c represented depends on the compiler.  this routine is more
c accurate if it is exactly represented.
      data pihi / 3.140625e0 /
      data pilo / 9.676535897 9323846e-4 /
      data pi2 / 1.5707963267 94896619e0 /
      data pirec / 0.3183098861 8379067e0 /
      data pi2rec / 0.63661977236 7581343e0 /
      data ntsn, xsml, xwarn, xmax / 0, 3*0.0 /
c
      if (ntsn.ne.0) go to 10
      ntsn = inits (sincs, 10, 0.1*r1mach(3))
c
      xsml = sqrt (2.0*r1mach(3))
      xmax = 1.0/r1mach(4)
      xwarn = sqrt(xmax)
c
 10   absx = abs(x)
      y = absx + pi2
      if (y.gt.xmax) call seteru (
     1  42hcos     no precision because abs(x) is big, 42, 2, 2)
      if (y.gt.xwarn) call seteru (
     1  54hcos     answer lt half precision because abs(x) is big,
     2  54, 1, 1)
c
      cos = 1.0
      if (absx.lt.xsml) return
c
      xn = aint (y*pirec+0.5)
      n2 = amod (xn, 2.0) + 0.5
      xn = xn - 0.5
      f = (absx-xn*pihi) - xn*pilo
c
      cos = f + f*csevl (2.0*(f*pi2rec)**2 - 1.0, sincs, ntsn)
      if (n2.ne.0) cos = -cos
      if (abs(cos).gt.1.0) cos = sign (1.0, cos)
c
      return
      end