From mboxrd@z Thu Jan  1 00:00:00 1970
X-Spam-Checker-Version: SpamAssassin 3.4.4 (2020-01-24) on polar.synack.me
X-Spam-Level: 
X-Spam-Status: No, score=-0.9 required=5.0 tests=BAYES_00,FORGED_GMAIL_RCVD,
	FREEMAIL_FROM autolearn=no autolearn_force=no version=3.4.4
X-Google-Thread: 103376,2b8382e0387f53c5
X-Google-NewGroupId: yes
X-Google-Attributes: gida07f3367d7,domainid0,public,usenet
X-Google-Language: ENGLISH,ASCII-7-bit
Path: 
 g2news2.google.com!postnews.google.com!y16g2000yqk.googlegroups.com!not-for-mail
From: marius63 <amado.alves@gmail.com>
Newsgroups: comp.lang.ada
Subject: Re: Eigenvalues to find roots of polynomials
Date: Tue, 26 Jul 2011 08:40:45 -0700 (PDT)
Organization: http://groups.google.com
Message-ID: 
 <8d424bae-53cb-4418-8ec1-1fbf3ec260e3@y16g2000yqk.googlegroups.com>
References: 
 <3ec148c9-a033-40a1-bc32-1aa52867478a@e35g2000yqc.googlegroups.com>
 <nospam-8E8EBD.02215326072011@news.aioe.org>
NNTP-Posting-Host: 188.82.29.171
Mime-Version: 1.0
Content-Type: text/plain; charset=ISO-8859-1
X-Trace: posting.google.com 1311694963 4575 127.0.0.1 (26 Jul 2011 15:42:43
 GMT)
X-Complaints-To: groups-abuse@google.com
NNTP-Posting-Date: Tue, 26 Jul 2011 15:42:43 +0000 (UTC)
Complaints-To: groups-abuse@google.com
Injection-Info: y16g2000yqk.googlegroups.com; posting-host=188.82.29.171;
 posting-account=3cDqWgoAAAAZXc8D3pDqwa77IryJ2nnY
User-Agent: G2/1.0
X-Google-Web-Client: true
X-Google-Header-Order: LECRUANKH
X-HTTP-UserAgent: Mozilla/5.0 (Macintosh; U; PPC Mac OS X 10_4_11; en)
 AppleWebKit/533.19.4 (KHTML,
 like Gecko) Version/4.1.3 Safari/533.19.4,gzip(gfe)
Xref: g2news2.google.com comp.lang.ada:21335
Date: 2011-07-26T08:40:45-07:00
List-Id: <comp.lang.ada>

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;