Function CABS(A1, A2:real):real; var X,Y:real; begin X:=Abs(A1); Y:=Abs(A2); If X = 0 Then CABS:=Y Else If Y = 0 Then CABS:=X Else If X > Y Then CABS:=X * Sqrt(1 + Sqrt(Y / X)) Else CABS:=Y * Sqrt(1 + Sqrt(X / Y)); End; Function CDIV1(A1, A2, B1, B2:real):real; var R,DEN:real; begin If Abs(B1) >= Abs(B2) Then begin R:=B2 / B1; DEN:=B1 + R * B2; CDIV1:=(A1 + A2 * R) / DEN; end Else begin R:=B1 / B2; DEN:=B2 + R * B1; CDIV1:=(A1 * R + A2) / DEN; end; end; Function CDIV2(A1, A2, B1, B2:real):real; var R,DEN:real; begin If Abs(B1) >= Abs(B2) Then begin R:=B2 / B1; DEN:=B1 + R * B2; CDIV2:=(A2 - A1 * R) / DEN; end Else begin R:=B1 / B2; DEN:=B2 + R * B1; CDIV2:=(A2 * R - A1) / DEN; end; end; Function CSQR1(X, Y:real):real; var U,W,V,R:real; begin If (X = 0) And (Y = 0) Then U:=0 Else begin If Abs(X) >= Abs(Y) Then W:=Sqrt(Abs(X)) * Sqrt(0.5 * (1 + Sqrt(1 + Sqrt(Abs(Y / X))))) Else begin R:=Abs(X / Y); W:=Sqrt(Abs(Y)) * Sqrt(0.5 * (R + Sqrt(1 + Sqrt(R)))); end; If X >= 0 Then begin U:=W; V:=Y / (2 * U); end Else begin If Y >= 0 Then V:=W Else V:=-W; U:=Y / (2 * V); end; end; CSQR1:=U; end; Function CSQR2(X, Y:real):real; var V,W,U,R:real; begin If (X = 0) And (Y = 0) Then V:=0 Else begin If Abs(X) >= Abs(Y) Then W:=Sqrt(Abs(X)) * Sqrt(0.5 * (1 + Sqrt(1 + Sqrt(Abs(Y / X))))) Else begin R:=Abs(X / Y); W:=Sqrt(Abs(Y)) * Sqrt(0.5 * (R + Sqrt(1 + Sqrt(R)))); end; If X >= 0 Then begin U:=W; V:=Y / (2 * U); end Else begin If Y >= 0 Then V:=W Else V:=-W; U:=Y / (2 * V); end; end; CSQR2:=V; end; Procedure LAGUER(A:matrx2;M:integer;var X:array of real; EPS:real;POLISH:boolean); const EPSS = 0.6e-7; MAXIT = 100; var ZERO,B,D,F,G,H:array[0..2] of real; G2,SQ,GP,GM,DX,X1:array[0..2] of real; ITER,J:integer; DXOLD,ERQ,ABX,DUM,DUM1,DUM2,CDX:real; begin ZERO[1]:=0; ZERO[2]:=0; DXOLD:=CABS(X[1], X[2]); For ITER:=1 To MAXIT do begin B[1]:=A[1, M + 1]; B[2]:=A[2, M + 1]; ERQ:=CABS(B[1], X[2]); D[1]:=ZERO[1]; D[2]:=ZERO[2]; F[1]:=ZERO[1]; F[2]:=ZERO[2]; ABX:=CABS(X[1], X[2]); For J:=M DownTo 1 do begin DUM:=X[1] * F[1] - X[2] * F[2] + D[1]; F[2]:=X[2] * F[1] + X[1] * F[2] + D[2]; F[1]:=DUM; DUM:=X[1] * D[1] - X[2] * D[2] + B[1]; D[2]:=X[2] * D[1] + X[1] * D[2] + B[2]; D[1]:=DUM; DUM:=X[1] * B[1] - X[2] * B[2] + A[1, J]; B[2]:=X[2] * B[1] + X[1] * B[2] + A[2, J]; B[1]:=DUM; ERQ:=CABS(B[1], B[2]) + ABX * ERQ; End; ERQ:=EPSS * ERQ; If CABS(B[1], B[2]) <= ERQ Then Exit Else begin G[1]:=CDIV1(D[1], D[2], B[1], B[2]); G[2]:=CDIV2(D[1], D[2], B[1], B[2]); G2[1]:=G[1] * G[1] - G[2] * G[2]; G2[2]:=2 * G[1] * G[2]; H[1]:=G2[1] - 2 * CDIV1(F[1], F[2], B[1], B[2]); H[2]:=G2[2] - 2 * CDIV2(F[1], F[2], B[1], B[2]); DUM1:=(M - 1) * (M * H[1] - G2[1]); DUM2:=(M - 1) * (M * H[2] - G2[2]); SQ[1]:=CSQR1(DUM1, DUM2); SQ[2]:=CSQR2(DUM1, DUM2); GP[1]:=G[1] + SQ[1]; GP[2]:=G[2] + SQ[2]; GM[1]:=G[1] - SQ[1]; GM[2]:=G[2] - SQ[2]; If CABS(GP[1], GP[2]) < CABS(GM[1], GM[2]) Then begin GP[1]:=GM[1]; GP[2]:=GM[2]; end; DX[1]:=CDIV1(M, 0, GP[1], GP[2]); DX[2]:=CDIV2(M, 0, GP[1], GP[2]); end; X1[1]:=X[1] - DX[1]; X1[2]:=X[2] - DX[2]; If (X[1] = X1[1]) And (X[2] = X1[2]) Then Exit; X[1]:=X1[1]; X[2]:=X1[2]; CDX:=CABS(DX[1], DX[2]); DXOLD:=CDX; If Not POLISH Then begin If CDX <= EPS * CABS(X[1], X[2]) Then begin Exit; end; end; end; ShowMessage('too many iterations'); End;