NUMAL Section 3.6.1

BEGIN SECTION : 3.6.1 (December, 1978)

AUTHORS:  TH.J. DEKKER AND TH.H.P. REYMER

CONTRIBUTOR:  TH.H.P. REYMER

INSTITUTE:  UNIVERSITY OF AMSTERDAM

RECEIVED:  770427;

BRIEF DESCRIPTION:

    THIS SECTION CONTAINS TWO PROCEDURES:
    A) ZERPOL   CALCULATES  ALL  ROOTS   OF  A   POLYNOMIAL  WITH  REAL
       COEFFICIENTS BY MEANS OF LAGUERRE'S METHOD;
    B) BOUNDS   CALCULATES  UPPERBOUNDS FOR THE ABSOLUTE ERROR IN GIVEN
       APPROXIMATED ZEROS OF A POLYNOMIAL WITH REAL COEFFICIENTS;

KEYWORDS:

    ZEROS;
    REAL POLYNOMIAL;
    LAGUERRE'S METHOD;
    ERROR BOUNDS;


SUBSECTION: ZERPOL

CALLING SEQUENCE:

    THE HEADING OF THE PROCEDURE READS:
    "INTEGER" "PROCEDURE" ZERPOL(N, A, EM, RE, IM, D);
    "VALUE" N; "INTEGER" N; "ARRAY" A, EM, RE, IM, D;
    "CODE" 34501 ;

    THE MEANING OF THE FORMAL PARAMETERS IS:
    N:      <ARITHMETIC EXPRESSION>;
            ENTRY:  THE DEGREE OF THE POLYNOMIAL;
    A:      <ARRAY IDENTIFIER>;
            "ARRAY" A[0 : N];
            ENTRY:  THE COEFFICIENTS OF THE POLYNOMIAL, IN  SUCH  A WAY
                    THAT
                    P(Z) = (..(A[N]*Z + A[N-1])*Z +..+ A[1])*Z + A[0];

    EM:     <ARRAY IDENTIFIER>;
            "ARRAY" EM[0 : 4];
            ENTRY:  EM[0]: MACHINE PRECISION;
                    EM[1]: THE MAXIMAL NUMBER OF ITERATIONS ALLOWED FOR
                           EACH ZERO;
            EXIT:   EM[2]: FAIL INDICATION;
                           0  SUCCESSFUL CALL;
                           1  UPON ENTRY DEGREE N <= 0;
                           2  UPON ENTRY LEADING COEFFICIENT A[N] = 0;
                           3  NUMBER OF ITERATIONS EXCEEDED EM[1];
                    EM[3]: NUMBER OF NEW STARTS IN THE LAST ITERATION;
                    EM[4]: TOTAL NUMBER OF ITERATIONS PERFORMED;
            FOR THE CD CYBER 70 SYSTEM SUITABLE VALUES ARE:
            EM[0]:= "-14;
            EM[1]:= 40; IF, UPON EXIT, EM[2] = 3  AND EM[3] < 5 THEN IT
                        MAY BE  USEFUL TO  START  AGAIN  WITH  A HIGHER
                        VALUE OF EM[1];
    RE, IM: <ARRAY IDENTIFIERS>;
            "ARRAY" RE, IM[1 : N];
            EXIT:   THE  REAL  AND  IMAGINARY PARTS OF THE ZEROS OF THE
                    POLYNOMIAL; THE  MEMBERS  OF  EACH  NONREAL COMPLEX
                    CONJUGATE PAIR ARE CONSECUTIVE;
    D:      <ARRAY IDENTIFIER>;
            "ARRAY" D[0 : N];
            EXIT:   IF THE CALL IS UNSUCCESSFUL AND ONLY N-K ZEROS HAVE
                    BEEN   FOUND,  THEN   D[0 :  K]   CONTAINS   THE
                    COEFFICIENTS   OF   THE   (DEFLATED)    POLYNOMIAL;
                    MOREOVER,  THEN THE ZEROS  FOUND  ARE  DELIVERED IN
                    RE, IM[ K + 1 : N], WHEREAS THE REMAINING  PARTS
                    OF RE AND IM CONTAIN NO INFORMATION;

    ZERPOL:= THE NUMBER, K, OF ZEROS NOT FOUND;

PROCEDURES USED:

    DWARF   = CP30003;
    GIANT   = CP30004;
    COMABS  = CP34340;
    COMSQRT = CP34343;

REQUIRED CENTRAL MEMORY:

    TOTAL SIZE OF LOCAL ARRAYS IS N + 16 REAL LOCATIONS;

RUNNING TIME:

    ROUGHLY PROPORTIONAL TO N**2;

METHOD AND PERFORMANCE:

    THE  PROCEDURE  USES  LAGUERRE'S  METHOD TO FIND ZEROS OF THE GIVEN
    POLYNOMIAL (SEE [2]);  WHEN A ZERO HAS BEEN FOUND, A COMPOSITE
    DEFLATION TECHNIQUE IS USED TO  OBTAIN A  NEW POLYNOMIAL OF LOWER
    DEGREE (SEE [1], [3], [4]); IF CONVERGENCE IS NOT APPARENT, SEVERAL
    RESTARTS, THE NUMBER OF WHICH DEPENDS ON  EM[1]  BUT HAS A MAXIMUM
    OF  5, ARE MADE IN THE NEIGHBOURHOOD OF THE ABSOLUTE LARGEST ZERO;
    THE  ACCURACY OF THE  CALCULATED  ZEROS  STRONGLY  DEPENDS  ON  THE
    POLYNOMIAL; A ROUGH INDICATION FOR THE ERROR IN A CALCULATED ZERO Z
    FOLLOWS FROM P(Z) / DP(Z), WHERE P DENOTES THE GIVEN POLYNOMIAL AND
    DP ITS  FIRST DERIVATIVE (SEE E.G. [5]); TO FIND A  TRUE UPPERBOUND
    FOR THESE ERRORS, ONE CAN USE PROCEDURE BOUNDS (SEE NEXT
    SUBSECTION); FOR  A MORE  DETAILED DESCRIPTION OF THE PROCEDURE AND
    TEST RESULTS SEE [4];

REFERENCES:

    [1] D.A. ADAMS, A STOPPING CRITERION FOR POLYNOMIAL ROOT FINDING,
        CACM 10, NO 10, PP. 655-658, OCTOBER 1967;
    [2] T.J. DEKKER, NEWTON-LAGUERRE ITERATION,
        MATHEMATISCH CENTRUM MR82, 1966;
    [3] G. PETERS AND J.H. WILKINSON, PRACTICAL PROBLEMS ARISING IN THE
        SOLUTION OF POLYNOMIAL EQUATIONS,
        J. INST. MATHS APPLICS 1971, NO. 8, PP. 16-35;
    [4] TH.H.P. REYMER, BEREKENING VAN NULPUNTEN VAN REELE POLYNOMEN EN
        FOUTGRENZEN VOOR DEZE NULPUNTEN,
        DOCTORAAL SCRIPTIE UVA, APRIL 1977;
    [5] J.H. WILKINSON, ROUNDING ERRORS IN ALGEBRAIC PROCESSES,
        PRENTICE HALL, 1963;

EXAMPLE OF USE:

    FOR AN EXAMPLE OF USE SEE PROCEDURE BOUNDS (NEXT SUBSECTION);


SUBSECTION: BOUNDS

CALLING SEQUENCE:

    THE HEADING OF THE PROCEDURE READS:
    "PROCEDURE" BOUNDS(N,A,RE,IM,RELE,ABSE,RECENTRE,IMCENTRE,BOUND);
    "VALUE" N, RELE, ABSE; "INTEGER" N; "REAL" RELE, ABSE;
    "ARRAY" A, RE, IM, RECENTRE, IMCENTRE, BOUND;
    "CODE" 34502 ;

    THE MEANING OF THE FORMAL PARAMETERS IS:
    N:        <ARITHMETIC EXPRESSION>;
              ENTRY:  DEGREE OF THE POLYNOMIAL;
    A:        <ARRAY IDENTIFIER>;
              "ARRAY" A[0 : N];
              ENTRY:  THE  COEFFICIENTS  OF  THE  POLYNOMIAL  OF  WHICH
                      RE[J] + I * IM[J]  ARE  THE  APPROXIMATED  ZEROS,
                      IN SUCH A WAY THAT
                      P(Z) = (..(A[N]*Z + A[N-1])*Z +..+A[1])*Z + A[0]
    RE, IM:   <ARRAY IDENTIFIERS>;
              "ARRAY" RE, IM[1 : N];
              ENTRY:  REAL AND IMAGINARY PARTS OF APPROXIMATED ZEROS OF
                      A  POLYNOMIAL,  SUCH  THAT  THE  MEMBERS  OF EACH
                      NONREAL COMPLEX CONJUGATE PAIR ARE CONSECUTIVE;
              EXIT:   A PERMUTATION OF THE INPUT DATA;
    RELE:     <ARITHMETIC EXPRESSION>;
              ENTRY:  RELATIVE ERROR IN THE  NON-VANISHING COEFFICIENTS
                      A[J] OF THE GIVEN POLYNOMIAL;
    ABSE:     <ARITHMETIC EXPRESSION>;
              ENTRY:  ABSOLUTE ERROR IN THE VANISHING COEFFICIENTS A[J]
                      GIVEN  POLYNOMIAL;  IF  THERE  ARE  NO  VANISHING
                      COEFFICIENTS, ABSE SHOULD BE ZERO;
    RECENTRE, IMCENTRE: <ARRAY IDENTIFIERS>;
              "ARRAY" RECENTRE, IMCENTRE[1 : N];
              EXIT:   REAL  AND IMAGINARY PARTS OF THE CENTERS OF DISKS
                      IN WHICH  SOME  NUMBER OF ZEROS OF THE POLYNOMIAL
                      GIVEN BY A ARE SITUATED;
                      THE NUMBER OF IDENTICAL CENTERS  DENOTES
                      THE NUMBER OF ZEROS IN THAT DISK;
    BOUND:    <ARRAY IDENTIFIER>;
              "ARRAY" BOUND[1 : N];
              EXIT:   RADIUS  OF THE  DISKS  WHOSE  CENTERS  ARE  GIVEN
                      CORRESPONDINGLY IN RECENTRE AND IMCENTRE;

PROCEDURES USED:

    ARREB = CP30002;
    GIANT = CP30004;

REQUIRED CENTRAL MEMORY:

    TOTAL SIZE OF LOCAL ARRAYS IS AT MOST 7 * N REAL LOCATIONS;

RUNNING TIME:

    APPROXIMATELY OF ORDER N**2;

METHOD AND PERFORMANCE:

    FROM  THE  APPROXIMATED  ZEROS  A  POLYNOMIAL IS  RECONSTRUCTED AND
    COMPARED  WITH  THE  GIVEN  POLYNOMIAL; SUBSEQUENTLY, THE PROCEDURE
    CALCULATES  DISKS  SUCH THAT THE NUMBER OF GIVEN APPROXIMATED ZEROS
    WITHIN EACH DISK EQUALS THE NUMBER OF ZEROS OF THE GIVEN POLYNOMIAL
    WITHIN  THAT  DISK;  UPON  EXIT  EVERY  TWO NON-IDENTICAL DISKS ARE
    DISJOINT;
    FOR A MORE DETAILED DESCRIPTION SEE [1], [2];

REFERENCES:

    [1] G. PETERS AND J.H. WILKINSON, PRACTICAL PROBLEMS ARISING IN THE
        SOLUTION OF POLYNOMIAL EQUATIONS, J.INST.MATHS APPLICS 1971,
        NO 8, PP. 16-35;
    [2] TH.H.P. REYMER, BEREKENING VAN NULPUNTEN VAN REELE POLYNOMEN EN
        FOUTGRENZEN VOOR DEZE NULPUNTEN,
        DOCTORAAL SCRIPTIE UVA, APRIL 1977;

EXAMPLE OF USE:

"BEGIN" "INTEGER" I, J;
  "ARRAY" A, D[0:7], RE, IM[1:7], EM[0:4];

  A[7]:= 1; A[6]:= -3; A[5]:= -3; A[4]:= 25; A[3]:= -46;
  A[2]:= 38; A[1]:= -12; A[0]:= 0;
  EM[0]:= "-14; EM[1]:= 40;
  I:= ZERPOL(7, A, EM, RE, IM, D);
  OUTPUT(61,"(""("COEFFICIENTS OF POLYNOMIAL:")",//")");
  "FOR" J:=7 "STEP" -1 "UNTIL" 0 "DO"
  OUTPUT(61,"("-2ZD3B")",A[J]);
  OUTPUT(61,"("//,"("NUMBER NOT FOUND ZEROS ")",3ZD,/,
     "("FAIL INDICATION ")",3ZD,/,"("NUMBER NEW STARTS ")",3ZD,/,
     "("NUMBER OF ITERATIONS ")",3ZD,/")",I,EM[2],EM[3],EM[4]);
  OUTPUT(61,"("/,"("ZEROS: ")",/")");
  "FOR" J:= I+1 "STEP" 1 "UNTIL" 7 "DO" "IF" IM[J] = 0
  "THEN" OUTPUT(61,"("/,N")",RE[J])
  "ELSE" OUTPUT(61,"("/,2(N)")",RE[J],IM[J]);

  "IF" I = 0 "THEN"
  "BEGIN" "ARRAY" RECENTRE, IMCENTRE, BOUND[1:7];

    BOUNDS(7, A, RE, IM, 0, 0, RECENTRE, IMCENTRE, BOUND);
    OUTPUT(61,"("2/,"("REAL AND IMAG. PART OF CENTRE + RADIUS")",/")");
    "FOR" J:= 1 "STEP" 1 "UNTIL" 7 "DO"
    OUTPUT(61,"("/,3(N)")",RECENTRE[J],IMCENTRE[J],BOUND[J])
  "END"
"END"

RESULTS :

COEFFICIENTS OF POLYNOMIAL:

   1     -3     -3     25    -46     38    -12      0

NUMBER NOT FOUND ZEROS    0
FAIL INDICATION    0
NUMBER NEW STARTS    0
NUMBER OF ITERATIONS   11

ZEROS:

+2.0000000000000"+000
-3.0000000000000"+000
+1.0000000000000"+000  -1.0000000000000"+000
+1.0000000000000"+000  +1.0000000000000"+000
+1.0000000083024"+000
+9.9999999169752"-001
+0.0000000000000"+000

REAL AND IMAG. PART OF CENTRE + RADIUS

+2.0000000000000"+000  +0.0000000000000"+000  +1.3238111117716"-011
-3.0000000000000"+000  +0.0000000000000"+000  +3.8857510604494"-013
+1.0000000000000"+000  -1.0000000000000"+000  +4.0912729775463"-012
+1.0000000000000"+000  +1.0000000000000"+000  +4.0912729775463"-012
+9.9999999999998"-001  +0.0000000000000"+000  +2.2533888428865"-006
+9.9999999999998"-001  +0.0000000000000"+000  +2.2533888428865"-006
+0.0000000000000"+000  +0.0000000000000"+000  +0.0000000000000"+000

SOURCE TEXT(S) :

"CODE" 34501;

"CODE" 34502;