function r = quadroot(a,flag)
%
% Two ways to compute quadratic roots of the equation 
% a(1) x^2 + a(2) x + a(3) = 0.
% We assume roots are both real and a(1)~=0.
%
% The mode flag = 0 corresponds to a naive implementation 
% of the quadratic root formulas so that one of the 
% roots will be computed inaccurately if b^2 >> 4ac. 
% The mode flag = 1 corresponds to a more stable implementation. 
% To see this, let a = (1, 1e10, 1) and run the code
% for both modes to see the difference.
%
% However, if one or both of b^2 and 4ac overflow/underflow, 
% this code could still fail even though the exact roots themselves
% should not overflow/underflow. 
%
% To avoid this problem, try to scale the input data so 
% the largest component in a will not be too large or too small.
%
% flag = 0: unstable way.
% flag = 1: more stable way.
%
% Written by Ming Gu for MA221, Fall 2008
%
r = zeros(2,1);
if (flag == 0)
   temp = sqrt(a(2)*a(2)-4*a(1)*a(3));
   r(1) = (-a(2)+temp)/(2*a(1));
   r(2) = (-a(2)-temp)/(2*a(1));
else
   temp = sqrt(a(2)*a(2)-4*a(1)*a(3));
   if (a(2) > 0)
       r(2) = (-a(2)-temp)/(2*a(1));
       r(1) = a(3)/r(2)/a(1);
   else
       r(1) = (-a(2)+temp)/(2*a(1));
       if (r(1) == 0)
           r(2) = 0;
       else
           r(2) = a(3)/r(1)/a(1);       
       end
   end
end