[cbrtresult,x]=cbrt(x);

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

if isempty(cbrtresult), cbrtresult=0; end;
if isempty(cbrt2), cbrt2=zeros(1,5); end;
if isempty(cbrtsq), cbrtsq=0; end;
if isempty(y), y=0; 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;
%***BEGIN PROLOGUE  CBRT
%***PURPOSE  Compute the cube root.
%***LIBRARY   SLATEC (FNLIB)
%***CATEGORY  C2
%***TYPE      SINGLE PRECISION (CBRT-S, DCBRT-D, CCBRT-C)
%***KEYWORDS  CUBE ROOT, ELEMENTARY FUNCTIONS, FNLIB, ROOTS
%***AUTHOR  Fullerton, W., (LANL)
%***DESCRIPTION
%
% CBRT(X) calculates the cube root of X.
%
%***REFERENCES  (NONE)
%***ROUTINES CALLED  R1MACH, R9PAK, R9UPAK
%***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  CBRT
if firstCall,   cbrt2(1)=[0.62996052494743658e0];  end;
if firstCall,   cbrt2(2)=[0.79370052598409974e0];  end;
if firstCall,   cbrt2(3)=[1.0e0];  end;
if firstCall,   cbrt2(4)=[1.25992104989487316e0];  end;
if firstCall,   cbrt2(5)=[1.58740105196819947e0];  end;
if firstCall,   niter=[0];  end;
firstCall=0;
%***FIRST EXECUTABLE STATEMENT  CBRT
if( niter==0 )
niter = fix(1.443.*log(-.106.*log(0.1.*r1mach(3))) + 1.);
end;
%
cbrtresult = 0.0;
if( x==0. )
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]=r9upak(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.
%
cbrtresult = .439581e0 + y.*(.928549e0+y.*(-.512653e0+y.*.144586e0));
%
for iter = 1 : niter;
cbrtsq = cbrtresult.*cbrtresult;
cbrtresult = cbrtresult +(y-cbrtresult.*cbrtsq)./(3.0.*cbrtsq);
end; iter = fix(niter+1);
%
[cbrtresult ,dumvar2,ixpnt]=r9pak(cbrt2(irem).*(abs(cbrtresult).*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 CBUNI