(15C) Bairstow's Method
02-25-2022, 01:20 AM
Post: #1
 Thomas Klemm Senior Member Posts: 1,814 Joined: Dec 2013
(15C) Bairstow's Method
Bairstow's Method

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 R2 - R7 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.

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-2)(x+3)(x^2-4x+13)$$

Solutions

$$x_1=1.5$$
$$x_2=2$$
$$x_3=-3$$
$$x_4=2+3i$$
$$x_5=2-3i$$

Attached File(s)
02-25-2022, 01:58 AM
Post: #2
 rprosperi Super Moderator Posts: 5,529 Joined: Dec 2013
RE: (15C) Bairstow's Method
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.

--Bob Prosperi
02-25-2022, 07:04 AM
Post: #3
 Thomas Klemm Senior Member Posts: 1,814 Joined: Dec 2013
RE: (15C) Bairstow's Method
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.
02-27-2022, 12:40 PM
Post: #4
 Albert Chan Senior Member Posts: 2,101 Joined: Jul 2018
RE: (15C) Bairstow's Method
(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² + b1*x + b0

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]
 « Next Oldest | Next Newest »

User(s) browsing this thread: 1 Guest(s)