[dcbrtresult,x]=dcbrt(x);

function [dcbrtresult,x]=dcbrt(x);
dcbrtresult=[];
persistent cbrt2 cbrtsq firstCall irem iter ixpnt n niter y z ; if isempty(firstCall),firstCall=1;end; 

;
if isempty(irem), irem=0; end;
if isempty(iter), iter=0; end;
if isempty(ixpnt), ixpnt=0; end;
if isempty(n), n=0; end;
if isempty(niter), niter=0; end;
if isempty(z), z=0; end;
%***BEGIN PROLOGUE  DCBRT
%***PURPOSE  Compute the cube root.
%***LIBRARY   SLATEC (FNLIB)
%***CATEGORY  C2
%***TYPE      doubleprecision (CBRT-S, DCBRT-D, CCBRT-C)
%***KEYWORDS  CUBE ROOT, ELEMENTARY FUNCTIONS, FNLIB, ROOTS
%***AUTHOR  Fullerton, W., (LANL)
%***DESCRIPTION
%
% DCBRT(X) calculates the doubleprecision cube root for
% doubleprecision argument X.
%
%***REFERENCES  (NONE)
%***ROUTINES CALLED  D1MACH, D9PAK, D9UPAK
%***REVISION HISTORY  (YYMMDD)
%   770601  DATE WRITTEN
%   890531  Changed all specific intrinsics to generic.  (WRB)
%   890531  REVISION DATE from Version 3.2
%   891214  Prologue converted to Version 4.0 format.  (BAB)
%***end PROLOGUE  DCBRT
if isempty(cbrt2), cbrt2=zeros(1,5); end;
if isempty(y), y=0; end;
if isempty(cbrtsq), cbrtsq=0; end;
if firstCall,   cbrt2(1)=[0.62996052494743658238360530363911d0];  end;
if firstCall,   cbrt2(2)=[0.79370052598409973737585281963615d0];  end;
if firstCall,   cbrt2(3)=[1.0d0];  end;
if firstCall,   cbrt2(4)=[1.25992104989487316476721060727823d0];  end;
if firstCall,   cbrt2(5)=[1.58740105196819947475170563927231d0];  end;
if firstCall,   niter=[0];  end;
firstCall=0;
%***FIRST EXECUTABLE STATEMENT  DCBRT
if( niter==0 )
niter = fix(1.443.*log(-.106.*log(0.1.*real(d1mach(3))))+ 1.0);
end;
%
dcbrtresult = 0.0d0;
if( x==0.0d0 )
csnil=dbstack(1); csnil=csnil(1).name(1)~='@';
if csnil&&~isempty(inputname(1)), assignin('caller','FUntemp',x); evalin('caller',[inputname(1),'=FUntemp;']); end
return;
end;
%
[dumvar1,y,n]=d9upak(abs(x),y,n);
ixpnt = fix(fix(n./3));
irem = fix(n - 3.*ixpnt + 3);
%
% THE APPROXIMATION BELOW IS A GENERALIZED CHEBYSHEV SERIES CONVERTED
% TO POLYNOMIAL FORM.  THE APPROX IS NEARLY BEST IN THE SENSE OF
% RELATIVE ERROR WITH 4.085 DIGITS ACCURACY.
%
z = y;
dcbrtresult = .439581e0 + z.*(.928549e0+z.*(-.512653e0+z.*.144586e0));
%
for iter = 1 : niter;
cbrtsq = dcbrtresult.*dcbrtresult;
dcbrtresult = dcbrtresult +(y-dcbrtresult.*cbrtsq)./(3.0d0.*cbrtsq);
end; iter = fix(niter+1);
%
[dcbrtresult ,dumvar2,ixpnt]=d9pak(cbrt2(irem).*(abs(dcbrtresult).*sign(x)),ixpnt);
%
csnil=dbstack(1); csnil=csnil(1).name(1)~='@';
if csnil&&~isempty(inputname(1)), assignin('caller','FUntemp',x); evalin('caller',[inputname(1),'=FUntemp;']); end
end
%DECK DCDOT