function r = quadroot2(a,flag)
%
% Two ways to compute quadratic roots.
% We assume roots are both real and a(1)~=0.
%
% This code is modified from quadroot.m to reduce 
% overflow/underflow possibilities. 
% It does not always work, however.
%
% flag = 0: unstable way.
% flag = 1: more stable way.
%
% Written by Ming Gu for MA221, Fall 2008
%
% The following step brings a to normal floating point
% range.
a = a/max(abs(a));
%
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