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 the library routines, enter the following command:

make -C numerics

To test the newly compiled library routines, enter the following command:

make test

If you choose to install the libraries into the existing GNAT run-time library, edit the PREFIX in numerics/makefile to match your own installation and enter the following command:

make -C numerics install

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 11-Aug-2010
Return home.