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