1 | ////////////////////////////////////////////////////////////////////////////// |
---|
2 | version="$Id$"; |
---|
3 | category="System and Control Theory"; |
---|
4 | info=" |
---|
5 | LIBRARY: jacobson.lib Algorithms for Smith and Jacobson Normal Form |
---|
6 | AUTHOR: Kristina Schindelar, Kristina.Schindelar@math.rwth-aachen.de, |
---|
7 | @* Viktor Levandovskyy, levandov@math.rwth-aachen.de |
---|
8 | |
---|
9 | OVERVIEW: |
---|
10 | We work over a ring R, that is an Euclidean principal ideal domain. |
---|
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. |
---|
29 | |
---|
30 | |
---|
31 | REFERENCES: |
---|
32 | @* [1] N. Jacobson, 'The theory of rings', AMS, 1943. |
---|
33 | @* [2] Manuel Avelino Insua Hermo, 'Varias perspectives sobre las bases de Groebner : |
---|
34 | @* Forma normal de Smith, Algorithme de Berlekamp y algebras de Leibniz'. |
---|
35 | @* PhD thesis, Universidad de Santiago de Compostela, 2005. |
---|
36 | @* [3] V. Levandovskyy, K. Schindelar 'Computing Jacobson normal form using Groebner bases', |
---|
37 | @* to appear in Journal of Symbolic Computation, 2010. |
---|
38 | |
---|
39 | PROCEDURES: |
---|
40 | smith(M[,eng1,eng2]); compute the Smith Normal Form of M over commutative ring |
---|
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 | |
---|
44 | SEE ALSO: control_lib |
---|
45 | |
---|
46 | KEYWORDS: Jacobson form; Jacobson normal form; Smith form; Smith normal form; matrix diagonalization |
---|
47 | |
---|
48 | "; |
---|
49 | |
---|
50 | LIB "poly.lib"; |
---|
51 | LIB "involut.lib"; // involution |
---|
52 | LIB "dmodapp.lib"; // engine |
---|
53 | LIB "qhmoduli.lib"; // Min |
---|
54 | LIB "latex.lib"; // tex output for printlevel=4 |
---|
55 | |
---|
56 | proc tstjacobson() |
---|
57 | { |
---|
58 | /* tests all procs for consistency */ |
---|
59 | example divideUnits; |
---|
60 | example smith; |
---|
61 | example jacobson; |
---|
62 | } |
---|
63 | |
---|
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 | } |
---|
88 | |
---|
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; |
---|
97 | |
---|
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 | } |
---|
141 | |
---|
142 | proc smith(matrix MA, list #) |
---|
143 | "USAGE: smith(M[, eng1, eng2]); M matrix, eng1 and eng2 are optional integers |
---|
144 | RETURN: matrix or list of matrices, depending on arguments |
---|
145 | ASSUME: Basering is a commutative polynomial ring in one variable |
---|
146 | PURPOSE: compute the Smith Normal Form of M with (optionally) transformation matrices |
---|
147 | THEORY: Groebner bases are used for the Smith form like in [2] and [3]. |
---|
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. |
---|
152 | @* U and V above are square unimodular (invertible) matrices. |
---|
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. |
---|
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 |
---|
160 | SEE ALSO: divideUnits, jacobson |
---|
161 | " |
---|
162 | { |
---|
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]; |
---|
184 | } |
---|
185 | else{BASIS=0;} |
---|
186 | |
---|
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 | { |
---|
214 | list rueckLI=diagonal_with_trafo(R,MA,BASIS); |
---|
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 | { |
---|
222 | matrix rueckm=diagonal_without_trafo(R,MA,BASIS); |
---|
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); |
---|
234 | print(s[2]); // non-normalized Smith form of m |
---|
235 | print(s[1]*m*s[3] - s[2]); // check U*M*V = D |
---|
236 | list t = divideUnits(s); |
---|
237 | print(t[2]); // the Smith form of m |
---|
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 | } |
---|
268 | |
---|
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; |
---|
297 | |
---|
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; |
---|
305 | matrix STDFIN; |
---|
306 | STDFIN=STD; |
---|
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); |
---|
316 | } |
---|
317 | else |
---|
318 | { |
---|
319 | STD_EX=EXR,transpose(STD); |
---|
320 | } |
---|
321 | dbprint(ppl,"Computing Groebner basis: start"); |
---|
322 | dbprint(ppl-1,STD); |
---|
323 | STD=engine(STD,BASIS); |
---|
324 | dbprint(ppl,"Computing Groebner basis: finished"); |
---|
325 | dbprint(ppl-1,STD); |
---|
326 | |
---|
327 | if(flag mod 2==1){s=ROW; st=COL;}else{s=COL; st=ROW;} |
---|
328 | dbprint(ppl,"Computing Groebner basis for transformation matrix: start"); |
---|
329 | dbprint(ppl-1,STD_EX); |
---|
330 | STD_EX=transpose(STD_EX); |
---|
331 | STD_EX=engine(STD_EX,BASIS); |
---|
332 | dbprint(ppl,"Computing Groebner basis for transformation matrix: finished"); |
---|
333 | dbprint(ppl-1,STD_EX); |
---|
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 |
---|
357 | |
---|
358 | if (flag mod 2 ==1) |
---|
359 | { |
---|
360 | TrafoL=transpose(LT)*TrafoL; |
---|
361 | } |
---|
362 | else |
---|
363 | { |
---|
364 | TrafoR=TrafoR*LT; |
---|
365 | } |
---|
366 | |
---|
367 | |
---|
368 | STDFIN=STD; |
---|
369 | /////////////////////////////////// check if the D-th row is finished /////////////////////////////////// |
---|
370 | COMPARE=insert(COMPARE,STDFIN); |
---|
371 | if(size(COMPARE)>=3) |
---|
372 | { |
---|
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); } |
---|
384 | |
---|
385 | |
---|
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 | { |
---|
394 | ff=0; |
---|
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; |
---|
407 | for(i=1; i<=fehlpos; i++) |
---|
408 | { |
---|
409 | SORT=gen(2); |
---|
410 | for (j=3;j<=ROW;j++) |
---|
411 | { |
---|
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) |
---|
424 | { |
---|
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 | |
---|
435 | for(i=1; i<=fehlpos; i++) |
---|
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 | } |
---|
444 | |
---|
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 | |
---|
454 | //Test |
---|
455 | if(TrafoLM*m*TrafoRM!=STDM){ return(0); } |
---|
456 | |
---|
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 | { |
---|
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; |
---|
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 | { |
---|
534 | int ppl = printlevel-voice+2; |
---|
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; |
---|
568 | STDFIN=STD; |
---|
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); |
---|
579 | STDFIN=STD; |
---|
580 | /////////////////////////////////// check if the D-th row is finished /////////////////////////////////// |
---|
581 | COMPARE=insert(COMPARE,STDFIN); |
---|
582 | if(size(COMPARE)>=3) |
---|
583 | { |
---|
584 | if(STD==COMPARE[3]){finish=1;} |
---|
585 | } |
---|
586 | ////////////////////////////////// change to the opposite module |
---|
587 | |
---|
588 | TSTD=transpose(STD); |
---|
589 | STD=TSTD; |
---|
590 | flag++; |
---|
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 | { |
---|
601 | ff=0; |
---|
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; |
---|
614 | for(i=1; i<=fehlpos; i++) |
---|
615 | { |
---|
616 | SORT=gen(2); |
---|
617 | for (j=3;j<=ROW;j++) |
---|
618 | { |
---|
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) |
---|
630 | { |
---|
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 | |
---|
641 | for(i=1; i<=fehlpos; i++) |
---|
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 | } |
---|
649 | |
---|
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 | |
---|
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 | // } |
---|
688 | |
---|
689 | proc jacobson(matrix MA, list #) |
---|
690 | "USAGE: jacobson(M, eng); M matrix, eng an optional int |
---|
691 | RETURN: list |
---|
692 | ASSUME: Basering is a (non-commutative) ring in two variables. |
---|
693 | PURPOSE: compute a weak Jacobson normal form of M over the basering |
---|
694 | THEORY: Groebner bases and involutions are used, following [3] |
---|
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 |
---|
697 | @* L[1], L[3] are square invertible polynomial (unimodular) matrices. |
---|
698 | @* Note, that M can be rectangular. |
---|
699 | @* The optional integer @code{eng2} determines the Groebner basis engine: |
---|
700 | @* 0 (default) ensures the use of 'slimgb' , otherwise 'std' is used. |
---|
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 |
---|
704 | SEE ALSO: divideUnits, smith |
---|
705 | " |
---|
706 | { |
---|
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; |
---|
739 | intvec v=0,1; |
---|
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; |
---|
747 | |
---|
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 | |
---|
759 | //back to the ring |
---|
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); |
---|
766 | matrix TRR=TR; |
---|
767 | matrix CON=divideByContent(MAA); |
---|
768 | |
---|
769 | list RUECK=CON*TL, CON*MAA, TRR; |
---|
770 | return(RUECK); |
---|
771 | } |
---|
772 | example |
---|
773 | { |
---|
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 |
---|
785 | list J = jacobson(m); |
---|
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 | |
---|
790 | //static |
---|
791 | proc triangle(matrix MA, int B) |
---|
792 | { |
---|
793 | int ppl = printlevel-voice+2; |
---|
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 | |
---|
836 | |
---|
837 | |
---|
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; |
---|
844 | string @s; |
---|
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;} |
---|
857 | |
---|
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++) |
---|
868 | { COM=COM[1..size(COM)],i; COM_EX=COM_EX[1..size(COM_EX)],i; } |
---|
869 | nz=size(STD_EX)-size(STD); |
---|
870 | |
---|
871 | //zero rows in the begining |
---|
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 | } |
---|
879 | } |
---|
880 | |
---|
881 | for(i=nz+1; i<=size(STD_EX); i++ ) |
---|
882 | {if( leadcoef(STD[i-nz])!=leadcoef(STD_EX[i]) ) {STD[i-nz]=leadcoef(STD_EX[i])*STD[i-nz];} |
---|
883 | } |
---|
884 | |
---|
885 | //assign the zero rows in STD_EX |
---|
886 | |
---|
887 | for (j=2; j<=nz; j++) |
---|
888 | { |
---|
889 | p=0; |
---|
890 | i=1; |
---|
891 | while(STD_EX[j-1][i]==0){i++;} |
---|
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 | } |
---|
900 | |
---|
901 | //assign the zero rows in STD |
---|
902 | for (j=2; j<=size(STD); j++) |
---|
903 | { |
---|
904 | i=size(transpose(STD)); |
---|
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 | } |
---|
916 | |
---|
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 | |
---|
939 | |
---|
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 |
---|
958 | |
---|
959 | if (flag mod 2 ==1) |
---|
960 | { |
---|
961 | TrafoL=transpose(LT)*TrafoL; |
---|
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 | } |
---|
968 | else |
---|
969 | { |
---|
970 | TrafoR=TrafoR*involution(LT,inv); |
---|
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); |
---|
976 | } |
---|
977 | |
---|
978 | ///////////////////////// check whether the alg terminated ///////////////// |
---|
979 | if(size(COMPARE)>=3) |
---|
980 | { |
---|
981 | if(STD==COMPARE[3]){finish=1;} |
---|
982 | } |
---|
983 | ////////////////////////////////// change to the opposite module |
---|
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"); |
---|
993 | } |
---|
994 | |
---|
995 | if (flag mod 2 ==0){ STD = involution(STD,inv);} else { STD = transpose(STD); } |
---|
996 | |
---|
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 | } |
---|
1011 | poly con=content(CON); |
---|
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() |
---|
1022 | { |
---|
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 | |
---|
1031 | //static proc Ex_RLC-circuit() |
---|
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 | |
---|
1080 | //static proc Ex_compare_output_with_maple_package_JanetOre() |
---|
1081 | { |
---|
1082 | ring w = 0,(x,d),Dp; |
---|
1083 | def W=nc_algebra(1,1); |
---|
1084 | setring W; |
---|
1085 | matrix m[3][3]=[d2,d+1,0],[d+1,0,d3-x2*d],[2d+1, d3+d2, d2]; |
---|
1086 | list J=jacobson(m,0); |
---|
1087 | print(J[1]*m*J[3]); |
---|
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; |
---|
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]; |
---|
1100 | list J=jacobson(m,0); |
---|
1101 | print(J[1]*m*J[3]); |
---|
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() |
---|
1109 | { |
---|
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); |
---|
1115 | print(J[1]*m*J[3]); |
---|
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); |
---|
1129 | print(J[1]*m*J[3]); |
---|
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); |
---|
1144 | print(J[1]*m*J[3]); |
---|
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 |
---|
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 |
---|
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 |
---|
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 |
---|
1164 | |
---|
1165 | // this one is quite nasty and time consumig |
---|
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 | |
---|
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]); |
---|
1176 | |
---|
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]); |
---|
1203 | print(J[1]*m*J[3]-J[2]); |
---|
1204 | |
---|
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 | |
---|
1215 | */ |
---|