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