double precision function dcbrt (x)
c june 1977 edition.   w. fullerton, c3, los alamos scientific lab.
      double precision x, cbrt2(5), y, cbrtsq,  d9pak, d1mach
      external d1mach, d9pak
c
      data cbrt2(1) / 0.6299605249 4743658238 3605303639 11 d0 /
      data cbrt2(2) / 0.7937005259 8409973737 5852819636 15 d0 /
      data cbrt2(3) / 1.0 d0 /
      data cbrt2(4) / 1.2599210498 9487316476 7210607278 23 d0 /
      data cbrt2(5) / 1.5874010519 6819947475 1705639272 31 d0 /
c
      data niter / 0 /
c
      if (niter.eq.0) niter = 1.443*alog(-.106*alog(0.1*sngl(d1mach(3)))
     1  ) + 1.0
c
      dcbrt = 0.d0
      if (x.eq.0.d0) return
c
      call d9upak (dabs(x), y, n)
      ixpnt = n/3
      irem = n - 3*ixpnt + 3
c
c the approximation below is a generalized chebyshev series converted
c to polynomial form.  the approx is nearly best in the sense of
c relative error with 4.085 digits accuracy.
c
      z = y
      dcbrt = .439581e0 + z*(.928549e0 + z*(-.512653e0 + z*.144586e0))
c
      do 10 iter=1,niter
        cbrtsq = dcbrt*dcbrt
        dcbrt = dcbrt + (y-dcbrt*cbrtsq)/(3.d0*cbrtsq)
 10   continue
c
      dcbrt = d9pak (cbrt2(irem)*dsign(dcbrt,x), ixpnt)
      return
c
      end