-- Generic_Complex_Polynomials with ARRAY_EXCEPTIONS, ELEMENTARY_FUNCTIONS_EXCEPTIONS; generic -- package USER_NAME is new Generic_Complex_Polynomials -- (Float, Complex, Complex_Vector); -- The type Complex typically comes from Complex_Types -- The type Complex_Vector typically comes from Complex_Arrays type REAL is digits <>; type COMPLEX is private; type COMPLEX_VECTOR is array (integer range <>) of COMPLEX; with function RE (X : COMPLEX) return REAL is <>; with function IM (X : COMPLEX) return REAL is <>; with procedure SET_RE (X : in out COMPLEX; RE : in REAL) is <>; with procedure SET_IM (X : in out COMPLEX; IM : in REAL) is <>; with function COMPOSE_FROM_CARTESIAN (RE, IM : REAL) return COMPLEX is <>; with function MODULUS (X : COMPLEX) return REAL is <>; with function "abs" (X : COMPLEX) return REAL is <>; with function ARGUMENT (X : COMPLEX) return REAL is <>; with function ARGUMENT (X : COMPLEX; CYCLE : REAL) return REAL is <>; with function COMPOSE_FROM_POLAR (MODULUS, ARGUMENT : REAL) return COMPLEX is <>; with function COMPOSE_FROM_POLAR (MODULUS, ARGUMENT, CYCLE : REAL) return COMPLEX is <>; with function "-" (RIGHT : COMPLEX) return COMPLEX is <>; with function CONJUGATE (X : COMPLEX) return COMPLEX is <>; with function "+" (LEFT, RIGHT : COMPLEX) return COMPLEX is <>; with function "-" (LEFT, RIGHT : COMPLEX) return COMPLEX is <>; with function "*" (LEFT, RIGHT : COMPLEX) return COMPLEX is <>; with function "/" (LEFT, RIGHT : COMPLEX) return COMPLEX is <>; with function "**" (LEFT : COMPLEX; RIGHT : INTEGER) return COMPLEX is <>; with function "+" (LEFT : REAL; RIGHT : COMPLEX) return COMPLEX is <>; with function "-" (LEFT : REAL; RIGHT : COMPLEX) return COMPLEX is <>; with function "*" (LEFT : REAL; RIGHT : COMPLEX) return COMPLEX is <>; with function "/" (LEFT : REAL; RIGHT : COMPLEX) return COMPLEX is <>; with function "/" (LEFT : COMPLEX; RIGHT : REAL) return COMPLEX is <>; package GENERIC_COMPLEX_POLYNOMIALS is -- COMPLEX POLYNOMIALS PACKAGE SPECIFICATION -- -- PURPOSE : TO PROVIDE POLYNOMIAL ARITHMETIC CAPABILITY IN ADA -- -- METHOD : THIS PACKAGE DEFINES THE TYPE COMPLEX_POLYNOMIAL -- -- THIS PACKAGE DEFINES THE OVERLOADED OPERATORS -- + - * / UNARY + UNARY - -- FOR THE COMPLEX_POLYNOMIAL TYPE -- -- THIS PACKAGE ALSO DEFINES THE OVERLOADED PROCEDURE -- AND FUNCTIONS: -- INTEGRAL for the COMPLEX_POLYNOMIAL type. -- DERIVITIVE for the COMPLEX_POLYNOMIAL type. -- EVALUATE for the COMPLEX_POLYNOMIAL type. -- ZERO_POLY for the COMPLEX_POLYNOMIAL type. -- MAKE_FROM_ROOTS for the COMPLEX_POLYNOMIAL type -- SOLVE_FOR_ROOTS for the COMPLEX_POLYNOMIAL type -- VECTOR_TO_POLYNOMIAL for COMPLEX_POLYNOMIAL type -- POLYNOMIAL_TO_VECTOR for COMPLEX_POLYNOMIAL type -- FORCE_POWER for COMPLEX_POLYNOMIAL type type COMPLEX_POLYNOMIAL is array ( INTEGER range <> ) of COMPLEX ; -- POLYNOMIALS ARE REPRESENTED BY THE COEFFICIENT VECTOR -- EXAMPLE: 3.7 + 2.5*X + 9.1*X**2 + ... + 4.2*X**(N-1) + 6.5*X**N -- IS CODED: -- C:COMPLEX_POLYNOMIAL (0..N) := (3.7,2.5,9.1, ... ,4.2,6.5); -- REPRESENTING THE GENERAL POLYNOMIAL: -- C(0) + C(1)*X + C(2)*X**2 + ... + C(N-1)*X**(N-1)+ C(N)*X**N -- -- -- highest power of return function "+" ( A , B : COMPLEX_POLYNOMIAL ) return COMPLEX_POLYNOMIAL ; -- MAX(A'LAST,B'LAST) function "+" ( A : COMPLEX_POLYNOMIAL ) return COMPLEX_POLYNOMIAL ; -- A'LAST function "-" ( A , B : COMPLEX_POLYNOMIAL ) return COMPLEX_POLYNOMIAL ; -- MAX(A'LAST,B'LAST) function "-" ( A : COMPLEX_POLYNOMIAL ) return COMPLEX_POLYNOMIAL ; -- A'LAST function "*" ( A , B : COMPLEX_POLYNOMIAL ) return COMPLEX_POLYNOMIAL ; -- A'LAST+B'LAST function "*" ( A : COMPLEX; B : COMPLEX_POLYNOMIAL ) return COMPLEX_POLYNOMIAL ; -- B'LAST function "*" ( A : COMPLEX_POLYNOMIAL; B : COMPLEX ) return COMPLEX_POLYNOMIAL ; -- A'LAST function "/" ( DIVIDEND , DIVISOR : COMPLEX_POLYNOMIAL ) return COMPLEX_POLYNOMIAL ; -- DIVIDEND'LAST-DIVISOR'LAST may be 0 function "/" ( DIVIDEND : COMPLEX_POLYNOMIAL; DIVISOR : COMPLEX ) return COMPLEX_POLYNOMIAL ; -- DIVIDEND'LAST function "rem" ( DIVIDEND , DIVISOR : COMPLEX_POLYNOMIAL)return COMPLEX_POLYNOMIAL ; -- DIVISOR'LAST-1 may be 0 function EVALUATE ( A : COMPLEX_POLYNOMIAL ; X : COMPLEX ) return COMPLEX ; function INTEGRAL ( A : COMPLEX_POLYNOMIAL ) return COMPLEX_POLYNOMIAL ; -- A'LAST+1 function DERIVITIVE ( A : COMPLEX_POLYNOMIAL ) return COMPLEX_POLYNOMIAL ; -- A'LAST-1 function ZERO_POLY ( N : INTEGER ) return COMPLEX_POLYNOMIAL ; -- N function MAKE_FROM_ROOTS ( V : COMPLEX_VECTOR ) return COMPLEX_POLYNOMIAL ; -- V'LENGTH function SOLVE_FOR_ROOTS ( POLY : COMPLEX_POLYNOMIAL ) return COMPLEX_VECTOR ; function VECTOR_TO_POLYNOMIAL ( V : COMPLEX_VECTOR ) return COMPLEX_POLYNOMIAL ; -- V'LENGTH-1 function POLYNOMIAL_TO_VECTOR ( P : COMPLEX_POLYNOMIAL ) return COMPLEX_VECTOR ; -- 1..P'LENGTH vector function FORCE_POWER ( P : COMPLEX_POLYNOMIAL; N : INTEGER ) return COMPLEX_POLYNOMIAL; -- N POLYNOMIAL_ERROR : exception; -- raised if P'FIRST /= 0 and invalid use end GENERIC_COMPLEX_POLYNOMIALS ; with TEXT_IO ; with GENERIC_ELEMENTARY_FUNCTIONS; package body GENERIC_COMPLEX_POLYNOMIALS is -- WRITTEN BY : JON SQUIRE , 27 FEB 1983 -- MODIFIED 11 JUNE 1987 -- MODIFIED 28 DEC 1993 for ISO standard packages -- AUGMENTED 16 JUNE 1994 for working with other packages -- MODIFIED 9 SEPT 1996 for similar structure to Complex_Arrays standard -- MODIFIED 24 SEPT to include Solve_For_Roots (don't need extra package) package Elementary_Functions is new Generic_Elementary_functions(Real); use Elementary_Functions; function MAX ( LEFT , RIGHT : INTEGER ) return INTEGER is begin if LEFT > RIGHT then return LEFT ; end if ; return RIGHT ; end MAX ; function "+" ( A , B : COMPLEX_POLYNOMIAL ) return COMPLEX_POLYNOMIAL is C : COMPLEX_POLYNOMIAL ( 0 .. MAX( A'LAST , B'LAST )) ; begin if A'FIRST /= 0 or B'FIRST /= 0 then raise POLYNOMIAL_ERROR; end if; if B'LAST <= A'LAST then C := A ; for I in B'RANGE loop C ( I ) := C ( I ) + B ( I ) ; end loop ; else C := B ; for I in A'RANGE loop C ( I ) := C ( I ) + A ( I ) ; end loop ; end if ; return C ; end "+" ; function "+" ( A : COMPLEX_POLYNOMIAL ) return COMPLEX_POLYNOMIAL is begin if A'FIRST /= 0 then raise POLYNOMIAL_ERROR; end if; return A ; end "+" ; function "-" ( A , B : COMPLEX_POLYNOMIAL ) return COMPLEX_POLYNOMIAL is C : COMPLEX_POLYNOMIAL ( 0 .. MAX( A'LAST , B'LAST )) ; begin if A'FIRST /= 0 or B'FIRST /= 0 then raise POLYNOMIAL_ERROR; end if; if B'LAST <= A'LAST then C := A ; for I in B'RANGE loop C ( I ) := C ( I ) - B ( I ) ; end loop ; else C := - B ; for I in A'RANGE loop C ( I ) := C ( I ) + A ( I ) ; end loop ; end if ; return C ; end "-" ; function "-" ( A : COMPLEX_POLYNOMIAL ) return COMPLEX_POLYNOMIAL is C : COMPLEX_POLYNOMIAL ( 0 .. A'LAST ) ; begin if A'FIRST /= 0 then raise POLYNOMIAL_ERROR; end if; for I in A'RANGE loop C ( I ) := - A ( I ) ; end loop ; return C ; end "-" ; function "*" ( A , B : COMPLEX_POLYNOMIAL ) return COMPLEX_POLYNOMIAL is C : COMPLEX_POLYNOMIAL ( 0 .. A'LAST + B'LAST ) ; begin if A'FIRST /= 0 or B'FIRST /= 0 then raise POLYNOMIAL_ERROR; end if; C := ZERO_POLY ( A'LAST + B'LAST ) ; for I in A'RANGE loop for J in B'RANGE loop C ( I + J ) := C ( I + J ) + A ( I ) * B ( J ) ; end loop ; end loop ; return C ; end "*" ; function "*" ( A : COMPLEX; B : COMPLEX_POLYNOMIAL ) return COMPLEX_POLYNOMIAL is C : COMPLEX_POLYNOMIAL ( B'RANGE ) ; begin if B'FIRST /= 0 then raise POLYNOMIAL_ERROR; end if; for J in B'RANGE loop C ( J ) := A * B ( J ) ; end loop ; return C ; end "*" ; function "*" ( A : COMPLEX_POLYNOMIAL; B : COMPLEX ) return COMPLEX_POLYNOMIAL is C : COMPLEX_POLYNOMIAL ( A'RANGE ) ; begin if A'FIRST /= 0 then raise POLYNOMIAL_ERROR; end if; for I in A'RANGE loop C ( I ) := A ( I ) * B ; end loop ; return C ; end "*" ; function "/" ( DIVIDEND , DIVISOR : COMPLEX_POLYNOMIAL ) return COMPLEX_POLYNOMIAL is QUOTIENT : COMPLEX_POLYNOMIAL ( 0 .. DIVIDEND'LAST - DIVISOR'LAST ) ; REMAINDER : COMPLEX_POLYNOMIAL ( DIVIDEND'RANGE ) ; REMAINDER_INDEX : INTEGER := DIVIDEND'LAST ; ZEROS : COMPLEX_POLYNOMIAL ( QUOTIENT'RANGE ) ; TRIAL : COMPLEX_POLYNOMIAL ( QUOTIENT'RANGE ) ; begin if DIVIDEND'FIRST /= 0 or DIVISOR'FIRST /= 0 or DIVISOR'LAST > DIVIDEND'LAST then raise POLYNOMIAL_ERROR; end if; if abs DIVISOR ( DIVISOR'LAST ) = 0.0 or abs DIVIDEND ( DIVIDEND'LAST ) = 0.0 then raise POLYNOMIAL_ERROR ; -- must be proper polynomials end if ; -- ZEROS := ZERO_POLY ( QUOTIENT'LAST ) ; REMAINDER := DIVIDEND ; for L in reverse 0 .. DIVIDEND'LAST - DIVISOR'LAST loop QUOTIENT ( L ) := REMAINDER ( REMAINDER_INDEX ) / DIVISOR ( DIVISOR'LAST ) ; TRIAL := ZEROS ; TRIAL ( L ) := QUOTIENT ( L ) ; REMAINDER := REMAINDER - TRIAL * DIVISOR ; REMAINDER_INDEX := REMAINDER_INDEX - 1 ; end loop ; return QUOTIENT ; end "/" ; function "/" ( DIVIDEND : COMPLEX_POLYNOMIAL; DIVISOR : COMPLEX ) return COMPLEX_POLYNOMIAL is C : COMPLEX_POLYNOMIAL ( DIVIDEND'RANGE ) ; begin if DIVIDEND'FIRST /= 0 then raise POLYNOMIAL_ERROR; end if; for I in DIVIDEND'RANGE loop C ( I ) := DIVIDEND ( I ) / DIVISOR ; end loop ; return C ; end "/"; function "rem" ( DIVIDEND , DIVISOR : COMPLEX_POLYNOMIAL ) return COMPLEX_POLYNOMIAL is QUOTIENT : COMPLEX_POLYNOMIAL ( 0 .. DIVIDEND'LAST - DIVISOR'LAST ) ; REMAINDER : COMPLEX_POLYNOMIAL ( 0 .. DIVIDEND'LAST ) ; REMAINDER_INDEX : INTEGER := DIVIDEND'LAST ; ZEROS : COMPLEX_POLYNOMIAL ( 0 .. QUOTIENT'LAST ) := ZERO_POLY ( QUOTIENT'LAST ) ; TRIAL : COMPLEX_POLYNOMIAL ( 0 .. QUOTIENT'LAST ) ; begin if DIVIDEND'FIRST /= 0 or DIVISOR'FIRST /= 0 or DIVIDEND'LAST < 1 then raise POLYNOMIAL_ERROR; end if; if abs DIVISOR ( DIVISOR'LAST ) = 0.0 then raise POLYNOMIAL_ERROR ; -- trying to divide by zero end if ; if DIVIDEND'LAST - DIVISOR'LAST < 0 then return DIVIDEND & COMPLEX_POLYNOMIAL'( DIVIDEND'LAST + 1 .. DIVISOR'LAST - 1 => COMPOSE_FROM_CARTESIAN(0.0,0.0) ) ; end if ; -- REMAINDER := DIVIDEND ; for L in reverse 0 .. DIVIDEND'LAST - DIVISOR'LAST loop QUOTIENT ( L ) := REMAINDER ( REMAINDER_INDEX ) / DIVISOR ( DIVISOR'LAST ) ; TRIAL := ZEROS ; TRIAL ( L ) := QUOTIENT ( L ) ; REMAINDER := REMAINDER - TRIAL * DIVISOR ; REMAINDER_INDEX := REMAINDER_INDEX - 1 ; end loop ; return REMAINDER ( 0 .. DIVISOR'LAST - 1 ) ; end "rem" ; function EVALUATE ( A : COMPLEX_POLYNOMIAL ; X : COMPLEX ) return COMPLEX is Y : COMPLEX ; begin if A'FIRST /= 0 or A'LAST < 0 then raise POLYNOMIAL_ERROR; end if; Y := COMPOSE_FROM_CARTESIAN(0.0,0.0) ; for I in reverse A'RANGE loop Y := Y * X + A ( I ) ; end loop ; return Y ; end EVALUATE ; function INTEGRAL ( A : COMPLEX_POLYNOMIAL ) return COMPLEX_POLYNOMIAL is I : COMPLEX_POLYNOMIAL ( 0 .. A'LAST + 1 ) ; begin if A'FIRST /= 0 or A'LAST < 0 then raise POLYNOMIAL_ERROR; end if; I ( 0 ) := COMPOSE_FROM_CARTESIAN(0.0,0.0) ; for J in 0 .. A'LAST loop I ( J + 1 ) := A ( J ) / COMPOSE_FROM_CARTESIAN( REAL(J+1),0.0) ; end loop ; return I ; end INTEGRAL ; function DERIVITIVE ( A : COMPLEX_POLYNOMIAL ) return COMPLEX_POLYNOMIAL is D : COMPLEX_POLYNOMIAL ( 0 .. A'LAST - 1 ) ; begin if A'FIRST /= 0 or A'LAST < 1 then raise POLYNOMIAL_ERROR; end if; for J in 0 .. A'LAST - 1 loop D ( J ) := A ( J + 1 ) * COMPOSE_FROM_CARTESIAN(REAL(J+1),0.0); end loop ; return D ; end DERIVITIVE ; function ZERO_POLY ( N : INTEGER ) return COMPLEX_POLYNOMIAL is C : COMPLEX_POLYNOMIAL ( 0 .. N ) ; begin if N < 0 then raise POLYNOMIAL_ERROR; end if; for I in 0 .. N loop C ( I ) := COMPOSE_FROM_CARTESIAN(0.0,0.0) ; end loop ; return C ; end ZERO_POLY ; function MAKE_FROM_ROOTS ( V : COMPLEX_VECTOR ) return COMPLEX_POLYNOMIAL is P : COMPLEX_POLYNOMIAL ( 0..V'LENGTH ) := ZERO_POLY(V'LENGTH); PR : COMPLEX_POLYNOMIAL(0..1) := (COMPOSE_FROM_CARTESIAN(0.0,0.0), COMPOSE_FROM_CARTESIAN(1.0,0.0)); begin if V'LENGTH < 1 then raise POLYNOMIAL_ERROR; end if; P(0) := COMPOSE_FROM_CARTESIAN(1.0,0.0); for I in V'RANGE loop PR(0) := -V(I) ; P := P(0..V'LENGTH-1) * PR; end loop ; return P ; -- will have range 0 .. V'LENGTH P(V'LENGTH)=1.0 end MAKE_FROM_ROOTS ; function SOLVE_FOR_ROOTS ( POLY : COMPLEX_POLYNOMIAL ) return COMPLEX_VECTOR is -- ALGORITHM 419 COLLECTED ALGORITHMS FROM ACM. -- ALGORITHM APPEARED IN COMM. ACM, VOL. 15, NO. 02, -- P. 097. from toms netlib.att.com -- converted to Ada using existing Complex_Types and -- Complex_Polynomials packages by Jon Squire 7/28/94 -- some goto's remain. some complex arithmetic is still in line -- from Fortran COMMON statements P, H, QP, QH, SH : complex_polynomial(0..50); S, T, PV : complex; ARE, MRE, ETA, INFIN, SMALNO, BASE : Real; NN : integer; -- set to DEGREE in CPOLY, then changed ZERO : Complex_Vector(0..POLY'last-1) := (others=>Compose_From_Cartesian(0.0,0.0)); FAIL : boolean := false; XX,YY,COSR,SINR,XXX,BND : Real; Z : complex; CONV : boolean; DEGREE : integer; -- internal subprograms used by Solve_For_Roots procedure POLYEV(NN : integer; S : complex; P : complex_polynomial; Q : in out complex_polynomial; PV : in out complex) is -- SUBROUTINE POLYEV(NN,SR,SI,PR,PI,QR,QI,PVR,PVI) -- EVALUATES A POLYNOMIAL P AT S BY THE HORNER RECURRENCE -- PLACING THE PARTIAL SUMS IN Q AND THE COMPUTED VALUE IN PV. begin Q(NN) := P(NN); -- order reversed, now 0..NN old 1 -> NN PV := Q(NN); for I in reverse 0..NN-1 loop PV := PV*S + P(I); Q(I) := PV; end loop; end POLYEV; procedure CALCT(BOOL : in out boolean) is -- SUBROUTINE CALCT(BOOL) -- COMPUTES T := -P(S)/H(S). -- BOOL - LOGICAL, SET TRUE IF H(S) IS ESSENTIALLY ZERO. HV : complex; begin --N := NN-1; -- N is coefficient of X, 1 in Ada -- EVALUATE H(S). --POLYEV(NN-1,S,H,QH,HV); -- not using constant term QH(NN) := H(NN); -- order reversed, now 0..NN old 1 -> NN HV := QH(NN); for I in reverse 1..NN-1 loop HV := HV*S + H(I); QH(I) := HV; end loop; -- special POLYEV that does not use constant term BOOL := (abs HV) <= ARE*10.0*(abs H(1)); if not BOOL then T := -PV/HV; return; end if; T := Compose_From_Cartesian(0.0,0.0); end CALCT; procedure NEXTH(BOOL : boolean) is -- SUBROUTINE NEXTH(BOOL) -- CALCULATES THE NEXT SHIFTED H POLYNOMIAL. -- BOOL - LOGICAL, IF .TRUE. H(S) IS ESSENTIALLY ZERO begin --N := NN-1; -- NN is now DEGREE and order reversed if not BOOL then for J in reverse 1..NN-1 loop -- 10 H(J) := T * QH(J+1) + QP(J); end loop; H(NN) := QP(NN); return; end if; -- IF H(S) IS ZERO REPLACE H WITH QH. for J in reverse 1..NN-1 loop -- 30 H(J) := QH(J+1); end loop; H(NN) := Compose_From_Cartesian(0.0,0.0); end NEXTH; -- H(1) our constant term, not set function ERREV(NN : integer; Q : complex_polynomial; MS : Real; MP : Real) return Real is -- DOUBLE PRECISION FUNCTION ERREV(NN,QR,QI,MS,MP,ARE,MRE) -- BOUNDS THE ERROR IN EVALUATING THE POLYNOMIAL BY THE HORNER -- RECURRENCE. -- QR,QI - THE PARTIAL SUMS -- MS -MODULUS OF THE POINT -- MP -MODULUS OF POLYNOMIAL VALUE -- ARE, MRE -ERROR BOUNDS ON COMPLEX ADDITION AND MULTIPLICATION E : Real; begin E := (abs Q(NN)) * MRE/(ARE+MRE); for I in reverse 0..NN-1 loop E := E*MS+abs Q(I); end loop; --10 E := E*(ARE+MRE)-MP*MRE; return E; end ERREV; function SCALE(NN : INTEGER; PT : complex_polynomial) return Real is HI, LO, MAX, MIN, X, SC : Real; L : integer; begin -- FIND LARGEST AND SMALLEST MODULI OF COEFFICIENTS. HI := SQRT(INFIN); LO := SMALNO/ETA; MAX := 0.0; MIN := INFIN; for I in 0..NN loop X := re(PT(I)); if X > MAX then MAX := X; end if; if X /= 0.0 and X= LO and MAX <= HI then return 1.0; end if; X := LO/MIN; if X <= 1.0 then SC := 1.0/(SQRT(MAX)*SQRT(MIN)); else SC := X; end if; if INFIN/SC > MAX then SC := 1.0; end if; L := integer(LOG(SC)/LOG(BASE)); return BASE**L; end SCALE; procedure CAUCHY(XOUT : out Real; NN : integer; P : in out complex_polynomial) is -- DOUBLE PRECISION FUNCTION CAUCHY(NN,PT,Q) -- CAUCHY COMPUTES A LOWER BOUND ON THE MODULI OF THE ZEROS OF A -- POLYNOMIAL - re P = old PT IS THE MODULUS OF THE COEFFICIENTS. -- old Q is now im P (actually scratch when leaving) X, XM, F, DX, DF : Real; begin P(0) := -P(0); -- order reversed -- COMPUTE UPPER ESTIMATE OF BOUND. -- N := NN-1; -- watch how used, NN is now degree, also order reversed X := EXP( (LOG(-re(P(0))) - LOG(re(P(NN))))/Real(NN) ); if re(P(1)) /= 0.0 then -- IF NEWTON STEP AT THE ORIGIN IS BETTER, USE IT. XM := -re(P(0))/re(P(1)); if XM> XM := X*0.1; F := re(P(NN)); for I in reverse 0..NN-1 loop F := F*XM+re(P(I)); end loop; if F > 0.0 then X := XM; goto L20; end if; DX := X; -- DO NEWTON ITERATION UNTIL X CONVERGES TO TWO DECIMAL PLACES. while abs(DX/X) > 0.005 loop set_im(P(NN),re(P(NN))); for I in reverse 0..NN-1 loop set_im(P(I), im(P(I+1))*X+re(P(I))); end loop; F := im(P(0)); DF := im(P(NN)); for I in reverse 1..NN-1 loop DF := DF*X+im(P(I)); end loop; DX := F/DF; X := X-DX; end loop; XOUT := X; end CAUCHY; procedure NOSHFT(L1 : integer) is -- COMPUTES THE DERIVATIVE POLYNOMIAL AS THE INITIAL H -- POLYNOMIAL AND COMPUTES L1 NO-SHIFT H POLYNOMIALS. XNI : Real; begin -- N = NN-1 NN is now degree, order is reversed for I in 1..NN loop XNI := Real(I); H(I) := XNI*P(I)/Real(NN); end loop; -- 10 for JJ in 1..L1 loop -- DO 50 if abs H(1) > ETA*10.0*(abs P(1)) then T := -P(0)/H(1); for J in 1..NN-1 loop H(J) := T*H(J+1)+P(J); end loop; H(NN) := P(NN); -- IF THE CONSTANT TERM IS ESSENTIALLY ZERO, SHIFT H COEFFICIENTS. else for J in 2..NN loop H(J-1) := H(J); end loop; H(NN) := Compose_From_Cartesian(0.0,0.0); end if; end loop; end NOSHFT; procedure VRSHFT(L3 : integer; Z : in out complex; CONV : out boolean) is -- SUBROUTINE VRSHFT(L3,ZR,ZI,CONV) -- CARRIES OUT THE THIRD STAGE ITERATION. -- L3 - LIMIT OF STEPS IN STAGE 3. -- ZR,ZI - ON ENTRY CONTAINS THE INITIAL ITERATE, IF THE -- ITERATION CONVERGES IT CONTAINS THE FINAL ITERATE ON EXIT. -- CONV - .TRUE. IF ITERATION CONVERGES MP,MS,OMP,RELSTP,R1,R2,TP : Real; B, BOOL : boolean; begin CONV := false; B := false; BOOL := false; S := Z; -- MAIN LOOP FOR STAGE THREE for I in 1..L3 loop -- 60 -- EVALUATE P AT S AND TEST FOR CONVERGENCE. POLYEV(NN,S,P,QP,PV); MP := abs PV; MS := abs S; if MP <= 20.0*ERREV(NN,QP,MS,MP) then -- POLYNOMIAL VALUE IS SMALLER IN VALUE THAN A BOUND ON THE ERROR -- IN EVALUATING P, TERMINATE THE ITERATION. CONV := true; Z := S; return; end if; if I /= 1 then --40 if not (B or MP=0.05) then -- 30 -- ITERATION HAS STALLED. PROBABLY A CLUSTER OF ZEROS. DO 5 FIXED -- SHIFT STEPS INTO THE CLUSTER TO FORCE ONE ZERO TO DOMINATE. TP := RELSTP; B := true; if RELSTP < ETA then TP := ETA; end if; R1 := SQRT(TP); R2 := re(S)*(1.0+R1)-im(S)*R1; S := Compose_From_Cartesian(re(S)*R1+im(S)*(1.0+R1),R2); POLYEV(NN,S,P,QP,PV); for J in 1..5 loop CALCT(BOOL); NEXTH(BOOL); end loop; OMP := INFIN; goto L50; -- EXIT IF POLYNOMIAL VALUE INCREASES SIGNIFICANTLY. end if; -- 30 if MP*0.1 > OMP then return; end if; end if; -- 40 OMP := MP; -- CALCULATE NEXT ITERATE. <> CALCT(BOOL); NEXTH(BOOL); CALCT(BOOL); if not BOOL then RELSTP := (abs T)/(abs S); S := S+T; end if; end loop; -- 60 end VRSHFT; procedure FXSHFT(L2 : integer; Z : in out complex; CONV : in out boolean) is -- SUBROUTINE FXSHFT(L2,ZR,ZI,CONV) -- COMPUTES L2 FIXED-SHIFT H POLYNOMIALS AND TESTS FOR CONVERGENCE. -- INITIATES A VARIABLE-SHIFT ITERATION AND RETURNS WITH THE -- APPROXIMATE ZERO IF SUCCESSFUL. -- L2 - LIMIT OF FIXED SHIFT STEPS -- ZR,ZI - APPROXIMATE ZERO IF CONV IS .TRUE. -- CONV - LOGICAL INDICATING CONVERGENCE OF STAGE 3 ITERATION OT, SVS : complex; TEST, PASD, BOOL : boolean; begin BOOL := false; --N := NN-1; -- NN now degree and polynomial reversed -- EVALUATE P AT S. POLYEV(NN,S,P,QP,PV); -- working on data from COMMON TEST := true; PASD := false; -- CALCULATE FIRST T := -P(S)/H(S). CALCT(BOOL); -- MAIN LOOP FOR ONE SECOND STAGE STEP. for J in 1..L2 loop --50 OT := T; -- COMPUTE NEXT H POLYNOMIAL AND NEW T. NEXTH(BOOL); CALCT(BOOL); Z := S+T; -- TEST FOR CONVERGENCE UNLESS STAGE 3 HAS FAILED ONCE OR THIS -- IS THE LAST H POLYNOMIAL . if not ( BOOL or (not TEST) or (J = L2)) then --50c if abs(T-OT) < 0.5*(abs Z) then --40 if PASD then --30 -- THE WEAK CONVERGENCE TEST HAS BEEN PASSED TWICE, START THE -- THIRD STAGE ITERATION, AFTER SAVING THE CURRENT H POLYNOMIAL -- AND SHIFT. for I in 1..NN loop -- 10 SH(I) := H(I); end loop; -- 10 SVS := S; VRSHFT(10,Z,CONV); if CONV then return; end if; -- THE ITERATION FAILED TO CONVERGE. TURN OFF TESTING AND RESTORE -- H,S,PV AND T. TEST := false; for I in 1..NN loop -- 20 H(I) := SH(I); end loop; -- 20 S := SVS; POLYEV(NN,S,P,QP,PV); CALCT(BOOL); goto L50; end if; --30 PASD := true; goto L50; end if; --40 PASD := false; end if; --50c <> null; end loop; -- 50 on L2 -- ATTEMPT AN ITERATION WITH FINAL H POLYNOMIAL FROM SECOND STAGE. VRSHFT(10,Z,CONV); end FXSHFT; begin -- INITIALIZATION OF CONSTANTS -- ETA THE MAXIMUM RELATIVE REPRESENTATION ERROR -- WHICH CAN BE DESCRIBED AS THE SMALLEST POSITIVE -- FLOATING-POINT NUMBER SUCH THAT 1.0D0 + ETA IS -- GREATER THAN 1.0D0. -- INFIN THE LARGEST FLOATING-POINT NUMBER -- SMALNO THE SMALLEST POSITIVE FLOATING-POINT NUMBER -- BASE THE BASE OF THE FLOATING-POINT NUMBER SYSTEM USED ETA := Real'epsilon; INFIN := Real'last; SMALNO := Real'small; BASE := 2.0; ARE := ETA; MRE := 2.0*sqrt(2.0)*ETA; XX := 0.70710678; YY := -XX; COSR := -0.060756474; SINR := 0.99756405; FAIL := false; DEGREE := POLY'last; T := Compose_From_Cartesian(0.0,0.0); S := Compose_From_Cartesian(0.0,0.0); NN := DEGREE; -- watch out how NN is used, polynomials backward -- ALGORITHM FAILS IF THE LEADING COEFFICIENT IS ZERO. if re(POLY(NN)) = 0.0 and im(POLY(NN)) = 0.0 then FAIL := true; -- raise ??? return ZERO; end if; -- MAKE A COPY OF THE COEFFICIENTS. for I in 0..NN loop P(I) := POLY(I); end loop; -- REMOVE THE ZEROS AT THE ORIGIN IF ANY. while re(P(0)) = 0.0 and im(P(0)) = 0.0 loop ZERO(DEGREE-NN+zero'first) := Compose_From_Cartesian(0.0,0.0); NN := NN-1; for I in 0..NN loop P(I) := P(I+1); -- adjust polynomial in COMMON end loop; end loop; -- MAKE A COPY OF THE COEFFICIENTS. for I in 0..NN loop SH(I) := Compose_From_Cartesian(abs P(I),0.0); -- was CMOD end loop; -- SCALE THE POLYNOMIAL. BND := SCALE (NN,SH); if BND /= 1.0 then for I in 0..NN loop P(I) := BND*P(I); end loop; end if; -- START THE ALGORITHM FOR ONE ZERO . <> if NN = 1 then -- CALCULATE THE FINAL ZERO AND RETURN. ZERO(DEGREE-1+zero'first) := -P(0)/P(1); return ZERO; end if; -- CALCULATE BND, A LOWER BOUND ON THE MODULUS OF THE ZEROS. for I in 0..NN loop SH(I) := Compose_From_Cartesian(abs P(I),0.0); end loop; CAUCHY(BND, NN, SH); -- OUTER LOOP TO CONTROL 2 MAJOR PASSES WITH DIFFERENT SEQUENCES -- OF SHIFTS. for CNT1 in 1..2 loop -- 100 -- FIRST STAGE CALCULATION, NO SHIFT. NOSHFT(5); -- INNER LOOP TO SELECT A SHIFT. for CNT2 in 1..9 loop -- 90 -- SHIFT IS CHOSEN WITH MODULUS BND AND AMPLITUDE ROTATED BY -- 94 DEGREES FROM THE PREVIOUS SHIFT XXX := COSR*XX-SINR*YY; YY := SINR*XX+COSR*YY; XX := XXX; S := Compose_From_Cartesian(BND*XX,BND*YY); -- SECOND STAGE CALCULATION, FIXED SHIFT. FXSHFT(10*CNT2,Z,CONV); if CONV then -- 80 -- THE SECOND STAGE JUMPS DIRECTLY TO THE THIRD STAGE ITERATION. -- IF SUCCESSFUL THE ZERO IS STORED AND THE POLYNOMIAL DEFLATED. ZERO(DEGREE-NN+zero'first) := Z; NN := NN-1; for I in 0..NN loop QP(I) := QP(I+1); -- adjust polynomial in COMMON end loop; for I in 0..NN loop P(I) := QP(I); end loop; goto L40; end if; -- 80 -- IF THE ITERATION IS UNSUCCESSFUL ANOTHER SHIFT IS CHOSEN. end loop; -- IF 9 SHIFTS FAIL, THE OUTER LOOP IS REPEATED WITH ANOTHER -- SEQUENCE OF SHIFTS. end loop; -- THE ZEROFINDER HAS FAILED ON TWO MAJOR PASSES. -- RETURN EMPTY HANDED. FAIL := true; -- raise ?? return ZERO; end SOLVE_FOR_ROOTS; function VECTOR_TO_POLYNOMIAL ( V : COMPLEX_VECTOR ) return COMPLEX_POLYNOMIAL is P : COMPLEX_POLYNOMIAL ( 0..V'LENGTH-1 ) ; begin if V'LENGTH < 1 then raise POLYNOMIAL_ERROR; end if; for I in V'RANGE loop P ( I - V'FIRST ) := V ( I ) ; end loop ; return P ; -- will have range 0 .. V'LENGTH-1 end VECTOR_TO_POLYNOMIAL ; function POLYNOMIAL_TO_VECTOR ( P : COMPLEX_POLYNOMIAL ) return COMPLEX_VECTOR is V : COMPLEX_VECTOR ( P'RANGE ) ; begin if P'FIRST /= 0 then raise POLYNOMIAL_ERROR; end if; for I in P'RANGE loop V ( I ) := P ( I ) ; end loop ; return V ; -- will have range needed by vector end POLYNOMIAL_TO_VECTOR ; function FORCE_POWER ( P : COMPLEX_POLYNOMIAL; N : INTEGER ) return COMPLEX_POLYNOMIAL is PR : COMPLEX_POLYNOMIAL(0..N); begin if P'FIRST /= 0 or N < 0 then raise POLYNOMIAL_ERROR; end if; if N <= P'LAST then PR := P(0..N); -- shorten, users responsibility if upper terms non zero else PR(P'RANGE) := P; PR(P'LAST+1 .. N) := (1..N-P'LAST=>COMPOSE_FROM_CARTESIAN(0.0,0.0)); end if; return PR; end FORCE_POWER; end GENERIC_COMPLEX_POLYNOMIALS ;