The Ins and Outs of Solving Quadratic Equations with
Floating-Point Arithmetic
Frédéric Goualard

To cite this version:
Frédéric Goualard. The Ins and Outs of Solving Quadratic Equations with Floating-Point Arithmetic.
2023. �hal-04116310�

HAL Id: hal-04116310
https://cnrs.hal.science/hal-04116310v1
Preprint submitted on 3 Jun 2023

HAL is a multi-disciplinary open accessL’archive ouverte pluridisciplinaire HAL, estarchive for the deposit and dissemination of sci- destinée au dépôt et à la diffusion de documents
entific research documents, whether they are pub- scientifiques de niveau recherche, publiés ou non,lished or not. The documents may come from émanant des établissements d’enseignement et deteaching and research institutions in France or recherche français ou étrangers, des laboratoires
abroad, or from public or private research centers. publics ou privés.
Copyright

Received: Added at production

Revised: Added at production

Accepted: Added at production

DOI: xxx/xxxx

SURVEY PAPER

The Ins and Outs of Solving Quadratic Equations
with Floating-Point Arithmetic
Frédéric Goualard
1 Nantes Université, École Centrale Nantes,

CNRS, LS2N, UMR 6004, Nantes, France
Correspondence*Frédéric Goualard. LS2N, 2, rue de la
Houssinière, BP 92208, F-44322 NANTESCEDEX 3. Email: frederic.goualard@ls2n.fr

Summary
Solving quadratic equations with radicals on a computer with floating-point arithmetic requires great care to handle correctly all possible parameters. Literature on
the subject glosses over the details, often considered as important but tedious topresent. As a consequence, most implementations are flawed in one way or another.
After having reviewed both the literature and the actual implementations in several
programming languages and applications, we present an algorithm inspired from an
exposition by Pat Sterbenz from 1974, adapted to take advantage of more recent researches in the field, which leads to a robust quadratic equation solver.
KEYWORDS:
Floating-point arithmetic, numerical errors, robustness, quadratic equation

1

INTRODUCTION

How many solutions to the quadratic equation:()
251 + 1𝑝(𝑥) = 𝑥2 + 1 + 2−52 𝑥 +
= 0?
(1)
253
Only one solution, as GNU GSL 2.7, Scilab 2023.0.0, and the Rust 1.69.0 mathematical library would have us believe? Tworeal solutions, as the Boost C++ libraries 1.82.0 and the Racket mathematical library v.8.8[cs] report? Or, maybe, two complex
solutions, as found by Octave 8.2.0, the Python Numpy library 1.22.3, and MATLAB R2023a? What about 2600 𝑝(𝑥) = 0? It
definitely should have the same solutions as 𝑝(𝑥) = 0, yet GNU GSL, the Rust mathematical library, the Racket mathematicallibrary and the Boost C++ libraries now consider that the equation has no real solution at all!
Equation (1) can be solved easily by hand, and there are indeed two real solutions very close to one another. On devising an up
to par quadratic equation solver, Forsythe 1 would say in 1966: “I venture to guess that no more than five quadratic solvers exist
anywhere that meet the general level of the specifications.” More than half a century later, it seems things have not changed, ormaybe for the worst.
The implementers of mathematical libraries are not the only ones to blame, however. Since Forsythe’s 1966 article, aptly titled“How do you solve a quadratic equation?”, several authors have written about the dangers lurking in the implementation of arobust quadratic equation solver, but very few, if any, have described in all the excruciating details the steps to take in order toavoid them. As Forsythe 1 put it, “They [The details] are extremely important to actual computing, but carry less general interest
than ideas just presented.” But, as the saying goes, the devil is in the details. For want of a precise exposition of a reference
algorithm to robustly handle all possible inputs, actual implementations seem all to be flawed in one way or another.
A robust algorithm to solve quadratic equations with radicals is the cornerstone of many more elaborate algorithms, like
Muller’s method 2 for example. As Forsythe quips, “School examples do factor with a frequency bewildering to anyone who
has done mathematics outside of school!”, but real life problems tend to exert a solver to the utmost, requiring it to deliverF. Goualard2

correct solutions for inputs large and small, and to handle gracefully exceptional values such as IEEE 754 3 Not-a-Numbers andinfinities, as well as quadratics degenerating into linear equations.
We present in Section 5 a complete algorithm to solve quadratic equations with radicals, for real solutions only—finding
complex solutions does not bring new problems, anyway. Contrary to our predecessors, no detail was considered too puny to
be exposed. Readers who only wish to implement a robust solver can jump straight to that section after having read Section 2
for the necessary notations. An implementation of the algorithm in the Julia language is given in Program 7; it is also available
on github. Most of the ideas underlying the algorithm were sketched by Sterbenz 4 in 1974 for a computer using hexadecimal
pre-IEEE 754 floating-point arithmetic, in a book that appears to have been long out-of-print. It is rather distressing to observethat many writers have tried since to reimplement, often less correctly, such an algorithm. In Section 3, we take a look at
algorithms presented in the literature since 1974, pointing out some of their flaws in the process; Section 4 does the same for
actual implementations in current software.

2

FLOATING-POINT ARITHMETICThe IEEE 754 standard 3 defines the representation and the properties of floating-point arithmetic on most modern computers.
We only present in this section the elements that are relevant to the understanding of the algorithms presented in Sections 3 to
5. A more in-depth exposition can be found by the interested reader in Muller et al.’s book. 5
The IEEE 754 standard defines a floating-point number 𝑣 (“float”, for short) as a binary1 value of the form:𝑣 = (−1)𝑠 × 𝜎 × 2𝐸 ,where:• 𝑠 ∈ {0, 1} is the sign bit;• 𝜎 = 𝑏0 .𝑏−1 𝑏−2 ⋯ 𝑏1−𝑝 (𝑏𝑖 ∈ {0, 1}), is the significand, with one bit for the integer part and 𝑝 − 1 bits for the fractionalpart 𝑓 ;• 𝐸 ∈ [𝐸min , 𝐸max ] is the exponent.
We note 𝑚 = (−1)𝑠 × 𝜎 the signed significand. Given 𝑣 a floating-point number, let sign (𝑣) be defined as:{−1 if 𝑣 < 0;sign (𝑣) =1 if 𝑣 ⩾ 0.

(2)

A floating-point number is stored as a binary string in memory, with 𝑠, 𝐸, and 𝑓 stored contiguously; the integer part of thesignificand is not stored as it is inferred from the value of the stored exponent. The size of the string determines the precision of
the floating-point format. Two formats of note are:
• The single precision format, stored as a 32-bit string;• The double precision format, stored as a 64-bit string.
To ensure continuity of the computation, the IEEE 754 standard defines special floating-point values:• −∞ and +∞, to obtain an affine extension of the real number system. When a result is too large in magnitude to be
represented in the floating-point format, it is replaced by an infinity. This is an overflow situation;√• When some computation has no meaning over the reals (e.g., −1, or 0∕0), its result is represented by an Not-a-Number(NaN).
Some configurations of the binary string used to store the exponent 𝐸 are reserved to encode these special values.
A normal float is a float whose integer part of the significand is “1”; a subnormal float (with an integer part equal to “0”)has an exponent equal to 𝐸min . An underflow occurs when the result of a computation is so small that it must be coded by asubnormal float. This is a situation we usually try to avoid as subnormals may lead to a loss of precision. For the normal floats, the1 The latest revisions of the standard define both decimal and binary floating-point arithmetic. The binary version is still the one used the most, though.

F. Goualard

3

3.75

3.5

3.25

3.0

2.5

2.75

2.25

0.125
0.25
0.375
0.5
0.625
0.75
0.875

0

1.0
1.125
1.25
1.375
1.5
1.625
1.75
1.875
2.0

-2.0
-1.875-1.75-1.625-1.5-1.375-1.25-1.125
-1.0

-2.5

-2.25

-2.75

-3.0

-3.25

-3.5

-3.75

-

-0.875-0.75-0.625-0.5-0.375-0.25-0.125multiplication and division by a power of 2 is errorless in the absence of underflow and overflow, as it only requires manipulating
their exponent; the division by a power of 2 that underflows may not be errorless as we cannot decrement the exponent past
𝐸min , and the binary point of the significand needs to be right-shifted in order to complete the operation, potentially losing bitsin the process.
Figure 1 shows the position on the real line of the floats from a format with 2 bits for the exponent and 4 bits for the significand.
As can be seen on the figure, the distance ulp (𝑣) from a float 𝑣 = 𝜎 × 2𝐸 to the next is not uniform, and equal to 21−𝑝 × 2𝐸 . The
quantity 21−𝑝 is called the epsilon “𝜀” of the format. It is sometimes defined as the distance from 1 to the smallest float greaterthan 1.

+Figure 1 Floating-point numbers on the real line (2 bits for the exponent and 4 bits for the significand). Normal number are inblue below the line, and subnormals are in magenta above it.

The set of floating-point numbers is not closed for arithmetic operations (the sum of two floats may not be a float, for example).
The results may, therefore, need rounding to be represented as floats. For the arithmetic operations, the IEEE 754 standard
mandates that the floating-point result be the closest float to the real result (correct rounding). For a real value 𝑤, we note f l (𝑤)its floating-point representation. Given two floating-point values 𝑣1 and 𝑣2 , we have:()⎧ ||f l 𝑣1 ⊤𝑣2 − (𝑣1 ⊤𝑣2 )|| ⩽ ulp(f l(𝑣1 ⊤𝑣2 )) , for ⊤ ∈ {+, −, ×, ÷};2|⎪|⎨ | (√ ) √ | ulp(f l(√𝑣 ))1⎪ |f l𝑣1 − 𝑣1 || ⩽
.
⎩ ||
2
|
thanks to the correct rounding. As a shorthand, we will write f l ⟨expr⟩ to express the fact that each step of the evaluation of thenumerical expression “expr” is rounded according to the floating-point format in use. Example:
⟨
⟩
( ( )
( ))
f l 𝑥2 − 𝑦2 ≡ f l f l 𝑥2 − f l 𝑦2 .
When solving quadratic equations, there are two phenomena that will be of importance because they may induce very large
errors in the resulting solutions, or even in the counting of the number of solutions:
Absorption. Summing or subtracting two floats requires equating their exponents. This is done by incrementing the exponentof the smallest value in magnitude, and right-shifting the binary point of its significand accordingly. As such, we may
lose some bits from the significand that cannot be represented anymore. Consider, for example the addition with 4-bit
significands 1.010 × 210 + 1.101 × 22 :1.010 ×210
1.010
×210
←→
+ 1.101 ×22+ 0.00000001101 ×210
The underlined red bits cannot be represented in the registers of the machine and are lost. The actual computation is then
1.010 × 210 + 0.000 × 210 . The larger value absorbed the smaller one;
Cancellation. When subtracting very close values, the only bits kept are the significands’ rightmost ones. If the operands ofthe subtraction originate from a previous computation, these bits may be the result of successive roundings. Consider, forexample the subtraction with 8-bit significands 1.0010111 × 20 − 1.0010010 × 20 , where the underlined red bits are theresults of roundings in previous computations:1.0010111 ×20− 1.0010010 ×20
0.0000101 ×20

F. Goualard

4

After renormalizing the result, we get a value 1.01 × 2−5 whose bits are all the uncertain ones. The subtraction of close
operands does not add new errors but it reveals and magnifies the past ones. The cancellation becomes a catastrophic
cancellation when most or all correct bits disappear in the subtraction.
The code and the examples presented in the next sections consider double precision only, even though the exposition of themain algorithm in Section 5 is done for an arbitrary size. Table 1 synthesizes the various parameters of the single and double
formats.
Table 1 Parameters of single and double precision IEEE 754 floating-point formats.

3

Name

𝑝

𝐸min

𝐸max

𝜀

Single precision
Double precision

24
53

-126
-1022

127
1023

2−23
2−52

A REVIEW OF THE LITERATURE

The method to solve quadratic equations of the form:
𝑎𝑥2 + 𝑏𝑥 + 𝑐 = 0,(𝑎, 𝑏, 𝑐) ∈ ℝ3

has been known for a very long time.According to Smith, 6 the ninth century indian mathematician Srı̄dhara presented it in its
current form in his book Ganita-Sara: First, compute the discriminant Δ:
Δ = 𝑏2 − 4𝑎𝑐(3)Then, depending on the sign of Δ, there are either no real solution, one, or two real solutions:⎧ Δ < 0 ∶ No real solution;⎪ Δ = 0∶ 𝑥 = − 𝑏 ;⎨2𝑎 √⎪ Δ > 0 ∶ 𝑥 = −𝑏± Δ .
1,2
⎩
2𝑎

(4)

Even though the naive algorithm based on Equation (4) was used without modification on early machines, 7 computer scientists
quickly understood that it was not a good idea to use an algorithm devised for real numbers to solve quadratics on computerswith floating-point arithmetic. As summarized by Metropolis (of Monte Carlo method fame) in a 1973 article, 8 there are mainly
three pitfalls:
1. The values 𝑎, 𝑏, and 𝑐 may be so small or so large that an underflow or an overflow arises when computing Δ, even though
the solutions are perfectly representable;
2. When 𝑏2 ≫ |4𝑎𝑐|, the formula for one of the solutions—depending on the sign of 𝑏—suffers from catastrophic
cancellation,2 leading to a very inaccurate result;
3. When 𝑏2 ≈ 4𝑎𝑐, the computation of the discriminant itself suffers from catastrophic cancellation. A large error in Δ maylead to inferring the wrong number of solutions.
Solutions to the first problem were already pretty well understood and practiced early on, even though they were not always
implemented systematically: to compute a complex division—a problem close to the computation of a discriminant—, Smith 9
proposed in 1962 a scaled formula in order to avoid underflows and overflows. Forsythe, 1 again, had presented in 1966 someexamples of equations that would lead to overflows or underflows, and he had suggested some scaling strategy to avoid them.
2 Or “fantastic cancellation,” as Metropolis puts it.

F. Goualard

5

Sadly, as many others after him, he did not go through with it in that article or in any of his later ones, alluding to the fact that
it was important but too wearisome to present.
In 1956, Muller 2 used an alternative formula—amusingly called the “citardauq, 10 ” for “quadratic” in reverse—originatingfrom Fagnano’s XVIIIth century work, 11 to avoid the second problem:√−𝑏 ± Δ
2𝑐
𝑥=
=
√ ,
2𝑎
−𝑏 ∓ Δ√
where “∓” is to be understood as the opposite sign to the one chosen for “±.” Since the sign between −𝑏 and Δ is different inboth formulas, it is possible to avoid catastrophic cancellation by using one formula for a solution and the other formula for theother solution, depending on the sign of 𝑏.
For the third problem, it seems that Kahan had already suggested to compute the discriminant with a precision twice the sizeof the working precision during a lecture at Stanford University 12 in 1966. In 2004, he introduced an algorithm 13 that avoided
the need for a larger precision by simulating it.
The major problems in the implementation of quadratic solvers were cleared up, at least theoretically, very early in the historyof modern computers. What has been lacking ever since is an exposition of all the steps to take in order to obtain a robust solver.
Authors have highlighted the necessity of handling all possible cases: large and small parameters, quadratics degenerating intolinear or constant equations, cancellation problems. . . but, despite our best efforts, we could not find a presentation summarizing
all these pitfalls and presenting their solutions in a clear algorithm. Details are always dismissed as trivial or tedious.
This absence of an accessible reference algorithm has led to some naive or flawed algorithms benefiting from undeserved
attention. Search on the web how to solve a quadratic equation and you will find over and over again the same naive algorithm
implementing the mathematical formulas cited above. Even more troubling, this is also the algorithm presented in many“educational settings”, such as the CASIO manual 14 or the Numworks calculator documentation. 15
Consider such a naive algorithm, as implemented in Program 1.
Program 1: Naive implementation in Julia of a quadratic equation solver.
function naive(a,b,c)a=Float64(a);b=Float64(b);c=Float64(c)delta = b*b - 4*a*c
if delta < 0return (NaN64,NaN64)
elseif delta == 0
return -b/(2*a)
else
x1 = (-b + sqrt(delta))/(2*a)
x2 = (-b - sqrt(delta))/(2*a)
if x1 > x2
return (x2,x1)
else
return (x1,x2)
end
end
end

According to Forsythe, 12 a satisfactory quadratic equation solver should handle correctly the cases where either 𝑎, 𝑏, or 𝑐 arezero, and for each representable solution 𝑥𝑖 , it should compute an approximation 𝑥̂𝑖 of it with a precision such that:|𝑥𝑖 − 𝑥̂𝑖 | 3|| ⩽ 𝜀.
|𝑥𝑖 |2| |In addition, Forsythe requires that extended precision be used only to compute the discriminant, if necessary.
Program 1 fails on all these requirements. If 𝑎 is zero, we get a division by zero. In addition, consider, for example, the equation:3
𝑝1 (𝑥) = 𝑥2 + 227 𝑥 + = 0.
4

F. Goualard6

It has two solutions:√⎧−227 + 254 −3≈ −5.59 × 10−9 ;⎪ 𝑥1 =
2√⎨−227 − 254 −3≈ −134217728.
⎪ 𝑥2 =
2⎩
Catastrophic cancellation prevents to compute 𝑥1 accurately, and we get with Program 1:{𝑥̂1 = −7.450580596923828 × 10−9 ;𝑥̂2 = −134217728.

Scale the equation to 2500 𝑝1 (𝑥) = 0, and Program 1 will return one solution. . . NaN. Scale the equation in the other directionto 2−1000 𝑝(𝑥) = 0, and there is only one solution again: -67108864, this time. Evidently, such lack of scale invariance is not
acceptable.
Some authors have been more diligent than others, however, and there are some works that still deserve some attention. The
rest of this section highlights some of these references in chronological order, pointing out where they are lacking. The Julia
code for all of the methods described below is available in the same Julia package on github as the final procedure in Section 5.

3.1

Nonweiler, 1968

In “Roots of low-order polynomial equations,” 16 Nonweiler presents procedures to solve quartic, cubic, and quadratic equations.
He explicitly states at the beginning of the quadratic solving procedure that it will fail on overflow when 𝑎 is zero “and in othercases”, which does not bode well for its overall robustness.
The polynomial:𝑝(𝑥) = 𝑎𝑥2 + 𝑏𝑥 + 𝑐 = 0
is rewritten as the monic polynomial:{𝑞(𝑥) = 𝑥 − 2𝑏 𝑥 + 𝑐 , with2′′

We can then compute Δ:−𝑏𝑏′ = 2𝑎.
𝑐 ′ = 𝑎𝑐Δ = 𝑏′2 − 𝑐 ′ .
If Δ is strictly negative, there are no real solutions. If Δ = 0, we return 𝑏′ ; otherwise:{√⎧Δ + 𝑏′ if 𝑏′ > 0;√
⎪ 𝑥1 =
𝑏′ − Δ if 𝑏′ ⩽ 0;⎨
⎪𝑥 = 𝑐 .
⎩ 2 𝑥
1The discriminant is evaluated as badly as in the naive procedure and is therefore plagued by the same cancellation problem. On
the other hand, the author avoids cancellation in computing the solutions by using the standard formula discerningly, according
to the sign of 𝑏, for one solution, and by using one of Viète’s formulas: 17{
𝑥1 + 𝑥2 = − 𝑎𝑏
𝑥1 𝑥2

=

𝑐
𝑎

for the other solution.
Apart from the transformation into a monic polynomial, no effort is made towards avoiding overflows and underflows.

3.2

Jenkins, 1975
18

Jenkins presented in 1975 a generic algorithm to find the roots of a polynomial, which is based on a quadratic solver. Thissolver handles separately the degenerate quadratics:
1. If 𝑎 = 0 and 𝑏 = 0, it returns (0, 0);
2. If 𝑎 = 0 and 𝑏 ≠ 0, it returns (− 𝑏𝑐 , 0);
3. If 𝑐 = 0, it returns (− 𝑎𝑏 , 0).

F. Goualard

7

In the general case, it tries to avoid overflows when computing Δ by first computing the quantity 𝑒:{
1 − 2𝑎× 2𝑐𝑏 ,
if |𝑏|⩾ |𝑐|;𝑏2𝑒=𝑏𝑏
× 2|𝑐| − sign(𝑐) ⋅ 𝑎, otherwise.
2√

When 𝑒 ⩾ 0 (there is no real solution otherwise), it computes 2Δ as:{√√𝑒 × |𝑏|,
if |𝑏|⩾ |𝑐|;
Δ
2
2
√
√
=
2
𝑒 × |𝑐|, otherwise.
Lastly, the solutions are computed as:√⎧Δ− 2𝑏 −sign(𝑏)⋅ 2⎪ 𝑥1 =
;
𝑎
⎪
{ 𝑐
⎨
𝑥1
if 𝑥1 ≠ 0,
⎪ 𝑥2 =𝑎⎪0 otherwise.
⎩in order to avoid cancellation.
The formulas used lead to inaccuracies for large or small parameters. Consider for example the equation 2−1073 (𝑥2 −𝑥−1) = 0,
which has the same two solutions as 𝑥2 −𝑥−1 = 0: 𝑥1 ≈ −0.6180339887498948 and 𝑥2 ≈ 1.618033988749895. This algorithm
will compute the solutions −0.5 and 1.5.

3.3

Franklin, 1977

In Chapter 1 of “Fundamental formulas of physics,” 19 edited by Menzel in 1977, Philip Franklin presents the solving of aquadratic equation as the very first formula. He does not consider the case 𝑎 = 0, which is explicitly excluded; in addition,
he only takes care of the possibility of cancellation in computing the solutions by using the citardauq when appropriate in themodified Fagnano’s formulas:√⎧−𝑏−sign(𝑏) 𝑏2 −4𝑎𝑐
𝑥
=
;
⎪ 1
2𝑎
⎨
−2𝑐
√
.
⎪ 𝑥2 =
𝑏+𝑠𝑖𝑔𝑛(𝑏) 𝑏2 −4𝑎𝑐⎩
Overflows and underflows are not handled. The case 𝑏2 ≈ 4𝑎𝑐, which could lead to an inaccurate computation of Δ, is notconsidered either.

3.4

Hamming, 1986

The method presented in “Numerical methods for scientists and engineers,” 20 authored by Hamming in 1986, does not offer
much that is new: Hamming draws the attention to the risk of cancellation when computing one of the solutions when 𝑏2 ≫ |4𝑎𝑐|.
His solutions is to compute the solution that is not affected by cancellation with the usual formula and to use Viète’s formulafor the other, as Nonweiler does:√{−𝑏−sign(𝑏) 𝑏2 −4𝑎𝑐
𝑥1 =
;
2𝑎
𝑥2 = 𝑎𝑥𝑐 .
1

What is interesting is that, though he is aware of all the other problems that can plague a quadratic equation solver, he dismisses
them with typical offhandedness: “This is not the complete answer on how to evaluate the formula; we still need to worry about(1) underflow, (2) overflow, and (3) 𝑏2 − 4𝑎𝑐 < 0, but these are not relevant here3 .”3.5

Young & Gregory, 1988

The procedure proposed by Young and Gregory in 1988 in their book “A survey of numerical mathematics” 21 is particularly wellthought out, but with peculiar features, and it uses floating-point arithmetic from a made-up computer. If one of the parameters
3 Our emphasis.

F. Goualard

8𝑎, 𝑏, or 𝑐 is “sufficiently small”—that is, smaller than some user-defined threshold for each parameter—, it is flushed to zero
before proceeding. Their procedure is quite convoluted, because, as they put it “[i]f we wish to develop a computer program
which can solve [the quadratic] for any given value of [𝑎], [𝑏], and [𝑐], then we must first analyze certain special cases.”
All special cases for either 𝑎, 𝑏, and 𝑐 being zero are taken into account separately. For the general case, the algorithm rewrites
the formulas to avoid overflows and underflows as much as possible. Cancellation in computing the solutions is handled as in
Section 3.4.
Overall, Young and Gregory’s algorithm handles correctly almost of the traps, except for the accurate computation of thediscriminant. The authors are aware of the necessity to compute it with a precision that is double the working one—and write
as much—but they do not offer any algorithm to do so on a computer without the proper format. In addition, being defined fora non IEEE 754 arithmetic, some care needs to be put into adapting the procedure for actual modern computers.

3.6

Baker, 1998

The 1998 article “You could learn a lot from a quadratic: I. Overloading considered harmful,” 22 by Baker, is particularly interesting as it starts by castigating the authors of the “Ada Language Reference Manual” from 1983 for presenting a flawed
procedure (the naive one) for solving quadratic equations. Baker proceeds by recalling the properties of floating-point arithmetic
and by presenting two procedures to avoid underflows and overflows. We will only consider the second one, as the first one isless efficient and precise.
Baker consciously avoids considering the cases where either 𝑎 or 𝑐 are zero. To avoid underflows and overflows, he takes
advantage of the representation of floating-point numbers to scale the parameters. His scaling is not error-free, however.
Given the quadratic polynomial:𝑝(𝑥) = 𝑎𝑥2 + 𝑏𝑥 + 𝑐,
we may express 𝑎 as:
𝑎 = 𝑚𝑎 2𝑟𝑎 2𝑒𝑎 ,with 𝑒𝑎 an even integer (𝑒𝑎 = 2𝑘𝑎 ) and 𝑟𝑎 ∈ {0, 1}.
Let us divide 𝑝(𝑥) by 𝑚𝑎 2𝑟𝑎 :𝑝(𝑥)𝑐𝑏= 2𝑒𝑎 𝑥2 +
𝑝𝑠 (𝑥) =
𝑥+
.
𝑟
𝑟
𝑎
𝑎
𝑚𝑎 2
𝑚𝑎 2
𝑚𝑎 2𝑟𝑎
⏟⏟⏟
⏟⏟⏟
𝑏1

𝑐1

As Baker puts it, “we can divide the equation through [...] without much risk of underflow or overflow”, which is not the samething as no risk at all: if either 𝑏 or 𝑐 are very large, for example, the division may be enough to trigger an overflow.
In the same way as for 𝑎, let 𝑐1 = 𝑚𝑐1 2𝑟𝑐1 2𝑒𝑐1 , with 𝑟𝑐1 ∈ {0, 1} and 𝑒𝑐1 = 2𝑘𝑐1 . Let us make a change of variable:
𝑥 = −𝑦 sign(𝑏1 )2𝑘𝑐1 −𝑘𝑎 .
We get:𝑞(𝑦) = 22𝑘𝑐1 𝑦2 − |𝑏1 |2𝑘𝑐1 −𝑘𝑎 𝑦 + 𝑚𝑐1 2𝑟𝑐1 22𝑘𝑐1 .
Divide 𝑞(𝑦) by 22𝑘𝑐1 :𝑞𝑠 (𝑦) =If 𝑏2 ⩾

√

𝑞(𝑦)
22𝑘𝑐1

= 𝑦2 − 2 |𝑏1 |2−𝑘𝑐1 −𝑘𝑎 −1 𝑦 + 𝑚𝑐1 2𝑟𝑐1 .
⏟⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏟⏟⏟⏟𝑏2|𝑐2 |, we have two solutions:
| |

𝑐2

√
⎧
⎪ 𝑦1 = 𝑏2 + 𝑏22 − 𝑐2.
⎨𝑐2
⎪ 𝑦2 = 𝑦1⎩
Note that 𝑐2 is obtained from a computation. Hence 𝑏22 − 𝑐2 may lead to some cancellation. Baker proposes to compute 𝑦1 as:
√
√
⎛
𝑐2 ⎞
√
⎜
√
𝑏2 ⎟
𝑦1 = 𝑏2 ⎜1 + 1 − ⎟ .
𝑏2 ⎟
⎜
⎝
⎠

F. Goualard

9

Eventually, we can resubstitute 𝑥 for 𝑦 to get the solutions of the original equation:( ){𝑥1 = − sign 𝑏1 2𝑘𝑐1 −𝑘𝑎 𝑦1 ;( ) 𝑘 −𝑘𝑥2 = − sign 𝑏1 2 𝑐1 𝑎 𝑦2 .
Baker avoids cancellation when computing solutions in the same way as Hamming in Section 3.4 and Young & Gregory in
Section 3.5. He does not take any step to accurately compute the discriminant, however, which leads to errors in the numberof solutions reported. Additionnally, his scaling procedure is slightly flawed when considering subnormal parameters: since
the largest exponent (1023, in double precision) is smaller than the absolute value of the smallest exponent (−1074, in double
precision), we may get an overflow when computing the exponent of 𝑏2 , as −𝑘𝑐1 − 𝑘𝑎 − 1 may be larger than 𝐸max . This is, forexample, the case with the equation:()
2−1073 𝑥2 − 𝑥 − 1 = 0,
which has the two solutions:√−1 ± 5𝑥1,2 =.
2
Baker’s procedure will return the two solutions 0 and ∞.

3.7

Higham, 2002

“Accuracy and stability of numerical algorithms,” 23 by Higham, is a great book on computing with floating-point arithmetic.
Concerning quadratic equations, the author has not much to say, however: he warns against the problem of computing thediscriminant accurately, pointing to Kahan’s work on the subject, and he suggests to use the method already seen several times
to avoid cancellation when computing the solutions:√{−𝑏−sign(𝑏) 𝑏2 −4𝑎𝑐
𝑥1 =
;
2𝑎
𝑥2 = 𝑎𝑥𝑐 .
1

As many authors, Higham evokes the problem of underflows and overflows without elaborating, stating that “[t]hese ideas
can be built into a general strategy [...] but the details are nontrivial.”

3.8

Nievergelt, 2003

True to its title, “How (not) to solve quadratic equations,” 24 by Nievergelt, explores various algorithms, solving with radicals
being only of them. In particular, Nievergelt notes that several softwares solve quadratics by computing the eigenvalues of thecompanion matrix, 25 which leads to inaccuracies, something we witnessed with Numpy, MATLAB and Octave for the example
in the introduction of this article.
To avoid cancellation when computing the solutions, Nievergelt suggests to use the modified Fagnano’s formulas:
−𝑐
⎧𝑥 =
√
;
⎪ 1 (𝑏∕2)+sign(𝑏)∗ (𝑏∕2)2 −4𝑎𝑐√⎨(𝑏∕2)+sign(𝑏) (𝑏∕2)2 −4𝑎𝑐
⎪ 𝑥2 =.
−𝑎⎩
The value 𝑏∕2 is used instead of 𝑏, which raises the possibility of some additional rounding error, should the value of 𝑏 be
subnormal.
He also suggests to consider the cases where either 𝑎, 𝑏, or 𝑐 are zero, without elaborating. Lastly, he advises to use a precisionthat is double the working one to avoid inaccuracies in computing the discriminant. Yet again, no method is proposed for systems
that do not offer such a precision.

3.9

Kahan, 2004

In the 2004 article “On the cost of floating-point computation without extra-precise arithmetic,” 13 Kahan presents an algorithm tosimulate double precision in order to compute the discriminant accurately. He also points out the necessity to handle exceptionalparameters and underflows/overflows in a quadratic solver, but without providing details.
A short MATLAB program to compute the solutions of a quadratic equation is provided in the article, without any provision
to handle degenerate equations, or underflows/overflows, as the main focus seems to be on the accurate computation of the

F. Goualard

10

discriminant, for which a very length MATLAB program is also provided. The necessary algorithm is actually quite short but
Kahan had to consider various environments, whose arithmetic had different properties.
In section 5, we will present a short algorithm to compute the discriminant according to Kahan’s method. It is interesting tonote that the properties of this algorithm have only been proved two years later by Boldo and others. 26,27,28

3.10

Einarsson, 2005

Despite its name, the book “Accuracy and reliability in scientific computing,” 29 by Einarsson, seems only concerned, when
considering the solving of a quadratic equation, by the potential cancellation when computing the solutions. His method is thesame as Menzel’s or Nievergelt’s, relying on the citardauq for one of the solutions. As the author says himself, “Even a simple
problem such as computing the roots of a quadratic equation needs great care.” Nevertheless, nothing more is said to that effect.

3.11

Press et al., 2007

The often-cited book “Numerical recipes: the art of scientific computing,” 30 by Press et al., has something to say about many
numerical problems, the solving of quadratic equations being one of them. Sadly, it is disappointingly light on the matter,
simply considering the cancellation problem when computing the solutions. Its solution is the one already seen several times:
use the naive formula to compute the solution unaffected by cancellation and obtain the other solution through one of Viète’s
formulas. Nothing is said about exceptional parameters, or underflows and overflows. The problem of accurately computing thediscriminant is not considered either.

3.12

Dahlquist & Björck, 2008

The book “Numerical Methods in Scientific Computing,” 31 by Dahlquist and Björck, presents the same method as Press et al.,with the same shortcomings. Apart from computing the solutions with the naive formula and one of Viète’s formulas, nothingis said about the computation the discriminant or the possibility of underflows and overflows.

3.13

McNamee & Pan, 2013

The book “Numerical methods for roots of polynomials,” 32 part II, by McNamee and Pan, devotes its section 12.4 to the errors
encountered when solving quadratic equations. The authors warn about the case where 𝑎 is zero or near zero, remarking thatit may happen often with some algorithms such as Muller’s method, which use a quadratic equation solver as one of theirsubroutine. The cancellation in computing solutions is addressed by using modified Fagnano’s formulas, as seen in Section 3.8,
for example.
The authors allude to the problem of underflows and overflows without addressing it in any way. Nothing is said about theaccurate computation of the discriminant and the handling of exceptional parameters.

3.14

Panchekha et al., 2015

The article “Automatically improving accuracy for floating point expressions,” 33 by Panchekha et al. is particular in that it uses
the Herbie tool to rewrite the naive formula in order to limit errors when evaluating it. For the purpose of demonstrating Herbie’s
abilities, the authors only consider the formula:√−𝑏 − 𝑏2 − 4𝑎𝑐𝑥=.
2𝑎
The tool is able to rewrite the expression to avoid both catastrophic cancellation and overflows for the double precision format:⎧ −𝑏+√4𝑎𝑐𝑏2 −4𝑎𝑐
if 𝑏 < 0;⎪2𝑎⎪ −𝑏−√𝑏2 −4𝑎𝑐−𝑏 − 𝑏2 − 4𝑎𝑐⇝⎨
if 0 ⩽ 𝑏 ⩽ 10127 ;2𝑎2𝑎⎪⎪ −𝑏 + 𝑐
if 𝑏 > 10127 .
⎩ 𝑎 𝑏√

F. Goualard

11

Though interesting, the article does not present a full-fledged algorithm to solve quadratic equations and some of the problems
evoked at the beginning of this section are not addressed.

3.15

Mastronardi & van Dooren, 2015

In “Revisiting the stability of computing the roots of a quadratic polynomial,” 34 Mastronardi and van Dooren present a scheme
to solve quadratic equations accurately. They exclude the case 𝑎 = 0, but consider the posssibilities that 𝑏 or 𝑐 be zero.
If 𝑐 = 0, we have:{;𝑥1 = −𝑏𝑎𝑥2 = 0.
If 𝑏 = 0, we have:

{

𝑥1 =

√

−𝑐
;
𝑎

𝑥2 = −𝑥1 .
The authors consider both complex and real roots, so they do not have to handle √differently the case 𝑎𝑐 > 0. On the other
hand, they do not consider the possibility of underflow or overflow when computing −𝑐∕𝑎.
If 𝑎, 𝑏, and 𝑐 are all different from 0, we start by computing the coefficients (1, 𝑏1 , 𝑐1 ) of the monic quadratic polynomial:𝑏𝑐𝑝(𝑥) = 𝑎𝑥2 + 𝑏𝑥 + 𝑐 ⇝ 𝑝𝑚 (𝑥) = 𝑥2 +
𝑥+
,
𝑎
𝑎
⏟⏟⏟
⏟⏟⏟
𝑏1

and we then make a change of variable:
𝑦=

𝑐1

√
−𝑥
, with 𝛼 = sign(𝑏1 ) |𝑐1 |,𝛼and consider the polynomial 𝑞(𝑦) = 𝑝𝑚 (−𝛼𝑦)∕𝛼 2 :𝑞(𝑦) = 𝑦2 − 2𝛽𝑦 + 𝑒,
with:
The roots of 𝑞(𝑦) are then:|𝑏 |𝛽 = √1 ,2 |𝑐1 |{

𝑒 = sign(𝑐1 ).

√𝑦1 = 𝛽 + 𝛽 2 − 𝑒;√𝑦2 = 𝛽 − 𝛽 2 − 𝑒.

They must be computed as follows:

}
√
⎧
2+1
𝛽
𝑦
=
𝛽
+
⎪ 1
if 𝑒 = −1;⎪ 𝑦2 = − 1𝑦1 √}⎪⎨ 𝑦1 = 𝛽 + (𝛽 + 1)(𝛽 − 1)
if 𝑒 = 1 ∧ 𝑏 ⩾ 1;
⎪ 𝑦 =−1𝑦1⎪ 2⎪ No real solution
if 𝑒 = 1 ∧ 𝑏 < 1.
⎩
The roots of the original polynomial are then computed as:{𝑥1 = −𝛼𝑦1 ;𝑥2 = −𝛼𝑦2 .
No special care is taken in computing the discriminant accurately. Note also that the transformations performed are noterror-free and, therefore, introduce more rounding errors in the process.

3.16

Beebe, 2017

The book “Mathematical-Function Computation Handbook: Programming Using the MathCW Portable Software Library,” 35
by Beebe, is a very large book that is a treasure trove of numerical algorithms. The entire chapter 16 is devoted to quadratic
equations.

F. Goualard

12

Beebe’s exposition is one of the few that explicitly states that the solver must handle all inputs, even NaNs, infinities and zeros,and that presents the formulas to use in those cases. He also presents the three major problems when solving a quadratic, notingthat we need extra precision for computing the discriminant. He gives an algorithm to simulate the double precision needed,which is akin to Kahan’s. 13
Given an equation of the form:
𝑎𝑥2 + 𝑏𝑥 + 𝑐 = 0 = 𝑚𝑎 2𝑒𝑎 𝑥2 + 𝑚𝑏 2𝑒𝑏 𝑥 + 𝑚𝑐 2𝑒𝑐 ,the algorithm starts by scaling all the coefficients by dividing them by 𝑒 = max(𝑒𝑎 , 𝑒𝑏 , 𝑒𝑐 ):𝑚𝑎 2𝑒𝑎 𝑥2 + 𝑚𝑏 2𝑒𝑏 𝑥 + 𝑚𝑐 2𝑒𝑐 = 0 ⇝ 𝑚𝑎 2𝑒𝑎 −𝑒 𝑥2 + 𝑚𝑏 2𝑒𝑏 −𝑒 𝑥 + 𝑚𝑐 2𝑒𝑐 −𝑒 = 0.
This reduces the possibilities of underflows and overflows without averting them entirely (The equation 2−1073 (𝑥2 − 𝑥 − 1) = 0
is not correctly solved, for example, and does not give the same answers as 𝑥2 − 𝑥 − 1 = 0). The solutions are then computedwith the modified Fagnano’s formulas.

4

A REVIEW OF THE IMPLEMENTATIONS IN SOFTWARE

An article may dismissively consider some of the problems that can plague a quadratic equation solver, referring the reader toother publications or to “common sense” for details. A piece of software does not have that luxury, and must implement eachand every strategy that is necessary to handle correctly all inputs.
Not all programming languages offer standard facilities to solve quadratic equations, maybe because it is considered so simple
that any programmer can implement his/her own when needed. A look at the information available on websites for programmers,such as StackOverflow, shows that the suggestions for the implementation of a quadratic solver usually range from the downrightflawed to the quite incomplete. Is the code implemented in “standard” libraries better?

4.1

C++ Boost libraries

The C++ Boost libraries are a huge set of code that is oftentimes the antechamber for code that will be included into a future
standard of C++. As such, it is usually of high quality.
Boost implements the boost::math::tools::quadratic_roots() function to solve quadratic equations with radicals.
The comment in the header of the source code for that function refers to a discussion on StackOverflow about a numericallystable method for solving quadratic equations.
The algorithm used considers separately the various cases where either 𝑎, 𝑏, or 𝑐 are zero. Nothing is done about the possibilityof underflow or overflow, however. To avoid cancellation in computing the solutions, it uses the modified Fagnano’s formulas.
The discriminant is accurately computed using Kahan’s algorithm, 13 with an fma instruction when available.
Overall, the algorithm used is quite good, apart for the absence of underflow/overflow management.

4.2

GNU Scientific Library

The GNU Scientific Library is a project with a long history starting in 1996 to provide the C language with a “standard” libraryfor numerical code. It offers the function gsl_poly_solve_quadratic() to solve quadratic equations with radicals.
The algorithm handles separately the case of a quadratic degenerating into a linear equation. Strangely, when 𝑎 and 𝑏 are zero,the function returns that the equation has no solution, without testing the value of 𝑐 (if 𝑐 is zero, it should report that any value
is a solution).
Nothing is done to avoid underflows and overflows. In addition, the discriminant is computed with the textbook formula(Δ = 𝑏2 − 4𝑎𝑐) with no extra precision. The documentation for the function acknowledges that by stating that “[the] errors
may cause a discrete change in the number of roots.”. Cancellation in computing the solutions is avoided by using the modified
Fagnano’s formulas, as the Boost libraries do.

F. Goualard

4.3

13

Numworks calculator

The Numworks calculator is a French scientific calculator extensively used in France high schools. It is also on the list of
permitted calculators for national exams in many other countries. It is an interesting device as it offers the possibility to program
it using the Python language, and its source code is fully available.
There are two ways to solve quadratic equations with the calculator: either use the module polynomial in the Python emulator,
or use the "Equations” application.
The polynomial module implements the naive algorithm from Program 1, plain and simple. It has, therefore, all its flaws.
The “Equations” application is embedded in an environment that provides it with symbolic facilities. As such, some underflows
and overflows may be averted through behind-the-scene scaling of the arguments. This is not always the case, however: for theequation 2−1073 (𝑥2 − 𝑥 − 1) = 0, it finds 0.5 as the only solution, while it has two:√
−1 ± 5.
2
According to the source code, the discriminant is computed without any extra precision, and the two solutions are obtainedwith the naive formulas. As a consequence, the application falls prey to all the traps laid on the path of quadratic equation solvers.
Though a fantastic tool in its own right, with an unusual policy regarding its source code, it is rather sad that a teaching device
such as the Numworks calculator should not offer any robust means to solve a “simple” quadratic equation.

4.4

Racket “math” library

The math module of the Racket language offers the function quadratic-solutions to solve quadratic equations.
The function uses an algorithm devised by Pavel Panchekha in 2021, and presented in an entry on his blog. As the
documentation states, the algorithm tries to handle cancellations and overflows correctly.
Degenerate cases are not handled correctly however, and the function may return incorrect results when either 𝑎, 𝑏, or 𝑐 are
zero (Issue #81, raised by the author). Cancellations in computing the solutions are avoided by using the modified Fagnano’s
formulas.
Panchekha did put some thought in the implementation of the code computing the discriminant to avoid losing accuracy andraising overflows. However, his algorithm relies on modified expressions to avoid overflows and underflows instead of scaling.
As a consequence, some accuracy may be lost for very large and very small parameters (for example, the√function reports thatthe equation 2−1073 (𝑥2 − 𝑥 − 1) = 0 has the two solutions −2∕3 and 1.5, while the true solutions are (−1 ± 5)∕2). There is also
an error in the computation of the discriminant, whereby some accuracy may be lost, misdiagnosing the number of solutionsof an equation (for example, the equation −312499999999𝑥2 + 707106781186𝑥 − 400000000000 = 0, taken from Nievergelt’s
article, 24 has the two very close solutions ≈ 1.131369396027 and ≈ 1.131372303775, while Racket’s algorithm considers thatit has no real solution (Issue #83, raised by the author).

4.5

Rust “roots” module

The roots module of the Rust language offers the function roots::find_roots_quadratic() to solve quadratic equations.
According to its source code, the algorithm used handles separately the linear case. The discriminant is computed with thenaive formula without any extra precision, and no provision is made to avoid underflows and overflows. On the other hand,cancellation is avoided in computing the solutions by using the modified Fagnano’s formulas.

4.6

Scilab

The primary way to solve a quadratic equation in Scilab seems to be with the roots() function, which, by default, computes
the eigenvalues of the companion matrix. There is, however, a report titled “Scilab is not naive” 36 on the Scilab website, which
presents an algorithm to solve quadratics with radicals. That algorithm seems to be used for quadratic equations in roots()when we pass the right parameter.
To avoid cancellation when computing the solutions, the algorithm uses the modified Fagnano’s formulas as seen already in
Section 3.8. The report does not consider the problem of accurately computing the discriminant, and restricts its handling of

F. Goualard

14

overflows and underflows to the computation of the discriminant only. Two different expressions are used to compute 𝑠 =in order to limit the possibility of overflow. Given 𝑏′ = 𝑏∕2:
• If |𝑏′ | > |𝑐|, we compute:{• If |𝑐| > |𝑏′ |, we compute:

{

√

Δ

𝑒 = 1 − 𝑏𝑎′ 𝑏𝑐′ ,√𝑠 = sign (𝑒) × |𝑏′ | |𝑒|;′𝑒 = 𝑏′ 𝑏𝑐 − 𝑎,√ √𝑠 = sign (𝑒) × |𝑐| |𝑒|.

Strangely, that strategy for averting—but not completely eliminating—overflows is neither included in the algorithm presented
in the report, nor present in the code of the C++ function FindQuadraticPolynomialRoots() of the Scilab source code on
github.

5

A COMPLETE AND ROBUST ALGORITHM

We have witnessed in recent articles and current implementations the consequences of the very traps we were warned againstin the literature from the sixties. A robust quadratic equation solver should:
1. Handle exceptional parameters (NaNs and infinities);
2. Handle degenerate cases where either 𝑎, 𝑏 or 𝑐 are zero;
3. Avoid cancellation when computing Δ with 𝑏2 ≈ 4𝑎𝑐;
4. Avoid underflows and overflows when the solutions are representable;
5. Avoid cancellation when computing the solutions when 𝑏2 ≫ |4𝑎𝑐|.
In his book “Floating-point computation” from 1974, Sterbenz had already considered in details four of the five problems for
a computer using hexadecimal pre-IEEE 754 arithmetic. The fifth problem, the accurate computation of Δ, he had dismissed by
advocating to perform the computation with twice the working precision, which may not always be possible. Our presentation
will address this problem by using Kahan’s algorithm from 2004. 13 We will also flesh out some of his remarks about underflows
or overflows affecting intermediate computations.
We will consider the five pitfalls in order, before giving the complete algorithm. Pitfall 5 will be considered together with
Pitfall 4.

5.1

Pitfall 1: Exceptional parameters

Handling exceptional parameters simply requires to detect them at the beginning of the function. Programming languages usuallyprovide methods to determine whether a float is an NaN or an infinity (e.g., isnan() and isinf(), in C and C++). We also
have to decide what we will return in that case. In our Julia implementation, we have chosen to return a different value thanwhen there are no real solutions (Δ < 0).

5.2

Pitfall 2: Degenerate cases

Table 2 shows all possible degenerate cases. When all three parameters are zero, all reals are solution. What should we returnthen? We have decided to return the pair (∞, −∞) in this situation. Since we will always return two solutions sorted in ascendingorder, this pair cannot be mistaken for the case where two solutions are not representable and their computation overflowed.
When there are no solutions, real or complex, we will return one NaN. When there are only two non-real solutions, we will
return a pair of NaNs. This is of course possible because the Julia language is flexible enough to permit it. With a more rigid
language, we would have to use a different convention.

F. Goualard

15

Table 2 Degenerate cases when solving a quadratic equation.
Case

Real solutions

1. 𝑎 = 0 ∧ 𝑏 = 0 ∧ 𝑐 = 0
2. 𝑎 = 0 ∧ 𝑏 = 0 ∧ 𝑐 ≠ 0
3. 𝑎 = 0 ∧ 𝑏 ≠ 0 ∧ 𝑐 = 0
4. 𝑎 = 0 ∧ 𝑏 ≠ 0 ∧ 𝑐 ≠ 0
5. 𝑎 ≠ 0 ∧ 𝑏 = 0 ∧ 𝑐 = 0
6. 𝑎𝑐 > 0 ∧ 𝑏 = 0ℝ∅{0}{− 𝑏𝑐 }
{0}
𝑥∈ℂ
√
√
𝑐
{− − 𝑎 , − 𝑎𝑐 }7. 𝑎𝑐 < 0 ∧ 𝑏 = 0
8. 𝑎 ≠ 0 ∧ 𝑏 ≠ 0 ∧ 𝑐 = 0{− 𝑎𝑏 , 0}

In Cases 4 and 8 in Table 2, the computation of −𝑐∕𝑏 and −𝑏∕𝑎 may underflow or overflow. That would mean that the solution
is not representable with the floating-point format used, and such an event is therefore unavoidable. The solution returned willthen be 0 or ±∞, depending on the event. In a private communication with Gardiner and Metropolis, 37 Householder suggested
to return in such a situation an unevaluated product as a pair of floats; We did not consider that solution as it would havecomplexified the kind of output managed by the function.
Case 7 is different: we may get an underflow or an overflow when computing −𝑐∕𝑎 while the square root would have broughtback the result into the representable domain. To avoid that, we will exploit our knowledge of the representation of IEEE 754
floating-point numbers. We have:√√√
𝑚𝑐 2𝐸𝑐𝑚𝑐 2𝐸𝑐 −𝐸𝑎
𝑐
− = −
−
=
𝑎
𝑚𝑎
𝑚𝑎 2𝐸𝑎Given 𝑀 and 𝐸 two integers (with 𝐸 ∈ {0, 1}) such that:𝐸𝑐 − 𝐸𝑎 = 2𝑀 + 𝐸
We may then write:√

𝑚 2𝐸𝑐 −𝐸𝑎
− 𝑐
=
𝑚𝑎

√

𝑚 22𝑀+𝐸
− 𝑐
= 2𝑀
𝑚𝑎

√
−

𝑚𝑐 2𝐸
𝑚𝑎

(5)

This is an error-free transformation.
The fraction below the square root cannot trigger any underflow or overflow anymore. In order to implement that expression,
the programming language must allow isolating the signed significand and the exponent of a floating-point number, which isthe case for many languages, particularly those that inherit most of their library from C. However, we have to be wary of theactual value of the exponent reported by the dedicated function. For subnormal numbers, the exponent() function in Julia willreport a value smaller than 𝐸min as it always considers the significand as normalized. As a consequence, the value of 𝑀 may
be outside the domain [𝐸min , 𝐸max ], which would lead 2𝑀 to underflow or overflow; in the case of an overflow, for example, if
the square root is less than 1, the result might still be representable. In order to avoid that situation, we set 𝑀 = 𝑀1 + 𝑀2 with
(𝑀1 , 𝑀2 ) ∈ [𝐸min , 𝐸max ]2 , and we compute Expression 5 as:
√
√
⎛
𝐸
𝑚
2
𝑚 2𝐸 ⎞𝑐2𝑀 −= 2𝑀1 ⎜2𝑀2 − 𝑐 ⎟
(6)
⎜
𝑚𝑎
𝑚𝑎 ⎟
⎝
⎠
using the function keep_exponent_in_check() in Program 2.

5.3

Pitfall 3: Cancellation in computing Δ

The computation of the discriminant Δ according to the textbook formula:Δ = 𝑏2 − 4𝑎𝑐

F. Goualard

16

Program 2: Ensuring that an exponent is always in the domain [𝐸min , 𝐸max ].
function keep_exponent_in_check(Kin)
if -1022 <= Kin <= 1023
return (Kin, 0)
elseif Kin < -1022
return (-1022, Kin+1022)
else
return (1023, Kin-1023)
end
end

entails the subtraction of two computed values, which leads to a possibility of catastrophic cancellation. Consider, for example,
our introductory example:
1
𝑥2 + (1 + 2−52 )𝑥 + ( + 2−53 ) = 0
4
We get:1
Δ = (1 + 2−52 )2 − 4( + 2−53 ) = 2−1044( )Unfortunately, 𝑏2 = 1 + 2−104 + 2−51 cannot be represented with double precision, and will be rounded to f l 𝑏2 = 1 + 2−51 .
Hence:
f l ⟨Δ⟩ = (1 + 2−51 ) − (1 + 2−51 ) = 0
That small rounding error leads us to misidentify the number of solutions.
In 2004, Kahan 13 proposed an algorithm to accurately compute the discriminant by taking advantage of the fact that theerror performed when multiplying two floats is itself a representable float. 38 It is therefore possible to collect those errors andto reinject them into the computation. Program 5 presents a version of Kahan’s algorithm that takes advantage of that when no
fma4 instruction is available. That version requires two auxiliary functions:• A function veltkamp_split(x) (Program 3) to separate a float x into two non-overlapping floats xhi and xlo (withx = xhi + xlo);• A function exactmult(x,y,pxy) (Program 4), which computes the error e = x × y − f l (x × y).

Program 3: Splitting a double precision float into two non-overlapping parts.
function veltkamp_split(x)gamma = 134217729*x # 134217729 = 2^((52/2+1))+1
delta = x - gamma
xhi = gamma + delta
xlo = x - xhi
return (xhi,xlo)
end

When an fma instruction is available, a shorter algorithm exists, shown in Program 6. Note also that, when 𝑏2 ≉ 4𝑎𝑐, the⟨⟩straight computation f l 𝑏2 − 4𝑎𝑐 is accurate enough.
̂ the value of the discriminant computed by Kahan’s algorithm, Boldo 27 proved that:
Given Δ,( )|
̂ || ⩽ 2 ulp Δ̂|Δ − Δ||provided no underflow or overflow occurs. Jeannerod et al. 28 refined that bound further.
We will see in the next section that we will always perform the computation of the discriminant in the absence of underflow
and overflow.
4 An fma instruction is an instruction on three values 𝑥, 𝑦, 𝑧 that can perform the operation 𝑥𝑦 + 𝑧 with only one rounding error.

F. Goualard

17

Program 4: Computation of the error in the multiplication of two floats.
function exactmult(x,y,pxy)(xhi,xlo) = veltkamp_split(x)(yhi,ylo) = veltkamp_split(y)t1 = -pxy + xhi*yhi
t2 = t1 + xhi*ylo
t3 = t2 + xlo*yhi
e = t3 + xlo*ylo
return e
end

Program 5: Accurate computation of the discriminant with Kahan’s method 13 without an fma instruction.
function kahan_discriminant(a,b,c)
d = b*b - 4*a*c
if 3*abs(d) >= b*b + 4*a*c # b^2 and 4ac are different enough?
return dend
p = b*bq = 4*a*c
dp = exactmult(b,b,p)dq = exactmult(4*a,c,q)d = (p-q) + (dp-dq)
return d
end

Program 6: Accurate computation of the discriminant with Kahan’s method 13 with an fma instruction.
function kahan_discriminant_fma(a,b,c)
d = b*b - 4*a*c
if 3*abs(d) >= b*b + 4*a*c # b^2 and 4ac are different enough?
return dend
p = b*bdp = fma(b,b,-p)
q = 4*a*c
dq = fma(4*a,c,-q)d = (p-q) + (dp-dq)
return d
end

5.4

Pitfalls 4 & 5: Underflows/overflows and cancellation in computing the solutions

When computing the discriminant, 𝑏2 , 4𝑎𝑐, or 𝑏2 − 4𝑎𝑐, even, when 𝑎𝑐 < 0, may underflow or overflow while the solutions are
perfectly representable. In order to avoid it, we take advantage of the fact that scaling by a factor does not change the solutions.
Starting from:𝑎𝑥2 + 𝑏𝑥 + 𝑐 = 0,
we use, once again, our knowledge of the representation of floats to express it as:
𝑚𝑎 2𝐸𝑎 𝑥2 + 𝑚𝑏 2𝐸𝑏 𝑥 + 𝑚𝑐 2𝐸𝑐 = 0.
The following scaled equation has the same solutions:()
2𝐸𝑎 −2𝐸𝑏 𝑚𝑎 2𝐸𝑎 𝑥2 + 𝑚𝑏 2𝐸𝑏 𝑥 + 𝑚𝑐 2𝐸𝑐 = 0.

F. Goualard18

We perform a change of variable:
𝑥 = 2𝐸𝑏 −𝐸𝑎 𝑦.
We get:()
2𝐸𝑎 −2𝐸𝑏 𝑚𝑎 2𝐸𝑎 𝑥2 + 𝑚𝑏 2𝐸𝑏 𝑥 + 𝑚𝑐 2𝐸𝑐 = 0()⇐⇒ 2𝐸𝑎 −2𝐸𝑏 𝑚𝑎 22𝐸𝑏 −𝐸𝑎 𝑦2 + 𝑚𝑏 22𝐸𝑏 −𝐸𝑎 𝑦 + 𝑚𝑐 2𝐸𝑐 = 0⇐⇒ 𝑚𝑎 𝑦2 + 𝑚𝑏 𝑦 + 𝑚𝑐 2𝐸𝑐 +𝐸𝑎 −2𝐸𝑏 = 0

Let us now consider the new equation:
𝑚𝑎 𝑦2 + 𝑚𝑏 𝑦 + 𝑚𝑐 2𝐸𝑐 +𝐸𝑎 −2𝐸𝑏 = 0,
where only the constant term may be the source of an underflow or overflow. Once we get the solutions for that equation, we
will perform a new change of variable from 𝑦 to 𝑥 to get the solutions sought.
When computing the discriminant, the term 4𝑚𝑎 𝑚𝑐 2𝐸𝑐 +𝐸𝑎 −2𝐸𝑏 is the only one that can underflow or overflow. Let us consider
the three separate cases:• If 4𝑚𝑎 𝑚𝑐 2𝐸𝑐 +𝐸𝑎 −2𝐸𝑏 ∈ [2𝐸min , 2𝐸max +1 ), we have neither underflow nor overflow;• If 4𝑚𝑎 𝑚𝑐 2𝐸𝑐 +𝐸𝑎 −2𝐸𝑏 < 2𝐸min , we have an underflow;• If 4𝑚𝑎 𝑚𝑐 2𝐸𝑐 +𝐸𝑎 −2𝐸𝑏 ⩾ 2𝐸max +1 , we have an overflow.
Knowing that ||𝑚𝑎 || and ||𝑚𝑐 || take their value in [21−𝑝 , 2), we deduce:• 𝐸𝑐 + 𝐸𝑎 − 2𝐸𝑏 ∈ [𝐸min + 2𝑝 − 4, 𝐸max − 3): neither underflow, nor overflow;• 𝐸𝑐 + 𝐸𝑎 − 2𝐸𝑏 < 𝐸min + 2𝑝 − 4: underflow when computing the discriminant;• 𝐸𝑐 + 𝐸𝑎 − 2𝐸𝑏 ⩾ 𝐸max − 3: overflow when computing the discriminant.
However, if we use the function veltkamp_split() from Program 3 to compute the discriminant, we have to take into
account the multiplication by 2⌊𝑝∕2⌋+1 + 1 it performs. To avoid an overflow, we then also require:(2⌊𝑝∕2⌋+1 + 1)𝑚𝑐 2𝐸𝑐 +𝐸𝑎 −2𝐸𝑏 < 2𝐸max +1 .
To simplify, we choose a tighter bound and replace 2⌊𝑝∕2⌋+1 + 1 by 2⌊𝑝∕2⌋+2 . We get:(2⌊𝑝∕2⌋+2 )𝑚𝑐 2𝐸𝑐 +𝐸𝑎 −2𝐸𝑏 < 2𝐸max +1 ,which gives:
𝐸𝑐 + 𝐸𝑎 − 2𝐸𝑏 < 𝐸max − 2 −as our new upper bound to avoid an overflow.

Case 𝐸𝑐 + 𝐸𝑎 − 2𝐸𝑏 ∈ [𝐸min + 2𝑝 − 4, 𝐸max − 2 −The discriminant of the equation:

⌊ ⌋
𝑝
2

⌊𝑝⌋
2

).

𝑚𝑎 𝑦2 + 𝑚𝑏 𝑦 + 𝑚𝑐 2𝐸𝑐 +𝐸𝑎 −2𝐸𝑏 = 0,⏟⏟⏟⏟⏟⏟⏟⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏟
𝑎′

𝑏′

𝑐′

can be computed without any risk of underflow or overflow using Kahan’s algorithm. If it is negative—no real solution—wereturn the pair (NaN,NaN); if it is zero, we return (−𝑏′ ∕(2𝑎′ )) × 2𝐸𝑏 −𝐸𝑎 . The range of 𝑎′ and 𝑏′ ensure that we cannot have any
spurious overflow or underflow from the division.
If the discriminant is positive, we have a risk of cancellation when computing one of the solutions. As seen in the previous
section, we can avoid it by using the modified Fagnano’s formulas:( ′) √
⎧
′
−𝑏
−
sign
𝑏
𝑏′2 − 4𝑎′ 𝑐 ′
⎪𝑦 =
⎪ 1
2𝑎′
⎨
2𝑐 ′
⎪𝑦 = −√
2⎪𝑏′ + sign (𝑏′ ) 𝑏′2 − 4𝑎′ 𝑐 ′⎩F. Goualard

19

Then, we get:{𝑥1 = 𝑦1 2𝐸𝑏 −𝐸𝑎.
𝑥2 = 𝑦2 2𝐸𝑏 −𝐸𝑎
As we have seen for Case 7 in Table 2, 𝐸𝑏 − 𝐸𝑎 may be outside the domain [𝐸min , 𝐸max ]. We use the same workaround as inEquation 6 and use the expressions:{𝑥1 = (𝑦1 2Ω )2Γ
,𝑥2 = (𝑦2 2Ω )2Γ
with Ω + Γ = 𝐸𝑏 − 𝐸𝑎 and (Ω, Γ) ∈ [𝐸min , 𝐸max ]2 . The same must also be done when Δ = 0.

Case 𝐸𝑐 + 𝐸𝑎 − 2𝐸𝑏 < 𝐸min + 2𝑝 − 4.
√
The computation of 4𝑎′ 𝑐 ′ underflows. Consequently, we have Δ ≈ 𝑏′ . We may compute the first solution as:𝑏′𝑦1 = − ′𝑎To get the second solution, we may reuse the scaling as in Equation (5) and express 𝑐 ′ as:𝑐 ′ = 𝑚𝑐 22𝑀+𝐸 = 22𝑀 𝑚𝑐 2𝐸
⏟⏟⏟

(7)

𝑐 ′′

with 𝑀 and 𝐸 two integers such that 𝐸𝑐 + 𝐸𝑎 − 2𝐸𝑏 = 2𝑀 + 𝐸 and 𝐸 ∈ {0, 1}. Using Viète’s formula, we obtain:𝑦2 =Once again, we get our original solutions by rescaling:
{

𝑐 ′′
𝑎′ 𝑦1

𝑥1 = 𝑦1 2𝐸𝑏 −𝐸𝑎
𝑥2 = 𝑦2 22𝑀+𝐸𝑏 −𝐸𝑎

Again, the exponent 𝐸𝑏 − 𝐸𝑎 and 2𝑀 + 𝐸𝑏 − 𝐸𝑎 may be outside [𝐸min , 𝐸max ]. We ensure the change of variable is done without
spurious underflow or overflow by splitting these exponents into two if necessary, as above.
⌊ ⌋Case 𝐸𝑐 + 𝐸𝑎 − 2𝐸𝑏 ⩾ 𝐸max − 2 − 𝑝 .
We cannot solve the equation:2

𝑎′ 𝑦2 + 𝑏′ 𝑦 + 𝑐 ′ = 0,
directly since the discriminant will overflow on computing 4𝑎′ 𝑐 ′ or when calling veltkamp_split() on 𝑐 ′ . If 𝑎𝑐 > 0, we have𝑏′2 − 4𝑎′ 𝑐 ′ < 0 and there are no real solutions. If 𝑎𝑐 < 0, the value of 4𝑎′ 𝑐 ′ is so large that we may ignore the contribution of𝑏′2 and 𝑏′ . Therefore:√√√| 𝑐′ || 𝑐 ′′ |± −4𝑎′ 𝑐 ′
𝑀
|
|
| |
𝑦1,2 ==±=±2||| 𝑎′ |
2𝑎′| 𝑎′ || |′′′′′The values of 𝑐 and 𝑐 ∕𝑎 may be computed without any risk of overflow, which was not necessarily the case for 𝑐 ′ and 𝑐 ′ ∕𝑎′ .
We rescale to get our original solutions:{𝑥1 = 𝑦1 2𝐸𝑏 −𝐸𝑎
𝑥2 = −𝑥1
As previously, the exponent 𝐸𝑏 − 𝐸𝑎 is split into two if necessary to ensure no spurious underflow or overflow.
Program 7 is the complete Julia function that embodies all the formulas from this section. We call it “sterbenz()” as itis largely inspired from Sterbenz’s notes with Kahan’s algorithm used to accurately compute the discriminant. The sign1()function that appears in its code implements the sign () function from Equation (2).

6

TESTS

Kahan 13 proposed to test quadratic equation solvers by considering equations of the form:
𝑀𝐹𝑛 𝑥2 − 2𝑀𝐹𝑛−1 𝑥 + 𝑀𝐹𝑛−2 = 0,

20

F. Goualard

Program 7: A robust quadratic equation solver using Sterbenz’s notes 4 , adapted to double precision IEEE 754 floats, and Kahan’salgorithm 13 to compute the discriminant.
function sterbenz(a,b,c)try
a=Float64(a);b=Float64(b);c=Float64(c)catch eprintln(stderr,"Parameters should be numbers!")
return nothingend
if isnan(a) || isnan(b) || isnan(c) ||isinf(a) || isinf(b) || isinf(c)return (NaN64)
end# Degenerate cases
if a == 0
if b == 0
if c == 0 # a==0, b==0, c==0
return (Inf64, -Inf64)
else# a==0, b==0, c!=0
return (NaN64)
end
else
if c == 0 # a==0, b!=0, c==0
return (0.0)
else# a==0, b!=0, c!=0
return (-c/b)
end
end
else
if b == 0
if c == 0 # a!=0, b==0, c==0
return (0.0)
else# a!=0, b==0, c!=0
if sign1(a) == sign1(c)return (NaN64, NaN64)
else
ea = exponent(a)ec = exponent(c)ecp = ec - ea
a2 = significand(a)dM = ecp & ~1 # dM = floor(ecp/2)*2
M = dM>>1# M = dM/2
E = ecp & 1# E = odd(ecp) ? 1 : 0
c3 = significand(c)*2.0^(E)
S = sqrt(-c3/a2)(M1,M2) = keep_exponent_in_check(M)x1 = (S*2.0^M2)*2.0^M1
return (-x1,x1)
end
end
else
if c == 0 # a!=0, b!=0, c==0
if sign1(a) == sign1(b)return (-b/a, 0.0)
else
return (0.0, -b/a)

end
else
# a!=0, b!=0, c!=0
ea = exponent(a)eb = exponent(b)ec = exponent(c)K = eb - ea
L = ea - 2*eb
ecp = ec + L
a2 = significand(a)b2 = significand(b)if ecp >= -920 && ecp < 995
c2 = significand(c)*2.0^ecp
delta = kahan_discriminant_fma(a2,b2,c2)
if delta < 0return (NaN64, NaN64)
end(K1,K2) = keep_exponent_in_check(K)if delta > 0
y1 = -(2*c2)/(b2+sign1(b)*sqrt(delta))
y2 = -(b2+sign1(b)*sqrt(delta))/(2*a2)
x1 = (y1*2.0^K2)*2.0^K1
x2 = (y2*2.0^K2)*2.0^K1
return (min(x1,x2),max(x1,x2))
end
return ((-b2/(2*a2))*2.0^K2)*2.0^K1
end
dM = ecp & ~1 # dM = floor(ecp/2)*2
M = dM>>1# M = dM/2
E = ecp & 1# E = odd(ecp) ? 1 : 0
c3 = significand(c)*2.0^(E)
S = sqrt(abs(c3/a2))
if ecp < -920
y1 = -b2/a2
y2 = c3/(a2*y1)(dMK1, dMK2) = keep_exponent_in_check(dM+K)(K1,K2) = keep_exponent_in_check(K)x1 = (y1*2.0^K2)*2.0^K1
x2 = (y2*2.0^dMK2)*2.0^(dMK1)
return (min(x1,x2),max(x1,x2))
end# ecp >= 995
if sign1(a) == sign1(c)return (NaN64, NaN64)
else(MK1,MK2) = keep_exponent_in_check(M+K)
x1 = (S*2.0^MK2)*2.0^(MK1)x2 = -x1
return (x2, x1)
end
end
end
end
end

where 𝐹𝑛 is the 𝑛-th element of the Fibonacci sequence (starting from 𝐹0 ), and 𝑀 = ⌊𝑅∕𝐹𝑛 ⌋, with 𝑅 an integer drawn uniformly
at random from [2𝑝−1 , 2𝑝 − 1]. The solutions are:√𝐹𝑛−1 ± (−1)𝑛.
𝐹𝑛
In order to consider only equations with parameters that are representable with double precision floats, we cannot exceed 𝐹76 .
To get real solutions, we also consider even values of 𝑛 only.
Table 3 summarizes the results for all the methods presented in Sections 3 to 5. Out of twenty algorithms, almost half of themeither compute the wrong number of solutions, or compute the wrong solutions. When a method computes only one result, itF. Goualard

21

Table 3 Results on the 38 randomized Fibonacci-based quadratic equations from 𝐹2 to 𝐹76 .
NameWrong Number of solutions

Wrong solutions

0
0
0
0
0
0
0
0
0
0
18
0
0
0
0
0
18
18
0
0

1
0
0
0
0
0
0
16
0
0
0
1
0
0
0
37
0
0
0
19

Baker
Beebe
Boost
Dahlquist/Björck
Einarsson
GSL
HighamJenkins
Kahan
Mastronardi/van Dooren
Naive
Nievergelt
Nonweiler
Numerical recipes
Panchekha (PLDI15)
Panchekha (racket)
Rust
ScilabSterbenz
Young/Gregory

is counted as a wrong solution only if the value does not correspond to any of the correct solutions. In our implementation, a
result 𝑥
̂ is considered correct with respect to a solution 𝑥 if:√()|̂| , |𝑥|
|𝑥
|
(8)
|
|̂ − 𝑥| ⩽ 𝜀 max |𝑥
The randomized Fibonacci-based tests advocated by Kahan 13 are difficult to solve correctly but not so challenging after all.
We have also defined a set of difficult quadratic equations that exercise the algorithms at their weakest points. All these tests,
and more, are present in our Julia QuadraticEquation package available on github. There are five categories:• Fifteen quadratics with no solution, not even complex ones (an NaN or an infinity is among the parameters);• Four quadratics with no real solution;• Seventeen degenerate quadratics (at least one parameter is zero);• One quadratic with one solution;• Eighteen quadratics with two real solutions.
Table 4 presents the results for all categories with the format x/y/z, where “x” is the number of correct results, “y” the numberof cases where the solver found the wrong number of solutions, and “z” the number of cases where the solutions were wrong.
If a solver finds the wrong number of solutions but at least one result is correct, it is counted in the first category of error but
not in the second. A result is considered correct with respect to a solution following the same definition as Equation (8). Notice
that, now, only Function sterbenz() keeps a perfect score.
Lastly, we generated one million quadratics randomly. The parameters are drawn among all finite floating-point numbers—
even subnormals—in such a way that the resulting quadratic has two real solutions. We computed the solutions with arbitrary
precision in order to be able to determine the size of the error for the results from each of our twenty methods. Table 5 reports
the largest relative error overall. Unfortunately, since many methods do not handle underflows and overflows correctly, mostresults are distressingly large (when the result was an infinite value, the error was not considered, however). We may note that

F. Goualard

22

Table 4 Correctness of the various algorithms on challenging quadratics
Name

No sol.

No real sol.

Degenerate

One sol.

Two sol.

Baker
Beebe
Boost
Dahlquist/Björck
Einarsson
GSL
HighamJenkins
Kahan
Mastronardi/van Dooren
Naive
Nievergelt
Nonweiler
Numerical recipes
Panchekha (PLDI15)
Panchekha (racket)
Rust
ScilabSterbenz
Young/Gregory

0/0/15
14/0/1
1/0/140/0/150/0/15
1/0/140/0/150/0/150/0/150/0/15
1/0/140/0/150/0/150/0/150/0/150/0/15
1/0/140/0/1515/0/01/0/14

4/0/0
4/0/0
4/0/0
4/0/0
4/0/0
4/0/0
4/0/0
4/0/0
4/0/0
1/0/3
4/0/0
1/0/3
4/0/0
4/0/0
4/0/0
4/0/0
4/0/0
4/0/0
4/0/0
1/0/3

2/0/15
6/10/1
13/3/2
3/10/11
4/10/5
13/3/2
3/10/11
6/10/1
3/10/7
6/2/9
5/9/11
0/0/17
7/1/9
4/10/5
5/10/6
6/10/6
12/2/4
7/9/1
17/0/014/2/1

0/1/0
0/1/0
0/1/1
0/1/1
0/1/1
0/1/1
0/1/1
0/1/0
0/1/1
0/1/1
0/1/1
0/1/0
1/0/0
0/1/1
0/1/1
0/1/1
0/1/1
1/0/0
1/0/0
0/1/1

16/0/214/0/413/0/5
14/0/414/0/413/0/5
14/0/417/0/1
13/0/5
14/0/4
10/2/7
13/1/5
13/3/414/0/413/0/515/0/3
13/2/4
16/2/1
18/0/013/1/5the algorithm used in the Racket language and the sterbenz() function exhibit both a good worst relative error and solve
correctly all quadratics. This also shows that testing randomly a procedure may not be the best way to ascertain its qualities as
it is still possible to devise tests that show flaws in Racket’s algorithm (See Tables 3 and 4, for example).
To better assess the precision achieved by the sterbenz() procedure, we also tracked the worst error for 200 000 000 random
quadratics. The largest relative error was around 1.76𝜀, a bit more than Forsythe’s wishes, 12 but quite satisfactory nonetheless.

7

CONCLUSION

Compare Program 1 and Program 7: one is 17 lines long and the other, 128 lines (if we count the ancillary functions keep_exponent_in_check() and kahan_discriminant_fma()); one must be used with great care and a restricted setof parameters—lest it return false results—, the other can take any parameter thrown at it and return precise solutions when they
can be represented as floating-point numbers.
Incredibly elaborate numerical software is built from layer upon layer of ever simpler code. The simplest code, such as aquadratic equation solver, should be flawless so that, as Forsythe 1 puts it, “[...]when a quadratic equation occurs in the midstof a complex and imperfectly understood computation, one can be sure that the quadratic equation solver can be relied upon to
do its part well and permit us to concentrate attention on the rest of the computation.”
The devising of Program 7 requires quite a good knowledge of the properties of IEEE 754 floating-point numbers. Should
we expect that from all the programmers who implement numerical algorithms into programming language libraries and applications? Probably not, and this is why those who devise these algorithms in the first place should not leave any detail out whenpresenting them.
It is rather disquieting to compare what has been available in the literature for more than half a century—some of the algorithms
published being quite good, if not perfect—with what is currently implemented in software. It is possible that the format chosen
for the exposition of their algorithm by experts in numerical methods is often not the most suited to appeal to those who have toF. Goualard

23

Table 5 Worst precision for 1 000 000 random quadratics with two real solutions.
Name
Baker
Beebe
Boost
Dahlquist/Björck
Einarsson
GSL
HighamJenkins
Kahan
Mastronardi/van Dooren
Naive
Nievergelt
Nonweiler
Numerical recipes
Panchekha (PLDI15)
Panchekha (racket)
Rust
ScilabSterbenz
Young/Gregory

Wrong solutions

Overflows

Worst relative precision

36
0
0
0
0
0
0
0
0
0
33640
0
79013
0
0
0
33640
33839
0
194605

11
93804
85388
24457424468324486924457426019
295388196961244574243559223242244869178856
0
109
71717
0
30596

14.43 × 𝜀
4.5 × 1015 × 𝜀
1.8 × 10161 × 𝜀
1.8 × 10161 × 𝜀
1.8 × 10161 × 𝜀
3.87 × 10302 × 𝜀
1.8 × 10161 × 𝜀
4.5 × 1015 × 𝜀
7.7 × 10302 × 𝜀
4.4 × 1015 × 𝜀
5.5 × 10296 × 𝜀
4.7 × 10161 × 𝜀
4.5 × 1015 × 𝜀
1.8 × 10161 × 𝜀
3.8 × 10302 × 𝜀
2.29 × 𝜀
3.1 × 10155 × 𝜀
8.4 × 10306 × 𝜀
1.52 × 𝜀
4.5 × 1015 × 𝜀

implement them. This rift has to be mended somehow if we expect the best numerical methods to drive out the bad ones fromnumerical software eventually.
This is not the end of the story, though. There is one problem that Program 7 cannot prevent, and it happens when theparameters of the quadratic equation are approximations themselves. This can occur because they are obtained from a previous
computation, or because they are provided by the user as decimal values or expressions. Consider for example the equation:
√
4
1
2
𝑥 +
𝑥+ =0
5
5
√
√
Its discriminant is clearly zero, and the only solution should be − 1∕5. Unfortunately, the parameters 4∕5 and 1∕5 have to berounded to be represented as floating-point numbers. Therefore, the actual equation that is being solved is not the one the user
intended, and even Program 7 will not get the expected solutions: it will return that there are no real solutions, which is true ofthe equation really considered:
⟨√ ⟩
( )
4
1
𝑥 + fl=0𝑥2 + f l
55
One solution to this problem may be to use interval arithmetic 39 to handle such cases by enclosing the real parameters into
intervals with floating-point bounds.
The code for all the algorithms presented in this article has been assembled into a Julia package, QuadraticEquation, whichis available on github. The package also contains challenging tests for the various solvers.

ACKNOWLEDGMENTS
Discussions with Dr. Christophe Jermann helped shape this article as it currently is. An invitation to present computer arithmetic
pitfalls to high school French mathematics teachers by Pr. Magali Hersant and Dr. Emmanuel Desmontils, co-directors of the

24

F. Goualard

Institute for Research on Teaching Mathematics (IREM des Pays de la Loire), was the initial impetus for the research presented
therein.

References
1. Forsythe GE. How do you solve a quadratic equation?. Technical Report CS-TR-66-40, Computer Science Department;
Stanford University: 1966.
2. Muller DE. A Method for Solving Algebraic Equations Using an Automatic Computer. Mathematical Tables and Other
Aids to Computation 1956; 10(56): 208–215. Publisher: American Mathematical Society.
3. IEEE Standard for Floating-Point Arithmetic. Tech. Rep. IEEE Std 754-2019 (Revision of IEEE 754-2008), Institute ofElectrical and Electronics Engineers; Piscataway, New Jersey, United States: 2019.
4. Sterbenz PH. Floating-point computation. Prentice-Hall . 1974.
5. Muller JM, Brisebarre N, Dinechin F, et al. Handbook of floating-point arithmetic. Birkhäuser . 2010.
6. Smith DE. General survey of the history of elementary mathematics. New York, NY: Dover . 1980.
7. Duncan FG, Huxtable DHR. The DEUCE Alphacode Translator. The Computer Journal 1960; 3(2): 98–107.
8. Metropolis NC. Analyzed Binary Computing. IEEE Transactions on Computers 1973; C-22(6): 573–576.
9. Smith RL. Algorithm 116: Complex division. Communications of the ACM 1962; 5(8): 435.
10. Goff GK. The Citardauq Formula. The Mathematics Teacher 1976; 69(7): 550–551. Publisher: National Council of Teachersof Mathematics.
11. Fagnano GCd. Produzioni matematiche del conte Giulio Carlo di Fagnano, marchese de’Toschi, e di Sant’Onorio, nobile
romano, e patrizio senogagliese . 1750.
12. Forsythe GE. What is a satisfactory quadratic equation solver?. Technical Report CS74, Computer Science Department;
Stanford University: 1967.
13. Kahan W. On the cost of floating-point computation without extra precise arithmetic. http://www.cs.berkeley.edu/~wkahan/
Qdrtcs.pdf; 2004.
14. CASIO . Computing with the scientific calculator. CASIO . 1986.
15. Numworks User Manual 20.0.0. 2023. https://www.numworks.com/resources/manual/.
16. Nonweiler TRF. Roots of low-order polynomial equations. Communications of the ACM 1968; 11(4): 269–270.
17. Vietæ F. Opera Mathematica. Bonaventuræ & Abrahami Elzeviriorum . 1646.
18. Jenkins MA. Algorithm 493: Zeros of a Real Polynomial [C2]. ACM Transactions on Mathematical Software 1975; 1(2):178–189.
19. Franklin P. Basic Mathematical Formulas. In: Menzel DH. , ed. Fundamental formulas of physics. 1. Dover. 1977 (pp.
1–106).
20. Hamming RW. Numerical methods for scientists and engineers. New York: Dover. 2nd ed. 1986.
21. Young DM, Gregory RT. A survey of numerical mathematics. New York: Dover Publications . 1988.
22. Baker HG. You could learn a lot from a quadratic: overloading considered harmful. ACM SIGPLAN Notices 1998; 33(1):30–38.

F. Goualard25

23. Higham NJ. Accuracy and stability of numerical algorithms. Philadelphia: Society for Industrial and Applied Mathematics.
2nd ed. 2002.
24. Nievergelt Y. How (Not) to Solve Quadratic Equations. The College Mathematics Journal 2003; 34(2): 90–104.
25. Edelman A, Murakami H. Polynomial roots from companion matrix eigenvalues. Mathematics of Computation 1995;
64(210): 763–776.
26. Boldo S, Daumas M, Kahan W, Melquiond G. Proof and certification for an accurate discriminant. Presentation at
SCAM’2006: 12th GAMM-IMACS International Symposium on Scientific Computing, Computer Arithmetic and Validated
Numerics; 2006.
27. Boldo S. Kahan’s Algorithm for a Correct Discriminant Computation at Last Formally Proven. IEEE Transactions on
Computers 2009; 58(2): 220–225.
28. Jeannerod CP, Louvet N, Muller JM. Further analysis of Kahan’s algorithm for the accurate computation of 2 × 2
determinants. Mathematics of Computation 2013; 82(284): 2245–2264.
29. Einarsson B. Accuracy and Reliability in Scientific Computing. Society for Industrial and Applied Mathematics . 2005.
30. Press WH, Teukolsky SA, Vetterling WT, Flannery BP. Numerical recipes: the art of scientific computing. Cambridge,
New York: Cambridge University Press. 3rd ed. 2007.
31. Dahlquist G, Björck o. Numerical Methods in Scientific Computing, Volume I.
Mathematics . 2008.

Society for Industrial and Applied

32. McNamee JM, Pan VY. Chapter 12 - Low-Degree Polynomials. In: . 16 of Numerical Methods for Roots of Polynomials Part II. Elsevier. 2013 (pp. 527–556).
33. Panchekha P, Sanchez-Stern A, Wilcox JR, Tatlock Z. Automatically improving accuracy for floating point expressions. In:
PLDI ’15. Association for Computing Machinery; 2015; New York, NY, USA: 1–11.
34. Mastronardi N, Dooren vP. Revisiting the stability of computing the roots of a quadratic polynomial. Electronic Transactions
on Numerical Analysis 2015; 44: 73–82.
35. Beebe NHF. The Mathematical-Function Computation Handbook: Programming Using the MathCW Portable Software
Library. Springer International Publishing . 2017.
36. Scilab . Scilab is not naive. Technical Report, Dassault systèmes; Vélizy-Villacoublay, France: 2010.
semanticscholar.org/paper/Scilab-Is-Not-Naive/c2e7b5c928eb18608d438e3b7a84142080d56164.

https://www.

37. Gardiner V, Metropolis N. A comprehensive approach to computer arithmetic. Research Report LA-4531, Los Alamos
Scientific Laboratory of the University of California; Los Alamos, New Mexico: 1970.
38. Dekker TJ. A floating-point technique for extending the available precision. Numerische Mathematik 1971; 18(3): 224–242.
39. Elishakoff I, Daphnis A. Exact enclosures of roots of interval quadratic equations by Sridhara’s and Fagnano’s or modified
Fagnano’s formulas. Applied Mathematics and Computation 2015; 271: 1024–1037.