(15C) Bairstow's Method

[Thomas Klemm]

Bairstow's Method

Links

Wikipedia
MathWorld

Algorithm

We are looking for a quadratic factor T(x) of the polynomial P(x). In general there will be a linear remainder R(x) after the division:

\(P(x) = Q(x) \cdot T(x) + R(x)\)

\(P(x) = a_n x^n + a_{n-1} x^{n-1} + \cdots + a_2 x^2 + a_1 x + a_0\)

\(T(x) = x^2 + p x + q\)

\(Q(x) = b_n x^{n-2} + b_{n-1} x^{n-3} + \cdots + b_4 x^2 + b_3 x + b_2\)

\(R(x) = b_1 ( x + p ) + b_0\)

If we set \(b_{n+1}\equiv b_{n+2}\equiv 0\), then \(\forall k\leq n\):

\(b_k = a_k - b_{k+1} p - b_{k+2} q\)

The two coefficients of the remainder R(x) depend on p and q:

\(b_0 = b_0(p, q)\)

\(b_1 = b_1(p, q)\)

To make T(x) a factor of P(x) the remainder R(x) must vanish. Let's assume that we already have p and q close to the solution. We want to know how to change these values to get a better solution. Therefore we set the first order Taylor series expansion to 0:

\(b_0(p + \delta p, q + \delta q) \approx b_0(p, q) + \frac{\partial b_0}{\partial p} \delta p + \frac{\partial b_0}{\partial q} \delta q = 0\)

\(b_1(p + \delta p, q + \delta q) \approx b_1(p, q) + \frac{\partial b_1}{\partial p} \delta p + \frac{\partial b_1}{\partial q} \delta q = 0 \)

It turns out that the partial derivatives of \(b_k\) with respect to p and q can be calculated with a similar recurrence formula:

\(c_{k+1} = - \frac{\partial b_k}{\partial p}\)

\(c_{k+2} = - \frac{\partial b_k}{\partial q}\)

\(c_k = b_k - c_{k+1} p - c_{k+2} q\)

Thus we have to solve:

\(\begin{matrix}
c_1 \delta p + c_2 \delta q &= b_0 \\
c_2 \delta p + c_3 \delta q &= b_1
\end{matrix}\)


Using Linear Regression

We have to solve a 2-dimensional linear equation:

\(\begin{bmatrix}
c_1 & c_2 \\
c_2 & c_3
\end{bmatrix}\cdot\begin{bmatrix}
\delta p \\
\delta q
\end{bmatrix}=\begin{bmatrix}
b_0 \\
b_1
\end{bmatrix}\)

But since the matrix is symmetric we can use the built-in function L.R. for linear regression.

We know that the following equation is solved for A and B:

\(\begin{bmatrix}
\sum x^2 & \sum x \\
\sum x & n
\end{bmatrix}\cdot\begin{bmatrix}
A \\
B
\end{bmatrix}=\begin{bmatrix}
\sum xy \\
\sum y
\end{bmatrix}\)

Cramer's Rule gives the solutions:

\(A=\frac{\begin{vmatrix}
\sum xy & \sum x \\
\sum y & n
\end{vmatrix}}{\begin{vmatrix}
\sum x^2 & \sum x \\
\sum x & n
\end{vmatrix}}=\frac{n\sum xy-\sum x\sum y}{n\sum x^2-(\sum x)^2}\)

\(B=\frac{\begin{vmatrix}
\sum x^2 & \sum xy \\
\sum x & \sum y
\end{vmatrix}}{\begin{vmatrix}
\sum x^2 & \sum x \\
\sum x & n
\end{vmatrix}}=\frac{\sum y\sum x^2-\sum x\sum xy}{n\sum x^2-(\sum x)^2}\)

For this to work we have to fill the registers R₂ - R₇ with the corresponding values of the equation:

\(\begin{matrix}
n & \rightarrow & R_2 \\
\sum x & \rightarrow & R_3 \\
\sum x^2 & \rightarrow & R_4 \\
\sum y & \rightarrow & R_5 \\
\sum y^2 & \rightarrow & R_6 \\
\sum xy & \rightarrow & R_7
\end{matrix}\)

This mapping leaves us with the following linear equation:

\(\begin{bmatrix}
R_4 & R_3 \\
R_3 & R_2
\end{bmatrix}\cdot\begin{bmatrix}
Y \\
X
\end{bmatrix}=\begin{bmatrix}
R_7 \\
R_5
\end{bmatrix}\)

Example:

\(\begin{bmatrix}
2 & -3 \\
-3 & 5
\end{bmatrix}\cdot\begin{bmatrix}
2 \\
-3
\end{bmatrix}=\begin{bmatrix}
13 \\
-21
\end{bmatrix}\)

CLEAR Σ
2 STO 4
-3 STO 3
5 STO 2
13 STO 7
-21 STO 5
L.R.

Y: 2
X: -3


Registers

\(\begin{matrix}
R_0 & p \\
R_1 & q \\
R_2 & c = c_j \\
R_3 & {c}′ = c_{j+1} \\
R_4 & {c}'' = c_{j+2} \\
R_5 & b = b_i \\
R_6 & \\
R_7 & {b}′ = b_{i+1} \\
R_8 & \text{index} = 9.fff \\
R_9 & a_n \\
R_{.0} & a_{n-1} \\
R_{.1} & a_{n-2} \\
R_{.2} & \cdots \\
R_{.3} & \cdots \\
\end{matrix}\)

Program

Code:

001 {    42 21 11 } f LBL A
002 {       42 32 } f ∑
003 {       45  8 } RCL 8
004 {       44 25 } STO I
005 {    42 21  0 } f LBL 0
006 {       45  5 } RCL 5
007 {       45  3 } RCL 3
008 {       44  4 } STO 4
009 {    45 20  1 } RCL × 1
010 {          30 } −
011 {       45  2 } RCL 2
012 {       44  3 } STO 3
013 {    45 20  0 } RCL × 0
014 {          30 } −
015 {       44  2 } STO 2
016 {       45 24 } RCL (i)
017 {       45  7 } RCL 7
018 {    45 20  1 } RCL × 1
019 {          30 } −
020 {       45  5 } RCL 5
021 {       44  7 } STO 7
022 {    45 20  0 } RCL × 0
023 {          30 } −
024 {       44  5 } STO 5
025 {    42  6 25 } f ISG I
026 {       22  0 } GTO 0
027 {       42 49 } f L.R.
028 {    44 40  0 } STO + 0
029 {          34 } x↔y
030 {    44 40  1 } STO + 1
031 {       43  1 } g →P
032 {       43 34 } g RND
033 {    43 30  0 } g TEST x≠0
034 {       22 11 } GTO A
035 {       45  8 } RCL 8
036 {           2 } 2
037 {          26 } EEX
038 {           3 } 3
039 {          16 } CHS
040 {          30 } −
041 {       44  8 } STO 8
042 {       44 25 } STO I
043 {           0 } 0
044 {          36 } ENTER
045 {    42 21  1 } f LBL 1
046 {       45 24 } RCL (i)
047 {       45  1 } RCL 1
048 {       43 33 } g R⬆
049 {          20 } ×
050 {          30 } −
051 {       45  0 } RCL 0
052 {       43 33 } g R⬆
053 {          20 } ×
054 {          30 } −
055 {       44 24 } STO (i)
056 {    42  6 25 } f ISG I
057 {       22  1 } GTO 1
058 {       43 32 } g RTN

Description

Initialization k = n
Lines 002 - 004
\((b, {b}′, c, {c}′, {c}'') \leftarrow (0, 0, 0, 0, 0)\)
\(I \leftarrow\) index

Partial Derivatives
Lines 005 - 026
At the end of the loop \(b = b_0\) but \(c = c_1\)!
That's why the second division is performed before the first.

Iteration k+1 → k
Lines 006 - 015
\((c, {c}′, {c}'') \leftarrow (b - c p - {c}′ q, c, {c}′)\)

Lines 016 - 024
\((b, {b}′) \leftarrow (a_k - b p - {b}′ q, b)\)

Solve the linear equation
Line 027
The trick using linear regression only works since the matrix is symmetric.
\(\begin{bmatrix}
c & {c}' \\
{c}' & {c}''
\end{bmatrix}\cdot\begin{bmatrix}
\delta p \\
\delta q
\end{bmatrix}=\begin{bmatrix}
b \\
{b}'
\end{bmatrix}\)

Find improved values for p and q
Lines 028 - 030
\(p \leftarrow p + \delta p\)
\(q \leftarrow q + \delta q\)

Stop Criterion
Lines 031 - 033
\(\left | (\delta p, \delta q) \right |< \varepsilon\)
\(\varepsilon = 10^{-n}\) if display is set to FIX n

Polynomial Division
Lines 035 - 057
We have to repeat the division here since we didn't keep \(b_k\) in lines 016 - 024.


Quadratic Solver

This program solves the quadratic equation: \(T(x)=x^2+px+q=0\)

Code:

059 {    42 21 12 } f LBL B
060 {       45  0 } RCL 0
061 {           2 } 2
062 {          16 } CHS
063 {          10 } ÷
064 {          36 } ENTER
065 {          36 } ENTER
066 {       43 11 } g x²
067 {    45 30  1 } RCL − 1
068 {          11 } √x̅
069 {          30 } −
070 {          34 } x↔y
071 {       43 36 } g LSTx
072 {          40 } +
073 {       43 32 } g RTN

Just be aware that this program can't find complex roots. Instead an Error 0 will be displayed.
However it's easy to find the complex solutions. Just use:

CHS
\(\sqrt{x}\)

The solutions then are: \(Y \pm iX\)


Example

\(P(x)=2x^5-9x^4+15x^3+65x^2-267x+234=0\)

Insert coefficients

Code:

2      STO 9
-9     STO .0
15     STO .1
65     STO .2
-267   STO .3
234    STO .4

Initialization

Code:

9.014  STO 8
1      STO 0
       STO 1

Alternatively use:

Code:

MATRIX 1

Run program

Code:

GSB A        -52.0000

RCL 0          1.5000
RCL 1         -4.5000

RCL 9          2.0000
RCL .0       -12.0000
RCL .1        42.0000
RCL .2       -52.0000

Code:

GSB B          1.5000
x<>y          -3.0000

Conclusion

\(2x^5-9x^4+15x^3+65x^2-267x+234=\)
\((x^2+1.5x-4.5)(2x^3-12x^2+42x-52)\)

Solutions

For \(x^2+1.5x-4.5=0\):
\(x_1=1.5\)
\(x_2=-3\)

Initialize guess

Code:

1      STO 0
       STO 1

Alternatively use:

Code:

MATRIX 1

Run program again

Code:

GSB A         -4.0000

RCL 0         -4.0000
RCL 1         13.0000

RCL 9          2.0000
RCL .0        -4.0000

Code:

GSB B          Error 0
←             -9.0000
CHS            9.0000
√x             3.0000
x<>y           2.0000

Code:

RCL .0        -4.0000
RCL 9          2.0000
÷             -2.0000
CHS            2.0000

Conclusion

\(2x^3-12x^2+42x-52=\)
\((x^2-4x+13)(2x-4)\)

Solutions

For \(x^2-4x+13=0\):
\(x_3=2+3i\)
\(x_4=2-3i\)

For \(2x-4=0\):
\(x_5=2\)

Summary

Factors

\(2x^5-9x^4+15x^3+65x^2-267x+234=\)
\((x^2+1.5x-4.5)(x^2-4x+13)(2x-4)=\)
\((x-1.5)(x+3)(x^2-4x+13)2(x-2)=\)
\((2x-3)(x+3)(x^2-4x+13)(x-2)\)

Solutions

\(x_1=1.5\)
\(x_2=-3\)
\(x_3=2+3i\)
\(x_4=2-3i\)
\(x_5=2\)

[rprosperi]

Thomas, hello!

It's good to see you're re-engaged here, with predictably impressive posts.

I have some Fortran stuff for Bairstow's method from the 70's in college, I've not seen much of this since then. I hope to take a closer look this weekend and maybe even find the old Fortran listings.

[Thomas Klemm]

It is primarily a copy of the corresponding page of its sibling HP-11C.
Some of the registers had to be adjusted and due to the RCL arithmetic even a few steps could be saved.
I used Torsten's excellent simulator to create the listing and a description in the attachment.

(02-25-2022 01:58 AM)rprosperi Wrote:  I have some Fortran stuff for Bairstow's method from the 70's in college, I've not seen much of this since then. I hope to take a closer look this weekend and maybe even find the old Fortran listings.

You may search for "Numerical Recipes in Fortran 77" and find in chapter 9.5 Roots of Polynomials on page 370 (page 400 of the PDF) a short description together with a Fortran program.

[Albert Chan]

(02-25-2022 01:20 AM)Thomas Klemm Wrote:  \(P(x) = Q(x) \cdot T(x) + R(x)\)

\(P(x) = a_n x^n + a_{n-1} x^{n-1} + \cdots + a_2 x^2 + a_1 x + a_0\)

\(T(x) = x^2 + p x + q\)

\(Q(x) = b_n x^{n-2} + b_{n-1} x^{n-3} + \cdots + b_4 x^2 + b_3 x + b_2\)

\(R(x) = b_1 ( x + p ) + b_0\)

quo(P*x²,T) = Q*x² + quo(R*x²,T) = Q*x² + b₁*x + b₀

This get all the b's. Same idea to get all the c's
HP Prime code (note: index 1 = head of list, index 0 = end of list)

Code:

#cas // P,T = polynomial, test quadratic
bairstow(P,T)
BEGIN
LOCAL b,c;
b := quo(extend(P,[0,0]),T);
c := quo(append(b,0),T);
c := linsolve([c[-2..-1],c[-1..0]],b[-1..0]);
RETURN T + [0,c[2],c[1]];
END;
#end

Example from Gerald's "Applied Numerical Analysis", p33

CAS> [1,1,1] // guess T = x^2+x+1
CAS> bairstow([1,-1.1,2.3,0.5,3.3],Ans)
[1,0.890309886867,1.06345302509]
[1,0.900024696323,1.10016130857]
[1,0.899999998585,1.09999999975]
[1,0.9,1.1]