John B. Matthews, M.D.
Return home.
Ada.Numerics.Generic_Real_Arrays.Generic_Roots
Ada.Numerics.Generic_Complex_Arrays.Generic_Roots
Download a gzip'd tar file containing the code.
View the generic roots implementation for complex arrays.
View the sample test output for complex roots.
This collection is derived from the Ada.Numerics packages found in GNU gcc 4.3.4 (FSF) and those found on Martin Dowie's Web Page
This collection includes the current library (numerics) bound to BLAS/LAPACK, as well as the original, partially completed generic array packages, modified to compile under GNAT 4.3 (oldnum):
AI95-0296: Ada.Numerics.Generic_Real_Arrays
AI95-0296: Ada.Numerics.Generic_Complex_Arrays
These have been extended with non-standard, experimental implementations of Generic_Roots for real and complex arrays:
AI95-0346: Ada.Numerics.Generic_Real_Arrays.Generic_Roots
AI95-0346: Ada.Numerics.Generic_Complex_Arrays.Generic_Roots
The implementation uses the Durand-Kerner-Weierstrass method:
Durand-Kerner method
The method of Weierstrass
Univariate Polynomial Root-Finding with Lower Computational Precision and Higher Convergence Rates
Installation
To build and install the libraries, edit the prefix in numerics/makefile to match your own installation and enter the following command:
make -C numerics
If all goes well, install the new library additions:
make -C numerics install
You can build the test code, too:
make test
Implementation
The specification and implementation of Ada.Numerics.Generic_Complex_Arrays.Generic_Roots are shown here. The implementation of Ada.Numerics.Generic_Real_Arrays.Generic_Roots uses an instance of this procedure for polynomials with real coefficients.
pragma License (Unrestricted);
generic
procedure Ada.Numerics.Generic_Complex_Arrays.Generic_Roots
(P : in Complex_Vector;
R : out Complex_Vector);
-- The array P defines a polynomial P = P[0] + P[1]*z + P[2]*z**2 + ...
-- + P[n]*z**n where n is P'Last. Constraint_Error is raised if P'First
-- is not zero and if P'Last is not greater than zero. Constraint_Error
-- is also raised if the bounds of the out parameter R are not one and
-- P'Last.
--
-- The components R[1] .. R[n] are the roots of the equation P=0
-- in some order.
pragma Pure (Generic_Roots);
pragma License (Modified_GPL);
------------------------------------------------------------------
--
-- Ada.Numerics.Generic_Complex_Arrays.Generic_Roots
--
-- Author: John B. Matthews, Gem City Software
-- Distribution: GMGPL
-- Last Modified: November 2007
--
------------------------------------------------------------------
--
-- Implementation notes:
-- This is an experimental implementation of AI-346.
--
-- This implementation uses the Durand-Kerner-Weierstrass method
-- to find the complex roots of a polynomial with complex
-- coefficients. The method requires a monic polynomial; some
-- error may occur when dividing by the leading coefficient.
--
-- Set Debug := True to enable diagnostic output.
--
with Ada.Text_IO;
procedure Ada.Numerics.Generic_Complex_Arrays.Generic_Roots (
P : in Complex_Vector;
R : out Complex_Vector) is
Debug : constant Boolean := False;
Epsilon : constant Real'Base := 1.0 / (10.0 ** Real'Digits);
Max_Count : constant Natural := 999;
-- Use Durand-Kerner-Weierstrass method
procedure OrderN (P : in Complex_Vector; R : out Complex_Vector);
-- Determine if the components of two vectors are unchanging
function Done (A, B : Complex_Vector) return Boolean;
-- Debug output
procedure Dump (A : Complex_Vector);
-- Evaluate the polynomial P at X by Horner's method
function Eval (P : Complex_Vector; X : Complex) return Complex;
function Done (A, B : Complex_Vector) return Boolean is
Unchanged : Boolean := True;
Change : Complex;
begin
for I in A'Range loop
Change := A (I) - B (I);
Unchanged := Unchanged and
Real'Base (abs (Re (Change))) < Epsilon and
Real'Base (abs (Im (Change))) < Epsilon;
end loop;
return Unchanged;
end Done;
procedure Dump (A : Complex_Vector) is
begin
for I in A'Range loop
Ada.Text_IO.Put (I'Img & ": ");
Ada.Text_IO.Put (Re (A (I))'Img);
Ada.Text_IO.Put (Im (A (I))'Img);
Ada.Text_IO.New_Line;
end loop;
end Dump;
function Eval (P : Complex_Vector; X : Complex) return Complex is
Result : Complex := P (P'Last);
begin
for I in reverse 0 .. P'Last - 1 loop
Result := Result * X + P (I);
end loop;
return Result;
end Eval;
procedure OrderN (P : in Complex_Vector; R : out Complex_Vector) is
Count : Natural := 1;
P0 : constant Complex_Vector (P'Range) := P / P (P'Last);
A0, A1 : Complex_Vector (R'Range);
One : constant Complex := Compose_From_Cartesian (1.0);
Result : Complex := Compose_From_Cartesian (0.4, 0.9);
begin
-- initialize A0
A0 (A0'First) := One;
for I in A0'First + 1 .. A0'Last loop
A0 (I) := A0 (I - 1) * Result;
end loop;
-- iterate
loop
for I in A0'Range loop
Result := One;
for J in A0'Range loop
if I /= J then
Result := Result * (A0 (I) - A0 (J));
end if;
end loop;
A1 (I) := A0 (I) - (Eval (P0, A0 (I)) / Result);
end loop;
Count := Count + 1;
exit when Count > Max_Count or Done (A0, A1);
A0 := A1;
end loop;
-- report results
R := A1;
if Debug then
Ada.Text_IO.Put_Line ("Iterations: " & Count'Img);
Dump (A1);
end if;
end OrderN;
begin
if P'First /= 0 or P'First = P'Last or
R'First /= 1 or P'Last /= R'Last then
raise Constraint_Error;
end if;
if P'Last > 1 then
OrderN (P, R);
end if;
end Ada.Numerics.Generic_Complex_Arrays.Generic_Roots;
Sample test output
The test driver (croot.adb) produces the following output using Long_Float complex arithmetic. A similar program (root.adb) exercises the real array version. On the test platform, Long_Float'Digits is 15. For each polynomial P and root r, the maximal value of abs (P(r) – 0) is shown. The first few entries are easily factored polynomials with integer coefficients. These are followed by a series of all one polynomials. Finally, there are a few examples with complex coefficients. On PowerPC (G4), the largest useful order appears to be about 38; on Intel (Pentium 4) & AMD (Athlon), extended precision gives correct results up to order 48. Beyond that, round-off errors cause the roots to converge to zero.
Ada.Numerics.Long_Complex_Arrays.Roots:
Poly: 1.00x^2 - 3.00x^1 + 2.00
Real: 1.00000000000000E+00
Real: 2.00000000000000E+00
Largest error: < 1.0E-15
Poly: 1.00x^2 - 2.00x^1 + 1.00
Real: 1.00000000000000E+00
Real: 1.00000000000000E+00
Largest error: < 1.0E-15
Poly: 1.00x^2 + 0.00x^1 + 4.00
Comp: 0.00000000000000E+00 +- 2.00000000000000E+00i
Largest error: < 1.0E-15
Poly: 4.00x^2 + 0.00x^1 + 1.00
Comp: 0.00000000000000E+00 +- 5.00000000000000E-01i
Largest error: < 1.0E-15
Poly: 1.00x^3 - 3.00x^2 + 2.00x^1 + 0.00
Real: 0.00000000000000E+00
Real: 1.00000000000000E+00
Real: 2.00000000000000E+00
Largest error: < 1.0E-15
Poly: 1.00x^4 + 0.00x^3 - 13.00x^2 + 0.00x^1 + 36.00
Real: -3.00000000000000E+00
Real: -2.00000000000000E+00
Real: 2.00000000000000E+00
Real: 3.00000000000000E+00
Largest error: < 1.0E-15
Poly: 1.00x^5 + 0.00x^4 - 13.00x^3 + 0.00x^2 + 36.00x^1 + 0.00
Real: -3.00000000000000E+00
Real: -2.00000000000000E+00
Real: 0.00000000000000E+00
Real: 2.00000000000000E+00
Real: 3.00000000000000E+00
Largest error: < 1.0E-15
Poly: 1.00x^2 + 1.00x^1 + 1.00
Comp: -5.00000000000000E-01 +- 8.66025403784439E-01i
Largest error: < 1.0E-15
Poly: 1.00x^3 + 1.00x^2 + 1.00x^1 + 1.00
Real: -1.00000000000000E+00
Comp: 0.00000000000000E+00 +- 1.00000000000000E+00i
Largest error: < 1.0E-15
Poly: 1.00x^4 + 1.00x^3 + 1.00x^2 + 1.00x^1 + 1.00
Comp: -8.09016994374948E-01 +- 5.87785252292473E-01i
Comp: 3.09016994374947E-01 +- 9.51056516295154E-01i
Largest error: < 1.0E-15
Poly: 1.00x^5 + 1.00x^4 + 1.00x^3 + 1.00x^2 + 1.00x^1 + 1.00
Real: -1.00000000000000E+00
Comp: -5.00000000000000E-01 +- 8.66025403784439E-01i
Comp: 5.00000000000000E-01 +- 8.66025403784439E-01i
Largest error: < 1.0E-15
Poly: 1.00x^6 + 1.00x^5 + 1.00x^4 + 1.00x^3 + 1.00x^2 + 1.00x^1 + 1.00
Comp: -9.00968867902419E-01 +- 4.33883739117558E-01i
Comp: -2.22520933956314E-01 +- 9.74927912181824E-01i
Comp: 6.23489801858734E-01 +- 7.81831482468030E-01i
Largest error: < 1.0E-15
Poly: 1.00x^7 + 1.00x^6 + 1.00x^5 + 1.00x^4 + 1.00x^3 + 1.00x^2 + 1.00x^1 + 1.00
Real: -1.00000000000000E+00
Comp: -7.07106781186548E-01 +- 7.07106781186548E-01i
Comp: 0.00000000000000E+00 +- 1.00000000000000E+00i
Comp: 7.07106781186548E-01 +- 7.07106781186548E-01i
Largest error: < 1.0E-15
Poly: 1.00x^8 + 1.00x^7 + 1.00x^6 + 1.00x^5 + 1.00x^4 + 1.00x^3 + 1.00x^2 + 1.00x^1 + 1.00
Comp: -9.39692620785908E-01 +- 3.42020143325669E-01i
Comp: -5.00000000000000E-01 +- 8.66025403784439E-01i
Comp: 1.73648177666930E-01 +- 9.84807753012208E-01i
Comp: 7.66044443118978E-01 +- 6.42787609686539E-01i
Largest error: < 1.0E-15
Poly: 1.00x^9 + 1.00x^8 + 1.00x^7 + 1.00x^6 + 1.00x^5 + 1.00x^4 + 1.00x^3 + 1.00x^2 + 1.00x^1 + 1.00
Real: -1.00000000000000E+00
Comp: -8.09016994374948E-01 +- 5.87785252292473E-01i
Comp: -3.09016994374947E-01 +- 9.51056516295154E-01i
Comp: 3.09016994374947E-01 +- 9.51056516295154E-01i
Comp: 8.09016994374948E-01 +- 5.87785252292473E-01i
Largest error: < 1.0E-15
Poly: 1.00x^10 + 1.00x^9 + 1.00x^8 + 1.00x^7 + 1.00x^6 + 1.00x^5 + 1.00x^4 + 1.00x^3 + 1.00x^2 + 1.00x^1 + 1.00
Comp: -9.59492973614497E-01 +- 2.81732556841430E-01i
Comp: -6.54860733945285E-01 +- 7.55749574354258E-01i
Comp: -1.42314838273285E-01 +- 9.89821441880933E-01i
Comp: 4.15415013001886E-01 +- 9.09631995354518E-01i
Comp: 8.41253532831181E-01 +- 5.40640817455598E-01i
Largest error: < 1.0E-15
Poly: 1.00x^11 + 1.00x^10 + 1.00x^9 + 1.00x^8 + 1.00x^7 + 1.00x^6 + 1.00x^5 + 1.00x^4 + 1.00x^3 + 1.00x^2 + 1.00x^1 + 1.00
Real: -1.00000000000000E+00
Comp: -8.66025403784439E-01 +- 5.00000000000000E-01i
Comp: -5.00000000000000E-01 +- 8.66025403784439E-01i
Comp: 0.00000000000000E+00 +- 1.00000000000000E+00i
Comp: 5.00000000000000E-01 +- 8.66025403784439E-01i
Comp: 8.66025403784439E-01 +- 5.00000000000000E-01i
Largest error: < 1.0E-15
Poly: 1.00x^12 + 1.00x^11 + 1.00x^10 + 1.00x^9 + 1.00x^8 + 1.00x^7 + 1.00x^6 + 1.00x^5 + 1.00x^4 + 1.00x^3 + 1.00x^2 + 1.00x^1 + 1.00
Comp: -9.70941817426052E-01 +- 2.39315664287558E-01i
Comp: -7.48510748171101E-01 +- 6.63122658240795E-01i
Comp: -3.54604887042536E-01 +- 9.35016242685415E-01i
Comp: 1.20536680255323E-01 +- 9.92708874098054E-01i
Comp: 5.68064746731156E-01 +- 8.22983865893656E-01i
Comp: 8.85456025653210E-01 +- 4.64723172043769E-01i
Largest error: = 1.14160444461578E-15
Poly: 1.00x^13 + 1.00x^12 + 1.00x^11 + 1.00x^10 + 1.00x^9 + 1.00x^8 + 1.00x^7 + 1.00x^6 + 1.00x^5 + 1.00x^4 + 1.00x^3 + 1.00x^2 + 1.00x^1 + 1.00
Real: -1.00000000000000E+00
Comp: -9.00968867902419E-01 +- 4.33883739117558E-01i
Comp: -6.23489801858734E-01 +- 7.81831482468030E-01i
Comp: -2.22520933956314E-01 +- 9.74927912181824E-01i
Comp: 2.22520933956314E-01 +- 9.74927912181824E-01i
Comp: 6.23489801858734E-01 +- 7.81831482468030E-01i
Comp: 9.00968867902419E-01 +- 4.33883739117558E-01i
Largest error: < 1.0E-15
Poly: 1.00x^14 + 1.00x^13 + 1.00x^12 + 1.00x^11 + 1.00x^10 + 1.00x^9 + 1.00x^8 + 1.00x^7 + 1.00x^6 + 1.00x^5 + 1.00x^4 + 1.00x^3 + 1.00x^2 + 1.00x^1 + 1.00
Comp: -9.78147600733806E-01 +- 2.07911690817759E-01i
Comp: -8.09016994374948E-01 +- 5.87785252292473E-01i
Comp: -5.00000000000000E-01 +- 8.66025403784439E-01i
Comp: -1.04528463267653E-01 +- 9.94521895368273E-01i
Comp: 3.09016994374947E-01 +- 9.51056516295154E-01i
Comp: 6.69130606358858E-01 +- 7.43144825477394E-01i
Comp: 9.13545457642601E-01 +- 4.06736643075800E-01i
Largest error: = 1.51837754008763E-15
Poly: 1.00x^15 + 1.00x^14 + 1.00x^13 + 1.00x^12 + 1.00x^11 + 1.00x^10 + 1.00x^9 + 1.00x^8 + 1.00x^7 + 1.00x^6 + 1.00x^5 + 1.00x^4 + 1.00x^3 + 1.00x^2 + 1.00x^1 + 1.00
Real: -1.00000000000000E+00
Comp: -9.23879532511287E-01 +- 3.82683432365090E-01i
Comp: -7.07106781186548E-01 +- 7.07106781186548E-01i
Comp: -3.82683432365090E-01 +- 9.23879532511287E-01i
Comp: 0.00000000000000E+00 +- 1.00000000000000E+00i
Comp: 3.82683432365090E-01 +- 9.23879532511287E-01i
Comp: 7.07106781186548E-01 +- 7.07106781186548E-01i
Comp: 9.23879532511287E-01 +- 3.82683432365090E-01i
Largest error: = 1.11022302462516E-15
Poly: 1.00x^37 + 1.00x^36 + 1.00x^35 + 1.00x^34 + 1.00x^33 + 1.00x^32 + 1.00x^31 + 1.00x^30 + 1.00x^29 + 1.00x^28 + 1.00x^27 + 1.00x^26 + 1.00x^25 + 1.00x^24 + 1.00x^23 + 1.00x^22 + 1.00x^21 + 1.00x^20 + 1.00x^19 + 1.00x^18 + 1.00x^17 + 1.00x^16 + 1.00x^15 + 1.00x^14 + 1.00x^13 + 1.00x^12 + 1.00x^11 + 1.00x^10 + 1.00x^9 + 1.00x^8 + 1.00x^7 + 1.00x^6 + 1.00x^5 + 1.00x^4 + 1.00x^3 + 1.00x^2 + 1.00x^1 + 1.00
Real: -1.00000000000000E+00
Comp: -9.86361303402722E-01 +- 1.64594590280734E-01i
Comp: -9.45817241700635E-01 +- 3.24699469204684E-01i
Comp: -8.79473751206489E-01 +- 4.75947393037074E-01i
Comp: -7.89140509396394E-01 +- 6.14212712689668E-01i
Comp: -6.77281571625741E-01 +- 7.35723910673132E-01i
Comp: -5.46948158122427E-01 +- 8.37166478262529E-01i
Comp: -4.01695424652970E-01 +- 9.15773326655057E-01i
Comp: -2.45485487140799E-01 +- 9.69400265939330E-01i
Comp: -8.25793454723323E-02 +- 9.96584493006670E-01i
Comp: 8.25793454723323E-02 +- 9.96584493006670E-01i
Comp: 2.45485487140799E-01 +- 9.69400265939330E-01i
Comp: 4.01695424652970E-01 +- 9.15773326655058E-01i
Comp: 5.46948158122427E-01 +- 8.37166478262529E-01i
Comp: 6.77281571625741E-01 +- 7.35723910673132E-01i
Comp: 7.89140509396394E-01 +- 6.14212712689668E-01i
Comp: 8.79473751206489E-01 +- 4.75947393037074E-01i
Comp: 9.45817241700635E-01 +- 3.24699469204684E-01i
Comp: 9.86361303402722E-01 +- 1.64594590280734E-01i
Largest error: = 1.09838798499774E-14
Poly: (1 + 3i)x^2 + (2 + 2i)x + (3 + i)
Comp: -8.00000000000000E-01 - 6.00000000000000E-01i
Comp: 0.00000000000000E+00 + 1.00000000000000E+00i
Largest error: < 1.0E-15
Poly: (1 + i)x^3 + (2 + i)x^2 + (3 + i)x + (4 + i)
Comp: -1.40135938330275E+00 + 2.88269653138075E-01i
Comp: -2.84985631785343E-01 - 1.30378640290474E+00i
Comp: 1.86345015088091E-01 + 1.51551674976666E+00i
Largest error: < 1.0E-15
Poly: (1 + i)x^3 + (2 + 2i)x^2 + (3 + 3i)x + (4 + 4i)
Real: -1.65062919143939E+00
Comp: -1.74685404280306E-01 - 1.54686888723140E+00i
Comp: -1.74685404280306E-01 + 1.54686888723140E+00i
Largest error: < 1.0E-15
Licenses
This is free software; you can redistribute it and/or modify it under terms of the GNU General Public License (GPL) as published by the Free Software Foundation, either version 2, or (at your option) any later version. It is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GPL for more details. A copy of the GPL, incorporated by reference, may be obtained from gnu.org.
Files in this collection that specify "Distribution: GPL" are distributed under terms of the GNU General Public License (GPL).
Files in this collection that specify "pragma License (Modified_GPL)" are distributed under the GNAT modified GPL:
As a special exception, if other files instantiate generics from this unit, or you link this unit with other files to produce an executable, this unit does not by itself cause the resulting executable to be covered by the GNU General Public License. This exception does not however invalidate any other reasons why the executable file might be covered by the GNU Public License.
Files in this collection that specify "pragma License (Unrestricted)" are distributed with no license restrictions.
Copyright 2007 John B. Matthews
Distributed without warranty under the terms of the GNU Public License.
Last updated 10-Jun-2009
Return home.