function cot (x) c may 1980 edition. w. fullerton, c3, los alamos scientific lab. dimension cotcs(8) external csevl, inits, r1mach c c series for cot on the interval 0. to 6.25000d-02 c with weighted error 3.76e-17 c log weighted error 16.42 c significant figures required 15.51 c decimal places required 16.88 c data cot cs( 1) / .2402591609 8295630e0 / data cot cs( 2) / -.0165330316 01500228e0 / data cot cs( 3) / -.0000429983 91931724e0 / data cot cs( 4) / -.0000001592 83223327e0 / data cot cs( 5) / -.0000000006 19109313e0 / data cot cs( 6) / -.0000000000 02430197e0 / data cot cs( 7) / -.0000000000 00009560e0 / data cot cs( 8) / -.0000000000 00000037e0 / c c pi2rec = 2/pi - 0.625 data pi2rec / .01161977236 75813430 e0 / data nterms, xmax, xsml, xmin, sqeps / 0, 4*0.0 / c if (nterms.ne.0) go to 10 nterms = inits (cotcs, 8, 0.1*r1mach(3)) xmax = 1.0/r1mach(4) xsml = sqrt (3.0*r1mach(3)) xmin = exp ( amax1(alog(r1mach(1)), -alog(r1mach(2))) + 0.01) sqeps = sqrt (r1mach(4)) c 10 y = abs(x) if (abs(x).lt.xmin) call seteru ( 1 48hcot abs(x) is zero or so small cot overflows, 48, 2, 2) if (y.gt.xmax) call seteru ( 1 42hcot no precision because abs(x) is big, 42, 3, 2) c c carefully compute y * (2/pi) = (aint(y) + rem(y)) * (.625 + pi2rec) c = aint(.625*y) + rem(.625*y) + y*pi2rec = aint(.625*y) + z c = aint(.625*y) + aint(z) + rem(z) c ainty = aint (y) yrem = y - ainty prodbg = 0.625*ainty ainty = aint (prodbg) y = (prodbg-ainty) + 0.625*yrem + y*pi2rec ainty2 = aint (y) ainty = ainty + ainty2 y = y - ainty2 c ifn = amod (ainty, 2.) if (ifn.eq.1) y = 1.0 - y c if (abs(x).gt.0.5 .and. y.lt.abs(x)*sqeps) call seteru ( 1 72hcot answer lt half precision, abs(x) too big or x near n* 2pi (n.ne.0), 72, 1, 1) c if (y.gt.0.25) go to 20 if (y.eq.0.0) call seteru (29hcot x is a multiple of pi, 1 29, 4, 2) cot = 1.0/y if (y.gt.xsml) cot = (0.5 + csevl (32.0*y*y-1., cotcs, nterms)) /y go to 40 c 20 if (y.gt.0.5) go to 30 cot = (0.5 + csevl (8.0*y*y-1., cotcs, nterms)) / (0.5*y) cot = (cot**2 - 1.0) * 0.5 / cot go to 40 c 30 cot = (0.5 + csevl (2.0*y*y-1., cotcs, nterms)) / (0.25*y) cot = (cot**2 - 1.0) * 0.5 / cot cot = (cot**2 - 1.0) * 0.5 / cot c 40 if (x.ne.0.) cot = sign (cot, x) if (ifn.eq.1) cot = -cot c return end