EXPLANATION FILE OF BAIRSTOW'S METHOD
=====================================
This program calculates the real or complex roots of a polynomial f(x).
The Bairstow's Method
A. Principles
----------
i=n n-i
Let be f(x) = Sum a x the given polynomial and (p,q) two numbers,
i=0 i
the division of f(x) by x^2 + p x + q gives:
i=n n-i i=n n-i
Sum a x = ( Sum b x ) (x^2 + px + q) +rx + s
i=0 i i=2 i
The coefficients, b , r and s are obtained from p, q by:
i
b = b = 0 and b = a - p b - q b for 2 <= i <= n + 2 (1)
0 1 i i-2 i-1 i-2
Then r = b and s = b + p b (2)
n+1 n+2 n+1
Coefficients b , r and s are functions of (p,q), we seek (p,q) such as
i
r(p,q) = s(p,q) = 0 (2 bis)
The Newton's method allows determining one root of system (2 bis).
By differential calculus, it can be linearized as follows:
| r + delta p dr/dp + delta q dr/dq = 0
| (3)
| s + delta p ds/dp + delta q ds/dq = 0
Here d sign means partial derivative sign.
So after testing (p,q), we'll test (p + delta p, q + delta q) where
(delta p, delta q) is a solution of linear system (3).
To calculate dr/dp, you just have to to calculate the db /dp terms by
recursion, derivating the formula (1):
db db
i-1 i-2
db / dp = b - p ----- - q -----
i i-1 dp dp
So if we put c = db /dp for 0 <= i <= n+2, the series c is defined
i i i
by the formulas:
c = c = 0, c = -b - p c - q c
0 1 i i-1 i-1 i-2
According to formulas (2): dr/dp = c and ds/dp = -q c
n+1 n
Beside, by derivating formulas (1), we can notice that db /dq = c
i i-1
for each i, so:
dr/dq = c and ds/dq = c + p c
n n+1 n
The linear system (3) can be easily solved by using Cramer's formulas:
2
d = c - c (c + b )
n+1 n n+2 n+1
c b - b c -b (q c + p c - b c
n n+2 n+1 n+1 n+1 n n+2 n+1
dp = ------------------ and dq = ----------------------------
d d
So we define a series of points (p, q), converging towards a point (p,q),
such as x^2 + p x + q is a divisor of f(x). Practically, we stop the
iterations when K=100 and we use the stopping criterion:
|delta p| + |delta q|
-------------------- < eps
|p| + |q|
with eps = 1E-3, for instance.
The equation x^2 + p x + q is then solved in the complex domain by usual
formulas. A first error approximation is given by |f(x0)/f'(x0)| according
to Newton's formula.
n n-i
The program continues with the new polynomial Sum b x , if the degree is
i=2 i
greater then 3.
B. Programming Technique
---------------------
The Bairstow program consists of several modules:
1) seek p and q as described above.
2) solve a 2nd degree equation (real and complex roots).
3) calculate |f/f'| by using a modified Hörner's method to deal with
complex numbers.
C. Example
-------
Let us seek all the real or complex roots of polynomial:
6 5 4 3 2
x - 127 x + 215 x + 28 x - 39 x + 20 x - 15 = 0
The program gives the following results:
roots error approximation
--------------------------------------------------------
0.03989619444 + i 0.44667717900 1E-11
0.03989619444 - i 0.44667717900 1E-11
0.5238337918 1E-6
-0.6457525530 1E-5
125,2821089 1E-8
1,760014124 1E-6
It is not likely that all decimals are correct.
Let us verify the first real root by applying the Hörner's rule and checking
the result error (not explained here):
x f(x) error
-------------------------------------
0.523833 -1.6E-4 5E-9
0.523834 -8.5E-5 5E-9
0.523835 -7.0E-6 5E-9
0.523836 7.1E-5 5E-9
So the first real root is r3 = 0.523835 +/- 1E-6.
Accuracy may be increased by using Newton's methode or this program with
eps = 1E-9 (instead of 1E-3).
For the other three real roots, we found:
-0.64575 +/- 1E-5, 125.2821089 +/- 1E-7 and 1.760013 +/- 1E-6.
It is more difficult to approximate the error for the complex roots.
Our estimation is: 0.0398962 +/- i 0.4466718 +/- 1E-7 .
From [BIBLI 06].
-----------------------------------------------
End of file bairstow.txt