function cbrt (x)
  c june 1977 edition.   w. fullerton, c3, los alamos scientific lab.
        dimension cbrt2(5)
        external  r1mach, r9pak
  c
        data cbrt2(1) / 0.6299605249 4743658e0 /
        data cbrt2(2) / 0.7937005259 8409974e0 /
        data cbrt2(3) / 1.0e0 /
        data cbrt2(4) / 1.2599210498 9487316e0 /
        data cbrt2(5) / 1.5874010519 6819947e0 /
  c
        data niter / 0 /
  c
        if (niter.eq.0) niter = 1.443*alog(-.106*alog(0.1*r1mach(3))) + 1.
  c
        cbrt = 0.0
        if (x.eq.0.) return
  c
        call r9upak (abs(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
        cbrt = .439581e0 + y*(.928549e0 + y*(-.512653e0 + y*.144586e0))
  c
        do 10 iter=1,niter
          cbrtsq = cbrt*cbrt
          cbrt = cbrt + (y-cbrt*cbrtsq)/(3.0*cbrtsq)
   10   continue
  c
        cbrt = r9pak (cbrt2(irem)*sign(cbrt,x), ixpnt)
        return
  c
        end