restart;maple_mode(1);cas_setup(0,0,0,1,0,1e-10,10,[1,50,0,25],0,0,0);#radians,pas de cmplx, pas de  Sqrt
gcd(6,3);#Attention, a priori pour lui c'est des polynomes, on a dela chance il les normalise bien.
gcd(x*(x+3),2*x);
igcd(4,6,8);
iquo(13,6);igcd(6,0);  
 igcdex(4,15,'a','b');#pour un couple de bezout:
 4*a+15*b;
delrows(matrix([[1,1],[2,3]]),2..2); #pour une colonne: delcols
 A:=diag([3,6,18,36]);B:=matrix(4,4,(i,j)->rand(7)-3):;normal(det(A)*det(B)-det(ismith(A*B)[2]));
 // ismith(A*B);
[U0,D0,V0]:=ismith(A*B):;D0-U0*A*B*V0;
/* -----------------------------------------------------------------------------------------------------*/
 minval:=proc(A)
 local m,u,v,i,j;
 m:=0;u:=0;v:=0;
 # NB: On retourne 0,0 si A est nulle
 for i from 1 to dim(A)[1] do
 for j from 1 to dim(A)[2] do
 if ((abs(A[i,j])>0) and ((m=0) or (abs(A[i,j])<m))) then m:=abs(A[i,j]);u:=i;v:=j; end if;
 od;
 od;
 u,v;end proc;
 A:=matrix([[6, 0, 7, 9], [6, 8, 6, 9], [3, 2, 4, 6], [3, 2, 3, 3]]);
minval(A);
 (i,j):=minval(A);
 U:=identity(4);
/* normalement il faudrait faire bouger l dans {1..4} sauf j, c'est plus simple de£ne pas mettre de test, et de corriger U[j,j] ensuite*/
 
 for l from 1 to 4 do U[j,l]:=-iquo(A[i,l],A[i,j]) od:U[j,j]:=1:U;
A*U; 
 trans:=proc(A,i,j)
 local n,U,V,l;
 n:=dim(A)[1];
 U:=identity(n);
 V:=identity(n);
 for l from 1 to n do
 V[l,i]:=-iquo(A[l,j],A[i,j]);
 od;
 //on corrige
 V[i,i]:=1;
 for l from 1 to n do
 U[j,l]:=-iquo(A[i,l],A[i,j]);
 od;
 //on corrige
 U[j,j]:=1;
 V*A*U;
 end proc;
 A:=matrix([[6, 1, 7, 9], [6, 8, 6, 9], [3, 1, 4, 6], [3, 2, 1, 3]]);
 trans(A,3,1);
 Zequiv:=proc(A)
 k:=0;
 n:=dim(A)[1];
 l:=[seq(0,k=1..n)];
 B:=A;
 (i,j):=minval(B);
 //minval retourne les coordonnees du minimum non nul,
 //ou un couple impossible (Ex ici (0,0) en mode maple) si B est nulle.
 while([i,j]<>[0,0]) and size(B)>1 do
 // attention, ne pas oublier de declarer les variables locales
 // dans trans pour  eviter les melanges
 B:=trans(B,i,j);
 if [i,j]=[minval(B)] then
 // on n'a que des zeros sur la ligne i et la colonne j sauf
 // en (i,j) on sauve donc B[i,j] et on raye ligne et colonne.
 k:=k+1; l[k]:=B[i,j];
 B:=delcols(B,j..j);
 B:=delrows(B,i..i);
 //le cas B de taille 1 pose probleme car on ne peut pas rayer
 //une colonne puis une ligne
 fi;
 // ASTUCE: voici comment assigner 2 valeurs d'un coup.
 (i,j):=minval(B);
 end do;
 // On a eventuellement change le signe du determinant
 diag([seq(l[k],k=1..n-1),B[1,1]]);
 end proc;
 A:=matrix([[2,2,2],[6,12,6],[6,4,6]]); #Exemple:
 Zequiv(A);
 A[3,3]:=12:A;
 Zequiv(A);
/*  On v\'erifie que A et Zequiv(A) ont bien meme forme de smith:*/
 ismith(A)[2], ismith(Zequiv(A))[2];
/*  Attention, on n'obtient pas forcement les diviseurs elementaires, par exemple:*/
 A:=matrix([[4,0,0],[0,6,0],[0,0,8]]);
 ismith(A)[2];
 Zequiv(A);
/*  ca sera forc\'ement le pgcd(d_1,...,d_n) car  l'algorithme reste sur la premiere ligne.*/
 f:=(i,j)->if (i-j)*(i-1)=0 then 1 else 0 fi;
matrix(3,3,f);A:=matrix(3,3,f)*A;Zequiv(A);
/* On recopie trans et Zequiv pour qu'ils ne travaillent que sur les  colonnes*/
 transC:=proc(A,i,j)
 //attention, si on ne met pas un local n, il modifiera la valeur de n dans ZequivC
 local n,U,l;
 n:=dim(A)[1];
 U:=identity(n);
 for l from 1 to n do
 U[j,l]:=-iquo(A[i,l],A[i,j]);
 od;
 //pour ne pas mettre de test dans la boucle precedente, on n'ecarte pas l=j, mais on le on corrige ici:
 U[j,j]:=1;
 A*U;
 end proc;
 ZequivC:=proc(A)
 k:=0;
 n:=dim(A)[1];
 l:=[seq(0,k=1..n)];
 B:=A;
 (i,j):=minval(B);
 while([i,j]<>[0,0]) and size(B)>1 do
 // attention, ne pas oublier de declarer les variables locales dans trans pour  eviter les melanges
 B:=transC(B,i,j);
 if [i,j]=[minval(B)] then
 k:=k+1; l[k]:=B[i,j];
 B:=delcols(B,j..j);
 B:=delrows(B,i..i);
 fi;
 // ASTUCE: voici comment assigner 2 valeurs d'un coup.
 (i,j):=minval(B);
 end do;
 // On a eventuellement change le signe du determinant
 diag([seq(l[k],k=1..n-1),B[1,1]]);
 end proc;
 elem:=proc(A)
 n:=dim(A)[1];
 d:=Zequiv(A);
 L:=[];
 for i from 1 to n-1 do
 T:=matrix(n+1-i,n+1-i,f);
 d:=ZequivC(T*d);
 L:=[op(L),d[1,1]];
 d:=delrows(delcols(d,1..1),1..1);
 od;
 [op(L),d[1,1]];
 end proc;
 A:=matrix([[2,2,2],[6,12,6],[6,4,12]]);
 elem(A);ismith(A)[2];
 elem(diag(4,6,16));ismith(diag(4,6,16))[2];