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