user
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 {\displaystyle ax^{2}+bx+c=0}}
.
Your program must correctly handle non-real roots, but it need not check that 

a
≠
0
{\displaystyle {\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 {\displaystyle a=1}}
, 

b
=
−

10

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

c
=
1
{\displaystyle {\displaystyle c=1}}
.
(For double-precision floats, set 

b
=
−

10

9
{\displaystyle {\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 {\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 {\displaystyle q={\sqrt {ac}}/b}}
 and 

f
=
1

/
2
+


1
−
4

q

2

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

and the two roots of the quardratic are: 

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

−
c

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

**Task**: do it better. This means that given 

a
=
1
{\displaystyle {\displaystyle a=1}}
, 

b
=
−

10

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

c
=
1
{\displaystyle {\displaystyle c=1}}
, both of the roots your program returns should be greater than 

10

−
11
{\displaystyle {\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 {\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.
assistant
(formerly Perl 6)

**Works with**: Rakudo version 2022.12

*Works with previous versions also but will return slightly less precise results.*

Raku has complex number handling built in.

```
for
[1, 2, 1],
[1, 2, 3],
[1, -2, 1],
[1, 0, -4],
[1, -10⁶, 1]
-> @coefficients {
    printf "Roots for %d, %d, %d\t=> (%s, %s)\n",
    |@coefficients, |quadroots(@coefficients);
}

sub quadroots (*[$a, $b, $c]) {
    ( -$b + $_ ) / (2 × $a),
    ( -$b - $_ ) / (2 × $a) 
    given
    ($b² - 4 × $a × $c ).Complex.sqrt.narrow
}

```

Output:

```
Roots for 1, 2, 1	=> (-1, -1)
Roots for 1, 2, 3	=> (-1+1.4142135623730951i, -1-1.4142135623730951i)
Roots for 1, -2, 1	=> (1, 1)
Roots for 1, 0, -4	=> (2, -2)
Roots for 1, -1000000, 1	=> (999999.999999, 1.00000761449337e-06)
```