A worked example: roots of a quadratic

Most people's first attempt at a solution would look something like this.

Here it is in action.

The checks that the supposed roots do at least make the polynomial very small (to around ${10}^{-16}$, the approximate "accuracy" of "double") are reassuring, but do not tell us the accuracy of the answers for the roots that we have obtained. Unfortunately not all runs are quite so convincing.

and

Let us try to get some rough estimates of the accuracy in this program. Firstly, we have seen that "double" is representented in binary to approximately 16 decimal places. Let us simplify this slightly and assume that exactly 16 decimal places are correct. (There are some cases when a double is not quite correct to this accuracy but we ignore these for now.) In the case when $x\ne 0$, to say we have $x$ correct to 16 decimal digits is to say the relative error is at most $5\times {10}^{-16}$. (The relative error of an estimate $\tilde{x}$ to a real value $x$ is the ratio of the absolute error $\left|x-\tilde{x}\right|$ by the actual value $\left|x\right|$. See Errors and precision in programming with floating point.)

Next we assume that $a,b,c$ are all given exactly and are all represented exactly by doubles. In other words there is no experimental error in the input. That is certainly the case with the coefficients in the three examples just given because they are plain integer values that are not too big.

The program computes values $-b$, $\sqrt{b^2-4ac}$, and $2a$. Of these, the most complicated is the second. But if $a,b,c$ are known exactly then the algorithms for $b^2$ and $4ac$ rounded to 16 decimal places should be accurate to 16 decimal places. So the relative error in $b^2$ and $4ac$ will be at most $5\times {10}^{-16}$. One could argue that if $b^2$ and $4ac$ are close together the relative error of the difference is much greater. That's true, but as it happens this issue is of only secondary importance for this problem, and the problems we are seeing are due to something else. We will return to this issue later, but for now will assume that $b^2$ and $4ac$ are rather different so that the severe decrease in relative error (called "catastrophic cancellation") in $b^2-4ac$ does not occur. For a similar reason, taking square roots of very small numbers can introduce large relative errors (because the derivative of $y=\sqrt{x}$ is large near $x=0$) but we again assume that this is not the case and analyse this properly later. Indeed in all three cases above $b^2$ is large compared to $4ac$.

We conclude that, except in these extreme cases, we might expect to get a relative error in $\sqrt{b^2-4ac}$ of approaching 16 decimal places, at least as small as ${10}^{-15}$ say. The same goes for $-b$ and $2a$.

Now look at the $\pm$ in $-b\pm \sqrt{b^2-4ac}$. In one case the signs of $-b$ and $\sqrt{b^2-4ac}$ are the same and in one case they are different. (Which way round this is depends on whether b is positive or negative.) When the signs are the same, the operation will cause no real problem, so

This doesn't quite say how good, but at least if $b^2$ and $4ac$ are not approximately equal we might hope to get close to 15 decimal places of accuracy. But

Looking at the second and third examples we can see this is the case, and one of our roots is probably very poor. It doesn't take much to check that it is the first of the two roots in all cases that is suspicious. Even when, in the last case, the polynomial was evaluated to $1$ not zero, the last root of approximately $-1000000000$ is rather accurate. But in any sort of terms the other root $0$ is probably rather inaccurate.

In any case, we have remaining a nagging worry that we do not know how many digits in the answer are good. For example, consider this run.

Rather than give exact analytic estimates of the errors (as a formula perhaps) we try to estimate the errors introduced at each stage. We will switch between absolute error and relative error as needed, and recall that they are related by

As noted, we believe we can trust the computer to give us full precision for double arithmetic of numbers known exactly. Full precision for double, at least for nonzero numbers that are represented by the format, is a relative error of about $5\times {10}^{-16}$. For simplicity we use this figure. If $x$ is so represented its approximate absolute error is $5\left|x\right|\times {10}^{-16}$.

What we do is take the last program apart and analyse the errors step by step, using the computer to do the calculations for us. I offer my program here and explain it in the notes following.

You see that we have broken the various calculations of long formulas into several steps. (This is generally speaking a very good thing. It takes the computer no longer to do several short steps than one long equivalent step, but the program is easier to read. In fact we may spot values or formulas that can be re-used and hence this may actually make the computation shorter.) For a quantity xxxx the quantity xxxx_err is an estimate of the maximum absolute error we expect. Some of the errors, such as in bsquared and fourac are computed just as indicated above using the presumed precision of double. In fact the computation $\left|x\right|\times 5\times {10}^{-16}$ is used so often it is programmed as a simple function. (This is also good programming style. As well as saving typing it means the number $5\times {10}^{-16}$ can be changed throughout in one single place.)

So you see how bsquared_err and fourac_err are computed. The discriminant is the difference of the values bsquared and fourac so the maximum error is the sum of the errors plus a possible additional round-off error. This gives discrim_err and takes us through the first section of the program.

The estimation of errors in the next section is a bit more tricky. Clearly we do rootdiscrim = sqrt(discrim); but what is the possible error here? Putting $\Delta$ for the discriminant and $\tilde{\Delta }$ for its approximation, we have $\tilde{\Delta }=\Delta +h$ where the error $h$ is in absolute value less than the error discrim_err. We can reasonably assume that $h$ is small compared to $\Delta$ and expanding a series,

So the relative error in $\sqrt{\Delta }$ is approximately $\frac{1}{2}\frac{h}{\Delta }$ and the absolute error is at most $\frac{1}{2}\frac{h}{\sqrt{\Delta }}$ plus a round-off term. (Or, done another way, you can observe that $\sqrt{x+h}\approx \sqrt{x}+h\frac{d\sqrt{x}}{dx}$. From this you get the same estimate of the error.) This explains the line rootdiscrim_err = 0.5 * discrim_err / rootdiscrim + basicerr(rootdiscrim);.

The absolute errors in a sum and difference are at most the sum of the errors plus round-off, hence top1_err and top2_err, and the absolute error in a quotient such as top1/2.0/a of a variable such as top1 by an exact quantity is the quotient of the error by the same quantity plus round-off. This explains the next section of code.

Finally, the number of available digits can be estimated from the relative error (for a nonzero quantity) by noting that if the relative error is within $5\times {10}^{-n}$ then $n$ digits are good. This leads to the logarithmic formula in the function numdigits(x, x_err). The final part of the program prints the answer and its estimated errors and expected number of digits accuracy.

We can run the program with the same test data, as follows.

Using the rather conservative estimates of errors, we see that the first case is accurate to most of the digits shown, and in the second and third cases only one of the roots is good to a reasonable number of digits. The answer in the fourth example is of intermediate quality.

The conclusion of the last section is that the "famous" formulas for the roots of a quadratic are good analytically, but sometimes poor numerically. What can we do about this? One option is to find a better algorithm. For cubics and higher, this is certainly the preferred approach, and we will discuss the bisection method later. The alternative is to look for alternative formulas, and this is what we do here.

In fact there are other useful, but not well-known, formulas for the roots of a quadratic. We have

and

So the roots are $2c/(-b\pm \sqrt{b^2-4ac})$. This formula is less well known because it is less useful analytically, but it turns out to be equally useful numerically. The new program computes the roots by both methods, and comparing the errors choose the best value.

The only part that needs to be explained further is the estimation of the error in the calculation x1alt = 2.0 * c / top2 and the similar one for x2alt. Again, $2c$ is a constant we suppose we know exactly. If we estimate a quantity $C/x$ by computing $C/\left(x+h\right)$ where $h$ is the absolute error in $x$ assumed small relative to $x$ then

So the relative error is $\frac{h}{x}$ and the absolute error is $\frac{hC}{x^2}$. This is the formula used.

Running the program with the same test data, we get

You will note a small change (a possible improvement but hardly important) in the first case, significant improvements in the second and third as the cancellation of digits is avoided in $-b+\sqrt{b^2-4ac}$, and little improvement in the fourth case which remains slightly problematic. The fourth case is less accurate because of cancellation of digits in $b^2-4ac$. To avoid the last problem a different formula or algorithm must be used.

Careful analysis of errors pays off! The additional confidence we have in the answers is significant. Using the computer to make much of the analysis also saves a lot of work. On the other hand, the calculation of the errors takes a significant amount of computer time, certainly much more than double the time required to just compute the formulas themselves. If this work was carried out inside a large loop the extra time required would be very noticable.

In cases when the "check" works to an accuracy about ${10}^{-16}$, this does not mean the answers are correct to the same relative accuracy. In the final case, the roots x1,x2 had p(x1),p(x2) of the order of ${10}^{-16}$ but the roots themselves are unlikely to be this accurate. Be suspicious of simple checks like this!

Note that the final program is by no means perfect. In terms of possible bugs, there are a number of places where division by zero can and sometimes does occur. C++ quietly records "infinity" in such cases which may or may not be appropriate. Also the program does not check that all intermediate values are represented in a reasonable way in the "double" format, or whether such values are stored to sufficient accuracy. Whether this matters for any particular input is not checked. In any case the program is not perfect for aesthetic reasons either: there is repeated code, and it may be better to compute the roots in an array. In any case it should never be necessary to compute all four of x1,x2,x1alt,x2alt and their errors. Two of these suffice, but it is necessary to work out which two.