```(* ----------------------------------------------------------------------------
* \$Id: LINALG.mi,v 1.1 1994/03/11 15:21:46 pesch Exp \$
* ----------------------------------------------------------------------------
* This file is part of MAS.
* ----------------------------------------------------------------------------
* Copyright (c) 1993 Universitaet Passau
* ----------------------------------------------------------------------------
* \$Log: LINALG.mi,v \$
* Revision 1.1  1994/03/11  15:21:46  pesch
* Counting real roots of multivariate polynomials, Diplomarbeit F. Lippold
*
* ----------------------------------------------------------------------------
*)

IMPLEMENTATION MODULE LINALG;
(* Linear algebra implementation module *)

(* Import lists and declarations *)

FROM MASSTOR IMPORT LIST, SIL, COMP, INV, ADV, LENGTH, FIRST, LIST1, RED;

FROM SACLIST IMPORT LELT, SUFFIX;

FROM SACI IMPORT ISUM, IQ, ISIGNF, INEG;

FROM MASBIOS IMPORT BLINES, SWRITE;

FROM LINALGI IMPORT IVSPROD, IMPROD;

FROM LINALGRN IMPORT MTRANS;

FROM MASI IMPORT IPROD;

CONST rcsidi = "\$Id: LINALG.mi,v 1.1 1994/03/11 15:21:46 pesch Exp \$";
CONST copyrighti = "Copyright (c) 1993 Universitaet Passau";

(* ---------------------------------------------------------------------- *)
(* Arbitrary domain linear algebra *)

PROCEDURE ADUM(D,n: LIST): LIST;
(* Arbitrary domain unit matrix.
n is an integer. The (n x n) unit matrix of domain D is returned. *)
VAR C,c,i,j,e: LIST;
BEGIN
C := SIL; i := 0;
WHILE i < n DO
c := SIL; i := i + 1; j := 0;
WHILE j < n DO
j := j + 1;
IF i = j THEN e := ADFI(D,1) ELSE e := ADFI(D,0) END;
c := COMP(e,c);
END;
c := INV(c); C := COMP(c,C);
END;
C := INV(C); RETURN(C);

PROCEDURE ADVSPROD(D,A,B: LIST): LIST;
(* Arbitrary domain vector scalar product.
A and B are vectors of the domain D. The arbitrary domain value
C = a1*b1 + ... + an*bn is returned. *)
VAR a,b,c,C: LIST;
BEGIN
WHILE A <> SIL DO
END;
RETURN(C);

PROCEDURE ADVSVPROD(A,b: LIST): LIST;
(* Arbitrary domain vector scalar vector product.
A is an arbitrary domain vector and b is a number of the same domain.
The arbitrary domain vector C = (a1*b, ..., an*b) is returned. *)
VAR a,c,C: LIST;
BEGIN
C := SIL;
WHILE A <> SIL DO
C := COMP(c,C);
END;
C := INV(C); RETURN(C);

PROCEDURE ADVVSUM(A,B: LIST): LIST;
(* Arbitrary domain vector vector sum.
A and B are arbitrary domain vectors. The arbitrary domain vector
C = (a1+b1, ..., an+bn) is returned. *)
VAR a,b,c,C: LIST;
BEGIN
C := SIL;
WHILE A <> SIL DO
C := COMP(c,C);
END;
C := INV(C); RETURN(C);

PROCEDURE ADSMPROD(A,b: LIST): LIST;
(* Arbitrary domain scalar and matrix product.
A is a arbitrary domain matrix. b is a arbitrary domain number.
The arbitrary domain matrix C = A * b is returned. *)
VAR a,c,C: LIST;
BEGIN
IF A = SIL THEN RETURN(A) END;
C := SIL;
WHILE A <> SIL DO
C := COMP(c,C);
END;
C := INV(C); RETURN(C);

PROCEDURE ADMSUM(A,B: LIST): LIST;
(* Arbitrary domain matrix sum.
A and B are arbitrary domain matrices. The arbitrary domain matrix
C = A + B is returned. *)
VAR a,b,c,C: LIST;
BEGIN
IF A = SIL THEN RETURN(B) END;
IF B = SIL THEN RETURN(A) END;
C := SIL;
WHILE A <> SIL DO
C := COMP(c,C);
END;
C := INV(C); RETURN(C);

PROCEDURE ADMPROD(D,A,B: LIST): LIST;
(* Arbitrary domain matrix product.
A and B are matrices of domain D. The matrix C = A * B of domain D is
returned, if the number of columns of A is equal to the number of rows
of B, otherwise the empty matrix is returned. *)
VAR H1,H2,C,C1,a,c,BP,BT,b: LIST;
BEGIN
C := SIL;
IF A = SIL THEN RETURN(C) END;
IF B = SIL THEN RETURN(C) END;
H1 := LENGTH(FIRST(A));
H2 := LENGTH(B);
IF H1 <> H2 THEN RETURN(C) END;
BT:=MTRANS(B);
WHILE A <> SIL DO ADV(A, a, A);
C1:=SIL; BP:=BT;
WHILE BP <> SIL DO ADV(BP,b,BP);
C1:=COMP(c,C1);
END;
C1 := INV(C1); C := COMP(C1,C);
END;
C := INV(C); RETURN(C);

(* Arbitrary domain vector write.
A is an arbitrary domain vector. A is written to the output stream. *)
VAR a : LIST;
BEGIN
SWRITE("(");
WHILE A <> SIL DO
IF A <> SIL THEN SWRITE(",") END;
END;
SWRITE(")");

(*Arbitrary domain matrix write.
A is an arbitrary domain  matrix. A is written to the output stream. *)
VAR a : LIST;
BEGIN
BLINES(0); SWRITE("(");
WHILE A <> SIL DO
IF A <> SIL THEN SWRITE(", "); BLINES(0) END;
END;
SWRITE(")"); BLINES(0);

PROCEDURE ADMTRACE(D,A: LIST): LIST;
(* Arbitrary domain matrix trace.
A is a matrix of domain D. The trace of A is returned. *)
VAR tr,i,a: LIST;
BEGIN
tr := ADFI(D,0); i := 0;
WHILE A <> SIL DO
ADV(A,a,A); i := i + 1;
END;
RETURN(tr);

PROCEDURE ADMPTRACE(D,A,B: LIST): LIST;
(* Arbitrary domain matrix product trace.
A and B are matrices of domain D. The trace of A*B is returned. *)
VAR tr,a,b,H1,H2,BT: LIST;
BEGIN
IF (A = SIL) OR (B = SIL) THEN RETURN(tr) END;
H1 := LENGTH(FIRST(A));
H2 := LENGTH(B);
IF H1 <> H2 THEN RETURN(tr) END;
BT:=MTRANS(B);
WHILE A <> SIL DO
END;
RETURN(tr);

PROCEDURE ADCHARPOL(D,Q: LIST): LIST;
(* Arbitrary domain characteristic polynomial.
Q is a p x p Matrix of domain D. The list al=(a(0),...,a(p)) is created
such that a(i) from D is the coefficient of X^(p-i) in det(XE-Q). *)
VAR L,l,A,B,sl,s,al,a,k,p: LIST;
BEGIN
IF Q = SIL THEN RETURN(al) END;
p := LENGTH(Q);
L := SIL; A := Q;
k := 1;
sl := LIST1(s);
al := SUFFIX(al,a);
WHILE k*k < p DO
k := k+1;
L := COMP(A,L);
sl := COMP(s,sl);
al := SUFFIX(al,a);
END;
B := A;
L := INV(L);
l := L;
WHILE k < p DO
k := k+1;
IF l = SIL THEN
l := L;
ELSE
l := RED(l);
END;
sl := COMP(s,sl);
al := SUFFIX(al,a);
END;
RETURN(al);

PROCEDURE ADSIG(D,Q: LIST): LIST;
(* Arbitrary domain signature.
Q is a symmetric p x p Matrix of domain D. The signature of Q ist returned.
ADCHARPOL is used. *)
VAR al,p,n,a,s,sp,sn,t: LIST;
BEGIN
p := 0; n := 0;
sp := ADSIGN(a); sn := sp; t := 1;
WHILE al <> SIL DO
ADV(al,a,al); t := -t; s := ADSIGN(a);
IF s <> 0 THEN
IF sp <> s   THEN p := p + 1; sp := s   END;
IF sn <> s*t THEN n := n + 1; sn := s*t END;
END;
END;
RETURN(p-n);

(* ---------------------------------------------------------------------- *)
(* Integer linear algebra *)

PROCEDURE IMTRACE(A: LIST): LIST;
(* Integral matrix trace.
A is an integral matrix. The trace of A is returned. *)
VAR tr,i,a: LIST;
BEGIN
tr := 0; i := 0;
WHILE A <> SIL DO
ADV(A,a,A); i := i + 1;
tr := ISUM(tr,LELT(a,i));
END;
RETURN(tr);
END IMTRACE;

PROCEDURE IMPTRACE(A,B: LIST): LIST;
(* Integral matrix product trace.
A and B are integral matrices. The trace of the matrix A*B is returned. *)
VAR tr,a,b,H1,H2,BT: LIST;
BEGIN
tr := 0;
IF (A = SIL) OR (B = SIL) THEN RETURN(tr) END;
H1 := LENGTH(FIRST(A));
H2 := LENGTH(B);
IF H1 <> H2 THEN RETURN(tr) END;
BT:=MTRANS(B);
WHILE A <> SIL DO
tr := ISUM(tr,IVSPROD(a,b));
END;
RETURN(tr);
END IMPTRACE;

PROCEDURE ICHARPOL(Q: LIST): LIST;
(* Integral matrix characteristic polynomial.
Q is an integral p x p Matrix. The list al = (a(0),...,a(p)) of integers
is created with a(i) is the coefficient of X^(p-i) in det(XE-Q). *)
VAR L,l,A,B,sl,s,al,a,k,p: LIST;
BEGIN
al := LIST1(1);
IF Q = SIL THEN RETURN(al) END;
p := LENGTH(Q);
L := SIL; A := Q;
k := 1;
s := IMTRACE(A);
sl := LIST1(s);
a := INEG(s);
al := SUFFIX(al,a);
WHILE k*k < p DO
k := k+1;
L := COMP(A,L);
A := IMPROD(A,Q);
s := IMTRACE(A);
sl := COMP(s,sl);
a := IQ(IVSPROD(sl,al),-k);
al := SUFFIX(al,a);
END;
B := A;
L := INV(L);
l := L;
WHILE k < p DO
k := k+1;
IF l = SIL THEN
B := IMPROD(B,A);
s := IMTRACE(B);
l := L;
ELSE
s := IMPTRACE(B,FIRST(l));
l := RED(l);
END;
sl := COMP(s,sl);
a := IQ(IVSPROD(sl,al),-k);
al := SUFFIX(al,a);
END;
RETURN(al);
END ICHARPOL;

PROCEDURE ISIG(Q: LIST): LIST;
(* Integral matrix signature.
Q is a symmetric integral p x p Matrix. The signature of Q ist returned.
ICHARPOL is used *)
VAR al,p,n,a,s,sp,sn,t: LIST;
BEGIN
p := 0; n := 0;
al := ICHARPOL(Q);
sp := ISIGNF(a); sn := sp; t := 1;
WHILE al <> SIL DO
ADV(al,a,al); t := -t; s := ISIGNF(a);
IF s <> 0 THEN
IF sp <> s   THEN p := p + 1; sp := s   END;
IF sn <> s*t THEN n := n + 1; sn := s*t END;
END;
END;
RETURN(p-n);
END ISIG;

PROCEDURE IMRTPROD(A,B: LIST): LIST;
(* Integral matrix right tensor product.
A and B are integral matrices. The matrix C is constructed by
replacing every entry a(i,j) of A by the matrix a(i,j)*B. *)
VAR C,c,a,ar,ars,Bs,b,br,brs: LIST;
BEGIN
C := SIL;
IF (A = SIL) OR (B = SIL) THEN RETURN(C) END;
WHILE A <> SIL DO
Bs := B;
WHILE Bs <> SIL DO
c := SIL;
ars := ar;
WHILE ars <> SIL DO
brs := br;
WHILE brs <> SIL DO