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;