////////////////////////////////////////////////////////////////////////////// version="$Id$"; category="System and Control Theory"; info=" LIBRARY: jacobson.lib Algorithms for Smith and Jacobson Normal Form AUTHOR: Kristina Schindelar, Kristina.Schindelar@math.rwth-aachen.de, @* Viktor Levandovskyy, levandov@math.rwth-aachen.de OVERVIEW: We work over a ring R, that is an Euclidean principal ideal domain. If R is commutative, we suppose R to be a polynomial ring in one variable. If R is non-commutative, we suppose R to have two variables, say x and d. We treat then the basering as the Ore localization of R with respect to the mult. closed set S = K[x] without 0. Thus, we treat basering as principal ideal ring with d a polynomial variable and x an invertible one. Note, that in computations no division by x will actually happen. Given a rectangular matrix M over R, one can compute unimodular (that is invertible) square matrices U and V, such that U*M*V=D is a diagonal matrix. Depending on the ring, the diagonal entries of D have certain properties. We call a square matrix D as before 'a weak Jacobson normal form of M'. It is known, that over the first rational Weyl algebra K(x), D can be further transformed into a diagonal matrix (1,1,...,1,f,0,..,0), where f is in K(x). We call such a form of D the strong Jacobson normal form. The existence of strong form in not guaranteed if one works with algebra, which is not rational Weyl algebra. REFERENCES: @* [1] N. Jacobson, 'The theory of rings', AMS, 1943. @* [2] Manuel Avelino Insua Hermo, 'Varias perspectives sobre las bases de Groebner : @* Forma normal de Smith, Algorithme de Berlekamp y algebras de Leibniz'. @* PhD thesis, Universidad de Santiago de Compostela, 2005. @* [3] V. Levandovskyy, K. Schindelar 'Computing Jacobson normal form using Groebner bases', @* to appear in Journal of Symbolic Computation, 2010. PROCEDURES: smith(M[,eng1,eng2]); compute the Smith Normal Form of M over commutative ring jacobson(M[,eng]); compute a weak Jacobson Normal Form of M over non-commutative ring divideUnits(L); create ones out of units in the output of smith or jacobson SEE ALSO: control_lib KEYWORDS: Jacobson form; Jacobson normal form; Smith form; Smith normal form; matrix diagonalization "; LIB "poly.lib"; LIB "involut.lib"; // involution LIB "dmodapp.lib"; // engine LIB "qhmoduli.lib"; // Min LIB "latex.lib"; // tex output for printlevel=4 proc tstjacobson() { /* tests all procs for consistency */ example divideUnits; example smith; example jacobson; } proc divideUnits(list L) "USAGE: divideUnits(L); list L RETURN: matrix or list of matrices ASSUME: L is an output of @code{smith} or a @code{jacobson} procedures, that is @* either L contains one rectangular matrix with elements only on the main diagonal @* or L consists of three matrices, where L[1] and L[3] are square invertible matrices @* while L[2] is a rectangular matrix with elements only on the main diagonal PURPOSE: divide out units from the diagonal and reflect this in transformation matrices EXAMPLE: example divideUnits; shows examples " { // check assumptions int s = size(L); if ( (s!=1) && (s!=3) ) { ERROR("The list has neither 1 nor 3 elements"); } // check sizes of matrices if (s==1) { list LL; LL[1] = L[1]; LL[2] = L[1]; L = LL; } // divide units out matrix M = L[2]; int ncM = ncols(M); int nrM = nrows(M); // matrix TM[nrM][nrM]; // m times m matrix matrix TM = matrix(freemodule(nrM)); int i; int nsize = Min(intvec(nrM,ncM)); poly p; number n; intvec lexp; for(i=1; i<=nsize; i++) { p = M[i,i]; lexp = leadexp(p); // TM[i,i] = 1; // check: is p a unit? if (p!=0) { if ( lexp == 0) { // hence p = n*1 n = leadcoef(p); TM[i,i] = 1/n; } } } int ppl = printlevel-voice+2; dbprint(ppl,"divideUnits: extra transformation matrix is: "); dbprint(ppl, TM); L[2] = TM*L[2]; if (s==3) { L[1] = TM*L[1]; return(L); } else { return(L[2]); } } example { "EXAMPLE:"; echo = 2; ring R=(0,m,M,L1,L2,m1,m2,g), D, lp; // two pendula example matrix P[3][4]=m1*L1*D^2,m2*L2*D^2,(M+m1+m2)*D^2,-1, m1*L1^2*D^2-m1*L1*g,0,m1*L1*D^2,0,0, m2*L2^2*D^2-m2*L2*g,m2*L2*D^2,0; list s=smith(P,1); // returns a list with 3 entries print(s[2]); // a diagonal form, close to the Smith form print(s[1]); // U, left transformation matrix list t = divideUnits(s); print(t[2]); // the Smith form of the matrix P print(t[1]); // U', modified left transformation matrix } proc smith(matrix MA, list #) "USAGE: smith(M[, eng1, eng2]); M matrix, eng1 and eng2 are optional integers RETURN: matrix or list of matrices, depending on arguments ASSUME: Basering is a commutative polynomial ring in one variable PURPOSE: compute the Smith Normal Form of M with (optionally) transformation matrices THEORY: Groebner bases are used for the Smith form like in [2] and [3]. NOTE: By default, just the Smith normal form of M is returned. @* If the optional integer @code{eng1} is non-zero, the list {U,D,V} is returned @* where U*M*V = D and the diagonal field entries of D are not normalized. @* The normalization of the latter can be done with the 'divideUnits' procedure. @* U and V above are square unimodular (invertible) matrices. @* Note, that the procedure works for a rectangular matrix M. @* @* The optional integer @code{eng2} determines the Groebner basis engine: @* 0 (default) ensures the use of 'slimgb' , otherwise 'std' is used. DISPLAY: If @code{printlevel}=1, progress debug messages will be printed, @* if @code{printlevel}>=2, all the debug messages will be printed. EXAMPLE: example smith; shows examples SEE ALSO: divideUnits, jacobson " { def R = basering; // check assume: R must be a polynomial ring in 1 variable setring R; if (nvars(R) >1 ) { ERROR("Basering must have exactly one variable"); } int eng = 0; int BASIS; if ( size(#)>0 ) { eng=1; if (typeof(#[1])=="int") { eng=int(#[1]); // zero can also happen } } if (size(#)==2) { BASIS=#[2]; } else{BASIS=0;} int ROW=ncols(MA); int COL=nrows(MA); //generate a module consisting of the columns of MA module m=MA[1]; int i; for(i=2;i<=ROW;i++) { m=m,MA[i]; } //if MA eqauls the zero matrix give back MA if ( (size(module(m))==0) and (size(transpose(module(m)))==0) ) { module L=freemodule(COL); matrix LM=L; L=freemodule(ROW); matrix RM=L; list RUECK=RM; RUECK=insert(RUECK,MA); RUECK=insert(RUECK,LM); return(RUECK); } if(eng==1) { list rueckLI=diagonal_with_trafo(R,MA,BASIS); list rueckLII=divisibility(rueckLI[2]); matrix CON=divideByContent(rueckLII[2]); list rueckL=CON*rueckLII[1]*rueckLI[1], CON*rueckLII[2], rueckLI[3]*rueckLII[3]; return(rueckL); } else { matrix rueckm=diagonal_without_trafo(R,MA,BASIS); list rueckL=divisibility(rueckm); matrix CON=divideByContent(rueckm); rueckm=CON*rueckL[2]; return(rueckm); } } example { "EXAMPLE:"; echo = 2; ring r = 0,x,Dp; matrix m[3][2]=x, x^4+x^2+21, x^4+x^2+x, x^3+x, 4*x^2+x, x; list s=smith(m,1); print(s[2]); // non-normalized Smith form of m print(s[1]*m*s[3] - s[2]); // check U*M*V = D list t = divideUnits(s); print(t[2]); // the Smith form of m } static proc diagonal_with_trafo( R, matrix MA, int B) { int ppl = printlevel-voice+2; int BASIS=B; int ROW=ncols(MA); int COL=nrows(MA); module m=MA[1]; int i,j,k; for(i=2;i<=ROW;i++) { m=m,MA[i]; } //add zero rows or columns //add zero rows or columns int adrow=0; for(i=1;i<=COL;i++) { k=0; for(j=1;j<=ROW;j++) { if(MA[i,j]!=0){k=1;} } if(k==0){adrow++;} } m=transpose(m); for(i=1;i<=adrow;i++){m=m,0;} m=transpose(m); list RINGLIST=ringlist(R); list o="C",0; list oo="lp",1; list ORD=o,oo; RINGLIST[3]=ORD; def r=ring(RINGLIST); setring r; //fix the required ordering map MAP=R,var(1); module m=MAP(m); int flag=1; ///////////////counts if the underlying ring is r (flag mod 2 == 1) or ro (flag mod 2 == 0) module TrafoL=freemodule(COL); module TrafoR=freemodule(ROW); module EXL=freemodule(COL); //because we start with transpose(m) module EXR=freemodule(ROW); option(redSB); option(redTail); module STD_EX; module Trafo; int s,st,p,ff; module LT,TSTD; module STD=transpose(m); int finish=0; int fehlpos; list pos; matrix END; matrix ENDSTD; matrix STDFIN; STDFIN=STD; list COMPARE=STDFIN; while(finish==0) { dbprint(ppl,"Going into the while cycle"); if(flag mod 2==1) { STD_EX=EXL,transpose(STD); } else { STD_EX=EXR,transpose(STD); } dbprint(ppl,"Computing Groebner basis: start"); dbprint(ppl-1,STD); STD=engine(STD,BASIS); dbprint(ppl,"Computing Groebner basis: finished"); dbprint(ppl-1,STD); if(flag mod 2==1){s=ROW; st=COL;}else{s=COL; st=ROW;} dbprint(ppl,"Computing Groebner basis for transformation matrix: start"); dbprint(ppl-1,STD_EX); STD_EX=transpose(STD_EX); STD_EX=engine(STD_EX,BASIS); dbprint(ppl,"Computing Groebner basis for transformation matrix: finished"); dbprint(ppl-1,STD_EX); //////// split STD_EX in STD and the transformation matrix STD_EX=transpose(STD_EX); STD=STD_EX[st+1]; LT=STD_EX[1]; ENDSTD=STD_EX; for(i=2; i<=ncols(ENDSTD); i++) { if (i<=st) { LT=LT,STD_EX[i]; } if (i>st+1) { STD=STD,STD_EX[i]; } } STD=transpose(STD); LT=transpose(LT); ////////////////////// compute the transformation matrices if (flag mod 2 ==1) { TrafoL=transpose(LT)*TrafoL; } else { TrafoR=TrafoR*LT; } STDFIN=STD; /////////////////////////////////// check if the D-th row is finished /////////////////////////////////// COMPARE=insert(COMPARE,STDFIN); if(size(COMPARE)>=3) { if(STD==COMPARE[3]){finish=1;} } ////////////////////////////////// change to the opposite module TSTD=transpose(STD); STD=TSTD; flag++; dbprint(ppl,"Finished one while cycle"); } if (flag mod 2!=0) { STD=transpose(STD); } //zero colums to the end matrix STDMM=STD; pos=list(); fehlpos=0; while( size(STD)+fehlpos-ncols(STDMM) < 0) { for(i=1; i<=ncols(STDMM); i++) { ff=0; for(j=1; j<=nrows(STDMM); j++) { if (STD[j,i]==0) { ff++; } } if(ff==nrows(STDMM)) { pos=insert(pos,i); fehlpos++; } } } int fehlposc=fehlpos; module SORT; for(i=1; i<=fehlpos; i++) { SORT=gen(2); for (j=3;j<=ROW;j++) { SORT=SORT,gen(j); } SORT=SORT,gen(1); STD=STD*SORT; TrafoR=TrafoR*SORT; } //zero rows to the end STDMM=transpose(STD); pos=list(); fehlpos=0; while( size(transpose(STD))+fehlpos-ncols(STDMM) < 0) { for(i=1; i<=ncols(STDMM); i++) { ff=0; for(j=1; j<=nrows(STDMM); j++) { if(transpose(STD)[j,i]==0){ ff++;}} if(ff==nrows(STDMM)) { pos=insert(pos,i); fehlpos++; } } } int fehlposr=fehlpos; for(i=1; i<=fehlpos; i++) { SORT=gen(2); for(j=3;j<=COL;j++){SORT=SORT,gen(j);} SORT=SORT,gen(1); SORT=transpose(SORT); STD=SORT*STD; TrafoL=SORT*TrafoL; } setring R; map MAPinv=r,var(1); module STD=MAPinv(STD); module TrafoL=MAPinv(TrafoL); matrix TrafoLM=TrafoL; module TrafoR=MAPinv(TrafoR); matrix TrafoRM=TrafoR; matrix STDM=STD; //Test if(TrafoLM*m*TrafoRM!=STDM){ return(0); } list RUECK=TrafoRM; RUECK=insert(RUECK,STDM); RUECK=insert(RUECK,TrafoLM); return(RUECK); } static proc divisibility(matrix M) { matrix STDM=M; int i,j; int ROW=nrows(M); int COL=ncols(M); module TrafoR=freemodule(COL); module TrafoL=freemodule(ROW); module SORT; matrix TrafoRM=TrafoR; matrix TrafoLM=TrafoL; list posdeg0; int posdeg=0; int act; int Sort=ROW; if(size(std(STDM))!=0) { while( size(transpose(STDM)[Sort])==0 ){Sort--;}} for(i=1;i<=Sort ;i++) { if(leadexp(STDM[i,i])==0){posdeg0=insert(posdeg0,i);posdeg++;} } //entries of degree 0 at the beginning for(i=1; i<=posdeg; i++) { act=posdeg0[i]; SORT=gen(act); for(j=1; j<=COL; j++){if(j!=act){SORT=SORT,gen(j);}} STDM=STDM*SORT; TrafoRM=TrafoRM*SORT; SORT=gen(act); for(j=1; j<=ROW; j++){if(j!=act){SORT=SORT,gen(j);}} STDM=transpose(SORT)*STDM; TrafoLM=transpose(SORT)*TrafoLM; } poly G; module UNITL=freemodule(ROW); matrix GCDL=UNITL; module UNITR=freemodule(COL); matrix GCDR=UNITR; for(i=posdeg+1; i<=Sort; i++) { for(j=i+1; j<=Sort; j++) { GCDL=UNITL; GCDR=UNITR; G=gcd(STDM[i,i],STDM[j,j]); ideal Z=STDM[i,i],STDM[j,j]; matrix T=lift(Z,G); GCDL[i,i]=T[1,1]; GCDL[i,j]=T[2,1]; GCDL[j,i]=-STDM[j,j]/G; GCDL[j,j]=STDM[i,i]/G; GCDR[i,j]=T[2,1]*STDM[j,j]/G; GCDR[j,j]=T[2,1]*STDM[j,j]/G-1; GCDR[j,i]=1; STDM=GCDL*STDM*GCDR; TrafoLM=GCDL*TrafoLM; TrafoRM=TrafoRM*GCDR; } } list RUECK=TrafoRM; RUECK=insert(RUECK,STDM); RUECK=insert(RUECK,TrafoLM); return(RUECK); } static proc diagonal_without_trafo( R, matrix MA, int B) { int ppl = printlevel-voice+2; int BASIS=B; int ROW=ncols(MA); int COL=nrows(MA); module m=MA[1]; int i; for(i=2;i<=ROW;i++) {m=m,MA[i];} list RINGLIST=ringlist(R); list o="C",0; list oo="lp",1; list ORD=o,oo; RINGLIST[3]=ORD; def r=ring(RINGLIST); setring r; //RICHTIGE ORDNUNG MACHEN map MAP=R,var(1); module m=MAP(m); int flag=1; ///////////////counts if the underlying ring is r (flag mod 2 == 1) or ro (flag mod 2 == 0) int act, j, ff; option(redSB); option(redTail); module STD=transpose(m); module TSTD; int finish=0; matrix STDFIN; STDFIN=STD; list COMPARE=STDFIN; while(finish==0) { dbprint(ppl,"Going into the while cycle"); dbprint(ppl,"Computing Groebner basis: start"); dbprint(ppl-1,STD); STD=engine(STD,BASIS); dbprint(ppl,"Computing Groebner basis: finished"); dbprint(ppl-1,STD); STDFIN=STD; /////////////////////////////////// check if the D-th row is finished /////////////////////////////////// COMPARE=insert(COMPARE,STDFIN); if(size(COMPARE)>=3) { if(STD==COMPARE[3]){finish=1;} } ////////////////////////////////// change to the opposite module TSTD=transpose(STD); STD=TSTD; flag++; dbprint(ppl,"Finished one while cycle"); } matrix STDMM=STD; list pos=list(); int fehlpos=0; while( size(STD)+fehlpos-ncols(STDMM) < 0) { for(i=1; i<=ncols(STDMM); i++) { ff=0; for(j=1; j<=nrows(STDMM); j++) { if (STD[j,i]==0) { ff++; } } if(ff==nrows(STDMM)) { pos=insert(pos,i); fehlpos++; } } } int fehlposc=fehlpos; module SORT; for(i=1; i<=fehlpos; i++) { SORT=gen(2); for (j=3;j<=ROW;j++) { SORT=SORT,gen(j); } SORT=SORT,gen(1); STD=STD*SORT; } //zero rows to the end STDMM=transpose(STD); pos=list(); fehlpos=0; while( size(transpose(STD))+fehlpos-ncols(STDMM) < 0) { for(i=1; i<=ncols(STDMM); i++) { ff=0; for(j=1; j<=nrows(STDMM); j++) { if(transpose(STD)[j,i]==0){ ff++;}} if(ff==nrows(STDMM)) { pos=insert(pos,i); fehlpos++; } } } int fehlposr=fehlpos; for(i=1; i<=fehlpos; i++) { SORT=gen(2); for(j=3;j<=COL;j++){SORT=SORT,gen(j);} SORT=SORT,gen(1); SORT=transpose(SORT); STD=SORT*STD; } //add zero rows or columns int adrow=COL-size(transpose(STD)); int adcol=ROW-size(STD); for(i=1;i<=adcol;i++){STD=STD,0;} STD=transpose(STD); for(i=1;i<=adrow;i++){STD=STD,0;} STD=transpose(STD); setring R; map MAPinv=r,var(1); module STD=MAPinv(STD); matrix STDM=STD; return(STDM); } // VL : engine replaced by the one from dmodapp.lib // cases are changed // static proc engine(module I, int i) // { // module J; // if (i==0) // { // J = std(I); // } // if (i==1) // { // J = groebner(I); // } // if (i==2) // { // J = slimgb(I); // } // return(J); // } proc jacobson(matrix MA, list #) "USAGE: jacobson(M, eng); M matrix, eng an optional int RETURN: list ASSUME: Basering is a (non-commutative) ring in two variables. PURPOSE: compute a weak Jacobson normal form of M over the basering THEORY: Groebner bases and involutions are used, following [3] NOTE: A list L of matrices {U,D,V} is returned. That is L[1]*M*L[3]=L[2], @* where L[2] is a diagonal matrix and @* L[1], L[3] are square invertible polynomial (unimodular) matrices. @* Note, that M can be rectangular. @* The optional integer @code{eng2} determines the Groebner basis engine: @* 0 (default) ensures the use of 'slimgb' , otherwise 'std' is used. DISPLAY: If @code{printlevel}=1, progress debug messages will be printed, @* if @code{printlevel}>=2, all the debug messages will be printed. EXAMPLE: example jacobson; shows examples SEE ALSO: divideUnits, smith " { def R = basering; // check assume: R must be a polynomial ring in 2 variables setring R; if ( (nvars(R) !=2 ) ) { ERROR("Basering must have exactly two variables"); } // check if MA is zero matrix and return corr. U,V if ( (size(module(MA))==0) and (size(transpose(module(MA)))==0) ) { int nr = nrows(MA); int nc = ncols(MA); matrix U = matrix(freemodule(nr)); matrix V = matrix(freemodule(nc)); list L; L[1]=U; L[2] = MA; L[3] = V; return(L); } int B=0; if ( size(#)>0 ) { B=1; if (typeof(#[1])=="int") { B=int(#[1]); // zero can also happen } } //change ring list RINGLIST=ringlist(R); list o="C",0; intvec v=0,1; list oo="a",v; v=1,1; list ooo="lp",v; list ORD=o,oo,ooo; RINGLIST[3]=ORD; def r=ring(RINGLIST); setring r; //fix the required ordering map MAP=R,var(1),var(2); matrix M=MAP(MA); module TrafoL, TrafoR; list TRIANGLE; TRIANGLE=triangle(M,B); TrafoL=TRIANGLE[1]; TrafoR=TRIANGLE[3]; module m=TRIANGLE[2]; //back to the ring setring R; map MAPR=r,var(1),var(2); module ma=MAPR(m); matrix MAA=ma; module TL=MAPR(TrafoL); module TR=MAPR(TrafoR); matrix TRR=TR; matrix CON=divideByContent(MAA); list RUECK=CON*TL, CON*MAA, TRR; return(RUECK); } example { "EXAMPLE:"; echo = 2; ring r = 0,(x,d),Dp; def R = nc_algebra(1,1); setring R; // the 1st Weyl algebra matrix m[2][2] = d,x,0,d; print(m); list J = jacobson(m); // returns a list with 3 entries print(J[2]); // a Jacobson Form D for m print(J[1]*m*J[3] - J[2]); // check that U*M*V = D /* now, let us do the same for the shift algebra */ ring r2 = 0,(x,s),Dp; def R2 = nc_algebra(1,s); setring R2; // the 1st shift algebra matrix m[2][2] = s,x,0,s; print(m); // matrix of the same for as above list J = jacobson(m); print(J[2]); // a Jacobson Form D, quite different from above print(J[1]*m*J[3] - J[2]); // check that U*M*V = D } //static proc triangle(matrix MA, int B) { int ppl = printlevel-voice+2; map inv=ncdetection(); int ROW=ncols(MA); int COL=nrows(MA); //generate a module consisting of the columns of MA module m=MA[1]; int i,j,s,st,p,k,ff,ex, nz, ii,nextp; for(i=2;i<=ROW;i++) { m=m,MA[i]; } int BASIS=B; //add zero rows or columns int adrow=0; for(i=1;i<=COL;i++) { k=0; for(j=1;j<=ROW;j++) { if(MA[i,j]!=0){k=1;} } if(k==0){adrow++;} } m=transpose(m); for(i=1;i<=adrow;i++){m=m,0;} m=transpose(m); int flag=1; ///////////////counts if the underlying ring is r (flag mod 2 == 1) or ro (flag mod 2 == 0) module TrafoL=freemodule(COL); module TrafoR=freemodule(ROW); module EXL=freemodule(COL); //because we start with transpose(m) module EXR=freemodule(ROW); option(redSB); option(redTail); module STD_EX,LT,TSTD, L, Trafo; module STD=transpose(m); int finish=0; list pos, COM, COM_EX; matrix END, ENDSTD, STDFIN; STDFIN=STD; list COMPARE=STDFIN; string @s; while(finish==0) { dbprint(ppl,"Going into the while cycle"); if(flag mod 2==1){STD_EX=EXL,transpose(STD); ex=2*COL;} else {STD_EX=EXR,transpose(STD); ex=2*ROW;} dbprint(ppl,"Computing Groebner basis: start"); dbprint(ppl-1,STD); STD=engine(STD,BASIS); dbprint(ppl,"Computing Groebner basis: finished"); dbprint(ppl-1,STD); if(flag mod 2==1){s=ROW; st=COL;}else{s=COL; st=ROW;} STD_EX=transpose(STD_EX); dbprint(ppl,"Computing Groebner basis for transformation matrix: start"); dbprint(ppl-1,STD_EX); STD_EX=engine(STD_EX,BASIS); dbprint(ppl,"Computing Groebner basis for transformation matrix: finished"); dbprint(ppl-1,STD_EX); COM=1; COM_EX=1; for(i=2; i<=size(STD); i++) { COM=COM[1..size(COM)],i; COM_EX=COM_EX[1..size(COM_EX)],i; } nz=size(STD_EX)-size(STD); //zero rows in the begining if(size(STD)!=size(STD_EX) ) { for(i=1; i<=size(STD_EX)-size(STD); i++) { COM_EX=0,COM_EX[1..size(COM_EX)]; } } for(i=nz+1; i<=size(STD_EX); i++ ) {if( leadcoef(STD[i-nz])!=leadcoef(STD_EX[i]) ) {STD[i-nz]=leadcoef(STD_EX[i])*STD[i-nz];} } //assign the zero rows in STD_EX for (j=2; j<=nz; j++) { p=0; i=1; while(STD_EX[j-1][i]==0){i++;} p=i-1; nextp=1; k=0; i=1; while(STD_EX[j][i]==0 and i<=p) { k++; i++;} if (k==p){ COM_EX[j]=-1; } } //assign the zero rows in STD for (j=2; j<=size(STD); j++) { i=size(transpose(STD)); while(STD[j-1][i]==0){i--;} p=i; i=size(transpose(STD[j])); while(STD[j][i]==0){i--;} if (i==p){ COM[j]=-1; } } for(j=1; j<=size(COM); j++) { if(COM[j]<0){COM=delete(COM,j);} } for(i=1; i<=size(COM_EX); i++) {ff=0; if(COM_EX[i]==0){ff=1;} else { for(j=1; j<=size(COM); j++) { if(COM_EX[i]==COM[j]){ff=1;} } } if(ff==0){COM_EX[i]=-1;} } L=STD_EX[1]; for(i=2; i<=size(COM_EX); i++) { if(COM_EX[i]!=-1){L=L,STD_EX[i];} } //////// split STD_EX in STD and the transformation matrix L=transpose(L); STD=L[st+1]; LT=L[1]; for(i=2; i<=s+st; i++) { if (i<=st) { LT=LT,L[i]; } if (i>st+1) { STD=STD,L[i]; } } STD=transpose(STD); STDFIN=matrix(STD); COMPARE=insert(COMPARE,STDFIN); LT=transpose(LT); ////////////////////// compute the transformation matrices if (flag mod 2 ==1) { TrafoL=transpose(LT)*TrafoL; dbprint(ppl-1,"Left transformation matrix:"); dbprint(ppl-1,TrafoL); @s = texobj("",matrix(TrafoL)); dbprint(ppl-2,"Left transformation matrix in TeX format:"); dbprint(ppl-2,@s); } else { TrafoR=TrafoR*involution(LT,inv); dbprint(ppl-1,"Right transformation matrix:"); dbprint(ppl-1,matrix(TrafoR)); @s = texobj("",TrafoR); dbprint(ppl-2,"Right transformation matrix in TeX format:"); dbprint(ppl-2,@s); } ///////////////////////// check whether the alg terminated ///////////////// if(size(COMPARE)>=3) { if(STD==COMPARE[3]){finish=1;} } ////////////////////////////////// change to the opposite module dbprint(ppl-1,"Matrix for the next step:"); dbprint(ppl-1,STD); dbprint(ppl-2,"Matrix for the next step in TeX format:"); @s = texobj("",matrix(STD)); dbprint(ppl-2,@s); TSTD=transpose(STD); STD=involution(TSTD,inv); flag++; dbprint(ppl,"Finished one while cycle"); } if (flag mod 2 ==0){ STD = involution(STD,inv);} else { STD = transpose(STD); } list REVERSE=TrafoL,STD,TrafoR; return(REVERSE); } static proc divideByContent(matrix M) { //find first entrie not equal to zero int i,k; k=1; vector CON; for(i=1;i<=ncols(M);i++) { if(leadcoef(M[i])!=0){CON=CON+leadcoef(M[i])*gen(k); k++;} } poly con=content(CON); matrix TL=1/con*freemodule(nrows(M)); return(TL); } /////interesting examples for smith//////////////// /* //static proc Ex_One_wheeled_bicycle() { ring RA=(0,m), D, lp; matrix bicycle[2][3]=(1+m)*D^2, D^2, 1, D^2, D^2-1, 0; list s=smith(bicycle, 1,0); print(s[2]); print(s[1]*bicycle*s[3]-s[2]); } //static proc Ex_RLC-circuit() { ring RA=(0,m,R1,R2,L,C), D, lp; matrix circuit[2][3]=D+1/(R1*C), 0, -1/(R1*C), 0, D+R2/L, -1/L; list s=smith(circuit, 1,0); print(s[2]); print(s[1]*circuit*s[3]-s[2]); } //static proc Ex_two_pendula() { ring RA=(0,m,M,L1,L2,m1,m2,g), D, lp; matrix pendula[3][4]=m1*L1*D^2,m2*L2*D^2,(M+m1+m2)*D^2,-1,m1*L1^2*D^2-m1*L1*g,0,m1*L1*D^2,0,0, m2*L2^2*D^2-m2*L2*g,m2*L2*D^2,0; list s=smith(pendula, 1,0); print(s[2]); print(s[1]*pendula*s[3]-s[2]); } //static proc Ex_linerized_satellite_in_a_circular_equatorial_orbit() { ring RA=(0,m,omega,r,a,b), D, lp; matrix satellite[4][6]= D,-1,0,0,0,0, -3*omega^2,D,0,-2*omega*r,-a/m,0, 0,0,D,-1,0,0, 0,2*omega/r,0,D,0,-b/(m*r); list s=smith(satellite, 1,0); print(s[2]); print(s[1]*satellite*s[3]-s[2]); } //static proc Ex_flexible_one_link_robot() { ring RA=(0,M11,M12,M13,M21,M22,M31,M33,K1,K2), D, lp; matrix robot[3][4]=M11*D^2,M12*D^2,M13*D^2,-1,M21*D^2,M22*D^2+K1,0,0,M31*D^2,0,M33*D^2+K2,0; list s=smith(robot, 1,0); print(s[2]); print(s[1]*robot*s[3]-s[2]); } */ /////interesting examples for jacobson//////////////// /* //static proc Ex_compare_output_with_maple_package_JanetOre() { ring w = 0,(x,d),Dp; def W=nc_algebra(1,1); setring W; matrix m[3][3]=[d2,d+1,0],[d+1,0,d3-x2*d],[2d+1, d3+d2, d2]; list J=jacobson(m,0); print(J[1]*m*J[3]); print(J[2]); print(J[1]); print(J[3]); print(J[1]*m*J[3]-J[2]); } // modification for shift algebra { ring w = 0,(x,s),Dp; def W=nc_algebra(1,s); setring W; matrix m[3][3]=[s^2,s+1,0],[s+1,0,s^3-x^2*s],[2*s+1, s^3+s^2, s^2]; list J=jacobson(m,0); print(J[1]*m*J[3]); print(J[2]); print(J[1]); print(J[3]); print(J[1]*m*J[3]-J[2]); } //static proc Ex_cyclic() { ring w = 0,(x,d),Dp; def W=nc_algebra(1,1); setring W; matrix m[3][3]=d,0,0,x*d+1,d,0,0,x*d,d; list J=jacobson(m,0); print(J[1]*m*J[3]); print(J[2]); print(J[1]); print(J[3]); print(J[1]*m*J[3]-J[2]); } // modification for shift algebra: gives a module with ann = s^2 { ring w = 0,(x,s),Dp; def W=nc_algebra(1,s); setring W; matrix m[3][3]=s,0,0,x*s+1,s,0,0,x*s,s; list J=jacobson(m,0); print(J[1]*m*J[3]); print(J[2]); print(J[1]); print(J[3]); print(J[1]*m*J[3]-J[2]); } // yet another modification for shift algebra: gives a module with ann = s // xs+1 is replaced by x*s+s { ring w = 0,(x,s),Dp; def W=nc_algebra(1,s); setring W; matrix m[3][3]=s,0,0,x*s+s,s,0,0,x*s,s; list J=jacobson(m,0); print(J[1]*m*J[3]); print(J[2]); print(J[1]); print(J[3]); print(J[1]*m*J[3]-J[2]); } // variations for above setup with shift algebra: // easy but instructive one: see the difference to Weyl case! matrix m[2][2]=s,x,0,s; print(m); // more interesting: matrix n[3][3]=s,x,0,0,s,x,s,0,x; // things blow up matrix m[2][3] = s,x^2*s,x^3*s,s*x^2,s*x+1,(x+1)^3; matrix m[2][3] = s,x^2*s,x^3*s,s*x^2,s*x+1,(x+s)^2; // variation matrix m[2][4] = s,x^2*s,x^3*s,s*x^2,s*x+1,(x+1)^3, (x+s)^2, x*s; // bug (matrix sizes matrix m[3][4] = s,x^2*s,x^3*s,s*x^2,s*x+1,(x+1)^3, (x+s)^2, x*s; // here the last row contains zeros // this one is quite nasty and time consumig matrix m[3][4] = s,x^2*s,x^3*s,s*x^2,s*x+1,(x+1)^3, (x+s)^2, x*s,x,x^2,x^3,s; // example from the paper: ring w = 0,(x,d),Dp; def W=nc_algebra(1,1); setring W; matrix m[2][2]=d^2-1,d+1,d^2+1,d-x; list J=jacobson(m,0); print(J[1]*m*J[3]); print(J[2]); print(J[1]); print(J[3]); print(J[1]*m*J[3]-J[2]); ring w2 = 0,(x,s),Dp; def W2=nc_algebra(1,s); setring W2; poly d = s; matrix m[2][2]=d^2-1,d+1,d^2+1,d-x; list J=jacobson(m,0); print(J[1]*m*J[3]); print(J[2]); print(J[1]); print(J[3]); print(J[1]*m*J[3]-J[2]); // here, both JNFs are cyclic // another example from the paper: ring w = 0,(x,d),Dp; def W=nc_algebra(1,1); setring W; matrix m[2][2]=-x*d+1, x^2*d, -d, x*d+1; list J=jacobson(m,0); print(J[1]*m*J[3]); print(J[2]); print(J[1]); print(J[3]); print(J[1]*m*J[3]-J[2]); ring w2 = 0,(x,s),Dp; def W2=nc_algebra(1,s); setring W2; poly d = s; matrix m[2][2]=-x*d+1, x^2*d, -d, x*d+1; list J=jacobson(m,0); print(J[1]*m*J[3]); print(J[2]); print(J[1]); print(J[3]); print(J[1]*m*J[3]-J[2]); // yet another example from the paper, also Middeke ring w = (0,y),(x,d),Dp; def W=nc_algebra(1,1); setring W; matrix m[2][2]=y^2*d^2+d+1, 1, x*d, x^2*d^2+d+y; list J=jacobson(m,0); print(J[1]*m*J[3]); print(J[2]); print(J[1]); print(J[3]); print(J[1]*m*J[3]-J[2]); */