module smithex; % Computation of the Smithex normal form of a matrix. %
% The function smithex computes the Smith normal form S of an n by m
% rectangular matrix of univariate polynomials in x.
%
% Specifically:
%
% - smithex(A,x) will return {S,P,Pinv} where S, P, and Pinv
% are such that P*S*Pinv = A.
symbolic procedure smithex(mat1,x);
begin
scalar A,Left,Right,tmp,isclear,g,L,R1,poly1,poly2,quo1,quo2,r,
lc,tquo,q,full_coeff_list,rule_list,input_mode;
integer i,j,n,m;
matrix_input_test(mat1);
input_mode := get(dmode!*,'dname);
if input_mode = 'modular
then rederr "ERROR: smithex does not work with modular on.";
all_integer_entries_test(mat1);
if input_mode neq 'arnum and input_mode neq 'rational
then on rational;
on combineexpt;
tmp := nest_input_smith(mat1,x);
A := car tmp;
full_coeff_list := cadr tmp;
n := car size_of_matrix(A); % No. of rows.
m := cadr size_of_matrix(A); % No. of columns.
Left := make_identity(n,n) ;
Right := make_identity(m,m);
for k:=1:min(n,m) do
<<
%
% Pivot selection from row k to column k.
%
i := k; while i <= n and getmat(A,i,k) = 0 do i := i+1;
j := k; while j <= m and getmat(A,k,j) = 0 do j := j+1;
if i > n and j > m then <<>>
else
<<
%
% Select smallest non-zero entry as pivot.
%
for l:=i+1:n do
<<
if getmat(A,l,k) = 0 then l := l+1
else if deg(getmat(A,l,k),x) < deg(getmat(A,i,k),x) then
i := l;
>>;
for l:=j+1:m do
<<
if getmat(A,k,l) = 0 then l := l+1
else if deg(getmat(A,k,l),x) < deg(getmat(A,k,j),x) then
j := l;
>>;
if i <= n and (j > m or deg(getmat(A,i,k),x) <
deg(getmat(A,k,j),x))
then
%
% Pivot is A(i,k), interchange row k,i if necessary.
%
<<
if i neq k then
<<
for l:=k:m do
<<
tmp := getmat(A,i,l);
setmat(A,i,l,getmat(A,k,l));
setmat(A,k,l,tmp);
>>;
for l:=1:n do
<<
tmp := getmat(Left,l,i);
setmat(Left,l,i,getmat(Left,l,k));
setmat(Left,l,k,tmp);
>>;
>>
>>
else
%
% Pivot is A(k,j), interchange column k,j if necessary.
%
<<
if j neq k then
<<
for l:=k:n do
<<
tmp := getmat(A,l,j);
setmat(A,l,j,getmat(A,l,k));
setmat(A,l,k,tmp);
>>;
for l:=1:m do
<<
tmp := getmat(Right,j,l);
setmat(Right,j,l,getmat(Right,k,l));
setmat(Right,k,l,tmp);
>>;
>>;
>>;
isclear := nil;
while not isclear do
%
% 0 out column k from k+1 to n.
%
<<
for i:=k+1:n do
<<
if getmat(A,i,k) = 0 then <<>>
else
<<
poly1 := getmat(A,k,k);
poly2 := getmat(A,i,k);
tmp := calc_exgcd(poly1,poly2,x);
g := car tmp;
L := cadr tmp;
R1 := caddr tmp;
quo1 := get_quo(poly1,g);
quo2 := get_quo(poly2,g);
for j:=k+1:m do
<<
tmp := {'plus,{'times,L,getmat(A,k,j)},{'times,R1,
getmat(A,i,j)}};
setmat(A,i,j,{'plus,{'times,quo1,getmat(A,i,j)},{'times,
{'minus,quo2},getmat(A,k,j)}});
setmat(A,k,j,tmp);
>>;
for j:=1:n do
<<
tmp := {'plus,{'times,quo1,getmat(Left,j,k)},
{'times,quo2,getmat(Left,j,i)}};
setmat(Left,j,i,{'plus,{'times,{'minus,R1},
getmat(Left,j,k)},{'times,L,getmat(Left,j,i)}});
setmat(Left,j,k,tmp);
>>;
setmat(A,k,k,g);
setmat(A,i,k,0);
>>;
>>;
isclear := t;
%
% 0 out row k from k+1 to m.
%
for i:=k+1:m do
<<
q := get_quo(getmat(A,k,i),getmat(A,k,k));
setmat(A,k,i,get_rem(getmat(A,k,i),getmat(A,k,k)));
for j:=1:m do
<<
setmat(Right,k,j,{'plus,getmat(Right,k,j),{'times,q,
getmat(Right,i,j)}});
>>;
>>;
for i:=k+1:m do
<<
if getmat(A,k,i) = 0 then <<>>
else
<<
poly1 := getmat(A,k,k);
poly2 := getmat(A,k,i);
tmp := calc_exgcd(poly1,poly2,x);
g := car tmp;
L := cadr tmp;
R1 := caddr tmp;
quo1 := get_quo(poly1,g);
quo2 := get_quo(poly2,g);
for j:=k+1:n do
<<
tmp := {'plus,{'times,L,getmat(A,j,k)},{'times,R1,
getmat(A,j,i)}};
setmat(A,j,i,{'plus,{'times,quo1,getmat(A,j,i)},{'times,
{'minus,quo2},getmat(A,j,k)}});
setmat(A,j,k,tmp);
>>;
for j:=1:m do
<<
tmp := {'plus,{'times,quo1,getmat(Right,k,j)},
{'times,quo2,getmat(Right,i,j)}};
setmat(Right,i,j,{'plus,{'times,{'minus,R1},
getmat(Right,k,j)},
{'times,L,getmat(Right,i,j)}});
setmat(Right,k,j,tmp);
>>;
setmat(A,k,k,g);
setmat(A,k,i,0);
isclear := nil;
>>;
>>;
>>;
>>;
>>;
r := 0;
%
% At this point, A is diagonal: some A(i,i) may be zero.
% Move non-zero's up also making all entries unit normal.
%
for i:=1:min(n,m) do
<<
if getmat(A,i,i) neq 0 then
<<
r := r+1;
% Watch out for integers giving lc = 0.
if lcof(getmat(A,i,i),x) = 0 then lc := getmat(A,i,i)
else lc := lcof(getmat(A,i,i),x);
setmat(A,r,r,{'quotient,getmat(A,i,i),lc});
if i = r then
<<
for j:=1:m do
<<
setmat(Right,i,j,{'times,getmat(Right,i,j),lc});
>>;
>>
else
<<
setmat(A,i,i,0);
for j:=1:n do
<<
tmp := getmat(Left,j,r);
setmat(Left,j,r,getmat(Left,j,i));
setmat(Left,j,i,tmp);
>>;
for j:=1:m do
<<
tmp := {'times,getmat(Right,i,j),lc};
setmat(Right,i,j,{'quotient,getmat(Right,r,j),lc});
setmat(Right,r,j,tmp);
>>;
>>;
>>;
>>;
%
% Now make A(i,i) | A(i+1,i+1) for 1 <= i < r.
%
for i:=1:r-1 do
<<
j:=i+1;
<<
while getmat(A,i,i) neq 1 and j <= r do
<<
poly1 := getmat(A,i,i);
poly2 := getmat(A,j,j);
tmp := calc_exgcd(poly1,poly2,x);
g := car tmp;
L := cadr tmp;
R1 := caddr tmp;
quo1 := get_quo(poly1,g);
quo2 := get_quo(poly2,g);
setmat(A,i,i,g);
setmat(A,j,j,{'times,quo1,getmat(A,j,j)});
for k:=1:n do
<<
tmp := {'plus,{'times,quo1,getmat(Left,k,i)},{'times,quo2,
getmat(Left,k,j)}};
setmat(Left,k,j,{'plus,{'times,{'minus,R1},
getmat(Left,k,i)},{'times,L,getmat(Left,k,j)}});
setmat(Left,k,i,tmp);
>>;
for k:=1:m do
<<
tquo := {'times,R1,quo2};
tmp := {'plus,{'times,{'plus,1,{'minus,tquo}},
getmat(Right,i,k)},{'times,tquo,getmat(Right,j,k)}};
setmat(Right,j,k,{'plus,{'minus,getmat(Right,i,k)},
getmat(Right,j,k)});
setmat(Right,i,k,tmp);
>>;
j := j+1;
>>;
>>;
>>;
%
% 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;
A := de_nest_mat(A);
Left := de_nest_mat(Left);
Right := de_nest_mat(Right);
clearrules rule_list;
%
% Return to original mode.
%
if input_mode neq 'rational and input_mode neq 'arnum then
<<
% onoff('nil,t) doesn't work so...
if input_mode = 'nil then off rational
else onoff(input_mode,t);
>>;
off combineexpt;
return {'list,A,Left,Right};
end;
flag ('(smithex),'opfn); % So it can be used from algebraic mode.
symbolic procedure get_coeffs_smith(poly,x);
%
% Gets all kernels in a poly.
%
% This is the version form smithex. The polynomials are
% known to be in x (smithex(matrix,x)) so this is removed
% from the output so as to not be nested in
% nest_input_smith.
%
begin
scalar ker_list_num,ker_list_den,new_list;
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);
if ker_list_num = nil then new_list := nil
else
<<
% Remove x.
for i:=1:length ker_list_num do
<<
if car ker_list_num = x then new_list := new_list
else new_list := car ker_list_num.new_list;
ker_list_num := cdr ker_list_num;
>>;
>>;
return new_list;
end;
symbolic procedure nest_input_smith(A,x);
%
% Takes a matrix and converts all elements into nested form.
% Also finds all coefficients and returns them in a list.
%
% With Smithex all polynomials are in x (input by user) so
% this is removed from the full_coeff_list (see get_coeffs),
% and is thus not nested.
%
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_smith(getmat(AA,i,j),x);
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;
switch int_test;
symbolic procedure all_integer_entries_test(mat1);
begin
on int_test;
for i:=1:car size_of_matrix(mat1) do
<<
for j:=1:cadr size_of_matrix(mat1) do
<<
if not numberp getmat(mat1,i,j) and !*int_test
then off int_test;
>>;
>>;
if !*int_test then prin2t
"*** WARNING: all matrix entries are integers.
If calculations in Z(the integers) are required, use smithex_int.";
% system "sleep 3";
end;
endmodule;
end;