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
.