% quadratic formula comparison

clear
format long

a=1;
b=2;
c=5;

rootsT=(-b+[1 -1]*sqrt(b^2-4*a*c))/(2*a); % traditional

rootsC=-2*c./(b+[1 -1]*sqrt(b^2-4*a*c)); % Citardauq

residualT=a*rootsT.^2+b*rootsT+c;
residualC=a*rootsC.^2+b*rootsC+c;

fprintf('\na=%g b=%g c=%g\n',a,b,c)
fprintf('Traditional: %0.15g %0.15g residuals=%0.15g %0.15g\n',rootsT(1),rootsT(2),residualT(1),residualT(2)) 
fprintf('Citardauq:   %0.15g %0.15g residuals=%0.15g %0.15g\n\n',rootsC(1),rootsC(2),residualC(1),residualC(2)) 

% looks the same--for large a, we're fine either way
% but consider small a: as a approaches 0, the solution is x=-c/b

a=1e-6; % try small values like 1e-6 or 1e-8
b=2;
c=5;

rootsT=(-b+[1 -1]*sqrt(b^2-4*a*c))/(2*a); % traditional
residualT=a*rootsT.^2+b*rootsT+c;

rootsC=-2*c./(b+[1 -1]*sqrt(b^2-4*a*c)); % Citardauq
residualC=a*rootsC.^2+b*rootsC+c;

fprintf('a=%g b=%g c=%g\n',a,b,c)
fprintf('Traditional: %0.15g %0.15g residuals=%0.15g %0.15g\n',rootsT(1),rootsT(2),residualT(1),residualT(2)) 
fprintf('Citardauq:   %0.15g %0.15g residuals=%0.15g %0.15g\n\n',rootsC(1),rootsC(2),residualC(1),residualC(2)) 


% typically one roots gives tiny residual error but other can be
% problematic (dividing by small difference if a is very small)

% better if a is very small: choose root x of Citardauq with same sign as b 
% (so adds instead of taking difference), then other root is c/(a*x)

root1=-2*c./(b+sign(b)*sqrt(b^2-4*a*c));
root2=c/(a*root1);
roots=[root1 root2];
residual=a*roots.^2+b*roots+c;
fprintf('Citardauq alt:  %0.15g %0.15g residuals=%0.15g %0.15g\n\n',roots(1),roots(2),residual(1),residual(2)) 

% another alternative if a is very small: choose root x of Citardauq with same sign as b 
% (so adds instead of taking difference), then use traditional formula for
% other root

root1=-2*c./(b+sign(b)*sqrt(b^2-4*a*c));
root2=(-b-sign(b)*sqrt(b^2-4*a*c))/(2*a);
roots=[root1 root2];
residual=a*roots.^2+b*roots+c;
fprintf('Citardauq alt2: %0.15g %0.15g residuals=%0.15g %0.15g\n\n',roots(1),roots(2),residual(1),residual(2)) 