code 34173; comment MCA 2405; procedure EQILBR(A, N, EM, D, INT); value N; integer N; array A, EM, D; integer array INT; begin integer I, IM, I1, P, Q, J, T, COUNT, EXPONENT, NI; real C, R, EPS, OMEGA, FACTOR; procedure MOVE(K); value K; integer K; begin real DI; NI:= Q - P; T:= T + 1; if K ^= I then begin ICHCOL(1, N, K, I, A); ICHROW(1, N, K, I, A); DI:= D[I]; D[I]:= D[K]; D[K]:= DI end end MOVE; FACTOR:= 1 / (2 * LN(2)); comment MORE GENERALLY: LN(BASE); EPS:= EM[0]; OMEGA:= 1 / EPS; T:= P:= 1; Q:= NI:= I:= N; COUNT:= (N + 1) * N // 2; for J:= 1 step 1 until N do begin D[J]:= 1; INT[J]:= 0 end; for I:= if I < Q then I + 1 else P while COUNT > 0 and NI > 0 do begin COUNT:= COUNT - 1; IM:= I - 1; I1:= I + 1; C:= SQRT(TAMMAT(P, IM, I, I, A, A) + TAMMAT(I1, Q, I, I, A, A)); R:= SQRT(MATTAM(P, IM, I, I, A, A) + MATTAM(I1, Q, I, I, A, A)); if C * OMEGA <= R * EPS then begin INT[T]:= I; MOVE(P); P:= P + 1 end else if R * OMEGA <= C * EPS then begin INT[T]:= -I; MOVE(Q); Q:= Q - 1 end else begin EXPONENT:= LN(R / C) * FACTOR; if ABS(EXPONENT) > 1 then begin NI:= Q - P; C:= 2 ** EXPONENT; R:= 1 / C; D[I]:= D[I] * C; for J:= 1 step 1 until IM, I1 step 1 until N do begin A[J,I]:= A[J,I] * C; A[I,J]:= A[I,J] * R end end else NI:= NI - 1 end end end EQILBR eop