This task has been clarified. Its programming examples are in need of review to ensure that they still fit the requirements of the task.
Write a program to find the roots of a quadratic equation, i.e., solve the equation 



a

x

2


+
b
x
+
c
=
0


{\displaystyle ax^{2}+bx+c=0}

.
Your program must correctly handle non-real roots, but it need not check that 



a
≠
0


{\displaystyle a\neq 0}

.

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 Computer Methods for Mathematical Computations, George Forsythe, Michael Malcolm, and Cleve Moler suggest trying the naive algorithm with 



a
=
1


{\displaystyle a=1}

, 



b
=
−

10

5




{\displaystyle b=-10^{5}}

, and 



c
=
1


{\displaystyle c=1}

.
(For double-precision floats, set 



b
=
−

10

9




{\displaystyle b=-10^{9}}

.)
Consider the following implementation in 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;

Output:
X1 = 1.00000E+06 X2 = 0.00000E+00
As we can see, the second root has lost all significant figures. The right answer is that X2 is about 




10

−
6




{\displaystyle 10^{-6}}

. The naive method is numerically unstable.

Suggested by Middlebrook (D-OA), a better numerical method: to define two parameters 



q
=


a
c



/

b


{\displaystyle q={\sqrt {ac}}/b}

 and 



f
=
1

/

2
+


1
−
4

q

2





/

2


{\displaystyle f=1/2+{\sqrt {1-4q^{2}}}/2}



and the two roots of the quardratic are: 






−
b

a


f


{\displaystyle {\frac {-b}{a}}f}

 and 






−
c


b
f





{\displaystyle {\frac {-c}{bf}}}




Task: do it better. This means that given 



a
=
1


{\displaystyle a=1}

, 



b
=
−

10

9




{\displaystyle b=-10^{9}}

, and 



c
=
1


{\displaystyle c=1}

, both of the roots your program returns should be greater than 




10

−
11




{\displaystyle 10^{-11}}

. Or, if your language can't do floating-point arithmetic any more precisely than single precision, your program should be able to handle 



b
=
−

10

6




{\displaystyle b=-10^{6}}

. Either way, show what your program gives as the roots of the quadratic in question. See page 9 of
"What Every Scientist Should Know About Floating-Point Arithmetic" for a possible algorithm.

/*REXX program finds the  roots  (which may be complex)  of a quadratic function.       */
parse arg  a b c .                               /*obtain the specified arguments: A B C*/
call quad  a,b,c                                 /*solve quadratic function via the sub.*/
r1= r1/1;  r2= r2/1;   a= a/1;  b= b/1;  c= c/1  /*normalize numbers to a new precision.*/
if r1j\=0  then r1=r1||left('+',r1j>0)(r1j/1)"i" /*Imaginary part? Handle complex number*/
if r2j\=0  then r2=r2||left('+',r2j>0)(r2j/1)"i" /*   "        "      "       "      "  */
              say '    a ='   a                  /*display the normalized value of   A. */
              say '    b ='   b                  /*   "     "       "       "    "   B. */
              say '    c ='   c                  /*   "     "       "       "    "   C. */
      say;    say 'root1 ='   r1                 /*   "     "       "       "   1st root*/
              say 'root2 ='   r2                 /*   "     "       "       "   2nd root*/
exit 0                                           /*stick a fork in it,  we're all done. */
/*──────────────────────────────────────────────────────────────────────────────────────*/
quad: parse arg aa,bb,cc;     numeric digits 200 /*obtain 3 args; use enough dec. digits*/
      $= sqrt(bb**2-4*aa*cc);       L= length($) /*compute  SQRT (which may be complex).*/
      r= 1 /(aa+aa);   ?= right($, 1)=='i'       /*compute reciprocal of 2*aa;  Complex?*/
      if ?  then do;  r1= -bb   *r;   r2=r1;          r1j= left($,L-1)*r;   r2j=-r1j;  end
            else do;  r1=(-bb+$)*r;   r2=(-bb-$)*r;   r1j= 0;               r2j= 0;    end
      return
/*──────────────────────────────────────────────────────────────────────────────────────*/
sqrt: procedure; parse arg x 1 ox; if x=0  then return 0; d= digits(); m.= 9; numeric form
      numeric digits 9; h= d+6; x=abs(x); parse value format(x,2,1,,0) 'E0' with g 'E' _ .
      g=g*.5'e'_%2;   do j=0  while h>9;      m.j=h;              h=h%2+1;       end /*j*/
                      do k=j+5  to 0  by -1;  numeric digits m.k; g=(g+x/g)*.5;  end /*k*/
      numeric digits d;         return (g/1)left('i', ox<0)     /*make complex if OX<0. */