function tan (x)
c may 1980 edition.   w. fullerton, c3, los alamos scientific lab.
      dimension tancs(11)
      external aint, csevl, inits, r1mach, sqrt
c
c series for tan        on the interval  0.          to  6.25000d-02
c                                        with weighted error   3.14e-18
c                                         log weighted error  17.50
c                               significant figures required  16.59
c                                    decimal places required  18.02
c
      data tan cs( 1) /    .2262793276 31293578e0 /
      data tan cs( 2) /    .0430179131 465489618e0 /
      data tan cs( 3) /    .0006854461 068256508e0 /
      data tan cs( 4) /    .0000110453 269475970e0 /
      data tan cs( 5) /    .0000001781 747790392e0 /
      data tan cs( 6) /    .0000000028 744968582e0 /
      data tan cs( 7) /    .0000000000 463748541e0 /
      data tan cs( 8) /    .0000000000 007481760e0 /
      data tan cs( 9) /    .0000000000 000120704e0 /
      data tan cs(10) /    .0000000000 000001947e0 /
      data tan cs(11) /    .0000000000 000000031e0 /
c
c pi2rec = 2/pi - 0.625
      data pi2rec / .01161977236 75813430 e0 /
      data nterms, xmax, xsml, sqeps /0, 0., 0., 0./
c
      if (nterms.ne.0) go to 10
      nterms = inits (tancs, 11, 0.1*r1mach(3))
      xmax = 1.0/r1mach(4)
      xsml = sqrt (3.0*r1mach(3))
      sqeps = sqrt (r1mach(4))
c
 10   y = abs(x)
      if (y.gt.xmax) call seteru (
     1  42htan     no precision because abs(x) is big, 42, 2, 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. - y
c
      if (1.-y.lt.abs(x)*sqeps) call seteru ( 
     $61htan: answer lt half prec, abs(x) big or x near pi/2 or 3*pi/2,
     2    69, 1, 1)
      if (y.eq.1.) call seteru (27htan     x is pi/2 or 3*pi/2, 27, 3,2)
c
      if (y.gt.0.25) go to 20
      tan = y
      if (y.gt.xsml) tan = y*(1.5 + csevl(32.*y**2-1., tancs, nterms))
      go to 40
c
 20   if (y.gt.0.5) go to 30
      tan = 0.5*y*(1.5 + csevl (8.*y**2-1., tancs, nterms))
      tan = 2.0*tan/(1.-tan**2)
      go to 40
c
 30   tan = 0.25*y*(1.5 + csevl (2.*y**2-1., tancs, nterms))
      tan = 2.0*tan/(1.-tan**2)
      tan = 2.0*tan/(1.-tan**2)
c
 40   if (x.ne.0.) tan = sign (tan, x)
      if (ifn.eq.1) tan = -tan
c
      return
      end