module froben; % Computation of the frobenius normal form of a matrix. %
% The function frobenius computes the Frobenius normal form F of a
% matrix A, the transformation matrix P and its inverse P^(-1).
%
% Specifically:
%
% - frobenius(A) will return {F,P,Pinv} where F, P, and Pinv
% are such that P*F*Pinv=A.
% Global description of the algorithm:
%
% For a given n by n matrix A over a field K, let L be the linear
% transformation of K^n induced by A. A polycyclic basis of K^n with
% respect to L is a basis of the following form:
% v1,L*v1,.,L^(d1-1)*v1,v2,L*v2,.,L^(d2-1)*v2,.,vr,L*vr,., L^(dr-1)*vr
% such that v1,L*v1,..,L^(d1-1)*v1,..,vi,L*vi,..,L^(di-1)*vi,L^di*vi
% are linearly dependent for i=1..r.
% It is easy to see that the matrix B of L with respect to a polycyclic
% basis is of the form plist_to_polycompanion(plist,x), where plist is
% a list of monic elements of K[x] of strictly increasing degree (for
% a description of plist_to_polycompanion see below).
% The computation of a polycyclic basis of K^n and the transformation
% matrix from A to B is performed in the function cyclic_vectors.
% Next we view K^n as a K[x]-module via x*v=B*v. Suppose that
% B=plist_to_polycompanion(plist,x), where plist=[p1,..,pr] and
% degree(pi)=di. Let G be the r by r upper triangular matrix such that
% G(i,j) satisfies:
% pj=G(1,j)+G(2,j)*x^d1+G(3,j)*x^d2+..+G(j,j)*x^d(j-1),
% where degree(G(j,j))=dj-d(j-1) and degree(G(i,j))<di-d(i-1) (d0=0).
% Let R be the K[x]-submodule of K[x]^r generated by the columns of G.
% Representants for the elements of the quotient module K[x]^r/R are
% the vectors [L1,L2,..,Lr] where degree(Li)<di-d(i-1). By taking the
% coefficients of the Li the quotient module is identified with K^n.
% The multiplication by x on the quotient module is identified with
% the multiplication by B on K^n.
% Next we compute the Smith normal form S of G. Say L*S*R=G. If R' is
% the K[x]-submodule of K[x]^r generated by the columns of S we get
% the following diagram:
%
% ~ ~ ~
% K^n <------- K[x]^r/R' -------> K[x]^r/R -------> K^n
% L
% | | | |
% | | | |
% |F |x |x |B
% | | | |
% | | | |
% \ / \ / \ / \ /
% ~ ~ ~
% K^n <------- K[x]^r/R' -------> K[x]^r/R -------> K^n
% L
%
% Here F is in Frobenius normal form and thus it is the Frobenius
% normal form of B (and thus of A). The computation of the Smith
% normal form of G is performed in the function cyclic_to_frobenius.
symbolic procedure frobenius(A);
begin
scalar AA,P,Pinv,ans,tmp,full_coeff_list,rule_list,input_mode;
matrix_input_test(A);
if (car size_of_matrix(A)) neq (cadr size_of_matrix(A))
then rederr "ERROR: expecting a square matrix. ";
input_mode := get(dmode!*,'dname);
%
% If modular or arnum are on then we keep them on else we want
% rational on.
%
if input_mode neq 'modular and input_mode neq 'arnum and
input_mode neq 'rational then on rational;
on combineexpt;
tmp := nest_input(A);
AA := car tmp;
full_coeff_list := cadr tmp;
tmp := frobeniusform(AA);
ans := car tmp;
P := cadr tmp;
Pinv := caddr tmp;
%
% Set up rule list for removing nests.
%
rule_list := {'co,2,{'~,'int}}=>'int when numberp(int);
for each elt in full_coeff_list do
<<
tmp := {'co,2,{'~,elt}}=>elt;
rule_list := append (tmp,rule_list);
>>;
%
% Remove nests.
%
let rule_list;
ans := de_nest_mat(ans);
P := de_nest_mat(P);
Pinv := de_nest_mat(Pinv);
clearrules rule_list;
%
% Return to original mode.
%
if input_mode neq 'modular and input_mode neq 'arnum and
input_mode neq 'rational then
<<
% onoff('nil,t) doesn't work so ...
if input_mode = 'nil then off rational
else onoff(input_mode,t);
>>;
off combineexpt;
return {'list,ans,P,Pinv};
end;
flag ('(frobenius),'opfn); % So it can be used from algebraic mode.
symbolic procedure frobeniusform(A);
begin
scalar ans,plist,tmp,P,Pinv,inv_fact,T1,Tinv,V,Vinv,x;
x := mkid('x,0);
tmp := cyclic_vectors(A,x);
plist := car tmp;
V := cadr tmp;
Vinv := caddr tmp;
tmp := cyclic_to_frobenius(plist,x);
inv_fact := car tmp;
T1 := cadr tmp;
Tinv := caddr tmp;
P:= reval {'times,V,T1};
Pinv:= reval {'times,Tinv,Vinv};
ans := invfact_to_frobenius(inv_fact,x);
return {ans,P,Pinv};
end;
symbolic procedure basis(n,i);
%
% Basis creates an element of the natural basis of a vector space.
%
begin
scalar VV;
VV := mkmatrix(1,n);
setmat(VV,1,i,1);
return VV;
end;
symbolic procedure calc_exgcd(poly1,poly2,x);
%
% Extended Euclidean algorithm for polynomials.
% poly1, and poly2 are polynomials in x.
% Returns gcd, s1, and t1 such that s1 * poly1 + t1 * poly2 = gcd,
% with degree(s1,x)<degree(poly2,x) and degree(t1,x)<degree(poly1,x).
%
begin
scalar gcd,c,c1,c2,d,d1,d2,q,r,r1,r2,s1,t1;
if poly1 = 0 and poly2 = 0 then return {0,0,0} else
<<
poly1 := reval poly1;
poly2 := reval poly2;
c := reval norm(poly1,x);
d := reval norm(poly2,x);
c1 := 1;
d1 := 0;
c2 := 0;
d2 := 1;
while reval d neq 0 do
<<
q := reval get_quo(c,d);
r := reval {'plus,c,{'minus,{'times,q,d}}};
r1 := reval {'plus,c1,{'minus,{'times,q,d1}}};
r2 := reval {'plus,c2,{'minus,{'times,q,d2}}};
c := reval d;
c1 := reval d1;
c2 := reval d2;
d := reval r;
d1 := reval r1;
d2 := reval r2;
>>;
gcd := reval norm(c,x);
s1 := reval{'quotient,c1,{'times,unit(poly1,x),unit(c,x)}};
t1 := reval{'quotient,c2,{'times,unit(poly2,x),unit(c,x)}};
return {gcd,s1,t1};
>>;
end;
symbolic procedure norm(poly,x);
begin
scalar normal;
if poly = 0 then normal := 0
else if lcof(poly,x) = 0 then normal := 1
else normal := reval{'quotient,poly,lcof(poly,x)};
return normal;
end;
symbolic procedure unit(poly,x);
begin
scalar unit1;
if poly = 0 then unit1 := 1
else if lcof(poly,x) = 0 then unit1 := poly
else unit1 := reval lcof(poly,x);
return unit1;
end;
symbolic procedure companion(poly,x);
%
% Takes as input a monic univariate polynomial in a variable x.
% Returns a companion matrix associated with the polynomial poly(x).
% If C := companion(p,x) and p is a0+a1*x+...+x^n (a univariate monic
% polynomial), them C(i,n) = -coeff(p,x,i-1), C(i,i-1) = 1 (i=2..n)
% and C(i,j) = 0 for all other i and j.
%
% Can be used independently from algebraic mode.
%
begin
scalar mat1;
integer n;
n := deg(poly,x);
if de_nest(reval coeffn(poly,x,n)) neq 1
then rederr {"ERROR: polynomial",poly," is not monic."};
mat1 := mkmatrix(n,n);
setmat(mat1,1,n,{'minus,coeffn(poly,x,0)});
for i:=2:n do
<<
setmat(mat1,i,i-1,1);
>>;
for j:=2:n do
<<
setmat(mat1,j,n,{'minus,coeffn(poly,x,j-1)});
>>;
return mat1;
end;
flag('(companion),'opfn); % So it can be used independently from
% algebraic mode.
symbolic procedure compute_g(r,dd,plist,x);
begin
scalar G,tmp,new_elt;
G := mkmatrix(r,r);
for j:=1:r do
<<
for i:=1:j-1 do
<<
new_elt := 0;
for k:=getmat(dd,1,i):getmat(dd,1,i+1)-1 do
<<
tmp := {'times,coeffn(nth(plist,j),x,k),{'expt,x,{'plus,k,
{'minus,getmat(dd,1,i)}}}};
new_elt := {'plus,new_elt,tmp};
>>;
setmat(G,i,j,new_elt);
>>;
new_elt := 0;
for k:=getmat(dd,1,j):getmat(dd,1,j+1) do
<<
tmp := {'times,coeffn(nth(plist,j),x,k),{'expt,x,{'plus,k,
{'minus,getmat(dd,1,j)}}}};
new_elt := {'plus,new_elt,tmp};
>>;
setmat(G,j,j,new_elt);
>>;
return G;
end;
symbolic procedure copy_mat(A);
%
% Creates a copy of the input and returns it;
%
begin
scalar C;
integer row_dim,col_dim;
matrix_input_test(A);
row_dim := car size_of_matrix(A);
col_dim := cadr size_of_matrix(A);
C := mkmatrix(row_dim,col_dim);
for i:=1:row_dim do
<<
for j:=1:col_dim do
<<
setmat(C,i,j,getmat(A,i,j));
>>;
>>;
return C;
end;
symbolic procedure cyclic_to_frobenius(plist,x);
%
% A matrix B=plist_to_polycompanion(plist,x) is transformed to its
% Frobenius normal form F.
% If F=diag(C1,..,Cr) where Ci is the companion matrix associated with
% pi, then cyclic_to_frobenius will return {p1,..,pr}.
% Let G be the matrix as described before. We compute the Smith normal
% form S of G. Then S=diag(p1,..,pr), where pi in K[x] such that pi
% pi divides p(i+1) for i=1..(r-1), and
% F=invfact_to_frobenius({p1,..,pr},x) is the frobenius normal form of
% B (for description of invfact_to_frobenius see invfact_to_frobenius)
% .
% Remark: to compute the smith normal form of G we car simplify G
% using the fact that G is upper triangular. Then we use a modified
% version of smithex.
begin
scalar dd,D,US,S,G,C,T1,Tinv,inv_fact,L,Linv,columnT,rowT,rr,q,
columnTinv,rowTinv,tmp,tmp1;
integer r,n;
r := length plist;
dd := mkmatrix(1,r+1);
for j:=1:r do
<<
setmat(dd,1,j+1,deg(nth(plist,j),x));
>>;
n:= getmat(dd,1,r+1);
%%%%%%%%%%%%%%%%%%%
% Compute matrix G.
%%%%%%%%%%%%%%%%%%%
G:=compute_g(r,dd,plist,x);
%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%
% Compute smith normal form of G.
%%%%%%%%%%%%%%%%%%%
tmp := uppersmith(G,x);
US:=car tmp;
L := cadr tmp;
Linv := caddr tmp;
tmp:=mysmith(US,L,Linv,x);
S:=car tmp;
L := cadr tmp;
Linv := caddr tmp;
%%%%%%%%%%%%%%%%%%%
D := mkmatrix(1,r);
for i:=1:r do
<<
setmat(D,1,i,deg(getmat(S,i,i),x));
>>;
%%%%%%%%%%%%%%%%%%%
% Compute transformation matrix.
%%%%%%%%%%%%%%%%%%%
C := mkmatrix(1,r);
T1 := mkmatrix(n,n);
columnT:=0;
for i:=1:r do
<<
for k:=1:r do
<<
setmat(C,1,k,getmat(L,k,i));
>>;
for j:=1:getmat(D,1,i) do
<<
columnT:=columnT+1;
for ii:=r step -1 until 1 do
<<
q:=get_quo(getmat(C,1,ii),getmat(G,ii,ii));
rr:=get_rem(getmat(C,1,ii),getmat(G,ii,ii));
setmat(C,1,ii,rr);
for jj:=1:(ii-1) do
<<
setmat(C,1,jj,reval {'plus,reval getmat(C,1,jj),{'times,
{'minus,q},reval getmat(G,jj,ii)}});
>>;
>>;
rowT:=0;
for ii:=1:r do
<<
tmp := reval{'plus,getmat(dd,1,ii+1),{'minus,
getmat(dd,1,ii)}};
for jj:=1:tmp do
<<
rowT:=rowT+1;
tmp1 := coeffn(getmat(C,1,ii),x,jj-1);
setmat(T1,rowT,columnT,tmp1);
>>;
>>;
for ii:=1:r do setmat(C,1,ii,{'times,getmat(C,1,ii),x});
>>;
>>;
%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%
% Compute inverse transformation matrix
%%%%%%%%%%%%%%%%%%%
<<
Tinv := mkmatrix(n,n);
columnTinv:=0;
for i:=1:r do
<<
for k:=1:r do setmat(C,1,k,getmat(Linv,k,i));
for j:=1:reval {'plus,getmat(dd,1,i+1),{'minus,
getmat(dd,1,i)}} do
<<
columnTinv:=columnTinv+1;
rowTinv:=0;
for ii:=1:r do
<<
setmat(C,1,ii,get_rem(getmat(C,1,ii),getmat(S,ii,ii)));
for jj:=1:reval getmat(D,1,ii) do
<<
rowTinv:=rowTinv+1;
setmat(Tinv,rowTinv,columnTinv,reval
coeffn(getmat(C,1,ii),x,jj-1));
>>;
>>;
for ii:=1:r do setmat(C,1,ii,{'times,getmat(C,1,ii),x});
>>;
>>;
>>;
%%%%%%%%%%%%%%%%%%%
inv_fact := {};
for i:=1:r do
<<
if getmat(D,1,i)>0 then
<<
inv_fact := append(inv_fact,{getmat(S,i,i)});
>>;
>>;
return {inv_fact,T1,Tinv};
end;
symbolic procedure cyclic_vectors(A,x);
%
% cyclic_vectors computes a polycyclic basis of K^n with respect to A.
% If this basis is (b1,..,bn)=
% (v1,A*v1,..,A^(d1-1)*v1,v2,A*v2,.,A*(d2-1)*v2,..,vr,A*vr,..,A^(dr-1)
% *vr) and a1*b1+..+a(d1+..+di)*b(d1+..+di)+A^di*vi=0 we set
% pi=a1+a2*x+..+a(d1+..+di)*x^(d1+..+di-1)+x^(d1+..di).
% cyclic_vectors returns {p1,..,pr}.
% The matrix of A on this basis (b1,..,bn) is
% plist_to_polycompanion([p1,..,pr],x).
%
begin
scalar V,Vinv,plist,U,Uinv,S,carrier,lincomb,VV,UU,SS,l,car,c,
tmp,ans,q,break;
integer n,r;
n := car size_of_matrix(A);
U := mkmatrix(n,n);
S := mkmatrix(n,n);
plist := {};
V := mkmatrix(n,n);
Vinv := mkmatrix(n,n);
carrier := mkvect(n);
lincomb := mkvect(n);
r := 0; % No. of elements already computed.
while r<n do
<<
%%%%%%%%%%%%%%%%%%%
% Start new cycle.
%%%%%%%%%%%%%%%%%%%
q := 1;
while getv(carrier,q) neq nil do q:=q+1; % Find car gap.
VV := basis(n,q);
%%%%%%%%%%%%%%%%%%%
break := nil; % Controls break out of loop.
while not break do
<<
UU := copy_mat(VV);
for i:=1:n do putv(lincomb,i,0);
% Always VV=UU+U*lincomb.
for i:=1:n do
<<
car:= getv(carrier,i);
if car neq nil and getmat(UU,1,i) neq 0 then
<<
c := {'quotient,getmat(UU,1,i),getmat(U,i,car)};
setmat(UU,1,i,0);
for j:=i+1:n do
<<
tmp := {'times,c,getmat(U,j,car)};
setmat(UU,1,j,reval {'plus,getmat(UU,1,j),{'minus,{'times,
c,getmat(U,j,car)}}});
>>;
putv(lincomb,car,c);
>>;
>>;
q := 1;
while q<=n and reval getmat(UU,1,q)=0 do
<<
q:=q+1;
>>;
if q<=n then
<<
% New element of basis.
r:=r+1;
putv(carrier,q,r); % This basis-element carries coordinates q.
% Always U=V*S.
for j:=q:n do setmat(U,j,r,getmat(UU,1,j));
for j:=1:n do setmat(V,j,r,getmat(VV,1,j));
for j:=1:r-1 do
<<
tmp:=getv(lincomb,j);
for l:=j+1:r-1 do tmp:={'plus,tmp,{'times,getmat(S,j,l),
getv(lincomb,l)}};
setmat(S,j,r,{'minus,tmp});
>>;
setmat(S,r,r,1);
% Compute A*VV.
for i:=1:n do
<<
tmp:=0;
for j:=1:n do
<<
tmp:=reval{'plus,tmp,reval{'times,reval getmat(A,i,j),
reval getmat(VV,1,j)}};
>>;
setmat(UU,1,i,tmp);
>>;
for i:=1: cadr size_of_matrix(UU) do
<<
setmat(VV,1,i,getmat(UU,1,i));
>>;
>>
else
<<
break := t;
>>;
>>;
%%%%%%%%%%%%%%%%%%%
% New cycle found
%%%%%%%%%%%%%%%%%%%
SS := mkmatrix(1,r);
for j:=1:r do
<<
tmp:=reval getv(lincomb,j);
for l:=j+1:r do
<<
tmp:=reval{'plus,tmp,{'times,reval getmat(S,j,l),
reval getv(lincomb,l)}};
>>;
setmat(SS,1,j,tmp);
>>;
ans := nil;
for j:=1:r do
<<
tmp := {'times,getmat(SS,1,r+1-j),{'expt,x,r-j}};
ans := reval{'plus,ans,tmp};
>>;
tmp := reval{'plus,{'expt,x,r},{'minus,ans}};
plist := append(plist,{tmp});
%%%%%%%%%%%%%%%%%%%
>>; % End while r<n.
Uinv := inv(U,carrier);
for i:=1:n do
<<
for j:=1:n do
<<
tmp:=0;
for l:=i:n do
<<
tmp:=reval {'plus,tmp,{'times,reval getmat(S,i,l),
reval getmat(Uinv,l,j)}};
>>;
setmat(Vinv,i,j,tmp);
>>;
>>;
return {plist,V,Vinv};
end;
symbolic procedure de_nest(input);
%
% Takes simple nested input and de-nests it.
%
begin
scalar output;
if atom input then output := input
else if car input neq 'co then output := input
else output := caddr input;
return output;
end;
symbolic procedure de_nest_mat(mat1);
%
% Removes nests from each elt of input matrix.
% Rules being applied from outside.
%
begin
integer row_dim,col_dim;
row_dim := car size_of_matrix(mat1);
col_dim := cadr size_of_matrix(mat1);
for i:=1:row_dim do
<<
for j:=1:col_dim do
<<
setmat(mat1,i,j,getmat(mat1,i,j));
>>;
>>;
return mat1;
end;
% Allow variable input.
put('diagi,'psopfn,'diag);
symbolic procedure diag(uu);
%
% Takes square or scalar matrix entries and creates a matrix with
% these matrices on the diagonal.
%
begin
scalar bigA,arg,input,u;
integer nargs,n,Aidx,stp,bigsize,smallsize;
u := car uu;
input := u;
bigsize:=0;
nargs:=length input;
for i:=1:nargs do
<<
arg:=car input;
% If scalar entry.
if algebraic length(arg) = 1 then bigsize:=bigsize+1
else
<<
bigsize:=bigsize+car size_of_matrix(arg);
>>;
input := cdr input;
>>;
bigA := mkmatrix(bigsize,bigsize);
Aidx:=1;
input := u;
for k:=1:nargs do
<<
arg:=car input;
% If scalar entry.
if algebraic length(arg) = 1 then
<<
setmat(bigA,Aidx,Aidx,arg);
Aidx:=Aidx+1;
input := cdr input;
>>
else
<<
smallsize:= car size_of_matrix(arg);
stp:=smallsize+Aidx-1;
for i:=Aidx:stp do
<<
for j:=Aidx:stp do
<<
arg:=car input;
% Find (i-Aidx+1)'th row.
arg := cdr arg;
<<
n:=1;
while n < (i-Aidx+1) do
<<
arg := cdr arg;
n:=n+1;
>>;
>>;
arg := car arg;
%
% Find (j-Aidx+1)'th column elt of i'th row.
%
<<
n:=1;
while n < (j-Aidx+1) do
<<
arg := cdr arg;
n:=n+1;
>>;
>>;
arg := car arg;
setmat(bigA,i,j,arg);
>>;
>>;
Aidx := Aidx+smallsize;
input := cdr input;
>>;
>>;
return biga;
end;
symbolic procedure get_coeffs(poly);
%
% Gets all kernels in a poly.
%
begin
scalar ker_list_num,ker_list_den;
ker_list_num := kernels !*q2f simp reval num poly;
ker_list_den := kernels !*q2f simp reval den poly;
ker_list_num := union(ker_list_num,ker_list_den);
return ker_list_num;
end;
symbolic procedure get_quo(poly1,poly2);
%
% Gets quotient of two polys.
%
begin
scalar quo1,input1,input2;
if input1 = 0 and input2 = 0 then return 0
else
<<
input1 := reval poly1;
input2 := reval poly2;
algebraic (quo1 := (input1-remainder(input1,input2))/input2);
quo1 := reval quo1;
return quo1;
>>;
end;
symbolic procedure get_rem(poly1,poly2);
%
% Gets remainder of two polys.
%
begin
scalar rem1,input1,input2;
input1 := reval poly1;
input2 := reval poly2;
algebraic (rem1 := remainder(input1,input2));
rem1 := reval rem1;
return rem1;
end;
symbolic procedure inv(U,carrier);
%
% inv computes the inverse of a permuted upper triangular matrix. The
% permutation is given by carrier.
%
begin
scalar Uinv,tmp;
integer n;
n:= car size_of_matrix(U);
Uinv := mkmatrix(n,n);
for i:=1:n do
<<
for j:=1:i-1 do
<<
tmp:=0;
for k:=j:i-1 do
<<
tmp := {'plus,tmp,{'times,getmat(U,i,getv(carrier,k)),
getmat(Uinv,getv(carrier,k),j)}};
>>;
setmat(Uinv,getv(carrier,i),j,{'quotient,{'minus,tmp},
getmat(U,i,getv(carrier,i))});
>>;
setmat(Uinv,getv(carrier,i),i,{'quotient,1,getmat(U,i,
getV(carrier,i))});
for j:=i+1:n do setmat(Uinv,getv(carrier,i),j,0);
>>;
return Uinv;
end;
symbolic procedure invfact_to_frobenius(inv_fact,x);
%
% For plist={p1,..,pr] where pi is a monic polynomial in x,
% invfact_to_frobenius(plist,x) makes a square matrix with diagonal
% blocks C1,..,Cr where Ci is the companion matrix to pi.
%
begin
scalar diag_mat,tmp;
integer num;
num := length inv_fact;
tmp:=for i:=1:num collect companion(nth(inv_fact,i),x);
diag_mat := reval{'diagi, tmp};
return diag_mat;
end;
symbolic procedure make_identity(row_dim,col_dim);
%
% Creates identity matrix.
%
begin
scalar A;
A := mkmatrix(row_dim,col_dim);
for i:=1:row_dim do
<<
for j:=1:col_dim do
<<
if i=j then setmat(A,i,i,1);
>>
>>;
return A;
end;
symbolic procedure matrix_input_test(A);
begin
if not eqcar(A,'mat) then rederr
{"ERROR: `",A,"' is non matrix input."}
else return A;
end;
symbolic procedure mysmith(US,L,Linv,x);
%
% The Smith normal form S of a matrix US is computed. L and Linv are
% also computed where L*S*R=US.
% For description of mysmith see smithex.
%
begin
scalar S,a,b,g,jj,s1,t1,tmp,isclear,q,lc,poly1,poly2,input1,input2;
integer n,r;
n:= car size_of_matrix(US);
S := copy_mat(US);
for k:=1:n do
<<
isclear := nil;
while not isclear do
<<
for i:= k+1:n do
<<
if getmat(S,i,k) = 0 then <<>>
else
<<
poly1 := getmat(S,k,k);
poly2 := getmat(S,i,k);
tmp := calc_exgcd(poly1,poly2,x);
g := car tmp;
s1 := cadr tmp ;
t1 := caddr tmp ;
a := get_quo(poly1,g);
b := get_quo(poly2,g);
for j:=k+1:n do
<<
input1 := getmat(S,k,j);
input2 := getmat(S,i,j);
tmp := {'plus,{'times,s1,input1},{'times,t1,input2}};
setmat(S,i,j,{'plus,{'times,a,input2},{'minus,{'times,b,
input1}}});
setmat(S,k,j,tmp);
>>;
for j:=1:n do
<<
tmp := reval{'plus,{'times,a,getmat(L,j,k)},{'times,b,
getmat(L,j,i)}};
setmat (L,j,i,reval{'plus,{'times,{'minus,t1},
getmat(L,j,k)},{'times,s1,getmat(L,j,i)}});
setmat (L,j,k,tmp);
>>;
for j:=1:n do
<<
tmp := reval{'plus,{'times,s1,getmat(Linv,k,j)},
{'times,t1,getmat(Linv,i,j)}};
setmat (Linv,i,j,reval{'plus,{'times,a,getmat(Linv,i,j)},
{'times,{'minus,b},getmat(Linv,k,j)}});
setmat (Linv,k,j,tmp);
>>;
setmat(S,k,k,g);
setmat(S,i,k,0);
>>;
>>;
isclear := t;
for i:=k+1:n do
<<
poly1:=getmat(S,k,i);
poly2:=getmat(S,k,k);
setmat(S,k,i,get_rem(poly1,poly2));
q := get_quo(poly1,poly2);
>>;
for i:=k+1:n do
<<
if getmat(S,k,i) = 0 then <<>>
else
<<
poly1:=getmat(S,k,k);
poly2:=getmat(S,k,i);
tmp := calc_exgcd(poly1,poly2,x);
g:= car tmp;
s1 := cadr tmp;
t1 := caddr tmp;
a:=get_quo(poly1,g);
b:=get_quo(poly2,g);
for j:=k+1:n do
<<
input1 := getmat(S,j,k);
input2 := getmat(S,j,i);
tmp := {'plus,{'times,s1,input1},{'times,t1,input2}};
setmat(S,j,i,{'plus,{'times,a,input2},{'minus,{'times,b,
input1}}});
setmat(S,j,k,tmp);
>>;
setmat(S,k,k,g);
setmat(S,k,i,0);
isclear := nil;
>>;
>>;
>>;
>>;
r:=0;
for i:=1:n do
<<
if getmat(S,i,i) neq 0 then
<<
r:=r+1;
% Watch out for integers giving lc = 0.
if off_mod_lcof(getmat(S,i,i),x) = 0 then lc := getmat(S,i,i)
else lc := off_mod_lcof(getmat(S,i,i),x);
setmat(S,r,r,{'quotient,getmat(S,i,i),lc});
if i neq r then
<<
setmat(S,i,i,0);
for j:=1:n do
<<
tmp := reval getmat(L,j,r);
setmat(L,j,r,reval getmat(L,i,j));
setmat(L,j,i,tmp);
>>;
for j:=1:n do
<<
tmp := reval getmat(Linv,r,j);
setmat(Linv,r,j,reval getmat(Linv,i,j));
setmat(Linv,i,j,tmp);
>>;
>>;
>>;
>>;
for i:=1:r-1 do
<<
jj:=i+1;
<<
while reval getmat(S,i,i) neq 1 and jj <= r do
<<
poly1:=reval getmat(S,i,i);
poly2:=reval getmat(S,jj,jj);
tmp := calc_exgcd(poly1,poly2,x);
g:= car tmp;
s1 := cadr tmp;
t1 := caddr tmp;
a:=get_quo(poly1,g);
b:=get_quo(poly2,g);
setmat(S,i,i,g);
setmat(S,jj,jj,{'times,a,poly2});
for k:=1:n do
<<
tmp := reval {'plus,{'times,a,getmat(L,k,i)},{'times,b,
getmat(L,k,jj)}};
setmat (L,k,jj,reval {'plus,{'times,{'minus,t1},
getmat(L,k,i)},{'times,s1,getmat(L,k,jj)}});
setmat (L,k,i,tmp);
>>;
for k:=1:n do
<<
tmp := reval {'plus,{'times,s1,getmat(Linv,i,k)},{'times,t1,
getmat(Linv,jj,k)}};
setmat(Linv,jj,k,reval {'plus,{'times,a,getmat(Linv,jj,k)},
{'times,{'minus,b},getmat(Linv,i,k)}});
setmat(Linv,i,k,tmp);
>>;
jj:=jj+1;
>>;
>>;
>>;
return {S,L,Linv};
end;
symbolic procedure nest_input(A);
%
% Takes a matrix and converts all elements into nested form.
% Also finds union of all coefficients in all elements and
% returns them in a list, along with the new matrix.
%
begin
scalar tmp,coeff_list,full_coeff_list,AA;
integer row_dim,col_dim;
full_coeff_list := nil;
coeff_list := nil;
AA := copy_mat(A);
row_dim := car size_of_matrix(AA);
col_dim := cadr size_of_matrix(AA);
for i := 1:row_dim do
<<
for j := 1:col_dim do
<<
coeff_list := get_coeffs(getmat(AA,i,j));
if coeff_list = nil then <<>>
else full_coeff_list := union(coeff_list,full_coeff_list);
for each elt in coeff_list do
<<
tmp := {'co,2,elt};
setmat(AA,i,j,algebraic (sub(elt=tmp,getmat(AA,i,j))));
>>;
>>;
>>;
return {AA,full_coeff_list};
end;
symbolic procedure off_mod_lcof(input,x);
begin
if !*modular then
<<
off modular;
input := lcof (input,x);
on modular;
>>
else input := lcof (input,x);
return input;
end;
symbolic procedure off_mod_reval(input);
%
% In certain cases it is required to reval with modular off,
% eg: when calculating degrees of polys.
%
begin
if !*modular then
<<
off modular;
input := reval input;
on modular;
>>
else input := reval input;
return input;
end;
flag('(off_mod_reval),'opfn); % So it can be used from
% algebraic mode.
symbolic procedure plist_to_polycompanion(plist,x);
%
% This is not used.
% It is here to help explain what's going on.
%
% If a=a0+a1*x+x^2, b=b0+b1*x+b2*x^2+x^3 and
% c=c0+c1*x+c2*x^2+c3*x^3+c4*x^4+x^5, then
% plist_to_polycompanion({a,b,c},x) yields
%
% [ 0 -a0 -b0 0 -c0 ]
% [ ]
% [ 1 -a1 -b1 0 -c1 ]
% [ ]
% [ 0 0 -b2 0 -c2 ]
% [ ]
% [ 0 0 0 0 -c3 ]
% [ ]
% [ 0 0 0 1 -c4 ]
%
begin
scalar d,A;
integer r,n;
r := length plist;
d := mkvect(r);
putv(d,0,0);
for i:=1:r do putv(d,i,deg(nth(plist,i),x));
n := getv(d,r);
A := mkmatrix(n,n);
for i:=1:r do
<<
for j:=getv(d,i-1)+2:getv(d,i) do setmat(A,j,j-1,1);
for j:=i:r do
<<
for k:=getv(d,i-1)+1:getv(d,i) do
<<
setmat(A,k,getv(d,j),{'minus,coeffn(nth(plist,j),x,k-1)});
>>;
>>;
>>;
return A;
end;
symbolic procedure size_of_matrix(A);
%
% Takes matrix and returns list {no. of rows, no. of columns).
%
begin
integer row_dim,col_dim;
matrix_input_test(A);
row_dim := -1 + length A;
col_dim := length cadr A;
return {row_dim,col_dim};
end;
symbolic procedure uppersmith(G,x);
%
% An upper triangular matrix G is simplified. Entry G(i,j) is reduced
% modulo gcd(G(i,i),G(j,j)). L and L^(-1) are also comnputed where
% L*G'*R=G, where G' is the reduced matrix.
%
begin
scalar US,L,Linv,g,s1,t1,q,r,tmp;
integer n;
n:= car size_of_matrix(G);
US:=copy_mat(G);
L := make_identity(n,n);
Linv := make_identity(n,n);
for j:=2:n do
<<
for i:=1:j-1 do
<<
tmp:=calc_exgcd(getmat(US,i,i),getmat(US,j,j),x);
g:= car tmp;
s1:= cadr tmp;
t1 := caddr tmp;
q := get_quo(getmat(US,i,j),g);
r := get_rem(getmat(US,i,j),g);
setmat(US,i,j,r);
for k:=1:i-1 do
<<
tmp := getmat(US,k,i);
setmat(US,k,j,{'plus,getmat(US,k,j),{'times,{'minus,q},s1,
getmat(US,k,i)}});
>>;
for k:=j+1:n do
<<
setmat(US,i,k,{'plus,getmat(US,i,k),{'times,{'minus,q},t1,
getmat(US,j,k)}});
>>;
for k:=1:i do
<<
setmat(L,k,j,{'plus,getmat(L,k,j),{'times,q,t1,
getmat(L,k,i)}});
>>;
setmat(Linv,i,j,{'times,{'minus,q},t1});
>>;
>>;
return {US,L,Linv};
end;
endmodule;
end;