Artifact c3e84e275c546f6aaa1ae5f54ff34de6daff9d30aafc5bc4e83e0987141a6733:
- Executable file
r37/packages/normform/smithex1.red
— part of check-in
[f2fda60abd]
at
2011-09-02 18:13:33
on branch master
— Some historical releases purely for archival purposes
git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/trunk/historical@1375 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 10285) [annotate] [blame] [check-ins using] [more...]
- Executable file
r38/packages/normform/smithex1.red
— part of check-in
[f2fda60abd]
at
2011-09-02 18:13:33
on branch master
— Some historical releases purely for archival purposes
git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/trunk/historical@1375 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 10285) [annotate] [blame] [check-ins using]
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;