module smithex1;
% %
%**********************************************************************%
% The function smithex_int computes the Smith normal form S of an n by
% m rectangular matrix of integers.
%
% Specifically:
%
% - smithex_int(A) will return {S,P,Pinv} where S, P, and Pinv
% are such that inverse(P)*A*P = S.
symbolic procedure smithex_int(B);
begin
scalar Left,Right,isclear,A;
integer n,m,i,j,k,l,tmp,g,ll,rr,int1,int2,quo1,quo2,r,sgn,rrquo,
q,input_mode;
matrix_input_test(B);
input_mode := get(dmode!*,'dname);
if input_mode = 'modular
then rederr "ERROR: smithex_int does not work with modular on.";
integer_entries_test(B);
A := copy_mat(B);
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 abs(getmat(A,l,k)) < abs(getmat(A,i,k)) then i := l;
>>;
for l:=j+1:m do
<<
if getmat(A,k,l) = 0 then l := l+1
else if abs(getmat(A,k,l)) < abs(getmat(A,k,j)) then j := l;
>>;
if i<=n and (j>m or abs(getmat(A,i,k))<abs(getmat(A,k,j))) 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
<<
int1 := getmat(A,k,k);
int2 := getmat(A,i,k);
tmp := (calc_exgcd_int(int1,int2));
g := car tmp;
ll := cadr tmp;
rr := caddr tmp;
quo1 := get_quo_int(getmat(A,k,k),g);
quo2 := get_quo_int(getmat(A,i,k),g);
%
% We have ll A(k,k)/g + rr A(i,k)/g = 1
%
% [ ll rr ] [ A[k,k] A[k,j] ] [ g ... ]
% [ ] [ ] = [ ]
% [ -quo2 quo1 ] [ A[i,k] A[i,j] ] [ 0 ... ]
%
% for j = k+1..m where note ll quo1 + rr quo2 = 1
%
for j:=k+1:m do
<<
tmp := ll*getmat(A,k,j)+rr*getmat(A,i,j);
setmat(A,i,j,quo1*getmat(A,i,j)-quo2*getmat(A,k,j));
setmat(A,k,j,tmp);
>>;
for j:=1:n do
<<
tmp := quo1*getmat(Left,j,k)+quo2*getmat(Left,j,i);
setmat(Left,j,i,-rr*getmat(Left,j,k)+ll*
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_int(getmat(A,k,i),getmat(A,k,k));
setmat(A,k,i,get_rem_int(getmat(A,k,i),getmat(A,k,k)));
for j:=1:m do
<<
setmat(Right,k,j,getmat(Right,k,j)+q*getmat(Right,i,j));
>>;
>>;
for i:=k+1:m do
<<
if getmat(A,k,i) = 0 then <<>>
else
<<
tmp := calc_exgcd_int( getmat(A,k,k),getmat(A,k,i));
g := car tmp;
ll := cadr tmp;
rr := caddr tmp;
quo1 := get_quo_int(getmat(A,k,k),g);
quo2 := get_quo_int(getmat(A,k,i),g);
for j:=k+1:n do
<<
tmp := ll*getmat(A,j,k) + rr*getmat(A,j,i);
setmat(A,j,i,quo1*getmat(A,j,i)-quo2*getmat(A,j,k));
setmat(A,j,k,tmp);
>>;
for j:=1:m do
<<
tmp := quo1*getmat(Right,k,j)+quo2*getmat(Right,i,j);
setmat(Right,i,j,-rr*getmat(Right,k,j)+ll*
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;
sgn := algebraic (sign(getmat(A,i,i)));
setmat(A,r,r,sgn*getmat(A,i,i));
if i = r then
<<
for j:=1:m do
<<
setmat(Right,i,j,getmat(Right,i,j)*sgn);
>>;
>>
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 := getmat(Right,i,j)*sgn;
setmat(Right,i,j,getmat(Right,r,j)*sgn);
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
<<
int1 := getmat(A,i,i);
int2 := getmat(A,j,j);
g := car (calc_exgcd_int(int1,int2));
ll := cadr (calc_exgcd_int(int1,int2));
rr := caddr (calc_exgcd_int(int1,int2));
quo1 := get_quo_int(getmat(A,i,i),g);
quo2 := get_quo_int(getmat(A,j,j),g);
setmat(A,i,i,g);
setmat(A,j,j,quo1*getmat(A,j,j));
for k:=1:n do
<<
tmp := quo1*getmat(Left,k,i)+quo2*getmat(Left,k,j);
setmat(Left,k,j,-rr*getmat(Left,k,i)+ll*
getmat(Left,k,j));
setmat(Left,k,i,tmp);
>>;
for k:=1:m do
<<
rrquo := rr*quo2;
tmp := (1-rrquo)*getmat(Right,i,k)+rrquo*
getmat(Right,j,k);
setmat(Right,j,k,-getmat(Right,i,k)+getmat(Right,j,k));
setmat(Right,i,k,tmp);
>>;
j := j+1;
>>;
>>;
>>;
return {'list,A,Left,Right};
end;
flag ('(smithex_int),'opfn); % So it can be used from algebraic mode.
symbolic procedure calc_exgcd_int(int1,int2);
begin
integer gcd,c,c1,c2,d,d1,d2,q,r,r1,r2,s1,t1;
if int1 = 0 and int2 = 0 then return {0,0,0}
else
<<
c := reval int1;
d := reval int2;
c1 := 1;
d1 := 0;
c2 := 0;
d2 := 1;
while d neq 0 do
<<
q := get_quo_int(c,d);
r := c - q*d;
r1 := c1 - q*d1;
r2 := c2 - q*d2;
c := d;
c1 := d1;
c2 := d2;
d := r;
d1 := r1;
d2 := r2;
>>;
gcd := abs(c);
s1 := c1;
t1 := c2;
if c < 0 then
<<
s1 := -s1;
t1 := -t1;
>>;
return {gcd,s1,t1};
>>;
end;
symbolic procedure get_quo_int(int1,int2);
begin
integer quo1,input1,input2;
input1 := reval int1;
input2 := reval int2;
if input1 = 0 and input2 = 0 then return 0
else
<<
if input1 < 0 and input2 < 0 then
<<
(input1 := abs(input1));
(input2 := abs(input2));
>>;
if (input1/input2) < 0 then
<<
quo1 :=ceiling(input1/input2);
>>
else
<<
quo1 :=floor(input1/input2);
>>;
return quo1;
>>;
end;
symbolic procedure get_rem_int(int1,int2);
begin
integer rem1,input1,input2;
input1 := reval int1;
input2 := reval int2;
rem1 := input1 - get_quo_int(input1,input2)*input2;
return rem1;
end;
symbolic procedure integer_entries_test(B);
begin
for i:=1:car size_of_matrix(B) do
<<
for j:=1:cadr size_of_matrix(B) do
<<
if not numberp getmat(B,i,j) then rederr
"ERROR: matrix contains non_integer entries. Try smithex. "
>>;
>>;
end;
endmodule;
end;