{{Clarified-review}}Write a program to find the roots of a quadratic equation, i.e., solve the equation <math>ax^2 + bx + c = 0</math>.
Your program must correctly handle non-real roots, but it need not check that <math>a \neq 0</math>.

The problem of solving a quadratic equation is a good example of how dangerous it can be to ignore the peculiarities of floating-point arithmetic.
The obvious way to implement the quadratic formula suffers catastrophic loss of accuracy when one of the roots to be found is much closer to 0 than the other.
In their classic textbook on numeric methods ''[http://www.pdas.com/fmm.htm Computer Methods for Mathematical Computations]'', George Forsythe, Michael Malcolm, and Cleve Moler suggest trying the naive algorithm with <math>a = 1</math>, <math>b = -10^5</math>, and <math>c = 1</math>.
(For double-precision floats, set <math>b = -10^9</math>.)
Consider the following implementation in [[Ada]]:
<lang ada>with Ada.Text_IO;                        use Ada.Text_IO;
with Ada.Numerics.Elementary_Functions;  use Ada.Numerics.Elementary_Functions;

procedure Quadratic_Equation is
   type Roots is array (1..2) of Float;
   function Solve (A, B, C : Float) return Roots is
      SD : constant Float := sqrt (B**2 - 4.0 * A * C);
      AA : constant Float := 2.0 * A;
   begin
      return ((- B + SD) / AA, (- B - SD) / AA);
   end Solve;

   R : constant Roots := Solve (1.0, -10.0E5, 1.0);
begin
   Put_Line ("X1 =" & Float'Image (R (1)) & " X2 =" & Float'Image (R (2)));
end Quadratic_Equation;</lang>
{{out}}
<pre>X1 = 1.00000E+06 X2 = 0.00000E+00</pre>
As we can see, the second root has lost all significant figures. The right answer is that <code>X2</code> is about <math>10^{-6}</math>. The naive method is numerically unstable.

Suggested by Middlebrook (D-OA), a better numerical method: to define two parameters <math> q = \sqrt{a c} / b </math> and <math> f = 1/2 + \sqrt{1 - 4 q^2} /2 </math>

and the two roots of the quardratic are: <math> \frac{-b}{a} f </math> and <math> \frac{-c}{b f} </math>


'''Task''': do it better. This means that given <math>a = 1</math>, <math>b = -10^9</math>, and <math>c = 1</math>, both of the roots your program returns should be greater than <math>10^{-11}</math>. Or, if your language can't do floating-point arithmetic any more precisely than single precision, your program should be able to handle <math>b = -10^6</math>. Either way, show what your program gives as the roots of the quadratic in question. See page 9 of
[https://web.archive.org/web/20080921074325/http://dlc.sun.com/pdf//800-7895/800-7895.pdf "What Every Scientist Should Know About Floating-Point Arithmetic"] for a possible algorithm.
