[2c0350] | 1 | ////////////////////////////////////////////////////////////////////////////// |
---|
[e67578d] | 2 | version="$Id: jacobson.lib,v 1.10 2009-03-05 20:36:11 levandov Exp $"; |
---|
[2c0350] | 3 | category="System and Control Theory"; |
---|
| 4 | info=" |
---|
| 5 | LIBRARY: jacobson.lib Algorithms for Smith and Jacobson Normal Form |
---|
[e67578d] | 6 | AUTHOR: Kristina Schindelar, Kristina.Schindelar@math.rwth-aachen.de, |
---|
| 7 | @* Viktor Levandovskyy, levandov@math.rwth-aachen.de |
---|
[2c0350] | 8 | |
---|
[e67578d] | 9 | THEORY: We work over a ring R, that is an Euclidean principal ideal domain. |
---|
| 10 | @* If R is commutative, we suppose R to be a polynomial ring in one variable. |
---|
| 11 | @* If R is non-commutative, we suppose R to have two variables, say x and d. |
---|
| 12 | @* We treat then the basering as the Ore localization of R |
---|
[7f3ad4] | 13 | @* with respect to the mult. closed set S = K[x] without 0. |
---|
[e67578d] | 14 | @* Thus, we treat basering as principal ideal ring with d a polynomial |
---|
| 15 | @* variable and x an invertible one. |
---|
[2c0350] | 16 | @* Note, that in computations no division by x will actually happen. |
---|
[e67578d] | 17 | @* |
---|
| 18 | @* Given a rectangular matrix M over R, one can compute unimodular (that is |
---|
| 19 | @* invertible) square matrices U and V, such that U*M*V=D is a diagonal matrix. |
---|
| 20 | @* Depending on the ring, the diagonal entries of D have certain properties. |
---|
[2c0350] | 21 | |
---|
[7f3ad4] | 22 | REFERENCES: |
---|
[e67578d] | 23 | @* (1) N. Jacobson, 'The theory of rings', AMS, 1943. |
---|
| 24 | @* (2) Manuel Avelino Insua Hermo, 'Varias perspectives sobre las bases de Groebner : |
---|
| 25 | @* Forma normal de Smith, Algorithme de Berlekamp y algebras de Leibniz'. |
---|
| 26 | @* PhD thesis, Universidad de Santiago de Compostela, 2005. |
---|
[2c0350] | 27 | |
---|
| 28 | |
---|
| 29 | PROCEDURES: |
---|
| 30 | smith(M[,eng1,eng2]); compute the Smith Normal Form of M over commutative ring |
---|
[e67578d] | 31 | jacobson(M[,eng]); compute a weak Jacobson Normal Form of M over non-commutative ring |
---|
| 32 | divideUnits(L); create ones out of units in the output of smith or jacobson |
---|
| 33 | |
---|
[2c0350] | 34 | |
---|
| 35 | SEE ALSO: control_lib |
---|
| 36 | "; |
---|
| 37 | |
---|
| 38 | LIB "poly.lib"; |
---|
[e67578d] | 39 | LIB "involut.lib"; // involution |
---|
| 40 | LIB "dmodapp.lib"; // engine |
---|
| 41 | LIB "qhmoduli.lib"; // Min |
---|
| 42 | |
---|
| 43 | proc divideUnits(list L) |
---|
| 44 | "USAGE: divideUnits(L); list L |
---|
| 45 | RETURN: matrix or list of matrices |
---|
| 46 | ASSUME: L is an output of @code{smith} or a @code{jacobson} procedures, that is |
---|
| 47 | @* either L contains one rectangular matrix with elements only on the main diagonal |
---|
| 48 | @* or L consists of three matrices, where L[1] and L[3] are square invertible matrices |
---|
| 49 | @* while L[2] is a rectangular matrix with elements only on the main diagonal |
---|
| 50 | PURPOSE: divide out units from the diagonal and reflect this in transformation matrices |
---|
| 51 | EXAMPLE: example divideUnits; shows examples |
---|
| 52 | " |
---|
| 53 | { |
---|
| 54 | // check assumptions |
---|
| 55 | int s = size(L); |
---|
| 56 | if ( (s!=1) && (s!=3) ) |
---|
| 57 | { |
---|
| 58 | ERROR("The list has neither 1 nor 3 elements"); |
---|
| 59 | } |
---|
| 60 | // check sizes of matrices |
---|
| 61 | |
---|
| 62 | if (s==1) |
---|
| 63 | { |
---|
| 64 | list LL; LL[1] = L[1]; LL[2] = L[1]; |
---|
| 65 | L = LL; |
---|
| 66 | } |
---|
| 67 | |
---|
| 68 | |
---|
| 69 | // divide units out |
---|
| 70 | matrix M = L[2]; |
---|
| 71 | int ncM = ncols(M); int nrM = nrows(M); |
---|
| 72 | // matrix TM[nrM][nrM]; // m times m matrix |
---|
| 73 | matrix TM = matrix(freemodule(nrM)); |
---|
| 74 | int i; int nsize = Min(intvec(nrM,ncM)); |
---|
| 75 | poly p; number n; intvec lexp; |
---|
| 76 | |
---|
| 77 | for(i=1; i<=nsize; i++) |
---|
| 78 | { |
---|
| 79 | p = M[i,i]; |
---|
| 80 | lexp = leadexp(p); |
---|
| 81 | // TM[i,i] = 1; |
---|
| 82 | // check: is p a unit? |
---|
| 83 | if (p!=0) |
---|
| 84 | { |
---|
| 85 | if ( lexp == 0) |
---|
| 86 | { |
---|
| 87 | // hence p = n*1 |
---|
| 88 | n = leadcoef(p); |
---|
| 89 | TM[i,i] = 1/n; |
---|
| 90 | } |
---|
| 91 | } |
---|
| 92 | } |
---|
| 93 | int ppl = printlevel-voice+2; |
---|
| 94 | dbprint(ppl,"divideUnits: extra transformation matrix is: "); |
---|
| 95 | dbprint(ppl, TM); |
---|
| 96 | L[2] = TM*L[2]; |
---|
| 97 | if (s==3) |
---|
| 98 | { |
---|
| 99 | L[1] = TM*L[1]; |
---|
| 100 | return(L); |
---|
| 101 | } |
---|
| 102 | else |
---|
| 103 | { |
---|
| 104 | return(L[2]); |
---|
| 105 | } |
---|
| 106 | } |
---|
| 107 | example |
---|
| 108 | { "EXAMPLE:"; echo = 2; |
---|
| 109 | ring R=(0,m,M,L1,L2,m1,m2,g), D, lp; // two pendula example |
---|
| 110 | matrix P[3][4]=m1*L1*D^2,m2*L2*D^2,(M+m1+m2)*D^2,-1, |
---|
| 111 | m1*L1^2*D^2-m1*L1*g,0,m1*L1*D^2,0,0, |
---|
| 112 | m2*L2^2*D^2-m2*L2*g,m2*L2*D^2,0; |
---|
| 113 | list s=smith(P,1); // returns a list with 3 entries |
---|
| 114 | print(s[2]); // a diagonal form, close to the Smith form |
---|
| 115 | print(s[1]); // U, left transformation matrix |
---|
| 116 | list t = divideUnits(s); |
---|
| 117 | print(t[2]); // the Smith form of the matrix P |
---|
| 118 | print(t[1]); // U', modified left transformation matrix |
---|
| 119 | } |
---|
[2c0350] | 120 | |
---|
| 121 | proc smith(matrix MA, list #) |
---|
| 122 | "USAGE: smith(M[, eng1, eng2]); M matrix, eng1 and eng2 are optional integers |
---|
[e67578d] | 123 | RETURN: matrix or list of matrices, depending on arguments |
---|
| 124 | ASSUME: Basering is a commutative polynomial ring in one variable |
---|
| 125 | PURPOSE: compute the Smith Normal Form of M with (optionally) transformation matrices |
---|
| 126 | THEORY: Groebner bases are used for the Smith form like in (2). |
---|
| 127 | NOTE: By default, just the Smith normal form of M is returned. |
---|
| 128 | @* If the optional integer @code{eng1} is non-zero, the list {U,D,V} is returned |
---|
| 129 | @* where U*M*V = D and the diagonal field entries of D are not normalized. |
---|
| 130 | @* The normalization of the latter can be done with the 'divideUnits' procedure. |
---|
| 131 | @* U and V above are square unimodular (invertible) matrices. |
---|
| 132 | @* Note, that the procedure works for a rectangular matrix M. |
---|
| 133 | @* |
---|
| 134 | @* The optional integer @code{eng2} determines the Groebner basis engine: |
---|
| 135 | @* 0 (default) ensures the use of 'slimgb' , otherwise 'std' is used. |
---|
[2c0350] | 136 | DISPLAY: If @code{printlevel}=1, progress debug messages will be printed, |
---|
| 137 | @* if @code{printlevel}>=2, all the debug messages will be printed. |
---|
| 138 | EXAMPLE: example smith; shows examples |
---|
[e67578d] | 139 | SEE ALSO: divideUnits, jacobson |
---|
[2c0350] | 140 | " |
---|
[7f3ad4] | 141 | { |
---|
[2c0350] | 142 | def R = basering; |
---|
| 143 | // check assume: R must be a polynomial ring in 1 variable |
---|
| 144 | setring R; |
---|
| 145 | if (nvars(R) >1 ) |
---|
| 146 | { |
---|
| 147 | ERROR("Basering must have exactly one variable"); |
---|
| 148 | } |
---|
| 149 | |
---|
| 150 | int eng = 0; |
---|
| 151 | int BASIS; |
---|
| 152 | if ( size(#)>0 ) |
---|
| 153 | { |
---|
| 154 | eng=1; |
---|
| 155 | if (typeof(#[1])=="int") |
---|
| 156 | { |
---|
| 157 | eng=int(#[1]); // zero can also happen |
---|
| 158 | } |
---|
| 159 | } |
---|
| 160 | if (size(#)==2) |
---|
| 161 | { |
---|
| 162 | BASIS=#[2]; |
---|
[7f3ad4] | 163 | } |
---|
[2c0350] | 164 | else{BASIS=0;} |
---|
[7f3ad4] | 165 | |
---|
[2c0350] | 166 | |
---|
| 167 | int ROW=ncols(MA); |
---|
| 168 | int COL=nrows(MA); |
---|
| 169 | |
---|
| 170 | //generate a module consisting of the columns of MA |
---|
| 171 | module m=MA[1]; |
---|
| 172 | int i; |
---|
| 173 | for(i=2;i<=ROW;i++) |
---|
| 174 | { |
---|
| 175 | m=m,MA[i]; |
---|
| 176 | } |
---|
| 177 | |
---|
| 178 | //if MA eqauls the zero matrix give back MA |
---|
| 179 | if ( (size(module(m))==0) and (size(transpose(module(m)))==0) ) |
---|
| 180 | { |
---|
| 181 | module L=freemodule(COL); |
---|
| 182 | matrix LM=L; |
---|
| 183 | L=freemodule(ROW); |
---|
| 184 | matrix RM=L; |
---|
| 185 | list RUECK=RM; |
---|
| 186 | RUECK=insert(RUECK,MA); |
---|
| 187 | RUECK=insert(RUECK,LM); |
---|
| 188 | return(RUECK); |
---|
| 189 | } |
---|
| 190 | |
---|
| 191 | if(eng==1) |
---|
| 192 | { |
---|
[7f3ad4] | 193 | list rueckLI=diagonal_with_trafo(R,MA,BASIS); |
---|
[2c0350] | 194 | list rueckLII=divisibility(rueckLI[2]); |
---|
| 195 | matrix CON=divideByContent(rueckLII[2]); |
---|
| 196 | list rueckL=CON*rueckLII[1]*rueckLI[1], CON*rueckLII[2], rueckLI[3]*rueckLII[3]; |
---|
| 197 | return(rueckL); |
---|
| 198 | } |
---|
| 199 | else |
---|
| 200 | { |
---|
[7f3ad4] | 201 | matrix rueckm=diagonal_without_trafo(R,MA,BASIS); |
---|
[2c0350] | 202 | list rueckL=divisibility(rueckm); |
---|
| 203 | matrix CON=divideByContent(rueckm); |
---|
| 204 | rueckm=CON*rueckL[2]; |
---|
| 205 | return(rueckm); |
---|
| 206 | } |
---|
| 207 | } |
---|
| 208 | example |
---|
| 209 | { "EXAMPLE:"; echo = 2; |
---|
| 210 | ring r = 0,x,Dp; |
---|
| 211 | matrix m[3][2]=x, x^4+x^2+21, x^4+x^2+x, x^3+x, 4*x^2+x, x; |
---|
| 212 | list s=smith(m,1); |
---|
[e67578d] | 213 | print(s[2]); // non-normalized Smith form of m |
---|
[2c0350] | 214 | print(s[1]*m*s[3] - s[2]); // check U*M*V = D |
---|
[e67578d] | 215 | list t = divideUnits(s); |
---|
| 216 | print(t[2]); // the Smith form of m |
---|
[2c0350] | 217 | } |
---|
| 218 | |
---|
| 219 | static proc diagonal_with_trafo( R, matrix MA, int B) |
---|
| 220 | { |
---|
| 221 | |
---|
| 222 | int ppl = printlevel-voice+2; |
---|
| 223 | |
---|
| 224 | int BASIS=B; |
---|
| 225 | int ROW=ncols(MA); |
---|
| 226 | int COL=nrows(MA); |
---|
| 227 | module m=MA[1]; |
---|
| 228 | int i,j,k; |
---|
| 229 | for(i=2;i<=ROW;i++) |
---|
| 230 | { |
---|
| 231 | m=m,MA[i]; |
---|
| 232 | } |
---|
| 233 | |
---|
| 234 | |
---|
| 235 | //add zero rows or columns |
---|
| 236 | //add zero rows or columns |
---|
| 237 | int adrow=0; |
---|
| 238 | for(i=1;i<=COL;i++) |
---|
| 239 | { |
---|
| 240 | k=0; |
---|
| 241 | for(j=1;j<=ROW;j++) |
---|
| 242 | { |
---|
| 243 | if(MA[i,j]!=0){k=1;} |
---|
| 244 | } |
---|
| 245 | if(k==0){adrow++;} |
---|
| 246 | } |
---|
[7f3ad4] | 247 | |
---|
[2c0350] | 248 | m=transpose(m); |
---|
| 249 | for(i=1;i<=adrow;i++){m=m,0;} |
---|
| 250 | m=transpose(m); |
---|
| 251 | |
---|
| 252 | list RINGLIST=ringlist(R); |
---|
| 253 | list o="C",0; |
---|
| 254 | list oo="lp",1; |
---|
| 255 | list ORD=o,oo; |
---|
| 256 | RINGLIST[3]=ORD; |
---|
| 257 | def r=ring(RINGLIST); |
---|
| 258 | setring r; |
---|
| 259 | //fix the required ordering |
---|
| 260 | map MAP=R,var(1); |
---|
| 261 | module m=MAP(m); |
---|
| 262 | |
---|
| 263 | int flag=1; ///////////////counts if the underlying ring is r (flag mod 2 == 1) or ro (flag mod 2 == 0) |
---|
| 264 | |
---|
| 265 | module TrafoL=freemodule(COL); |
---|
| 266 | module TrafoR=freemodule(ROW); |
---|
| 267 | module EXL=freemodule(COL); //because we start with transpose(m) |
---|
| 268 | module EXR=freemodule(ROW); |
---|
| 269 | |
---|
| 270 | option(redSB); |
---|
| 271 | option(redTail); |
---|
| 272 | module STD_EX; |
---|
| 273 | module Trafo; |
---|
| 274 | |
---|
| 275 | int s,st,p,ff; |
---|
[7f3ad4] | 276 | |
---|
[2c0350] | 277 | module LT,TSTD; |
---|
| 278 | module STD=transpose(m); |
---|
| 279 | int finish=0; |
---|
| 280 | int fehlpos; |
---|
| 281 | list pos; |
---|
| 282 | matrix END; |
---|
| 283 | matrix ENDSTD; |
---|
[7f3ad4] | 284 | matrix STDFIN; |
---|
| 285 | STDFIN=STD; |
---|
[2c0350] | 286 | list COMPARE=STDFIN; |
---|
| 287 | |
---|
| 288 | while(finish==0) |
---|
| 289 | { |
---|
| 290 | dbprint(ppl,"Going into the while cycle"); |
---|
| 291 | |
---|
| 292 | if(flag mod 2==1) |
---|
| 293 | { |
---|
| 294 | STD_EX=EXL,transpose(STD); |
---|
[7f3ad4] | 295 | } |
---|
| 296 | else |
---|
[2c0350] | 297 | { |
---|
| 298 | STD_EX=EXR,transpose(STD); |
---|
| 299 | } |
---|
[7f3ad4] | 300 | dbprint(ppl,"Computing Groebner basis: start"); |
---|
| 301 | dbprint(ppl-1,STD); |
---|
[2c0350] | 302 | STD=engine(STD,BASIS); |
---|
[7f3ad4] | 303 | dbprint(ppl,"Computing Groebner basis: finished"); |
---|
| 304 | dbprint(ppl-1,STD); |
---|
[2c0350] | 305 | |
---|
| 306 | if(flag mod 2==1){s=ROW; st=COL;}else{s=COL; st=ROW;} |
---|
[7f3ad4] | 307 | dbprint(ppl,"Computing Groebner basis for transformation matrix: start"); |
---|
| 308 | dbprint(ppl-1,STD_EX); |
---|
[2c0350] | 309 | STD_EX=transpose(STD_EX); |
---|
| 310 | STD_EX=engine(STD_EX,BASIS); |
---|
[7f3ad4] | 311 | dbprint(ppl,"Computing Groebner basis for transformation matrix: finished"); |
---|
| 312 | dbprint(ppl-1,STD_EX); |
---|
[2c0350] | 313 | |
---|
| 314 | //////// split STD_EX in STD and the transformation matrix |
---|
| 315 | STD_EX=transpose(STD_EX); |
---|
| 316 | STD=STD_EX[st+1]; |
---|
| 317 | LT=STD_EX[1]; |
---|
| 318 | |
---|
| 319 | ENDSTD=STD_EX; |
---|
| 320 | for(i=2; i<=ncols(ENDSTD); i++) |
---|
| 321 | { |
---|
| 322 | if (i<=st) |
---|
| 323 | { |
---|
| 324 | LT=LT,STD_EX[i]; |
---|
| 325 | } |
---|
| 326 | if (i>st+1) |
---|
| 327 | { |
---|
| 328 | STD=STD,STD_EX[i]; |
---|
| 329 | } |
---|
| 330 | } |
---|
| 331 | |
---|
| 332 | STD=transpose(STD); |
---|
| 333 | LT=transpose(LT); |
---|
| 334 | |
---|
| 335 | ////////////////////// compute the transformation matrices |
---|
[7f3ad4] | 336 | |
---|
[2c0350] | 337 | if (flag mod 2 ==1) |
---|
[7f3ad4] | 338 | { |
---|
[2c0350] | 339 | TrafoL=transpose(LT)*TrafoL; |
---|
| 340 | } |
---|
[7f3ad4] | 341 | else |
---|
| 342 | { |
---|
[2c0350] | 343 | TrafoR=TrafoR*LT; |
---|
| 344 | } |
---|
| 345 | |
---|
| 346 | |
---|
[7f3ad4] | 347 | STDFIN=STD; |
---|
[2c0350] | 348 | /////////////////////////////////// check if the D-th row is finished /////////////////////////////////// |
---|
| 349 | COMPARE=insert(COMPARE,STDFIN); |
---|
[7f3ad4] | 350 | if(size(COMPARE)>=3) |
---|
| 351 | { |
---|
[2c0350] | 352 | if(STD==COMPARE[3]){finish=1;} |
---|
| 353 | } |
---|
| 354 | ////////////////////////////////// change to the opposite module |
---|
| 355 | TSTD=transpose(STD); |
---|
| 356 | STD=TSTD; |
---|
| 357 | flag++; |
---|
| 358 | dbprint(ppl,"Finished one while cycle"); |
---|
| 359 | } |
---|
| 360 | |
---|
| 361 | |
---|
| 362 | if (flag mod 2!=0) { STD=transpose(STD); } |
---|
[7f3ad4] | 363 | |
---|
| 364 | |
---|
[2c0350] | 365 | //zero colums to the end |
---|
| 366 | matrix STDMM=STD; |
---|
| 367 | pos=list(); |
---|
| 368 | fehlpos=0; |
---|
| 369 | while( size(STD)+fehlpos-ncols(STDMM) < 0) |
---|
| 370 | { |
---|
| 371 | for(i=1; i<=ncols(STDMM); i++) |
---|
| 372 | { |
---|
[7f3ad4] | 373 | ff=0; |
---|
[2c0350] | 374 | for(j=1; j<=nrows(STDMM); j++) |
---|
| 375 | { |
---|
| 376 | if (STD[j,i]==0) { ff++; } |
---|
| 377 | } |
---|
| 378 | if(ff==nrows(STDMM)) |
---|
| 379 | { |
---|
| 380 | pos=insert(pos,i); fehlpos++; |
---|
| 381 | } |
---|
| 382 | } |
---|
| 383 | } |
---|
| 384 | int fehlposc=fehlpos; |
---|
| 385 | module SORT; |
---|
[7f3ad4] | 386 | for(i=1; i<=fehlpos; i++) |
---|
[2c0350] | 387 | { |
---|
| 388 | SORT=gen(2); |
---|
| 389 | for (j=3;j<=ROW;j++) |
---|
[7f3ad4] | 390 | { |
---|
[2c0350] | 391 | SORT=SORT,gen(j); |
---|
| 392 | } |
---|
| 393 | SORT=SORT,gen(1); |
---|
| 394 | STD=STD*SORT; |
---|
| 395 | TrafoR=TrafoR*SORT; |
---|
| 396 | } |
---|
| 397 | |
---|
| 398 | //zero rows to the end |
---|
| 399 | STDMM=transpose(STD); |
---|
| 400 | pos=list(); |
---|
| 401 | fehlpos=0; |
---|
| 402 | while( size(transpose(STD))+fehlpos-ncols(STDMM) < 0) |
---|
[7f3ad4] | 403 | { |
---|
[2c0350] | 404 | for(i=1; i<=ncols(STDMM); i++) |
---|
| 405 | { |
---|
| 406 | ff=0; for(j=1; j<=nrows(STDMM); j++) |
---|
| 407 | { |
---|
| 408 | if(transpose(STD)[j,i]==0){ ff++;}} |
---|
| 409 | if(ff==nrows(STDMM)) { pos=insert(pos,i); fehlpos++; } |
---|
| 410 | } |
---|
| 411 | } |
---|
| 412 | int fehlposr=fehlpos; |
---|
| 413 | |
---|
[7f3ad4] | 414 | for(i=1; i<=fehlpos; i++) |
---|
[2c0350] | 415 | { |
---|
| 416 | SORT=gen(2); |
---|
| 417 | for(j=3;j<=COL;j++){SORT=SORT,gen(j);} |
---|
| 418 | SORT=SORT,gen(1); |
---|
| 419 | SORT=transpose(SORT); |
---|
| 420 | STD=SORT*STD; |
---|
| 421 | TrafoL=SORT*TrafoL; |
---|
| 422 | } |
---|
[7f3ad4] | 423 | |
---|
[2c0350] | 424 | setring R; |
---|
| 425 | map MAPinv=r,var(1); |
---|
| 426 | module STD=MAPinv(STD); |
---|
| 427 | module TrafoL=MAPinv(TrafoL); |
---|
| 428 | matrix TrafoLM=TrafoL; |
---|
| 429 | module TrafoR=MAPinv(TrafoR); |
---|
| 430 | matrix TrafoRM=TrafoR; |
---|
| 431 | matrix STDM=STD; |
---|
| 432 | |
---|
[7f3ad4] | 433 | //Test |
---|
[2c0350] | 434 | if(TrafoLM*m*TrafoRM!=STDM){ return(0); } |
---|
[7f3ad4] | 435 | |
---|
[2c0350] | 436 | list RUECK=TrafoRM; |
---|
| 437 | RUECK=insert(RUECK,STDM); |
---|
| 438 | RUECK=insert(RUECK,TrafoLM); |
---|
| 439 | return(RUECK); |
---|
| 440 | } |
---|
| 441 | |
---|
| 442 | |
---|
| 443 | static proc divisibility(matrix M) |
---|
| 444 | { |
---|
| 445 | matrix STDM=M; |
---|
| 446 | int i,j; |
---|
| 447 | int ROW=nrows(M); |
---|
| 448 | int COL=ncols(M); |
---|
| 449 | module TrafoR=freemodule(COL); |
---|
| 450 | module TrafoL=freemodule(ROW); |
---|
| 451 | module SORT; |
---|
| 452 | matrix TrafoRM=TrafoR; |
---|
| 453 | matrix TrafoLM=TrafoL; |
---|
| 454 | list posdeg0; |
---|
| 455 | int posdeg=0; |
---|
| 456 | int act; |
---|
| 457 | int Sort=ROW; |
---|
| 458 | if(size(std(STDM))!=0) |
---|
| 459 | { while( size(transpose(STDM)[Sort])==0 ){Sort--;}} |
---|
| 460 | |
---|
| 461 | for(i=1;i<=Sort ;i++) |
---|
| 462 | { |
---|
| 463 | if(leadexp(STDM[i,i])==0){posdeg0=insert(posdeg0,i);posdeg++;} |
---|
| 464 | } |
---|
| 465 | //entries of degree 0 at the beginning |
---|
| 466 | for(i=1; i<=posdeg; i++) |
---|
| 467 | { |
---|
| 468 | act=posdeg0[i]; |
---|
| 469 | SORT=gen(act); |
---|
| 470 | for(j=1; j<=COL; j++){if(j!=act){SORT=SORT,gen(j);}} |
---|
| 471 | STDM=STDM*SORT; |
---|
| 472 | TrafoRM=TrafoRM*SORT; |
---|
| 473 | SORT=gen(act); |
---|
| 474 | for(j=1; j<=ROW; j++){if(j!=act){SORT=SORT,gen(j);}} |
---|
| 475 | STDM=transpose(SORT)*STDM; |
---|
| 476 | TrafoLM=transpose(SORT)*TrafoLM; |
---|
| 477 | } |
---|
| 478 | |
---|
| 479 | poly G; |
---|
| 480 | module UNITL=freemodule(ROW); |
---|
| 481 | matrix GCDL=UNITL; |
---|
| 482 | module UNITR=freemodule(COL); |
---|
| 483 | matrix GCDR=UNITR; |
---|
| 484 | for(i=posdeg+1; i<=Sort; i++) |
---|
| 485 | { |
---|
| 486 | for(j=i+1; j<=Sort; j++) |
---|
| 487 | { |
---|
[7f3ad4] | 488 | GCDL=UNITL; |
---|
| 489 | GCDR=UNITR; |
---|
| 490 | G=gcd(STDM[i,i],STDM[j,j]); |
---|
| 491 | ideal Z=STDM[i,i],STDM[j,j]; |
---|
| 492 | matrix T=lift(Z,G); |
---|
| 493 | GCDL[i,i]=T[1,1]; |
---|
| 494 | GCDL[i,j]=T[2,1]; |
---|
| 495 | GCDL[j,i]=-STDM[j,j]/G; |
---|
| 496 | GCDL[j,j]=STDM[i,i]/G; |
---|
| 497 | GCDR[i,j]=T[2,1]*STDM[j,j]/G; |
---|
| 498 | GCDR[j,j]=T[2,1]*STDM[j,j]/G-1; |
---|
| 499 | GCDR[j,i]=1; |
---|
| 500 | STDM=GCDL*STDM*GCDR; |
---|
| 501 | TrafoLM=GCDL*TrafoLM; |
---|
| 502 | TrafoRM=TrafoRM*GCDR; |
---|
[2c0350] | 503 | } |
---|
| 504 | } |
---|
| 505 | list RUECK=TrafoRM; |
---|
| 506 | RUECK=insert(RUECK,STDM); |
---|
| 507 | RUECK=insert(RUECK,TrafoLM); |
---|
| 508 | return(RUECK); |
---|
| 509 | } |
---|
| 510 | |
---|
| 511 | static proc diagonal_without_trafo( R, matrix MA, int B) |
---|
| 512 | { |
---|
[7f3ad4] | 513 | int ppl = printlevel-voice+2; |
---|
[2c0350] | 514 | |
---|
| 515 | int BASIS=B; |
---|
| 516 | int ROW=ncols(MA); |
---|
| 517 | int COL=nrows(MA); |
---|
| 518 | module m=MA[1]; |
---|
| 519 | int i; |
---|
| 520 | for(i=2;i<=ROW;i++) |
---|
| 521 | {m=m,MA[i];} |
---|
| 522 | |
---|
| 523 | |
---|
| 524 | list RINGLIST=ringlist(R); |
---|
| 525 | list o="C",0; |
---|
| 526 | list oo="lp",1; |
---|
| 527 | list ORD=o,oo; |
---|
| 528 | RINGLIST[3]=ORD; |
---|
| 529 | def r=ring(RINGLIST); |
---|
| 530 | setring r; |
---|
| 531 | //RICHTIGE ORDNUNG MACHEN |
---|
| 532 | map MAP=R,var(1); |
---|
| 533 | module m=MAP(m); |
---|
| 534 | |
---|
| 535 | int flag=1; ///////////////counts if the underlying ring is r (flag mod 2 == 1) or ro (flag mod 2 == 0) |
---|
| 536 | |
---|
| 537 | |
---|
| 538 | int act, j, ff; |
---|
| 539 | option(redSB); |
---|
| 540 | option(redTail); |
---|
| 541 | |
---|
| 542 | |
---|
| 543 | module STD=transpose(m); |
---|
| 544 | module TSTD; |
---|
| 545 | int finish=0; |
---|
| 546 | matrix STDFIN; |
---|
[7f3ad4] | 547 | STDFIN=STD; |
---|
[2c0350] | 548 | list COMPARE=STDFIN; |
---|
| 549 | |
---|
| 550 | while(finish==0) |
---|
| 551 | { |
---|
| 552 | dbprint(ppl,"Going into the while cycle"); |
---|
| 553 | dbprint(ppl,"Computing Groebner basis: start"); |
---|
| 554 | dbprint(ppl-1,STD); |
---|
| 555 | STD=engine(STD,BASIS); |
---|
| 556 | dbprint(ppl,"Computing Groebner basis: finished"); |
---|
| 557 | dbprint(ppl-1,STD); |
---|
[7f3ad4] | 558 | STDFIN=STD; |
---|
[2c0350] | 559 | /////////////////////////////////// check if the D-th row is finished /////////////////////////////////// |
---|
| 560 | COMPARE=insert(COMPARE,STDFIN); |
---|
[7f3ad4] | 561 | if(size(COMPARE)>=3) |
---|
| 562 | { |
---|
[2c0350] | 563 | if(STD==COMPARE[3]){finish=1;} |
---|
[7f3ad4] | 564 | } |
---|
| 565 | ////////////////////////////////// change to the opposite module |
---|
[2c0350] | 566 | |
---|
[7f3ad4] | 567 | TSTD=transpose(STD); |
---|
| 568 | STD=TSTD; |
---|
| 569 | flag++; |
---|
[2c0350] | 570 | dbprint(ppl,"Finished one while cycle"); |
---|
| 571 | } |
---|
| 572 | |
---|
| 573 | matrix STDMM=STD; |
---|
| 574 | list pos=list(); |
---|
| 575 | int fehlpos=0; |
---|
| 576 | while( size(STD)+fehlpos-ncols(STDMM) < 0) |
---|
| 577 | { |
---|
| 578 | for(i=1; i<=ncols(STDMM); i++) |
---|
| 579 | { |
---|
[7f3ad4] | 580 | ff=0; |
---|
[2c0350] | 581 | for(j=1; j<=nrows(STDMM); j++) |
---|
| 582 | { |
---|
| 583 | if (STD[j,i]==0) { ff++; } |
---|
| 584 | } |
---|
| 585 | if(ff==nrows(STDMM)) |
---|
| 586 | { |
---|
| 587 | pos=insert(pos,i); fehlpos++; |
---|
| 588 | } |
---|
| 589 | } |
---|
| 590 | } |
---|
| 591 | int fehlposc=fehlpos; |
---|
| 592 | module SORT; |
---|
[7f3ad4] | 593 | for(i=1; i<=fehlpos; i++) |
---|
[2c0350] | 594 | { |
---|
| 595 | SORT=gen(2); |
---|
| 596 | for (j=3;j<=ROW;j++) |
---|
[7f3ad4] | 597 | { |
---|
[2c0350] | 598 | SORT=SORT,gen(j); |
---|
| 599 | } |
---|
| 600 | SORT=SORT,gen(1); |
---|
| 601 | STD=STD*SORT; |
---|
| 602 | } |
---|
| 603 | |
---|
| 604 | //zero rows to the end |
---|
| 605 | STDMM=transpose(STD); |
---|
| 606 | pos=list(); |
---|
| 607 | fehlpos=0; |
---|
| 608 | while( size(transpose(STD))+fehlpos-ncols(STDMM) < 0) |
---|
[7f3ad4] | 609 | { |
---|
[2c0350] | 610 | for(i=1; i<=ncols(STDMM); i++) |
---|
| 611 | { |
---|
| 612 | ff=0; for(j=1; j<=nrows(STDMM); j++) |
---|
| 613 | { |
---|
| 614 | if(transpose(STD)[j,i]==0){ ff++;}} |
---|
| 615 | if(ff==nrows(STDMM)) { pos=insert(pos,i); fehlpos++; } |
---|
| 616 | } |
---|
| 617 | } |
---|
| 618 | int fehlposr=fehlpos; |
---|
| 619 | |
---|
[7f3ad4] | 620 | for(i=1; i<=fehlpos; i++) |
---|
[2c0350] | 621 | { |
---|
| 622 | SORT=gen(2); |
---|
| 623 | for(j=3;j<=COL;j++){SORT=SORT,gen(j);} |
---|
| 624 | SORT=SORT,gen(1); |
---|
| 625 | SORT=transpose(SORT); |
---|
| 626 | STD=SORT*STD; |
---|
| 627 | } |
---|
[7f3ad4] | 628 | |
---|
[2c0350] | 629 | //add zero rows or columns |
---|
| 630 | |
---|
| 631 | int adrow=COL-size(transpose(STD)); |
---|
| 632 | int adcol=ROW-size(STD); |
---|
| 633 | |
---|
| 634 | for(i=1;i<=adcol;i++){STD=STD,0;} |
---|
| 635 | STD=transpose(STD); |
---|
| 636 | for(i=1;i<=adrow;i++){STD=STD,0;} |
---|
| 637 | STD=transpose(STD); |
---|
| 638 | |
---|
| 639 | setring R; |
---|
| 640 | map MAPinv=r,var(1); |
---|
| 641 | module STD=MAPinv(STD); |
---|
| 642 | matrix STDM=STD; |
---|
| 643 | return(STDM); |
---|
| 644 | } |
---|
| 645 | |
---|
| 646 | |
---|
[e67578d] | 647 | // VL : engine replaced by the one from dmodapp.lib |
---|
| 648 | // cases are changed |
---|
| 649 | |
---|
| 650 | // static proc engine(module I, int i) |
---|
| 651 | // { |
---|
| 652 | // module J; |
---|
| 653 | // if (i==0) |
---|
| 654 | // { |
---|
| 655 | // J = std(I); |
---|
| 656 | // } |
---|
| 657 | // if (i==1) |
---|
| 658 | // { |
---|
| 659 | // J = groebner(I); |
---|
| 660 | // } |
---|
| 661 | // if (i==2) |
---|
| 662 | // { |
---|
| 663 | // J = slimgb(I); |
---|
| 664 | // } |
---|
| 665 | // return(J); |
---|
| 666 | // } |
---|
[2c0350] | 667 | |
---|
| 668 | proc jacobson(matrix MA, list #) |
---|
| 669 | "USAGE: jacobson(M, eng); M matrix, eng an optional int |
---|
[7f3ad4] | 670 | RETURN: list |
---|
[e67578d] | 671 | ASSUME: Basering is a (non-commutative) ring in two variables. |
---|
| 672 | PURPOSE: compute a weak Jacobson Normal Form of M over the basering |
---|
| 673 | THEORY: Groebner bases and involutions are used, generalizing an idea of (2) |
---|
[2c0350] | 674 | NOTE: A list L of matrices {U,D,V} is returned. That is L[1]*M*L[3]=L[2], where |
---|
| 675 | @* L[2] is a diagonal matrix and L[1], L[3] square invertible (unimodular) matrices. |
---|
| 676 | @* Note, that M can be rectangular. |
---|
[e67578d] | 677 | @* The optional integer @code{eng2} determines the Groebner basis engine: |
---|
| 678 | @* 0 (default) ensures the use of 'slimgb' , otherwise 'std' is used. |
---|
[2c0350] | 679 | DISPLAY: If @code{printlevel}=1, progress debug messages will be printed, |
---|
| 680 | @* if @code{printlevel}>=2, all the debug messages will be printed. |
---|
| 681 | EXAMPLE: example jacobson; shows examples |
---|
[e67578d] | 682 | SEE ALSO: divideUnits, smith |
---|
[2c0350] | 683 | " |
---|
[7f3ad4] | 684 | { |
---|
[2c0350] | 685 | def R = basering; |
---|
| 686 | // check assume: R must be a polynomial ring in 2 variables |
---|
| 687 | setring R; |
---|
| 688 | if ( (nvars(R) !=2 ) ) |
---|
| 689 | { |
---|
| 690 | ERROR("Basering must have exactly two variables"); |
---|
| 691 | } |
---|
| 692 | |
---|
| 693 | // check if MA is zero matrix and return corr. U,V |
---|
| 694 | if ( (size(module(MA))==0) and (size(transpose(module(MA)))==0) ) |
---|
| 695 | { |
---|
| 696 | int nr = nrows(MA); |
---|
| 697 | int nc = ncols(MA); |
---|
| 698 | matrix U = matrix(freemodule(nr)); |
---|
| 699 | matrix V = matrix(freemodule(nc)); |
---|
| 700 | list L; L[1]=U; L[2] = MA; L[3] = V; |
---|
| 701 | return(L); |
---|
| 702 | } |
---|
| 703 | |
---|
| 704 | int B=0; |
---|
| 705 | if ( size(#)>0 ) |
---|
| 706 | { |
---|
| 707 | B=1; |
---|
| 708 | if (typeof(#[1])=="int") |
---|
| 709 | { |
---|
| 710 | B=int(#[1]); // zero can also happen |
---|
| 711 | } |
---|
| 712 | } |
---|
| 713 | |
---|
| 714 | //change ring |
---|
| 715 | list RINGLIST=ringlist(R); |
---|
| 716 | list o="C",0; |
---|
[fec31e1] | 717 | intvec v=0,1; |
---|
[2c0350] | 718 | list oo="a",v; |
---|
| 719 | v=1,1; |
---|
| 720 | list ooo="lp",v; |
---|
| 721 | list ORD=o,oo,ooo; |
---|
| 722 | RINGLIST[3]=ORD; |
---|
| 723 | def r=ring(RINGLIST); |
---|
| 724 | setring r; |
---|
[7f3ad4] | 725 | |
---|
[2c0350] | 726 | //fix the required ordering |
---|
| 727 | map MAP=R,var(1),var(2); |
---|
| 728 | matrix M=MAP(MA); |
---|
| 729 | |
---|
| 730 | module TrafoL, TrafoR; |
---|
| 731 | list TRIANGLE; |
---|
| 732 | TRIANGLE=triangle(M,B); |
---|
| 733 | TrafoL=TRIANGLE[1]; |
---|
| 734 | TrafoR=TRIANGLE[3]; |
---|
| 735 | module m=TRIANGLE[2]; |
---|
| 736 | |
---|
[7f3ad4] | 737 | //back to the ring |
---|
[2c0350] | 738 | setring R; |
---|
| 739 | map MAPR=r,var(1),var(2); |
---|
| 740 | module ma=MAPR(m); |
---|
| 741 | matrix MAA=ma; |
---|
| 742 | module TL=MAPR(TrafoL); |
---|
| 743 | module TR=MAPR(TrafoR); |
---|
[fec31e1] | 744 | matrix TRR=TR; |
---|
[2c0350] | 745 | matrix CON=divideByContent(MAA); |
---|
[7f3ad4] | 746 | |
---|
[fec31e1] | 747 | list RUECK=CON*TL, CON*MAA, TRR; |
---|
[2c0350] | 748 | return(RUECK); |
---|
| 749 | } |
---|
| 750 | example |
---|
[7f3ad4] | 751 | { |
---|
[2c0350] | 752 | "EXAMPLE:"; echo = 2; |
---|
| 753 | ring r = 0,(x,d),Dp; |
---|
| 754 | def R = nc_algebra(1,1); setring R; // the 1st Weyl algebra |
---|
| 755 | matrix m[2][2] = d,x,0,d; print(m); |
---|
| 756 | list J = jacobson(m); // returns a list with 3 entries |
---|
| 757 | print(J[2]); // a Jacobson Form D for m |
---|
| 758 | print(J[1]*m*J[3] - J[2]); // check that U*M*V = D |
---|
| 759 | /* now, let us do the same for the shift algebra */ |
---|
| 760 | ring r2 = 0,(x,s),Dp; |
---|
| 761 | def R2 = nc_algebra(1,s); setring R2; // the 1st shift algebra |
---|
| 762 | matrix m[2][2] = s,x,0,s; print(m); // matrix of the same for as above |
---|
[7f3ad4] | 763 | list J = jacobson(m); |
---|
[2c0350] | 764 | print(J[2]); // a Jacobson Form D, quite different from above |
---|
| 765 | print(J[1]*m*J[3] - J[2]); // check that U*M*V = D |
---|
| 766 | } |
---|
| 767 | |
---|
| 768 | static proc triangle( matrix MA, int B) |
---|
| 769 | { |
---|
[7f3ad4] | 770 | int ppl = printlevel-voice+2; |
---|
[2c0350] | 771 | |
---|
| 772 | map inv=ncdetection(); |
---|
| 773 | int ROW=ncols(MA); |
---|
| 774 | int COL=nrows(MA); |
---|
| 775 | |
---|
| 776 | //generate a module consisting of the columns of MA |
---|
| 777 | module m=MA[1]; |
---|
| 778 | int i,j,s,st,p,k,ff,ex, nz, ii,nextp; |
---|
| 779 | for(i=2;i<=ROW;i++) |
---|
| 780 | { |
---|
| 781 | m=m,MA[i]; |
---|
| 782 | } |
---|
| 783 | int BASIS=B; |
---|
| 784 | |
---|
| 785 | //add zero rows or columns |
---|
| 786 | int adrow=0; |
---|
| 787 | for(i=1;i<=COL;i++) |
---|
| 788 | { |
---|
| 789 | k=0; |
---|
| 790 | for(j=1;j<=ROW;j++) |
---|
| 791 | { |
---|
| 792 | if(MA[i,j]!=0){k=1;} |
---|
| 793 | } |
---|
| 794 | if(k==0){adrow++;} |
---|
| 795 | } |
---|
| 796 | |
---|
| 797 | m=transpose(m); |
---|
| 798 | for(i=1;i<=adrow;i++){m=m,0;} |
---|
| 799 | m=transpose(m); |
---|
| 800 | |
---|
| 801 | |
---|
| 802 | int flag=1; ///////////////counts if the underlying ring is r (flag mod 2 == 1) or ro (flag mod 2 == 0) |
---|
| 803 | |
---|
| 804 | module TrafoL=freemodule(COL); |
---|
| 805 | module TrafoR=freemodule(ROW); |
---|
| 806 | module EXL=freemodule(COL); //because we start with transpose(m) |
---|
| 807 | module EXR=freemodule(ROW); |
---|
| 808 | |
---|
| 809 | option(redSB); |
---|
| 810 | option(redTail); |
---|
| 811 | module STD_EX,LT,TSTD, L, Trafo; |
---|
| 812 | |
---|
[7f3ad4] | 813 | |
---|
| 814 | |
---|
[2c0350] | 815 | module STD=transpose(m); |
---|
| 816 | int finish=0; |
---|
| 817 | list pos, COM, COM_EX; |
---|
| 818 | matrix END, ENDSTD, STDFIN; |
---|
| 819 | STDFIN=STD; |
---|
| 820 | list COMPARE=STDFIN; |
---|
| 821 | |
---|
| 822 | |
---|
| 823 | while(finish==0) |
---|
| 824 | { |
---|
| 825 | dbprint(ppl,"Going into the while cycle"); |
---|
| 826 | if(flag mod 2==1){STD_EX=EXL,transpose(STD); ex=2*COL;} else {STD_EX=EXR,transpose(STD); ex=2*ROW;} |
---|
| 827 | |
---|
| 828 | dbprint(ppl,"Computing Groebner basis: start"); |
---|
| 829 | dbprint(ppl-1,STD); |
---|
| 830 | STD=engine(STD,BASIS); |
---|
| 831 | dbprint(ppl,"Computing Groebner basis: finished"); |
---|
| 832 | dbprint(ppl-1,STD); |
---|
| 833 | if(flag mod 2==1){s=ROW; st=COL;}else{s=COL; st=ROW;} |
---|
[7f3ad4] | 834 | |
---|
[2c0350] | 835 | STD_EX=transpose(STD_EX); |
---|
| 836 | dbprint(ppl,"Computing Groebner basis for transformation matrix: start"); |
---|
| 837 | dbprint(ppl-1,STD_EX); |
---|
| 838 | STD_EX=engine(STD_EX,BASIS); |
---|
| 839 | dbprint(ppl,"Computing Groebner basis for transformation matrix: finished"); |
---|
| 840 | dbprint(ppl-1,STD_EX); |
---|
| 841 | |
---|
| 842 | COM=1; |
---|
| 843 | COM_EX=1; |
---|
| 844 | for(i=2; i<=size(STD); i++) |
---|
[7f3ad4] | 845 | { COM=COM[1..size(COM)],i; COM_EX=COM_EX[1..size(COM_EX)],i; } |
---|
| 846 | nz=size(STD_EX)-size(STD); |
---|
[2c0350] | 847 | |
---|
[7f3ad4] | 848 | //zero rows in the begining |
---|
[2c0350] | 849 | |
---|
| 850 | if(size(STD)!=size(STD_EX) ) |
---|
| 851 | { |
---|
| 852 | for(i=1; i<=size(STD_EX)-size(STD); i++) |
---|
| 853 | { |
---|
| 854 | COM_EX=0,COM_EX[1..size(COM_EX)]; |
---|
| 855 | } |
---|
[7f3ad4] | 856 | } |
---|
| 857 | |
---|
| 858 | |
---|
| 859 | |
---|
[2c0350] | 860 | |
---|
| 861 | for(i=nz+1; i<=size(STD_EX); i++ ) |
---|
[7f3ad4] | 862 | {if( leadcoef(STD[i-nz])!=leadcoef(STD_EX[i]) ) {STD[i-nz]=leadcoef(STD_EX[i])*STD[i-nz];} |
---|
[2c0350] | 863 | } |
---|
| 864 | |
---|
[7f3ad4] | 865 | //assign the zero rows in STD_EX |
---|
[2c0350] | 866 | |
---|
| 867 | for (j=2; j<=nz; j++) |
---|
| 868 | { |
---|
[7f3ad4] | 869 | p=0; |
---|
[2c0350] | 870 | i=1; |
---|
| 871 | while(STD_EX[j-1][i]==0){i++;}; |
---|
| 872 | p=i-1; |
---|
| 873 | nextp=1; |
---|
| 874 | k=0; |
---|
| 875 | i=1; |
---|
| 876 | while(STD_EX[j][i]==0 and i<=p) |
---|
| 877 | { k++; i++;} |
---|
| 878 | if (k==p){ COM_EX[j]=-1; } |
---|
| 879 | } |
---|
[7f3ad4] | 880 | |
---|
[2c0350] | 881 | //assign the zero rows in STD |
---|
| 882 | for (j=2; j<=size(STD); j++) |
---|
| 883 | { |
---|
[7f3ad4] | 884 | i=size(transpose(STD)); |
---|
[2c0350] | 885 | while(STD[j-1][i]==0){i--;} |
---|
| 886 | p=i; |
---|
| 887 | i=size(transpose(STD[j])); |
---|
| 888 | while(STD[j][i]==0){i--;} |
---|
| 889 | if (i==p){ COM[j]=-1; } |
---|
| 890 | } |
---|
| 891 | |
---|
| 892 | for(j=1; j<=size(COM); j++) |
---|
| 893 | { |
---|
| 894 | if(COM[j]<0){COM=delete(COM,j);} |
---|
| 895 | } |
---|
[7f3ad4] | 896 | |
---|
[2c0350] | 897 | for(i=1; i<=size(COM_EX); i++) |
---|
| 898 | {ff=0; |
---|
| 899 | if(COM_EX[i]==0){ff=1;} |
---|
| 900 | else |
---|
| 901 | { for(j=1; j<=size(COM); j++) |
---|
| 902 | { if(COM_EX[i]==COM[j]){ff=1;} } |
---|
| 903 | } |
---|
| 904 | if(ff==0){COM_EX[i]=-1;} |
---|
| 905 | } |
---|
| 906 | |
---|
| 907 | L=STD_EX[1]; |
---|
| 908 | for(i=2; i<=size(COM_EX); i++) |
---|
| 909 | { |
---|
| 910 | if(COM_EX[i]!=-1){L=L,STD_EX[i];} |
---|
| 911 | } |
---|
| 912 | |
---|
| 913 | //////// split STD_EX in STD and the transformation matrix |
---|
| 914 | |
---|
| 915 | L=transpose(L); |
---|
| 916 | STD=L[st+1]; |
---|
| 917 | LT=L[1]; |
---|
| 918 | |
---|
[7f3ad4] | 919 | |
---|
[2c0350] | 920 | for(i=2; i<=s+st; i++) |
---|
| 921 | { |
---|
| 922 | if (i<=st) |
---|
| 923 | { |
---|
| 924 | LT=LT,L[i]; |
---|
| 925 | } |
---|
| 926 | if (i>st+1) |
---|
| 927 | { |
---|
| 928 | STD=STD,L[i]; |
---|
| 929 | } |
---|
| 930 | } |
---|
| 931 | |
---|
| 932 | STD=transpose(STD); |
---|
| 933 | STDFIN=matrix(STD); |
---|
| 934 | COMPARE=insert(COMPARE,STDFIN); |
---|
| 935 | LT=transpose(LT); |
---|
| 936 | |
---|
| 937 | ////////////////////// compute the transformation matrices |
---|
[7f3ad4] | 938 | |
---|
[2c0350] | 939 | if (flag mod 2 ==1) |
---|
[7f3ad4] | 940 | { |
---|
[2c0350] | 941 | TrafoL=transpose(LT)*TrafoL; |
---|
| 942 | } |
---|
[7f3ad4] | 943 | else |
---|
| 944 | { |
---|
[2c0350] | 945 | TrafoR=TrafoR*involution(LT,inv); |
---|
| 946 | } |
---|
| 947 | |
---|
| 948 | |
---|
| 949 | ///////////////////////// check whether the alg termined ///////////////// |
---|
| 950 | if(size(COMPARE)>=3) |
---|
[7f3ad4] | 951 | { |
---|
[2c0350] | 952 | if(STD==COMPARE[3]){finish=1;} |
---|
[7f3ad4] | 953 | } |
---|
[2c0350] | 954 | ////////////////////////////////// change to the opposite module |
---|
| 955 | TSTD=transpose(STD); |
---|
| 956 | STD=involution(TSTD,inv); |
---|
| 957 | flag++; |
---|
| 958 | dbprint(ppl,"Finished one while cycle"); |
---|
| 959 | } |
---|
| 960 | |
---|
[7f3ad4] | 961 | if (flag mod 2 ==0){ STD = involution(STD,inv);} else { STD = transpose(STD); } |
---|
| 962 | |
---|
[2c0350] | 963 | list REVERSE=TrafoL,STD,TrafoR; |
---|
| 964 | return(REVERSE); |
---|
| 965 | } |
---|
| 966 | |
---|
| 967 | static proc divideByContent(matrix M) |
---|
| 968 | { |
---|
| 969 | //find first entrie not equal to zero |
---|
| 970 | int i,k; |
---|
| 971 | k=1; |
---|
| 972 | vector CON; |
---|
| 973 | for(i=1;i<=ncols(M);i++) |
---|
| 974 | { |
---|
| 975 | if(leadcoef(M[i])!=0){CON=CON+leadcoef(M[i])*gen(k); k++;} |
---|
| 976 | } |
---|
[7f3ad4] | 977 | poly con=content(CON); |
---|
[2c0350] | 978 | matrix TL=1/con*freemodule(nrows(M)); |
---|
| 979 | return(TL); |
---|
| 980 | } |
---|
| 981 | |
---|
| 982 | |
---|
| 983 | /////interesting examples for smith//////////////// |
---|
| 984 | |
---|
| 985 | /* |
---|
| 986 | |
---|
| 987 | //static proc Ex_One_wheeled_bicycle() |
---|
[7f3ad4] | 988 | { |
---|
[2c0350] | 989 | ring RA=(0,m), D, lp; |
---|
| 990 | matrix bicycle[2][3]=(1+m)*D^2, D^2, 1, D^2, D^2-1, 0; |
---|
| 991 | list s=smith(bicycle, 1,0); |
---|
| 992 | print(s[2]); |
---|
| 993 | print(s[1]*bicycle*s[3]-s[2]); |
---|
| 994 | } |
---|
| 995 | |
---|
| 996 | |
---|
[7f3ad4] | 997 | //static proc Ex_RLC-circuit() |
---|
[2c0350] | 998 | { |
---|
| 999 | ring RA=(0,m,R1,R2,L,C), D, lp; |
---|
| 1000 | matrix circuit[2][3]=D+1/(R1*C), 0, -1/(R1*C), 0, D+R2/L, -1/L; |
---|
| 1001 | list s=smith(circuit, 1,0); |
---|
| 1002 | print(s[2]); |
---|
| 1003 | print(s[1]*circuit*s[3]-s[2]); |
---|
| 1004 | } |
---|
| 1005 | |
---|
| 1006 | |
---|
| 1007 | //static proc Ex_two_pendula() |
---|
| 1008 | { |
---|
| 1009 | ring RA=(0,m,M,L1,L2,m1,m2,g), D, lp; |
---|
| 1010 | 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, |
---|
| 1011 | m2*L2^2*D^2-m2*L2*g,m2*L2*D^2,0; |
---|
| 1012 | list s=smith(pendula, 1,0); |
---|
| 1013 | print(s[2]); |
---|
| 1014 | print(s[1]*pendula*s[3]-s[2]); |
---|
| 1015 | } |
---|
| 1016 | |
---|
| 1017 | //static proc Ex_linerized_satellite_in_a_circular_equatorial_orbit() |
---|
| 1018 | { |
---|
| 1019 | ring RA=(0,m,omega,r,a,b), D, lp; |
---|
| 1020 | matrix satellite[4][6]= |
---|
| 1021 | D,-1,0,0,0,0, |
---|
| 1022 | -3*omega^2,D,0,-2*omega*r,-a/m,0, |
---|
| 1023 | 0,0,D,-1,0,0, |
---|
| 1024 | 0,2*omega/r,0,D,0,-b/(m*r); |
---|
| 1025 | list s=smith(satellite, 1,0); |
---|
| 1026 | print(s[2]); |
---|
| 1027 | print(s[1]*satellite*s[3]-s[2]); |
---|
| 1028 | } |
---|
| 1029 | |
---|
| 1030 | //static proc Ex_flexible_one_link_robot() |
---|
| 1031 | { |
---|
| 1032 | ring RA=(0,M11,M12,M13,M21,M22,M31,M33,K1,K2), D, lp; |
---|
| 1033 | 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; |
---|
| 1034 | list s=smith(robot, 1,0); |
---|
| 1035 | print(s[2]); |
---|
| 1036 | print(s[1]*robot*s[3]-s[2]); |
---|
| 1037 | } |
---|
| 1038 | |
---|
| 1039 | */ |
---|
| 1040 | |
---|
| 1041 | |
---|
| 1042 | /////interesting examples for jacobson//////////////// |
---|
| 1043 | |
---|
| 1044 | /* |
---|
| 1045 | |
---|
[7f3ad4] | 1046 | //static proc Ex_compare_output_with_maple_package_JanetOre() |
---|
| 1047 | { |
---|
[2c0350] | 1048 | ring w = 0,(x,d),Dp; |
---|
| 1049 | def W=nc_algebra(1,1); |
---|
| 1050 | setring W; |
---|
[7f3ad4] | 1051 | matrix m[3][3]=[d2,d+1,0],[d+1,0,d3-x2*d],[2d+1, d3+d2, d2]; |
---|
[2c0350] | 1052 | list J=jacobson(m,0); |
---|
[7f3ad4] | 1053 | print(J[1]*m*J[3]); |
---|
[2c0350] | 1054 | print(J[2]); |
---|
| 1055 | print(J[1]); |
---|
| 1056 | print(J[3]); |
---|
| 1057 | print(J[1]*m*J[3]-J[2]); |
---|
| 1058 | } |
---|
| 1059 | |
---|
| 1060 | // modification for shift algebra |
---|
| 1061 | { |
---|
| 1062 | ring w = 0,(x,s),Dp; |
---|
| 1063 | def W=nc_algebra(1,s); |
---|
| 1064 | setring W; |
---|
[7f3ad4] | 1065 | 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]; |
---|
[2c0350] | 1066 | list J=jacobson(m,0); |
---|
[7f3ad4] | 1067 | print(J[1]*m*J[3]); |
---|
[2c0350] | 1068 | print(J[2]); |
---|
| 1069 | print(J[1]); |
---|
| 1070 | print(J[3]); |
---|
| 1071 | print(J[1]*m*J[3]-J[2]); |
---|
| 1072 | } |
---|
| 1073 | |
---|
| 1074 | //static proc Ex_cyclic() |
---|
[7f3ad4] | 1075 | { |
---|
[2c0350] | 1076 | ring w = 0,(x,d),Dp; |
---|
| 1077 | def W=nc_algebra(1,1); |
---|
| 1078 | setring W; |
---|
| 1079 | matrix m[3][3]=d,0,0,x*d+1,d,0,0,x*d,d; |
---|
| 1080 | list J=jacobson(m,0); |
---|
[7f3ad4] | 1081 | print(J[1]*m*J[3]); |
---|
[2c0350] | 1082 | print(J[2]); |
---|
| 1083 | print(J[1]); |
---|
| 1084 | print(J[3]); |
---|
| 1085 | print(J[1]*m*J[3]-J[2]); |
---|
| 1086 | } |
---|
| 1087 | |
---|
| 1088 | // modification for shift algebra: gives a module with ann = s^2 |
---|
| 1089 | { |
---|
| 1090 | ring w = 0,(x,s),Dp; |
---|
| 1091 | def W=nc_algebra(1,s); |
---|
| 1092 | setring W; |
---|
| 1093 | matrix m[3][3]=s,0,0,x*s+1,s,0,0,x*s,s; |
---|
| 1094 | list J=jacobson(m,0); |
---|
[7f3ad4] | 1095 | print(J[1]*m*J[3]); |
---|
[2c0350] | 1096 | print(J[2]); |
---|
| 1097 | print(J[1]); |
---|
| 1098 | print(J[3]); |
---|
| 1099 | print(J[1]*m*J[3]-J[2]); |
---|
| 1100 | } |
---|
| 1101 | |
---|
| 1102 | // yet another modification for shift algebra: gives a module with ann = s |
---|
| 1103 | // xs+1 is replaced by x*s+s |
---|
| 1104 | { |
---|
| 1105 | ring w = 0,(x,s),Dp; |
---|
| 1106 | def W=nc_algebra(1,s); |
---|
| 1107 | setring W; |
---|
| 1108 | matrix m[3][3]=s,0,0,x*s+s,s,0,0,x*s,s; |
---|
| 1109 | list J=jacobson(m,0); |
---|
[7f3ad4] | 1110 | print(J[1]*m*J[3]); |
---|
[2c0350] | 1111 | print(J[2]); |
---|
| 1112 | print(J[1]); |
---|
| 1113 | print(J[3]); |
---|
| 1114 | print(J[1]*m*J[3]-J[2]); |
---|
| 1115 | } |
---|
| 1116 | |
---|
| 1117 | // variations for above setup with shift algebra: |
---|
| 1118 | |
---|
| 1119 | // easy but instructive one: see the difference to Weyl case! |
---|
| 1120 | matrix m[2][2]=s,x,0,s; print(m); |
---|
| 1121 | |
---|
| 1122 | // more interesting: |
---|
| 1123 | matrix n[3][3]=s,x,0,0,s,x,s,0,x; |
---|
| 1124 | |
---|
| 1125 | // things blow up |
---|
| 1126 | 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; |
---|
| 1127 | |
---|
| 1128 | // this one is quite nasty and time consuming |
---|
| 1129 | 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; |
---|
| 1130 | |
---|
| 1131 | */ |
---|
[fec31e1] | 1132 | |
---|