
  function z = MyCubeRoot(V)
% V satisfies 1/8<=V<=1
% Computes the cube root of V with relative error <= 10^-14

% Use linear interpolation to estimate the cube root of a number V
% that satisfies 1/8<=V<=1
L = (3+4*V)/7;
W = L;
H = V/(L*L);;
its=0;
% Newton update : xNew = x - (x^3-V)/(3x^2) = (x +  x  + V/x^2)/3
% So think of x = L = W and H = V/x^2. The boxes always have volume V
% and we get the new L (and W) by taking the average of  the current 
% box dimensions. Then we set the height H = V/L^2
M = max([L W H]);
m = min([L W H]);
while M-m > 10^-14*M
    L = (L+W+H)/3;
    W =  L;
    H = V/(L*W);
    M = max([L W H]);
    m = min([L W H]);
end
z = (L+W+H)/3;

% -1 for initial guesses not satisfying L*W*H = V