Eigenvalues to find roots of polynomials

[marius63 <amado.alves@gmail.com>]

I am trying to use Eigenvalues methods to find the roots of
polynomials... and failing miserably!
My math ability is not stellar. I used the information at

(1)
http://mathworld.wolfram.com/PolynomialRoots.html (equation 12)

and

(2)
http://en.wikipedia.org/wiki/Root_finding
(One possibility is to form the companion matrix of the polynomial.
Since the eigenvalues of this matrix coincide with the roots of the
polynomial, one can use any eigenvalue algorithm to find the roots of
the polynomial.)
http://en.wikipedia.org/wiki/Companion_matrix

I tried both schemes with the test data in (1). The output is first
the matrix, for verification, then the (wrong) roots (they should be
-1, 1, 2) and the result of applying them to the polynomial (and
yielding miserably non-zero values of course). The source code
follows. What am I missing? Thanks a lot.

$ ./find_roots 1 -2 -1 2
 5.00000E-01 1.00000E+00-5.00000E-01
 1.00000E+00 0.00000E+00 0.00000E+00
 0.00000E+00 1.00000E+00 0.00000E+00
-7.6156E-01=> 4.2003E-01
 4.1249E+00=>-1.6015E+01
 6.3667E-01=> 5.9465E-01

$ ./find_roots_wikipedia -2 -1 2
 0.00000000000000E+00,  0.00000000000000E+00, -2.00000000000000E+00,
 1.00000000000000E+00,  0.00000000000000E+00,  1.00000000000000E+00,
 0.00000000000000E+00,  1.00000000000000E+00,  2.00000000000000E+00,
-1.17008648662603E+00 => -3.38499151386795E-02
 6.88892182534018E-01 =>  1.96496630422799E+00
 2.48119430409202E+00 =>  3.00688836109107E+01


--  Trying to compute the roots of a polynomial using the method
--  described in eq. (12) of http://mathworld.wolfram.com/PolynomialRoots.html
--  (consulted July 2011)

with Ada.Command_Line;
with Ada.Numerics.Generic_Real_Arrays;
with Ada.Text_IO;
with Ada.Numerics.Generic_Complex_Types;
with Ada.Numerics.Generic_Complex_Arrays;

procedure Find_Roots is
   type Real is digits 5;
   package Real_Arrays is new Ada.Numerics.Generic_Real_Arrays (Real);
   package Complex_Types is new Ada.Numerics.Generic_Complex_Types
(Real);
   package Complex_Arrays is new Ada.Numerics.Generic_Complex_Arrays
(Real_Arrays, Complex_Types);
   N, M : Natural;
begin
   N := Ada.Command_Line.Argument_Count;
   --  N = number of arguments = degree of the polynomial.
   --  The input consists of all coficients in descending degree,
   --  e.g. for polynomial
   --     x^3 - 2x^2 - x + 2
   --  the input is
   --     1 -2 -1 2
   --  Roots of the above (for testing):
   --     -1, 1, 2

   declare
      C : Complex_Arrays.Complex_Matrix (1 .. N - 1, 1 .. N - 1) :=
(others => (others => (0.0, 0.0)));
      --  Using a complex matrix because real matrices for eigenvalues
must be symmetric.
      --  Below we convert all real values to complex (with null Im
part).
      --  The complex matrix must be Hamiltonian;
      --  I thought a non-symmetric converted to complex would
      --  yield a non-Hamiltonian matrix, but I tried and it
worked :-)... hmm... NOT :-(
      A0 : Real := Real'Value (Ada.Command_Line.Argument (N));
      use Ada.Text_IO;
   begin
      for I in 1 .. N - 1 loop
         C (1, I).Re := - Real'Value (Ada.Command_Line.Argument (N -
I)) / A0;
         if I < N - 1 then
            C (I + 1, I).Re := 1.0;
         end if;
      end loop;

      --  inspect matrix
      for i in c'range(1) loop
         for j in c'range(2) loop
            put (c(i,j).re'img);
         end loop;
         new_line;
      end loop;

      declare
         V : Real_Arrays.Real_Vector := Complex_Arrays.Eigenvalues
(C);
         use Real_Arrays;
         Z : Real;
         use Ada.Text_IO;
      begin
         for I in V'Range loop
            V (I) := 1.0 / V (I);
            Put (Real'Image (V (I)));
            Z := 0.0;
            for J in 1 .. N - 1 loop
               Z := Z + C (1, J).Re * V (I) ** J;
            end loop;
            Put ("=>" & Real'Image (Z));
            New_Line;
         end loop;
      end;
   end;
end;

--  Build with:
--  gnatmake find_roots -largs /Developer/SDKs/MacOSX10.4u.sdk/usr/lib/
libblas.dylib /Developer/SDKs/MacOSX10.4u.sdk/usr/lib/liblapack.dylib


--  Given an Nth degree polynomial A(N)X**N + ... + A(1)X + A(0),
--  the roots can be found by finding the eigenvalues L(I) of the
matrix
--     _                                                           _
--    |   -A(1)/A(0)   -A(2)/A(0)   -A(3)/A(0)   ...   -A(N)/A(0)   |
--    |        1            0            0       ...        0       |
--    |        0            1            0       ...        0       |
--    |       ...          ...           1       ...        0       |
--    |_       0            0            0       ...        0      _|
--
--  and taking R(i) = 1/L(I). This method can be computationally
expensive,
--  but is fairly robust at finding close and multiple roots.
--    (wolphram)


--  This program computes the roots of a polynomial using the method
--  described in http://en.wikipedia.org/wiki/Companion_matrix

--  This method is for monic polynomials, i.e. with coficient of the
term
--  of greater degree equal to 1, so the input is for the remainder
terms,
--  in descending order e.g. for polynomial
--     x^3 - 2x^2 - x + 2
--  the input is
--     -2 -1 2

with Ada.Command_Line;
with Ada.Numerics.Generic_Real_Arrays;
with Ada.Text_IO;
with Ada.Numerics.Generic_Complex_Types;
with Ada.Numerics.Generic_Complex_Arrays;

procedure Find_Roots_Wikipedia is
   type Real is digits 15;
   package Real_Arrays is new Ada.Numerics.Generic_Real_Arrays (Real);
   package Complex_Types is new Ada.Numerics.Generic_Complex_Types
(Real);
   package Complex_Arrays is new Ada.Numerics.Generic_Complex_Arrays
(Real_Arrays, Complex_Types);
   N : Natural;
begin
   N := Ada.Command_Line.Argument_Count;
   declare
      C : Complex_Arrays.Complex_Matrix (1 .. N, 1 .. N) := (others =>
(others => (0.0, 0.0)));
      --  Using a complex matrix because real matrices for eigenvalues
must be symmetric.
      --  Below we convert all real values to complex (with null Im
part).
      --  (Still, the complex matrix must be Hamiltonian...)
      use Ada.Text_IO;
   begin
      for I in 1 .. N loop
         C (I, N).Re := - Real'Value (Ada.Command_Line.Argument (N - I
+ 1));
         if I > 1 then
            C (I, I - 1).Re := 1.0;
         end if;
      end loop;

      --  inspect matrix
      for i in c'range(1) loop
         for j in c'range(2) loop
            put (c(i,j).re'img & ", ");
         end loop;
         new_line;
      end loop;

      declare
         V : Real_Arrays.Real_Vector := Complex_Arrays.Eigenvalues
(C);
         use Real_Arrays;
         Z : Real;
         use Ada.Text_IO;
      begin
         for I in V'Range loop
            Put (Real'Image (V (I)));
            Z := V (I) ** N;
            for J in 1 .. N loop
               Z := Z + C (N, J).Re * V (I) ** (J - 1);
            end loop;
            Put (" => " & Real'Image (Z));
            New_Line;
         end loop;
      end;
   end;
end;

--  Build with:
--  gnatmake find_roots_wikipedia -largs /Developer/SDKs/
MacOSX10.4u.sdk/usr/lib/libblas.dylib /Developer/SDKs/MacOSX10.4u.sdk/
usr/lib/liblapack.dylib