function atan (x)
c jan 1978 edition.   w. fullerton, c3, los alamos scientific lab.
      dimension atancs(9), tanp8(3), conpi8(4), pi8(4)
      external csevl, inits, r1mach, sqrt
c
c series for atan       on the interval  0.          to  4.00000d-02
c                                        with weighted error   1.00e-17
c                                         log weighted error  17.00
c                               significant figures required  16.38
c                                    decimal places required  17.48
c
      data atancs( 1) /    .4869011034 9241406e0 /
      data atancs( 2) /   -.0065108316 36717464e0 /
      data atancs( 3) /    .0000383458 28265245e0 /
      data atancs( 4) /   -.0000002687 22128762e0 /
      data atancs( 5) /    .0000000020 50093098e0 /
      data atancs( 6) /   -.0000000000 16450717e0 /
      data atancs( 7) /    .0000000000 00136509e0 /
      data atancs( 8) /   -.0000000000 00001160e0 /
      data atancs( 9) /    .0000000000 00000010e0 /
c
c xbndn = tan((2*n-1)*pi/16.0)
      data xbnd1 / +.1989123673 79658006 e+0                     /
      data xbnd2 / +.6681786379 19298919 e+0                     /
      data xbnd3 / +1.496605762 66548901 e+0                     /
      data xbnd4 / +5.027339492 12584810 e+0                     /
c
c tanp8(n) = tan(n*pi/8.)
      data tanp8(1) / .4142135623 73095048 e+0                     /
      data tanp8(2) / 1.0 e0 /
      data tanp8(3) / 2.414213562 37309504 e+0                     /
c
c conpi8(n) + pi8(n) = n*pi/8.0
      data conpi8(1) / 0.375 e0 /
      data conpi8(2) / 0.75 e0 /
      data conpi8(3) / 1.125 e0 /
      data conpi8(4) / 1.5   e0 /
c
      data pi8   (  1) / +.1769908169 87241548 e-1                     /
      data pi8   (  2) / +.3539816339 74483096 e-1                     /
      data pi8   (  3) / +.5309724509 61724644 e-1                     /
      data pi8   (  4) / 0.0707963267 948966192 e0 /
c
      data nterms, sqeps, xbig / 0, 2*0.0 /
c
      if (nterms.ne.0) go to 10
      nterms = inits (atancs, 9, 0.1*r1mach(3))
      sqeps = sqrt (6.0*r1mach(3))
      xbig = 1.0/r1mach(3)
c
 10   y = abs(x)
      if (y.gt.xbnd1) go to 20
c
      atan = x
      if (y.gt.sqeps) atan = x*(0.75+csevl(50.*y**2-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.0 + y*tanp8(n))
      atan = sign (conpi8(n) + (pi8(n) + t*(0.75 +
     1  csevl(50.0*t**2-1.0, atancs, nterms)) ), x)
      return
c
 30   atan = conpi8(4) + pi8(4)
      if (y.lt.xbig) atan = conpi8(4) + (pi8(4) - (0.75 +
     1  csevl (50./y**2-1.0, atancs, nterms))/y )
      atan = sign (atan, x)
      return
c
      end