double precision function datan (x)
c jan 1978 edition.   w. fullerton, c3, los alamos scientific lab.
      double precision x, atancs(16), tanp8(3), conpi8(4), pi8(4),
     1  sqeps, xbig, t, xbnd1, xbnd2, xbnd3, xbnd4, y, dcsevl,
     2  d1mach, dsqrt
      external d1mach, dcsevl, dsqrt, initds
c
c series for atan       on the interval  0.          to  4.00000e-02
c                                        with weighted error   4.83e-32
c                                         log weighted error  31.32
c                               significant figures required  30.70
c                                    decimal places required  31.92
c
      data atancs(  1) / +.4869011034 9241406474 6369159028 91 d+0     /
      data atancs(  2) / -.6510831636 7174641818 8697949459 74 d-2     /
      data atancs(  3) / +.3834582826 5245177653 5699924304 56 d-4     /
      data atancs(  4) / -.2687221287 6223146539 5954105187 88 d-6     /
      data atancs(  5) / +.2050093098 5824269846 6365146866 88 d-8     /
      data atancs(  6) / -.1645071739 5484269455 7341352853 48 d-10    /
      data atancs(  7) / +.1365097527 4390773423 8135284844 28 d-12    /
      data atancs(  8) / -.1160177959 1998246322 8913098346 66 d-14    /
      data atancs(  9) / +.1003833394 3866273835 7976574026 66 d-16    /
      data atancs( 10) / -.8807274715 2163859327 0736960000 00 d-19    /
      data atancs( 11) / +.7813632100 5661722180 5802666666 66 d-21    /
      data atancs( 12) / -.6995453514 8267456086 6133333333 33 d-23    /
      data atancs( 13) / +.6310590571 3702136004 2666666666 66 d-25    /
      data atancs( 14) / -.5729607537 0213874346 6666666666 66 d-27    /
      data atancs( 15) / +.5227479628 0602282666 6666666666 66 d-29    /
      data atancs( 16) / -.4832790391 1608320000 0000000000 00 d-31    /
c
c xbndn = tan((2*n-1)*pi/16.0)
      data xbnd1  / +.1989123673 7965800691 1597622644 67 d+0      /
      data xbnd2 / +.6681786379 1929891999 7757686523 08 d+0     /
      data xbnd3 / +1.496605762 6654890176 0113513494 24 d+0     /
      data xbnd4 / +5.027339492 1258481045 1497507106 40 d+0     /
c
c tanp8(n) = tan(n*pi/8.0)
      data tanp8 (  1) / +.4142135623 7309504880 1688724209 69 d+0     /
      data tanp8 (  2) / +1.0d0                                        /
      data tanp8 (  3) / +2.414213562 3730950488 0168872420 96 d+0     /
c
c conpi8(n) + pi8(n) = n*pi/8.0
      data conpi8(1) / 0.375 d0 /
      data conpi8(2) / 0.75  d0 /
      data conpi8(3) / 1.125 d0 /
      data conpi8(4) / 1.5   d0 /
c
      data pi8   (  1) / +.1769908169 8724154807 8304229099 37 d-1     /
      data pi8   (  2) / +.3539816339 7448309615 6608458198 75 d-1     /
      data pi8   (  3) / +.5309724509 6172464423 4912687298 13 d-1     /
      data pi8   (  4) / +.7079632679 4896619231 3216916397 51 d-1     /
c
      data nterms, sqeps, xbig / 0, 2*0.0d0 /
c
      if (nterms.ne.0) go to 10
      nterms = initds (atancs, 16, sngl(0.1d0*d1mach(3)))
      sqeps = dsqrt (6.0d0*d1mach(3))
      xbig = 1.0d0/d1mach(3)
c
 10   y = dabs(x)
      if (y.gt.xbnd1) go to 20
c
      datan = x
      if (y.gt.sqeps) datan = x*(0.75d0 + dcsevl (50.d0*y*y-1.0d0,
     1  atancs, nterms))
      return
c
 20   if (y.gt.xbnd4) go to 30
c
      n = 1
      if (y.gt.xbnd2) n = 2
      if (y.gt.xbnd3) n = 3
c
      t = (y - tanp8(n)) / (1.0d0 + y*tanp8(n))
      datan = dsign (conpi8(n) + (pi8(n) + t*(0.75d0 +
     1  dcsevl (50.0d0*t*t-1.0d0, atancs, nterms)) ), x)
      return
c
 30   datan = conpi8(4) + pi8(4)
      if (y.lt.xbig) datan = conpi8(4) + (pi8(4) - (0.75d0 +
     1  dcsevl (50.0d0/y**2-1.0d0, atancs, nterms))/y )
      datan = dsign (datan, x)
      return
c
      end