function c = cuberoot(a1,a2,a3)
%

Q = (a1*a1 - 3*a2)/9;
R = (2*a1^3 - 9*a1*a2 + 27*a3)/54;

x1 = 0;  x2 = 0;  x3 = 0;
if (Q^3-R^2) >= 0
  theta = acos(R/sqrt(Q));
  x1 = -2*sqrt(Q)*cos(theta/3) - a1/3;
  x2 = -2*sqrt(Q)*cos((theta+2*pi)/3) - a1/3;
  x3 = -2*sqrt(Q)*cos((theta+4*pi)/3) - a1/3;
end

if (R^2-Q^3) > 0
  x1 = -sign(R)*((sqrt(R^2-Q^3) + abs(R))^(1/3) + ...
                 Q/(sqrt(R^2-Q^3) + abs(R))^(1/3)) - a1/3;
end

c = [x1 x2 x3]';
