Approximating polynomial roots without complex arithmetics

Question

Suppose that we have polynomial of degree greater than two
$P(x)=a_{n}x^n+a_{n-1}x^{n-1}+_\cdots+a_{1}x+a_{0}$
with real coefficients and we want to approximate all roots of this polynomial
without complex arithmetic
It seems to be good idea to get quadratic factors of $P(x)$
$P(x)=Q(x)(x^2-px+q)+R(p,q)x+S(p,q)$
Now we want that after each iteration $R(p,q)$ and $S(p,q)$ be closer to zero
until given tolerance $\varepsilon$
Now we need to solve system of nonliear equations $$ R(p,q)=0\\ S(p,q)=0\\ $$
and how to solve this system and present solution in such way it is easy to write code for it

I would like to present the code I wrote

uses crt;
type TElem = double;
     TArray = array of TElem;
const EPS = 1e-12;
procedure Horner(deg:integer;a:TArray;z:TElem;var b:TArray;var alpha:TElem;var beta:TElem);
var k:integer;
begin
  alpha := a[deg];
  beta := 0;
  for k := deg - 1 downto 0 do
  begin
    b[k]:=alpha;
    beta := alpha + z * beta;
    alpha := a[k] + z * alpha;
  end;
end;
procedure Bairstow(deg:integer;a:TArray;var b:TArray;r0,s0:TElem;var r:TElem;var s:TElem);
var i:integer;
    d,dr,ds:TElem;
    c:TArray;
begin
  SetLength(c,deg+1);
  repeat
    b[deg] := a[deg];
    b[deg-1] := a[deg-1]+r0b[deg];
    c[deg] := 0;
    c[deg-1] := b[deg];
    for i := deg - 2 downto 0 do
    begin
      b[i] := a[i] + r0  b[i + 1] + s0 * b[i + 2];
      c[i] := b[i + 1] + r0 * c[i + 1] + s0 * c[i + 2];
    end;
    d := (c[0] * c[2] - sqr(c[1]));
    dr := (c[1] * b[1] - c[2] * b[0])/d;
    ds := (c[1] * b[0] - c[0] * b[1])/d;
 r0 := r0 + dr;
 s0 := s0 + ds;
 until ((abs(dr)<EPS)and(abs(ds) < EPS));
 r := r0;
 s := s0;
end;
procedure PolyRoots(deg:integer;a:TArray;p0,q0:TElem;var xr:TArray;var xi:TArray);
var i:integer;
    b,c:TArray;
    p,q,d:TElem;
    stop:boolean;
begin
  SetLength(b,deg+1);
  SetLength(c,deg+1);
  for i := 0 to deg do
      c[i] := a[i]/a[deg];
  for i := 0 to deg-1 do
  begin
    xr[i] := 0;
    xi[i] := 0;
  end;
  {Modyfikacja zaproponowana przez uzytkownika Lutz}
  if odd(deg) then
  begin
    d := 0;
    for i := 0 to deg - 1 do
        if abs(c[i]) > d then
           d := abs(c[i]);
    d := d + 1;
    repeat
      Horner(deg,c,d,b,p,q);
      xr[deg-1] := d - p/q;
      stop := (abs(xr[deg-1] - d) < EPS);
      d := xr[deg-1];
    until stop;
    deg := deg - 1;
    for i := 0 to deg do
        c[i] := b[i];

  end;
  while (deg >= 2) do
  begin
    Bairstow(deg,c,b,p0,q0,p,q);
    d := sqr(p) + 4q;
    if (d < 0) then
    begin
      xr[deg-1] := 0.5  p;
      xr[deg-2] := 0.5 * p;
      xi[deg-1] := -0.5 * sqrt(-d);
      xi[deg-2] := 0.5 * sqrt(-d);
    end
    else
    begin
      xr[deg-1] := 0.5 * (p - sqrt(d));
      xr[deg-2] := 0.5 * (p + sqrt(d));
    end;
    deg := deg - 2;
    for i := 0 to deg do
        c[i] := b[i+2];
    end;
    {
    if (deg = 1) then
       xr[0] := -c[0]/c[1];
    }
    end;
var
    i,deg:integer;
    a,xr,xi:TArray;
    p0,q0:TElem;
    esc : char;
begin
  ClrScr;
  WriteLn('Przyblizone rozwiazywanie rownan wielomianowych metoda Bairstowa');
  repeat
    WriteLn('Podaj stopien wielomianu');
    ReadLn(deg);
    SetLength(a,deg+1);
    SetLength(xr,deg);
    SetLength(xi,deg);
    WriteLn('Podaj poczatkowe wartosci wspolczynnikow trojmianu');
    ReadLn(p0,q0);
    WriteLn('Podaj wspolczynniki wielomianu');
    for i := deg downto 0 do
    begin
      Write('P[',i,']=');
      ReadLn(a[i]);
    end;
    PolyRoots(deg,a,p0,q0,xr,xi);
    WriteLn('Przyblizone wartosci pierwiastkow wielomianu to: ');
    for i := 0 to deg - 1 do
       if xi[i] < 0 then
          WriteLn('x[',i,']=',xr[i]:1:12,xi[i]:1:12,'i')
       else
          WriteLn('x[',i,']=',xr[i]:1:12,'+',xi[i]:1:12,'i');
    esc := ReadKey;
    until esc = #27;
 end.

I used Pascal to present the code which I have written so far because it produces clean code which is easy to read and understand even for those who are not programmers

This code still needs to be improved because of drawbacks and pitfalls of Newton's method which is used here

Suppose we want to stay with Newton's method Are there strategies for picking good initial guess or maybe there is modification of Newton's method which guaranties convergence

Answer 1

The (Lin-)Bairstow method is what you ask for, the Newton method applied to the task of splitting of a quadratic factor from a polynomial

The factorization task to solve is $$ P(x)=Q(x)·N(x),~~~~~N(x)=x^2-px+q. $$ With an approximation for the second factor, the linearization can be written as $$ P(x)=Q(x)·N(x)+Q(x)·ΔN(x)+ΔQ(x)·N(x) $$ The aim is to compute the update $ΔN$. This polynomial does not change under reduction of the whole equation and all operations in it modulo $N$ $$ [P(x)\bmod N(x)] = \bigl([Q(x)\bmod N(x)]·ΔN(x)\bigr)\bmod N(x) $$ This last system is now effectively a $2\times 2$ linear system for the coefficients $p,q$ in $ΔN(x)=-Δp·x+Δq$.

Repeat until the changes are small, or any other convergence test for Newton's method (tests for at least halving of the update length from one step to the next, tests for quadratic convergence,...)

In practical tests this iteration converges globally to some random factor for even degree polynomials. In odd-degree polynomials divergence can happen in a considerable part of the initial values, indicating a linear factor with root $x=q/p$. Either use that as a single root or split off a real root with some other method for polynomials with real coefficients, so that the next steps can proceed from even degree to even degree.