(15C) Bairstow's Method - Printable Version +- HP Forums (https://www.hpmuseum.org/forum) +-- Forum: HP Software Libraries (/forum-10.html) +--- Forum: General Software Library (/forum-13.html) +--- Thread: (15C) Bairstow's Method (/thread-18079.html) (15C) Bairstow's Method - Thomas Klemm - 02-25-2022 01:20 AM 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 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. 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-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$$ RE: (15C) Bairstow's Method - rprosperi - 02-25-2022 01:58 AM 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. RE: (15C) Bairstow's Method - Thomas Klemm - 02-25-2022 07:04 AM 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. RE: (15C) Bairstow's Method - Albert Chan - 02-27-2022 12:40 PM (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]