% This example solves a quadratic equation using the standard quadratic
% formula.  The smallest root is completely inaccurate.  A cure is
% illustrated.   The script produces comments as it proceeds.

FlapStartup
dr1 = rand(1);
dr2 = 10^(-FlapRlevel)*rand(1);
db = -(dr1 + dr2);
dc = dr1*dr2;
b = flap(db);
c = flap(dc);
rts = roots([1, b.d, c.d]);
frl = FlapRlevel;

PrintAry([b.d c.d], frl, '\nThe quadratic x^2 + b*x + c, where b and c are\n\n')
disp(' ')
disp('has the roots')
format long e
disp(' ')
disp(rts')
format short e
pause

disc = sqrt(b^2 - 4*c);
PrintAry(disc.d, frl, ...
[   '\nThe following calculations use the quadratic formula'... 
    'to compute the larger root\n\n   disc = sqrt(b^2 - 4*c) ='])
num1 =  (-b + disc);
PrintAry(num1.d, frl, ...
     '\n\n    num1 = (-b + disc) =')
r1 = num1/2;
PrintAry(r1.d, frl, ...
     '\n\n    r1 = num1/2 =')
disp(' ')
disp('Note that the root is accurate up to the rounding level.')
pause

num2 =  (-b - disc);
PrintAry(num2.d, frl, ...
[   '\nOn the other hand, when we calculate the second root, we get' ...)
    '\n\n    num2 =  (-b - disc) ='])
r2 = num2/2;
PrintAry(r2.d, frl, ...
     '\n\n    r2 = num2/2 = ')
disp(' ')
disp('Here the root is entirely inaccurate.')
disp(' ')
disp('The cancellation in computing num2 foretokens the inaccuracy.')
pause

disp('However, if we use the relation c = r1*r2 and compute')
disp('r2 = c/r1, we get')
disp(' ')
r2 = c/r1;
PrintAry(r2.d, FlapRlevel)
disp(' ')
disp('The root is computed to almost full accuracy')