code 34501;
integer procedure ZERPOL(N, A, EM, RE, IM, D);
value N; integer N; array A, EM, RE, IM, D;
begin integer I, TOTIT, IT, FAIL, START, UP, MAX, GIEX, ITMAX;
real X, Y, NEWF, OLDF, MAXRAD, AE, TOL, H1, H2, LN2;
array F[0 : 5], TRIES[1 : 10];
boolean procedure FUNCTION;
begin integer K, M1, M2;
real P, Q, QSQRT, F01, F02, F03, F11, F12, F13,
F21, F22, F23, STOP;
IT:= IT + 1;
P:= 2 * X; Q:= -(X * X + Y * Y); QSQRT:= SQRT(-Q);
F01:= F11:= F21:= D[0]; F02:= F12:= F22:= 0;
M1:= N - 4; M2:= N - 2;
STOP:= ABS(F01) * 0.8;
for K:= 1 step 1 until M1 do
begin F03:= F02; F02:= F01; F01:= D[K] + P * F02 + Q * F03;
F13:= F12; F12:= F11; F11:= F01 + P * F12 + Q * F13;
F23:= F22; F22:= F21; F21:= F11 + P * F22 + Q * F23;
STOP:= QSQRT * STOP + ABS(F01)
end;
if M1 < 0 then M1:= 0;
for K:= M1 + 1 step 1 until M2 do
begin F03:= F02; F02:= F01; F01:= D[K] + P * F02 + Q * F03;
F13:= F12; F12:= F11; F11:= F01 + P * F12 + Q * F13;
STOP:= QSQRT * STOP + ABS(F01)
end;
if N = 3 then F21:= 0;
F03:= F02; F02:= F01; F01:= D[N - 1] + P * F02 + Q * F03;
F[0]:= D[N] + X * F01 + Q * F02;
F[1]:= Y * F01;
F[2]:= F01 - 2 * F12 * Y * Y;
F[3]:= 2 * Y * (- X * F12 + F11);
F[4]:= 2 * (- X * F12 + F11) - 8 * Y * Y * (- X * F22 + F21);
F[5]:= Y * (6 * F12 - 8 * Y * Y * F22);
STOP:= QSQRT * (QSQRT * STOP + ABS(F01)) + ABS(F[0]);
NEWF:= F02:= COMABS(F[0], F[1]);
FUNCTION:= F02 < (2 * ABS(X * F01) - 8 * (ABS(F[0]) + ABS(F01)
* QSQRT) + 10 * STOP) * TOL * (1 + TOL) ** (4 * N + 3)
boolean procedure CONTROL;
if IT > ITMAX then
begin TOTIT:= TOTIT + IT; FAIL:= 3; goto EXIT end
else if IT = 0 then
begin integer I, H; real H1, SIDE;
MAXRAD:= 0; MAX:= (GIEX - LN(ABS(D[0])) / LN2) / N;
for I:= 1 step 1 until N do
begin H1:= if D[I] = 0 then 0
else EXP(LN(ABS(D[I] / D[0])) / I);
if H1 > MAXRAD then MAXRAD:= H1
end;
for I:= 1 step 1 until N - 1 do
if D[I] ^= 0 then
begin H:= (GIEX - LN(ABS(D[I])) / LN2) / (N - I);
if H < MAX then MAX:= H
end;
MAX:= MAX * LN2 / LN(N);
SIDE:= - D[1] / D[0];
SIDE:= if ABS(SIDE) < TOL then 0 else SIGN(SIDE);
if SIDE = 0 then
begin TRIES[7]:= TRIES[2]:= MAXRAD; TRIES[9]:= -MAXRAD;
TRIES[6]:= TRIES[4]:= TRIES[3]:= MAXRAD / SQRT(2);
TRIES[5]:= -TRIES[3]; TRIES[10]:= TRIES[8]:= TRIES[1]:= 0
end else
begin TRIES[8]:= TRIES[4]:= MAXRAD/ SQRT(2);
TRIES[1]:= SIDE * MAXRAD; TRIES[3]:= TRIES[4] * SIDE;
TRIES[6]:= MAXRAD; TRIES[7]:= -TRIES[3];
TRIES[9]:= -TRIES[1]; TRIES[2]:= TRIES[5]:= TRIES[10]:= 0
end;
if COMABS(X, Y) > 2 * MAXRAD then X:= Y:= 0;
CONTROL:= false
end else
begin if IT > 1 & NEWF >= OLDF then
begin UP:= UP+ 1;
if UP = 5 & START < 5 then
begin START:= START + 1; UP:= 0; X:= TRIES[2 * START - 1];
Y:= TRIES[2 * START]; CONTROL:= false
end else CONTROL:= true
end else CONTROL:= true
procedure DEFLATION;
if X = 0 & Y = 0 then N:= N - 1 else
begin integer I, SPLIT; real H1, H2;
array B[0 : N - 1];
if Y = 0 then
begin N:= N - 1; B[N]:= -D[N + 1] / X;
for I:= 1 step 1 until N do
B[N - I]:= (B[N - I + 1] - D[N - I + 1]) / X;
for I:= 1 step 1 until N do
D[I]:= D[I] + D[I - 1] * X
end else
begin H1:= - 2 * X; H2:= X * X + Y * Y;
N:= N - 2;
B[N]:= D[N + 2] / H2; B[N - 1]:= (D[N + 1] - H1 * B[N]) / H2;
for I:= 2 step 1 until N do
B[N - I]:= (D[N - I + 2] - H1 * B[N - I + 1] - B[N - I + 2])/H2;
D[1]:= D[1] - H1 * D[0];
for I:= 2 step 1 until N do
D[I]:= D[I] - H1 * D[I-1] - H2 * D[I-2]
end;
SPLIT:= N;
H2:= ABS(D[N] - B[N]) / (ABS(D[N]) + ABS(B[N]));
for I:= N - 1 step -1 until 0 do
begin H1:= ABS(D[I]) + ABS(B[I]);
if H1 > TOL then
begin H1:= ABS(D[I] - B[I]) / H1;
if H1 < H2 then begin H2:= H1; SPLIT:= I end
end
end;
for I := SPLIT + 1 step 1 until N do D[I]:= B[I];
D[SPLIT]:= (D[SPLIT] + B[SPLIT]) / 2
end OF DEFLATION;
procedure LAGUERRE;
begin integer M;
real S1RE, S1IM, S2RE, S2IM, DX, DY, H1, H2, H3, H4, H5, H6;
if ABS(F[0]) > ABS(F[1]) then
begin H1:= F[0]; H6:= F[1] / H1; H2:= F[2] + H6 * F[3];
H3:= F[3] - H6 * F[2]; H4:= F[4] + H6 * F[5];
H5:= F[5] - H6 * F[4]; H6:= H6 * F[1] + H1
end else
begin H1:= F[1]; H6:= F[0] / H1; H2:= H6 * F[2] + F[3];
H3:= H6 * F[3] - F[2]; H4:= H6 * F[4] + F[5];
H5:= H6 * F[5] - F[4]; H6:= H6 * F[0] + F[1]
S1RE:= H2 / H6; S1IM:= H3 / H6;
H2:= S1RE * S1RE - S1IM * S1IM; H3:= 2 * S1RE * S1IM;
S2RE:= H2 - H4 / H6; S2IM:= H3 - H5 / H6;
H1:= S2RE * S2RE + S2IM * S2IM;
H1:= if H1 ^= 0 then (S2RE * H2 + S2IM * H3) / H1 else 1;
M:= if H1 >= N - 1 then (if N > 1 then N - 1 else 1)
else if H1 > 1 then H1 else 1;
H1:= (N - M) / M;
COMSQRT(H1 * (N * S2RE - H2), H1 * (N * S2IM - H3), H2, H3);
if S1RE * H2 + S1IM * H3 < 0 then
begin H2:= - H2; H3:= - H3 end;
H2:= S1RE + H2; H3:= S1IM + H3;
H1:= H2 * H2 + H3 * H3;
if H1 = 0 then begin DX:= -N; DY:= N end else
begin DX:= - N * H2 / H1; DY:= N * H3 / H1 end;
H1:= ABS(X) * TOL + AE; H2:= ABS(Y) * TOL+ AE;
if ABS(DX) < H1 & ABS(DY) < H2 then
begin DX:= if DX = 0 then H1 else SIGN(DX) * H1;
DY:= if DY = 0 then H2 else SIGN(DY) * H2
end;
X:= X + DX; Y:= Y + DY;
if COMABS(X, Y) > 2 * MAXRAD then
begin H1:= if ABS(X) > ABS(Y) then ABS(X) else ABS(Y);
H2:= LN(H1) / LN2 + 1 - MAX;
if H2 > 0 then
begin H2:= 2 ** H2; X:= X / H2; Y:= Y / H2 end
end
end OF LAGUERRE;
TOTIT:= IT:= FAIL:= UP:= START:= 0; LN2:= LN(2);
NEWF:= GIANT; AE:= DWARF; GIEX:= LN(NEWF) / LN2 - 40;
TOL:= EM[0]; ITMAX:= EM[1];
for I:= 0 step 1 until N do D[I]:= A[N-I];
if N <= 0 then
begin FAIL:= 1; goto EXIT end
else if D[0] = 0 then
begin FAIL:= 2; goto EXIT end;
for I:= 1 while D[N] = 0 & N > 0 do
begin RE[N]:= IM[N]:= 0; N:= N - 1 end;
for I:= 1 while N > 2 do
begin if CONTROL then LAGUERRE;
OLDF:= NEWF;
if FUNCTION then
begin if Y ^= 0 & ABS(Y) < .1 then
begin real H; H:= Y; Y:= 0;
if ^ FUNCTION then Y:= H
end;
RE[N]:= X; IM[N]:= Y;
if Y ^= 0 then begin RE[N - 1]:= X; IM[N - 1]:= -Y end;
DEFLATION; TOTIT:= TOTIT + IT; UP:= START:= IT:= 0
end
end;
if N = 1 then begin RE[1]:= - D[1] / D[0]; IM[1]:= 0 end
else
begin real H1, H2;
H1:= - 0.5 * D[1] / D[0]; H2:= H1 * H1 - D[2] / D[0];
if H2 >= 0 then
begin RE[2]:= if H1 < 0 then H1 - SQRT(H2)
else H1 + SQRT(H2);
RE[1]:= D[2] / (D[0] * RE[2]);
IM[2]:= IM[1]:= 0
end else
begin RE[2]:= RE[1]:= H1;
IM[2]:= SQRT(-H2); IM[1]:= -IM[2]
end
end; N:= 0;
EXIT: EM[2]:= FAIL; EM[3]:= START; EM[4]:= TOTIT;
for I:= (N-1) "DIV" 2 step -1 until 0 do
begin TOL := D[I]; D[I]:= D[N-I]; D[N-I]:= TOL
end;
ZERPOL:= N
end OF ZERPOL;
eop