Solve the real quadratic equation \(c-2bz+az^2=0\) for real or complex roots.

[Thomas Klemm]

Abstract
This article refers to "Example 6: The Smaller Root of a Quadratic" in the "Advanced Functions Handbook" for the HP-15C (p. 191).
It attempts to explain how the \(\Sigma-\) key can be used to guarantee nine correct significant digits on a 10-digit calculator.
It is recommended that you are familiar with both programs "A" and "B" that can be found in "Example 6 Continued" (pp.205).
Both programs "A" and "B" are also attached to this article as listing and state-file for the Nonpareil emulator.

How to avoid destruction cancellation
There are two cases during the calculation of the solutions where cancellation may happen. One is when a solution is close to 0. This can be avoided by calculating the root with the bigger absolute value first and then use Viëta's formula \(x_1\cdot x_2=\frac{c}{a}\) to calculate the 2nd solution.

Quote:The method uses \(d=b^2-a\cdot c\).

If \(d<0\), then the roots are a complex conjugate pair
\[\frac{b\pm i\sqrt{-d}}{a}\]
If \(d\geq 0\), then the roots are real numbers \(x\) and \(y\) calculated by
\[
\begin{align}
s&=b+\sqrt{d}\operatorname{sgn}{b} \\
x&=\frac{s}{a} \\
y&=\left\{
\begin{matrix}
\frac{c}{s} & \quad \text{if}\ s\ne0 \\
0 & \quad \text{if}\ s=0 \\
\end{matrix}\right. \\
\end{align}
\]
The \(s\) calculation avoids destruction cancellation.

However there's another problem when calculating \(d=b^2-a\cdot c\). We may encounter cancellation when \(b^2\) and \(a\cdot c\) are close.

Quote:The program exploits a tricky property of the \(\Sigma-\) and \(\Sigma+\) keys whereby certain calculations can be carried out to 13 signficant digits before rounded back to 10.

How to calculate \(b^2 - a \cdot c\) with \(\Sigma xy\)

Let's have a look at the back label of the HP-15C:

R₂: \(n\)
R₃: \(\Sigma x\)
R₄: \(\Sigma x^2\)
R₅: \(\Sigma y\)
R₆: \(\Sigma y^2\)
R₇: \(\Sigma xy\)

This code snipet is from program "B":

Code:

047 -    45  8  RCL 8           ; b
048 -    43 11  x²              ; b²
049 -    44  7  STO 7           ; b² → Σxy
050 -    45 25  RCL I           ; a
051 -    45  9  RCL 9           ; c
052 -    43 49  Σ-              ; Σxy = b² - a⋅c
053 -    45  7  RCL 7           ; Σxy

Example
a = 11713
b = 735246
c = 46152709

Code:

b²  = 735246²        = 540586680516 ≈ 5.405866805⋅10¹¹
a⋅c = 11713⋅46152709 = 540586680517 ≈ 5.405866805⋅10¹¹
b² - a⋅c             =           -1 ≠ 0

When using the \(\Sigma-\) key the subtraction is executed with 13 digits. But keep in mind that we can store only 10 significant digits of b²:

Code:

735246
x²
STO 7
11713
ENTER
46152709
Σ-
RCL 7

The result is -17 as we're calculating:

Code:


  540586680500 (only 10 significant digits)
 -540586680517
           -17

Transformation of the coordinate system using Horner's method

We want to shift the coordinate system so that \(x_0\) becomes the new origin.
Horner's method can be used to calculate the coefficients of the polynomial in the new coordinates:

\[
\begin{matrix}
& a & -2b & c \\
x_0: & a & a\cdot x_0-2b & (a\cdot x_0-2b)\cdot x_0 + c \\
x_0: & a & 2a\cdot x_0-2b & \\
x_0: & a & & \\
\end{matrix}
\]

Now we can calculate the coefficients a', b' and c' of the quadratic polynomial in the new coordinates:

\[
\begin{align}
a' & = a \\
-2b' & = 2a\cdot x_0-2b \\
b' & = b - a\cdot x_0 \\
c' & = (a\cdot x_0-2b)\cdot x_0 + c \\
& = c -b\cdot x_0 -(b -a\cdot x_0)\cdot x_0 \\
& = c -b\cdot x_0 -b'\cdot x_0 \\
\end{align}
\]

It turnes out that all the operations are of the form: \(u - v\cdot x_0\)

We want to shift the coordinate system so that b' becomes 0:
\[
\begin{align}
b' & =b-a\cdot x_0 = 0 \\
x_0 & =\frac{b}{a} \\
\end{align}
\]

But in doing so we may encounter cancellation. To avoid it we use only the first 3 significant digits of \(x_0\). This is done by using the rounding function RND with the scientific display SCI 2. Together with the 10 digits of the other multiplicand this ensures that we are within the limits of 13 internal digits. Using Σ- avoids rounding the product to 10 digits before the subtraction is executed.

Thus we shift the coordinate system multiple times in steps of 3 digits. When we finally reach \(b=0\) the calculation becomes: \(d=-a\cdot c\)

This is the section of program "B" where the coefficients of the polynomial in the new coordinates are calculated:

Code:

012 -    45  8  RCL 8           ; b
013 -    44  7  STO 7           ; b → Σxy
014 - 45,10,25  RCL/ I          ; b/a
015 -    43 34  RND             ; x₀
016 -    45 25  RCL I           ; a
017 -    43 49  Σ-              ; Σxy = b - a⋅x₀
018 -    45  9  RCL 9           ; c
019 - 42, 4, 7  x⇄ 7            ; b'
020 -       34  x⇄y             ; x₀
021 -    45  8  RCL 8           ; b
022 -    43 49  Σ-              ; Σxy = c - b⋅x₀
023 -       33  R↓              ; x₀
024 -    43 49  Σ-              ; Σxy = c - b⋅x₀ - b'⋅x₀
025 -    45  7  RCL 7           ; c'

The rest of program "B" (lines 045-082) is similar to how the roots are calculated in program "A".
I hope that with this article I could shed some light on how the "Cadillac" Quadratic Solver works under the hood.