double precision function dtan (x)
c may 1980 edition.   w. fullerton, c3, los alamos scientific lab.
      double precision x, tancs(19), ainty, ainty2, pi2rec, sqeps,
     1  xmax, xsml, y, yrem, prodbg,  dint, dcsevl, d1mach, dsqrt
      external d1mach, dcsevl, dint, dsqrt, initds
c
c series for tan        on the interval  0.          to  6.25000e-02
c                                        with weighted error   1.44e-32
c                                         log weighted error  31.84
c                               significant figures required  30.92
c                                    decimal places required  32.48
c
      data tan cs(  1) / +.2262793276 3129357846 5786365317 52 d+0     /
      data tan cs(  2) / +.4301791314 6548961775 5834107480 67 d-1     /
      data tan cs(  3) / +.6854461068 2565088756 9294736234 61 d-3     /
      data tan cs(  4) / +.1104532694 7597098383 5788493696 96 d-4     /
      data tan cs(  5) / +.1781747790 3926312943 2385125889 40 d-6     /
      data tan cs(  6) / +.2874496858 2365265947 5296468324 71 d-8     /
      data tan cs(  7) / +.4637485419 5902995494 1374782343 63 d-10    /
      data tan cs(  8) / +.7481760904 1556138502 3416333082 15 d-12    /
      data tan cs(  9) / +.1207049700 2957544801 6445169478 24 d-13    /
      data tan cs( 10) / +.1947361081 2823019305 5138585845 33 d-15    /
      data tan cs( 11) / +.3141722487 4732446504 6145860266 66 d-17    /
      data tan cs( 12) / +.5068613255 5800153941 9048917333 33 d-19    /
      data tan cs( 13) / +.8177310515 9836540043 9799466666 66 d-21    /
      data tan cs( 14) / +.1319264341 2147384408 9514666666 66 d-22    /
      data tan cs( 15) / +.2128399549 7042377309 8666666666 66 d-24    /
      data tan cs( 16) / +.3433796019 2345945292 8000000000 00 d-26    /
      data tan cs( 17) / +.5539822212 1173811200 0000000000 00 d-28    /
      data tan cs( 18) / +.8937522779 4352810666 6666666666 66 d-30    /
      data tan cs( 19) / +.1441911137 1369130666 6666666666 66 d-31    /
c
c pi2rec = 2/pi - 0.625
      data pi2rec / .01161977236 7581343075 5350534900 57 d0 /
      data nterms, xmax, xsml, sqeps / 0, 3*0.d0 /
c
      if (nterms.ne.0) go to 10
      nterms = initds (tancs, 19, 0.1*sngl(d1mach(3)))
      xmax = 1.0d0/d1mach(4)
      xsml = dsqrt (3.0d0*d1mach(3))
      sqeps = dsqrt (d1mach(4))
c
 10   y = dabs(x)
      if (y.gt.xmax) call seteru (
     1  43hdtan    no precision because dabs(x) is big, 43, 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 = dint (y)
      yrem = y - ainty
      prodbg = 0.625d0*ainty
      ainty = dint (prodbg)
      y = (prodbg-ainty) + 0.625d0*yrem + pi2rec*y
      ainty2 = dint (y)
      ainty = ainty + ainty2
      y = y - ainty2
c
      ifn = dmod (ainty, 2.0d0)
      if (ifn.eq.1) y = 1.0d0 - y
      if (1.0d0-y.lt.dabs(x)*sqeps) call seteru ( 70hdtan    answer lt h
     1alf precision, dabs(x) big or x near pi/2 or 3*pi/2, 70, 1, 1)
      if (y.eq.1.d0) call seteru (27hdtan    x is pi/2 or 3*pi/2, 27, 3,
     1  2)
c
      if (y.gt.0.25d0) go to 20
      dtan = y
      if (y.gt.xsml) dtan = y*(1.5d0 + dcsevl (32.d0*y*y-1.d0, tancs,
     1  nterms) )
      go to 40
c
 20   if (y.gt.0.5d0) go to 30
      dtan = 0.5d0*y*(1.5d0 + dcsevl (8.d0*y*y-1.d0, tancs, nterms) )
      dtan = 2.0d0*dtan / (1.d0-dtan*dtan)
      go to 40
c
 30   dtan = 0.25d0*y*(1.5d0 + dcsevl (2.d0*y*y-1.d0, tancs, nterms) )
      dtan = 2.0d0*dtan / (1.0d0-dtan*dtan)
      dtan = 2.0d0*dtan / (1.0d0-dtan*dtan)
c
 40   if (x.ne.0.d0) dtan = dsign (dtan, x)
      if (ifn.eq.1) dtan = -dtan
c
      return
      end