Re: Eigenvalues to find roots of polynomials

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

And the winner is...

> <http://home.roadrunner.com/~jbmatthews/misc/groots.html>

John B. Matthews Ada...Generic_Roots works like a charm.
It passed the simple test case, and solved my quintic problem x^5 -
4178x + 4177

$ ./find_roots_mathews 1 0 0 0 -4178 4177
Index =>  0, Position =>  6, Value => 4177
Index =>  1, Position =>  5, Value => -4178
Index =>  2, Position =>  4, Value => 0
Index =>  3, Position =>  3, Value => 0
Index =>  4, Position =>  2, Value => 0
Index =>  5, Position =>  1, Value => 1
Root 1 => (Re=> 1.00000000000000E+00, Im=> 0.00000000000000E+00)
Root 2 => (Re=>-2.47582467311450E-01, Im=>-8.05881611194044E+00)
Root 3 => (Re=> 7.76752669636581E+00, Im=>-1.40129846432482E-45)
Root 4 => (Re=>-2.47582467311450E-01, Im=> 8.05881611194044E+00)
Root 5 => (Re=>-8.27236176174291E+00, Im=> 0.00000000000000E+00)

I know from previous mathematical derivation (the polynomial comes
from extracting a variable in a summation), that the solution is a non-
negative scalar different from one, so Root 3 is my number, with the
non-zero but very small value of the imaginary part certainly only a
residue of computation done inside Matthews's magic machine. Thanks a
very great lot. Source code follows. Sorry for the mispelling,
"Mathews".

--  Trying to compute the roots of a polynomial using John Mathews
procedure.

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;
with Ada.Numerics.Generic_Complex_Arrays.Generic_Roots;

procedure Find_Roots_Mathews 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);
   procedure Complex_Roots is new Complex_Arrays.Generic_Roots;
   M, N : Natural;
begin
   M := Ada.Command_Line.Argument_Count;
   N := M - 1;
   --  N = number of arguments = degree of the polynomial + 1.
   --  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
      P : Complex_Arrays.Complex_Vector (0 .. N);
      R : Complex_Arrays.Complex_Vector (1 .. N);
      use Ada.Text_IO;
   begin
      for I in 0 .. N loop
         Put_Line ("Index => " & I'img &
                   ", Position => " & Integer'Image(M - I) &
                   ", Value => " & Ada.Command_Line.Argument (M - I));
         P (I) := (Re => Real'Value (Ada.Command_Line.Argument (M -
I)), Im => 0.0);
      end loop;
      Complex_Roots (P, R);
      for I in R'Range loop
         Put_Line ("Root" & I'img &
                   " => (Re=>" & R(I).Re'img &
                   ", Im=>" & R(I).Im'img & ")");
      end loop;
   end;
end;