function  [x] = cubic (a,b,c,d)

% [x] = cubic (a,b,c,d)
%
%   Gives the roots of the cubic equation
%         ax^3 + bx^2 + cx + d = 0    (a <> 0 !!)
%   by Nickalls's method: R. W. D. Nickalls, ``A New Approach to
%   solving the cubic: Cardan's solution revealed,''
%   The Mathematical Gazette, 77(480)354-359, 1993.
%   dicknickalls@compuserve.com

%  Herman Bruyninckx 10 DEC 1996, 19 MAY 1999
%  Herman.Bruyninckx@mech.kuleuven.ac.be 
%  Dept. Mechanical Eng., Div. PMA, Katholieke Universiteit Leuven, Belgium
%  <http://www.mech.kuleuven.ac.be/~bruyninc>
%
%  Modified by Andrew Stein (University of Michigan) to 
%  handle vectors of coefficients
%
% THIS SOFTWARE COMES WITHOUT GUARANTEE.

sa = size(a); sb = size(b); sc = size(c); sd = size(d);
if any ( sa~=sb | sb~=sc | sc~=sd | sd~=sa)
    error('all vectors must be of equal size')
end
if all(sa~=1)
    error('function only accepts vectors')
end
if any(abs(a)<eps) 
  error('Coefficient of highest power must not be zero!\n'); 
end;

x = NaN * ones(length(a),3);

xN = -b/3./a;
yN = d + xN .* (c + xN .* (b + a.*xN));

two_a    = 2*a;
delta_sq = (b.*b-3.*a.*c)./(9*a.*a);
h_sq     = two_a .* two_a .* delta_sq.^3;
dis      = yN.*yN - h_sq;
pow      = 1/3;

ii = find(dis>=eps);
  % one real root:
  dis_sqrt = sqrt(dis(ii));
  r_p  = yN(ii) - dis_sqrt;
  r_q  = yN(ii) + dis_sqrt;
  p    = -sign(r_p) .* ( sign(r_p).*r_p./two_a(ii) ).^pow;
  q    = -sign(r_q) .* ( sign(r_q).*r_q./two_a(ii) ).^pow;
  x(ii,1) = xN(ii) + p + q;
  x(ii,2) = xN(ii) + p .* exp(2*pi*i/3) + q * exp(-2*pi*i/3);
  x(ii,3) = conj(x(ii,2));
ii = find(dis < -eps);
  % three distinct real roots:
  theta = acos(-yN(ii)./sqrt(h_sq(ii)))/3;
  delta = sqrt(delta_sq(ii));
  two_d = 2*delta;
  twop3 = 2*pi/3;
  x(ii,1) = xN(ii) + two_d.*cos(theta);
  x(ii,2) = xN(ii) + two_d.*cos(twop3-theta);
  x(ii,3) = xN(ii) + two_d.*cos(twop3+theta);
ii = find(dis<eps & dis >= -eps);
  % three real roots (two or three equal):
  delta = (yN(ii)./two_a(ii)).^pow;
  x(ii,1) = xN(ii) + delta; 
  x(ii,2) = xN(ii) + delta; 
  x(ii,3) = xN(ii) - 2*delta;