Quadratic Solver

Program


001- 42,21,11   LBL A
002-       36   ENTER
003-   43  33   R↑
004-       10   ÷
005-   43  33   R↑
006-   43  36   LSTx
007-       10   ÷
008-        2   2
009-       16   CHS
010-       10   ÷
011-       36   ENTER
012-       36   ENTER
013-   43  11   x²
014-   43  33   R↑
015-       30   -
016-       11   √x
017-       30   -
018-       34   x<>y
019-   43  36   LSTx
020-       40   +
021-   43  32   RTN

Example

4 x2 + 7 x + 3 = 0 
4
ENTER
7
ENTER
3
GSB A

Y:      -1.0000
X:      -0.7500

Bairstow's Method

Algorithm

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

P(x) = Q(x) · T(x) + R(x)
P(x) = an xn + an-1 xn-1 + … + a2 x2 + a1 x + a0
T(x) = x2 + p x + q
Q(x) = bn xn-2 + bn-1 xn-3 + … + b4 x2 + b3 x + b2
R(x) = b1 ( x + p ) + b0

If we set bn+1 ≡ bn+2 ≡ 0, then ∀ k ≤ n:

bk = ak - bk+1 p - bk+2 q

The two coefficients of the remainder R(x) depend on p and q:
b0 = b0( p, q )
b1 = b1( 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. Therfore we set the first order Taylor series expansion to 0:

b0( p + δp, q + δq ) ≈ b0( p, q ) + ∂b0/∂p δp + ∂b0/∂q δq = 0
b1( p + δp, q + δq ) ≈ b1( p, q ) + ∂b1/∂p δp + ∂b1/∂q δq = 0

It turns out that the partial derivatives of bk with respect to p and q can be calculated with a similar formula:
ck+1 = - ∂bk/∂p
ck+2 = - ∂bk/∂q

ck = bk - ck+1 p - ck+2 q

Thus we have to solve:
c1 δp + c2 δq = b0
c2 δp + c3 δq = b1

Registers


R0 :  c  =  cj
R1 :  c′ =  cj+1
R2 :  c″ =  cj+2
R3 :  b  =  bi
R4 :  
R5 :  b′ =  bi+1
R6 :  index (9.xxx)
R7 :  p
R8 :  q
R9 :  an
R.0:  an-1
R.1:  an-2
R.2:  …
R.3:  …

Program


022 - 42,21,12   LBL B                             053 - 44,40, 7   STO + 7            
023 -   42  32   CLEAR Σ                           054 -       34   x<>y         
024 -    45  6   RCL 6                             055 - 44,40, 8   STO + 8            
025 -   44  25   STO I                             056 -   42  26   →P            
026 - 42,21, 0   LBL 0                             057 -   43  34   RND                
027 -    45  3   RCL 3                             058 -   43  30   x≠0             
028 -    45  1   RCL 1                             059 -   22  12   GTO B              
029 -    44  2   STO 2                             060 -    45  6   RCL 6              
030 -    45  8   RCL 8                             061 -        2   2                  
031 -       20   ×                                 062 -       26   EEX                
032 -       30   -                                 063 -        3   3                  
033 -    45  0   RCL 0                             064 -       16   CHS                
034 -    44  1   STO 1                             065 -       30   -                  
035 -    45  7   RCL 7                             066 -    44  6   STO 6              
036 -       20   ×                                 067 -   44  25   STO I              
037 -       30   -                                 068 -        0   0                  
038 -    44  0   STO 0                             069 -       36   ENTER              
039 -   45  24   RCL (i)                           070 - 42,21, 1   LBL 1              
040 -    45  5   RCL 5                             071 -   45  24   RCL (i)            
041 -    45  8   RCL 8                             072 -    45  8   RCL 8              
042 -       20   ×                                 073 -   43  33   R↑            
043 -       30   -                                 074 -       20   ×            
044 -    45  3   RCL 3                             075 -       30   -                  
045 -    44  5   STO 5                             076 -    45  7   RCL 7              
046 -    45  7   RCL 7                             077 -   43  33   R↑            
047 -       20   ×                                 078 -       20   ×            
048 -       30   -                                 079 -       30   -                  
049 -    44  3   STO 3                             080 -   44  24   STO (i)            
050 -    42  6   ISG                               081 -    42  6   ISG                
051 -    22  0   GTO 0                             082 -    22  1   GTO 1              
052 -   42  49   L.R.                              083 -   43  32   RTN

Description

Initialization k = n

Lines 023 - 025
( b, b′, c, c′, c″) ← (0, 0, 0, 0, 0)
I ← index

Polynominal division

Lines 026 - 051
At the end of the loop b = b0 but c = c1!
That's why the second division is performed before the first.

Iteration k+1 → k

Lines 027 - 038
( c, c′, c″) ← (b - c p - c′ q, c, c′)

Lines 039 - 049
( b, b′) ← (ak - b p - b′ q, b)

Solve the linear equation

c δp + c′ δq = b
c′ δp + c″ δq = b′

Line 052
The trick using linear regression only works since the matrix is symmetric.

Find improved values for p and q

Lines 053 - 055
p ← p + δp
q ← q + δq

Stop Criterion

Lines 056 - 058
|(δp, δq)| < ε
ε = 10-n if display is set to FIX n

Polynominal Division

Lines 060 - 082
We have to repeat the division here since we didn't keep bk in lines 039 - 049.

Example

P(x) = x5 - 15 x4 + 85 x3 - 225 x2 + 274 x - 120 = 0 
1      STO 9
-15    STO .0
85     STO .1
-225   STO .2
274    STO .3
-120   STO .4
9.014  STO 6
1      STO 7
       STO 8
GSB B         -59.999999

RCL 7          -3.000000
RCL 8           2.000000

RCL 9           1.000000
RCL .0        -12.000000
RCL .1         47.000000
RCL .2        -59.999999
x5 - 15 x4 + 85 x3 - 225 x2 + 274 x - 120 =
(x2 - 3 x + 2)(x3 - 12 x2 + 47 x - 60) 
1      STO 7
       STO 8
GSB B          -5.000000

RCL 7          -7.000000
RCL 8          12.000000

RCL 9           1.000000
RCL .0         -5.000000
(x2 - 3 x + 2)(x3 - 12 x2 + 47 x - 60) =
(x2 - 3 x + 2)(x2 - 7 x + 12)(x - 5)

Now use the quadratic solver to find:
(x2 - 3 x + 2) = (x - 1)(x -2)
(x2 - 7 x + 12) = (x - 3)(x - 4)

Thus the factorisation of P(x) is:
(x - 1)(x -2)(x - 3)(x - 4)(x - 5) = 0

Solutions:
x1 = 1
x2 = 2
x3 = 3
x4 = 4
x5 = 5

Short & Sweet Math Challenges #6 [LONG]
Just for the record, the solutions are: 

  M = 9


     x3 - x2 - 8x + 4 = 0

        x = 3.14133611565 (0.00025653794)

  M = 20


     6x3 - 16x2 - 15x + 19 = 0

        x = 3.14159104610  (0.00000160749)

Just in case it might interest you, here is an assortment of results for larger values of M:


  M = 49

     19x3 - 47x2 - 30x - 31 = 0

        x = 3.14159235658  (2.97E-07)


  M = 99

     33x3 - 92x2 - 65x + 89 = 0

        x = 3.14159264441 (9.18E-09)


  M = 199

     37x3 - 114x2 - 36x + 91 = 0

        x = 3.14159265383  (2.39E-10 )

Links

Bairstow's Method
Numerical Recipes in Fortran
9.5 Roots of Polynomials