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