function cbrt (x)
c june 1977 edition.   w. fullerton, c3, los alamos scientific lab.
      dimension cbrt2(5)
      external alog, 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