A polynomial of degree n, P(x), has a factor of (x - r) if and only if P(r) = 0. Then r is said to be a root of the polynomial. Such a polynomial has a factorization of the form:
where Q(x) is a polynomial of degree n - 1. If P(x) has real coefficients and if x is a complex root of P(x), then the complex conjugate of x is also a root of P(x). A polynomial can have all real roots, or a mix of real roots and pairs of complex conjugate roots.P(x) = (x - r)Q(x) = 0
which may have pairs of complex conjugate roots:x2 - 2·a·x + a2 + b2
The advantage of this is that complex arithmetic is avoided. Polynomial deflation is used to simplify the problem.x = a ± ib
Zeno 1.2 is an interpreter for the Zeno programming language. It is an easy to learn and is suitable for educational purposes.
const TOLERANCE : real := 0.000000000001
const DIM : int := 11 % at least poly order + 2
type rvector : array[DIM] of real
type complex : record
re : real
im : real
end record
type cvector : array[DIM] of complex
program
var coeff : rvector
var roots : cvector
var r : complex
var num_roots, poly_order, i : int
% X^4 - 16 = 0
poly_order := 4
coeff[1] := -16
coeff[2] := 0
coeff[3] := 0
coeff[4] := 0
coeff[5] := 1
put "\nthe polynomial:"
print_polynomial( coeff, poly_order )
if Lin_Bairstow( coeff, roots, poly_order ) then
put "has roots:"
for i := 1...poly_order do
r.re := roots[i].re
r.im := roots[i].im
put cstr( r )
end for
else
put "could not solve"
end if
% X^5 - 4·X^4 + 1·X^3 + 4·X^2 + 10·X - 12 = 0
poly_order := 5
coeff[1] := -12
coeff[2] := 10
coeff[3] := 4
coeff[4] := 1
coeff[5] := -4
coeff[6] := 1
put "\nthe polynomial:"
print_polynomial( coeff, poly_order )
if Lin_Bairstow( coeff, roots, poly_order ) then
put "has roots:"
for i := 1...poly_order do
r.re := roots[i].re
r.im := roots[i].im
put cstr( r )
end for
else
put "could not solve"
end if
end program
%
% The Lin-Bairstow method improves on the deflation method to solve
% for the complex roots of the following polynomial:
%
% y = coeff(1) + coeff(2) X + coeff(3) X^2 +...+ coeff(n+1) X^n
%
% parameters:
%
% coeff must be an array with at least order+1 elements.
% roots output array of roots
% order order of polynomial
%
function Lin_Bairstow( var coeff : rvector,
var roots : cvector,
order : int ) : boolean
const SMALL : real := 0.00000001
var a, b, c, d : rvector
var alfa1, alfa2 : real
var beta1, beta2 : real
var delta1, delta2, delta3 : real
var i, j, k : int
var count : int
var n : int
var n1 : int
var n2 : int
n := order
n1 := n + 1
n2 := n + 2
% is the coefficient of the highest term zero?
if abs( coeff[1] ) < SMALL then
return false
end if
for i := 0...n do
a[n1-i+1] := coeff[i+1]
end for
% is highest coeff not close to 1?
if abs( a[2] - 1 ) > SMALL then
% adjust coefficients because a(1) != 1
for i := 3...n2 do
a[i] := a[i] / a[2]
end for
a[2] := 1
end if
% initialize root counter
count := 1
%
% start the main Lin-Bairstow iteration loop
% initialize the counter and guesses for the
% coefficients of quadratic factor:
%
% p(x) = x^2 + alfa1 * x + beta1
%
repeat
% these are just guesses
alfa1 := 2
beta1 := 1
repeat
b[1] := 0
d[1] := 0
b[2] := 1
d[2] := 1
j := 2
k := 1
for i := 3...n2 do
b[i] := a[i] - alfa1 * b[j] - beta1 * b[k]
d[i] := b[i] - alfa1 * d[j] - beta1 * d[k]
j := j + 1
k := k + 1
end for
j := n
k := n - 1
delta1 := d[j]^2 - ( d[n+1] - b[n+1] ) * d[k]
if delta1 = 0 then
return false % cannot divide by zero
end if
alfa2 := ( b[n+1] * d[j] - b[n1+1] * d[k] ) / delta1
beta2 := ( b[n1+1] * d[j] - ( d[n+1] - b[n+1] ) * b[n+1] ) / delta1
alfa1 := alfa1 + alfa2
beta1 := beta1 + beta2
until ( abs( alfa2 ) < TOLERANCE ) and ( abs( beta2 ) < TOLERANCE )
delta1 := alfa1^2 - 4 * beta1
if delta1 < 0 then % imaginary roots
delta2 := sqrt( abs( delta1 ) ) / 2
delta3 := -alfa1 / 2
roots[count].re := delta3
roots[count].im := +delta2
roots[count+1].re := delta3
roots[count+1].im := -delta2
else % roots are real
delta2 := sqrt( delta1 )
roots[count].re := ( delta2 - alfa1 ) / 2
roots[count].im := 0
roots[count+1].re := ( delta2 + alfa1 ) / ( -2 )
roots[count+1].im := 0
end if
% update root counter
count := count + 2
% reduce polynomial order
n := n - 2
n1 := n1 - 2
n2 := n2 - 2
% for n >= 2 calculate coefficients of the new polynomial
if n >= 2 then
for i := 2...n2 do
a[i] := b[i]
end for
end if
until n < 2
% obtain last single real root
if n = 1 then
roots[count].re := -b[3]
roots[count].im := 0
end if
return true
end function
% display polynomial in readable format
procedure print_polynomial( var coeff : rvector, order : int )
var i : int
var dot : char := chr( 183 )
for decreasing i := order+1...1 do
continue when coeff[i] = 0.0
if i = order+1 then
if coeff[i] ~= 1 then
put coeff[i], dot...
end if
put "X^", i-1...
elsif( i > 2 ) then
if coeff[i] < 0 then
put " - "...
else
put " + "...
end if
put abs( coeff[i] ), dot, "X^", i-1...
elsif i = 2 then
if coeff[i] < 0 then
put " - "...
else
put " + "...
end if
put abs( coeff[i] ), dot, "X"...
elsif i = 1 then
if coeff[i] < 0 then
put " - "...
else
put " + "...
end if
put abs( coeff[i] )...
end if
end for
put " = 0"
end procedure
% returns the absolute value of the argument
function abs( x : real ) : real
if x < 0.0 then
x := -x
end if
return x
end function
% returns a complex number as a string
function cstr( var c : complex ) : string
var s : string
var re, im : real
re := c.re
im := c.im
if im >= 0.0 then
s := frealstr(re,14,9 ) & " +" & frealstr(im,12,9 ) & "j"
else
im := -im
s := frealstr(re,14,9 ) & " -" & frealstr(im,12,9 ) & "j"
end if
return s
end function
the polynomial: X^4 - 16 = 0 has roots: 2.000000000 + 0.000000000j -2.000000000 + 0.000000000j 0.000000000 + 2.000000000j 0.000000000 - 2.000000000j the polynomial: X^5 - 4·X^4 + 1·X^3 + 4·X^2 + 10·X - 12 = 0 has roots: -1.000000000 + 1.000000000j -1.000000000 - 1.000000000j 2.000000000 + 0.000000000j 1.000000000 + 0.000000000j 3.000000000 + 0.000000000j