source:git/Singular/LIB/jacobson.lib@a470eb

spielwiese
Last change on this file since a470eb was a470eb, checked in by Viktor Levandovskyy <levandov@…>, 15 years ago
*levandov: minor bugfixes and docu git-svn-id: file:///usr/local/Singular/svn/trunk@11229 2c84dea3-7e68-4137-9b89-c4e89433aadc
• Property mode set to `100755`
File size: 21.9 KB
Line
1//////////////////////////////////////////////////////////////////////////////
2version="\$Id: jacobson.lib,v 1.5 2008-12-09 16:50:21 levandov Exp \$";
3category="System and Control Theory";
4info="
5LIBRARY: jacobson.lib     Algorithms for Smith and Jacobson Normal Form
6AUTHOR: Kristina Schindelar,     Kristina.Schindelar@math.rwth-aachen.de
7
8THEORY: We work over a ring R, that is a principal ideal domain.
9@*  If R is commutative, we suppose R to be a polynomial ring in one variable.
10@*  If R is non-commutative, we suppose R to be in two variables, say x and d. We treat
11@* then the Ore localization of R with respect to the mult.closed set S = K[x] without 0.
12@*  Given a rectangular matrix M over R, one can compute unimodular matrices U, V
13@*  such that U*M*V=D is a diagonal matrix.
14@*  Depending on the ring, the diagonal entries of D have certain properties.
15
16PROCEDURES:
17smith(R,M[,eng1,eng2]);  compute the Smith Normal Form of M over commutative ring
18jacobson(R,M[,eng]);     compute a weak Jacobson Normal Form of M over non-commutative ring, i.e. a diagonal matrix
19
21";
22  LIB "poly.lib";
23  LIB "involut.lib";
24
25proc smith(R, matrix MA, list #)
26"USAGE:  smith(R, M[, eng1, eng2]);  R ring, M matrix, eng1 and eng2 are optional int
27RETURN: list
28NOTE: By default, the diagonal matrix, the Smith normal form, is returned.
29@* R stays for the ring, where computations will take place.
30@* If optional integer eng1 is nonzero, smith returns
31@* a list of matrices L such that L[1]*M*L[3]=L[2] with L[2]_11|L[2]_22|...|L[2]_nn and
32@* L[1], L[3] unimodular matrices.
33@* eng2 determines the engine, that computes the Groebner basis. By default eng2 equals zero.
34@* If optional integer eng2 = 0 than std is used to caculate a Groebner basis ,
35@* If eng2 = 1 than groebner is used to caculate a Groebner basis
36@* If eng2 = 2 than slimgb is used to caculate a Groebner basis
37@* If @code{printlevel}=1, progress debug messages will be printed,
38@* if @code{printlevel}>=2, all the debug messages will be printed.
39EXAMPLE: example smith; shows examples
40"
41{
42  int eng = 0;
43  int BASIS;
44  if ( size(#)>0 )
45  {
46    eng=1;
47    if (typeof(#[1])=="int")
48    {
49      eng=int(#[1]); // zero can also happen
50    }
51  }
52  if (size(#)==2)
53  {
54  BASIS=#[2];
55  }
56  else{BASIS=0;}
57
58
59  int ROW=ncols(MA);
60  int COL=nrows(MA);
61
62  //generate a module consisting of the columns of MA
63  module m=MA[1];
64  int i;
65  for(i=2;i<=ROW;i++)
66  {
67    m=m,MA[i];
68  }
69
70  //if MA eqauls the zero matrix give back MA
71  if ( (size(module(m))==0) and (size(transpose(module(m)))==0) )
72  {
73    module L=freemodule(COL);
74    matrix LM=L;
75    L=freemodule(ROW);
76    matrix RM=L;
77    list RUECK=RM;
78    RUECK=insert(RUECK,MA);
79    RUECK=insert(RUECK,LM);
80    return(RUECK);
81  }
82
83  if(eng==1)
84  {
85    list rueckLI=diagonal_with_trafo(R,MA,BASIS);
86    list rueckLII=divisibility(rueckLI[2]);
87    matrix CON=divideByContent(rueckLII[2]);
88    list rueckL=CON*rueckLII[1]*rueckLI[1], CON*rueckLII[2], rueckLI[3]*rueckLII[3];
89    return(rueckL);
90  }
91  else
92  {
93    matrix rueckm=diagonal_without_trafo(R,MA,BASIS);
94    list rueckL=divisibility(rueckm);
95    matrix CON=divideByContent(rueckm);
96    rueckm=CON*rueckL[2];
97    return(rueckm);
98  }
99}
100example
101{ "EXAMPLE:"; echo = 2;
102  ring r = 0,x,Dp;
103  matrix m[3][2]=x, x^4+x^2+21, x^4+x^2+x, x^3+x, 4*x^2+x, x;
104  list s=smith(r,m,1);
105  print(s[2]);
106  print(s[1]*m*s[3]);
107}
108
109static proc diagonal_with_trafo( R, matrix MA, int B)
110{
111
112int ppl = printlevel-voice+2;
113
114  int BASIS=B;
115  int ROW=ncols(MA);
116  int COL=nrows(MA);
117  module m=MA[1];
118  int i,j,k;
119  for(i=2;i<=ROW;i++)
120  {
121    m=m,MA[i];
122  }
123
124
125  //add zero rows or columns
126  //add zero rows or columns
128  for(i=1;i<=COL;i++)
129  {
130  k=0;
131  for(j=1;j<=ROW;j++)
132  {
133  if(MA[i,j]!=0){k=1;}
134  }
136  }
137
138    m=transpose(m);
140    m=transpose(m);
141
142  list RINGLIST=ringlist(R);
143  list o="C",0;
144  list oo="lp",1;
145  list ORD=o,oo;
146  RINGLIST[3]=ORD;
147  def r=ring(RINGLIST);
148  setring r;
149  //fix the required ordering
150  map MAP=R,var(1);
151  module m=MAP(m);
152
153  int flag=1;  ///////////////counts if the underlying ring is r (flag mod 2 == 1) or ro (flag mod 2 == 0)
154
155    module TrafoL=freemodule(COL);
156    module TrafoR=freemodule(ROW);
158    module EXR=freemodule(ROW);
159
160    option(redSB);
161    option(redTail);
162    module STD_EX;
163    module Trafo;
164
165    int s,st,p,ff;
166
167    module LT,TSTD;
168    module STD=transpose(m);
169    int finish=0;
170    int fehlpos;
171    list pos;
172    matrix END;
173    matrix ENDSTD;
174    matrix STDFIN;
175    STDFIN=STD;
176    list COMPARE=STDFIN;
177
178    while(finish==0)
179    {
180      dbprint(ppl,"Going into the while cycle");
181
182      if(flag mod 2==1)
183      {
184        STD_EX=EXL,transpose(STD);
185      }
186      else
187      {
188        STD_EX=EXR,transpose(STD);
189      }
190        dbprint(ppl,"Computing Groebner basis: start");
191        dbprint(ppl-1,STD);
192      STD=engine(STD,BASIS);
193        dbprint(ppl,"Computing Groebner basis: finished");
194        dbprint(ppl-1,STD);
195
196      if(flag mod 2==1){s=ROW; st=COL;}else{s=COL; st=ROW;}
197        dbprint(ppl,"Computing Groebner basis for transformation matrix: start");
198        dbprint(ppl-1,STD_EX);
199      STD_EX=transpose(STD_EX);
200      STD_EX=engine(STD_EX,BASIS);
201        dbprint(ppl,"Computing Groebner basis for transformation matrix: finished");
202        dbprint(ppl-1,STD_EX);
203
204      //////// split STD_EX in STD and the transformation matrix
205      STD_EX=transpose(STD_EX);
206      STD=STD_EX[st+1];
207      LT=STD_EX[1];
208
209      ENDSTD=STD_EX;
210      for(i=2; i<=ncols(ENDSTD); i++)
211      {
212        if (i<=st)
213        {
214         LT=LT,STD_EX[i];
215        }
216        if (i>st+1)
217        {
218          STD=STD,STD_EX[i];
219        }
220      }
221
222      STD=transpose(STD);
223      LT=transpose(LT);
224
225      ////////////////////// compute the transformation matrices
226
227      if (flag mod 2 ==1)
228      {
229        TrafoL=transpose(LT)*TrafoL;
230      }
231      else
232      {
233        TrafoR=TrafoR*LT;
234      }
235
236
237      STDFIN=STD;
238      /////////////////////////////////// check if the D-th row is finished ///////////////////////////////////
239      COMPARE=insert(COMPARE,STDFIN);
240        if(size(COMPARE)>=3)
241       {
242        if(STD==COMPARE[3]){finish=1;}
243       }
244////////////////////////////////// change to the opposite module
245       TSTD=transpose(STD);
246       STD=TSTD;
247       flag++;
248      dbprint(ppl,"Finished one while cycle");
249    }
250
251
252   if (flag mod 2!=0) { STD=transpose(STD); }
253
254
255   //zero colums to the end
256    matrix STDMM=STD;
257    pos=list();
258    fehlpos=0;
259    while( size(STD)+fehlpos-ncols(STDMM) < 0)
260    {
261      for(i=1; i<=ncols(STDMM); i++)
262      {
263        ff=0;
264        for(j=1; j<=nrows(STDMM); j++)
265        {
266          if (STD[j,i]==0) { ff++; }
267        }
268        if(ff==nrows(STDMM))
269        {
270          pos=insert(pos,i); fehlpos++;
271        }
272      }
273    }
274    int fehlposc=fehlpos;
275    module SORT;
276    for(i=1; i<=fehlpos; i++)
277    {
278      SORT=gen(2);
279      for (j=3;j<=ROW;j++)
280      {
281        SORT=SORT,gen(j);
282      }
283      SORT=SORT,gen(1);
284      STD=STD*SORT;
285      TrafoR=TrafoR*SORT;
286    }
287
288    //zero rows to the end
289    STDMM=transpose(STD);
290    pos=list();
291    fehlpos=0;
292    while( size(transpose(STD))+fehlpos-ncols(STDMM) < 0)
293    {
294      for(i=1; i<=ncols(STDMM); i++)
295      {
296        ff=0; for(j=1; j<=nrows(STDMM); j++)
297        {
298           if(transpose(STD)[j,i]==0){ ff++;}}
299           if(ff==nrows(STDMM)) { pos=insert(pos,i); fehlpos++; }
300      }
301    }
302    int fehlposr=fehlpos;
303
304    for(i=1; i<=fehlpos; i++)
305    {
306      SORT=gen(2);
307      for(j=3;j<=COL;j++){SORT=SORT,gen(j);}
308      SORT=SORT,gen(1);
309      SORT=transpose(SORT);
310      STD=SORT*STD;
311      TrafoL=SORT*TrafoL;
312    }
313
314    setring R;
315    map MAPinv=r,var(1);
316    module STD=MAPinv(STD);
317    module TrafoL=MAPinv(TrafoL);
318    matrix TrafoLM=TrafoL;
319    module TrafoR=MAPinv(TrafoR);
320    matrix TrafoRM=TrafoR;
321    matrix STDM=STD;
322
323    //Test
324    if(TrafoLM*m*TrafoRM!=STDM){ return(0); }
325
326    list RUECK=TrafoRM;
327    RUECK=insert(RUECK,STDM);
328    RUECK=insert(RUECK,TrafoLM);
329    return(RUECK);
330}
331
332
333static proc divisibility(matrix M)
334   {
335    matrix STDM=M;
336    int i,j;
337    int ROW=nrows(M);
338    int COL=ncols(M);
339    module TrafoR=freemodule(COL);
340    module TrafoL=freemodule(ROW);
341    module SORT;
342    matrix TrafoRM=TrafoR;
343    matrix TrafoLM=TrafoL;
344    list posdeg0;
345    int posdeg=0;
346    int act;
347    int Sort=ROW;
348    if(size(std(STDM))!=0)
349    { while( size(transpose(STDM)[Sort])==0 ){Sort--;}}
350
351    for(i=1;i<=Sort ;i++)
352    {
354    }
355    //entries of degree 0 at the beginning
356    for(i=1; i<=posdeg; i++)
357    {
358      act=posdeg0[i];
359      SORT=gen(act);
360      for(j=1; j<=COL; j++){if(j!=act){SORT=SORT,gen(j);}}
361      STDM=STDM*SORT;
362      TrafoRM=TrafoRM*SORT;
363      SORT=gen(act);
364      for(j=1; j<=ROW; j++){if(j!=act){SORT=SORT,gen(j);}}
365      STDM=transpose(SORT)*STDM;
366      TrafoLM=transpose(SORT)*TrafoLM;
367    }
368
369    poly G;
370    module UNITL=freemodule(ROW);
371    matrix GCDL=UNITL;
372    module UNITR=freemodule(COL);
373    matrix GCDR=UNITR;
374    for(i=posdeg+1; i<=Sort; i++)
375    {
376      for(j=i+1; j<=Sort; j++)
377      {
378        GCDL=UNITL;
379        GCDR=UNITR;
380        G=gcd(STDM[i,i],STDM[j,j]);
381        ideal Z=STDM[i,i],STDM[j,j];
382        matrix T=lift(Z,G);
383        GCDL[i,i]=T[1,1];
384        GCDL[i,j]=T[2,1];
385        GCDL[j,i]=-STDM[j,j]/G;
386        GCDL[j,j]=STDM[i,i]/G;
387        GCDR[i,j]=T[2,1]*STDM[j,j]/G;
388        GCDR[j,j]=T[2,1]*STDM[j,j]/G-1;
389        GCDR[j,i]=1;
390        STDM=GCDL*STDM*GCDR;
391        TrafoLM=GCDL*TrafoLM;
392        TrafoRM=TrafoRM*GCDR;
393      }
394    }
395    list RUECK=TrafoRM;
396    RUECK=insert(RUECK,STDM);
397    RUECK=insert(RUECK,TrafoLM);
398    return(RUECK);
399}
400
401static proc diagonal_without_trafo( R, matrix MA, int B)
402{
403  int ppl = printlevel-voice+2;
404
405  int BASIS=B;
406  int ROW=ncols(MA);
407  int COL=nrows(MA);
408  module m=MA[1];
409  int i;
410  for(i=2;i<=ROW;i++)
411  {m=m,MA[i];}
412
413
414  list RINGLIST=ringlist(R);
415  list o="C",0;
416  list oo="lp",1;
417  list ORD=o,oo;
418  RINGLIST[3]=ORD;
419  def r=ring(RINGLIST);
420  setring r;
421  //RICHTIGE ORDNUNG MACHEN
422  map MAP=R,var(1);
423  module m=MAP(m);
424
425  int flag=1; ///////////////counts if the underlying ring is r (flag mod 2 == 1) or ro (flag mod 2 == 0)
426
427
428    int act, j, ff;
429    option(redSB);
430    option(redTail);
431
432
433    module STD=transpose(m);
434    module TSTD;
435    int finish=0;
436    matrix STDFIN;
437    STDFIN=STD;
438    list COMPARE=STDFIN;
439
440    while(finish==0)
441    {
442      dbprint(ppl,"Going into the while cycle");
443      dbprint(ppl,"Computing Groebner basis: start");
444      dbprint(ppl-1,STD);
445      STD=engine(STD,BASIS);
446      dbprint(ppl,"Computing Groebner basis: finished");
447      dbprint(ppl-1,STD);
448      STDFIN=STD;
449      /////////////////////////////////// check if the D-th row is finished ///////////////////////////////////
450      COMPARE=insert(COMPARE,STDFIN);
451        if(size(COMPARE)>=3)
452       {
453        if(STD==COMPARE[3]){finish=1;}
454       }
455        ////////////////////////////////// change to the opposite module
456
457        TSTD=transpose(STD);
458        STD=TSTD;
459        flag++;
460        dbprint(ppl,"Finished one while cycle");
461    }
462
463    matrix STDMM=STD;
464    list pos=list();
465    int fehlpos=0;
466    while( size(STD)+fehlpos-ncols(STDMM) < 0)
467    {
468      for(i=1; i<=ncols(STDMM); i++)
469      {
470        ff=0;
471        for(j=1; j<=nrows(STDMM); j++)
472        {
473          if (STD[j,i]==0) { ff++; }
474        }
475        if(ff==nrows(STDMM))
476        {
477          pos=insert(pos,i); fehlpos++;
478        }
479      }
480    }
481   int fehlposc=fehlpos;
482    module SORT;
483    for(i=1; i<=fehlpos; i++)
484    {
485      SORT=gen(2);
486      for (j=3;j<=ROW;j++)
487      {
488        SORT=SORT,gen(j);
489      }
490      SORT=SORT,gen(1);
491      STD=STD*SORT;
492    }
493
494    //zero rows to the end
495    STDMM=transpose(STD);
496    pos=list();
497    fehlpos=0;
498    while( size(transpose(STD))+fehlpos-ncols(STDMM) < 0)
499    {
500      for(i=1; i<=ncols(STDMM); i++)
501      {
502        ff=0; for(j=1; j<=nrows(STDMM); j++)
503        {
504           if(transpose(STD)[j,i]==0){ ff++;}}
505           if(ff==nrows(STDMM)) { pos=insert(pos,i); fehlpos++; }
506      }
507    }
508    int fehlposr=fehlpos;
509
510    for(i=1; i<=fehlpos; i++)
511    {
512      SORT=gen(2);
513      for(j=3;j<=COL;j++){SORT=SORT,gen(j);}
514      SORT=SORT,gen(1);
515      SORT=transpose(SORT);
516      STD=SORT*STD;
517    }
518
519   //add zero rows or columns
520
523
525    STD=transpose(STD);
527    STD=transpose(STD);
528
529    setring R;
530    map MAPinv=r,var(1);
531    module STD=MAPinv(STD);
532    matrix STDM=STD;
533    return(STDM);
534}
535
536
537
538static proc engine(module I, int i)
539{
540  module J;
541  if (i==0)
542  {
543    J = std(I);
544  }
545  if (i==1)
546  {
547    J = groebner(I);
548  }
549  if (i==2)
550  {
551    J = slimgb(I);
552  }
553  return(J);
554}
555
556proc jacobson(R, matrix MA, list #)
557 "USAGE:  jacobson(R, matrix MA, eng);  R ring, M matrix, eng an optional int
558RETURN: list
559NOTE: A list of matrices L such that L[1]*M*L[3]=L[2] such that L[2] is a diagonal matrix and
560@* L[1], L[3] unimodular matrices.
561@* R stays for the ring, where computations will take place.
562@* eng determines the engine, that computes the Groebner basis. By default eng equals zero.
563@* If eng = 0 than std is used to caculate a Groebner basis
564@* If eng = 1 than groebner is used to caculate a Groebner basis
565@* If eng = 2 than slimgb is used to caculate a Groebner basis
566@* If @code{printlevel}=1, progress debug messages will be printed,
567@* if @code{printlevel}>=2, all the debug messages will be printed.
568EXAMPLE: example jacobson; shows examples
569"
570{
571  int B=0;
572  if ( size(#)>0 )
573  {
574    B=1;
575    if (typeof(#[1])=="int")
576    {
577    B=int(#[1]); // zero can also happen
578    }
579  }
580
581  //change ring
582  list RINGLIST=ringlist(R);
583  list o="C",0;
584  intvec v=1,0;
585  list oo="a",v;
586  v=1,1;
587  list ooo="lp",v;
588  list ORD=o,oo,ooo;
589  RINGLIST[3]=ORD;
590  def r=ring(RINGLIST);
591  setring r;
592
593  //fix the required ordering
594  map MAP=R,var(1),var(2);
595  matrix M=MAP(MA);
596
597  module TrafoL, TrafoR;
598  list TRIANGLE;
599  TRIANGLE=triangle(M,B);
600  TrafoL=TRIANGLE[1];
601  TrafoR=TRIANGLE[3];
602  module m=TRIANGLE[2];
603
604  //back to the ring
605  setring R;
606  map MAPR=r,var(1),var(2);
607  module ma=MAPR(m);
608  matrix MAA=ma;
609  module TL=MAPR(TrafoL);
610  module TR=MAPR(TrafoR);
611  matrix CON=divideByContent(MAA);
612
613list RUECK=CON*TL, CON*MAA, TR;
614return(RUECK);
615}
616example
617{ "EXAMPLE:"; echo = 2;
618  ring r = 0,(x,d),Dp;
619  def R=nc_algebra(1,1); // the Weyl algebra
620  setring R;
621  matrix m[2][2]=d,x,0,d; print(m);
622  list J=jacobson(R,m); // returns a list with 3 entries
623  print(J[2]); // a Jacobson Form D
624  print(J[1]*m*J[3]); // check that U*M*V = D
625}
626
627
628
629static proc triangle( matrix MA, int B)
630{
631 int ppl = printlevel-voice+2;
632
633  map inv=ncdetection();
634  int ROW=ncols(MA);
635  int COL=nrows(MA);
636
637  //generate a module consisting of the columns of MA
638  module m=MA[1];
639  int i,j,s,st,p,k,ff,ex, nz, ii,nextp;
640  for(i=2;i<=ROW;i++)
641  {
642    m=m,MA[i];
643  }
644  int BASIS=B;
645
646  //add zero rows or columns
648  for(i=1;i<=COL;i++)
649  {
650  k=0;
651  for(j=1;j<=ROW;j++)
652  {
653  if(MA[i,j]!=0){k=1;}
654  }
656  }
657
658    m=transpose(m);
660    m=transpose(m);
661
662
663    int flag=1;  ///////////////counts if the underlying ring is r (flag mod 2 == 1) or ro (flag mod 2 == 0)
664
665    module TrafoL=freemodule(COL);
666    module TrafoR=freemodule(ROW);
668    module EXR=freemodule(ROW);
669
670    option(redSB);
671    option(redTail);
672    module STD_EX,LT,TSTD, L, Trafo;
673
674
675
676    module STD=transpose(m);
677    int finish=0;
678    list pos, COM, COM_EX;
679    matrix END, ENDSTD, STDFIN;
680    STDFIN=STD;
681    list COMPARE=STDFIN;
682
683
684  while(finish==0)
685    {
686      dbprint(ppl,"Going into the while cycle");
687      if(flag mod 2==1){STD_EX=EXL,transpose(STD); ex=2*COL;} else {STD_EX=EXR,transpose(STD); ex=2*ROW;}
688
689      dbprint(ppl,"Computing Groebner basis: start");
690      dbprint(ppl-1,STD);
691      STD=engine(STD,BASIS);
692      dbprint(ppl,"Computing Groebner basis: finished");
693      dbprint(ppl-1,STD);
694      if(flag mod 2==1){s=ROW; st=COL;}else{s=COL; st=ROW;}
695
696      STD_EX=transpose(STD_EX);
697      dbprint(ppl,"Computing Groebner basis for transformation matrix: start");
698      dbprint(ppl-1,STD_EX);
699      STD_EX=engine(STD_EX,BASIS);
700      dbprint(ppl,"Computing Groebner basis for transformation matrix: finished");
701      dbprint(ppl-1,STD_EX);
702
703      COM=1;
704      COM_EX=1;
705      for(i=2; i<=size(STD); i++)
706         { COM=COM[1..size(COM)],i; COM_EX=COM_EX[1..size(COM_EX)],i; }
707      nz=size(STD_EX)-size(STD);
708
709      //zero rows in the begining
710
711       if(size(STD)!=size(STD_EX) )
712       {
713         for(i=1; i<=size(STD_EX)-size(STD); i++)
714         {
715           COM_EX=0,COM_EX[1..size(COM_EX)];
716         }
717       }
718
719
720
721
722       for(i=nz+1; i<=size(STD_EX); i++ )
724        }
725
726        //assign the zero rows in STD_EX
727
728       for (j=2; j<=nz; j++)
729       {
730        p=0;
731        i=1;
732        while(STD_EX[j-1][i]==0){i++;};
733        p=i-1;
734        nextp=1;
735        k=0;
736        i=1;
737        while(STD_EX[j][i]==0 and i<=p)
738        { k++; i++;}
739        if (k==p){ COM_EX[j]=-1; }
740       }
741
742      //assign the zero rows in STD
743      for (j=2; j<=size(STD); j++)
744       {
745        i=size(transpose(STD));
746        while(STD[j-1][i]==0){i--;}
747        p=i;
748        i=size(transpose(STD[j]));
749        while(STD[j][i]==0){i--;}
750        if (i==p){ COM[j]=-1; }
751       }
752
753       for(j=1; j<=size(COM); j++)
754        {
755        if(COM[j]<0){COM=delete(COM,j);}
756        }
757
758      for(i=1; i<=size(COM_EX); i++)
759        {ff=0;
760         if(COM_EX[i]==0){ff=1;}
761         else
762          {  for(j=1; j<=size(COM); j++)
763               { if(COM_EX[i]==COM[j]){ff=1;} }
764          }
765         if(ff==0){COM_EX[i]=-1;}
766        }
767
768      L=STD_EX[1];
769      for(i=2; i<=size(COM_EX); i++)
770       {
771        if(COM_EX[i]!=-1){L=L,STD_EX[i];}
772       }
773
774      //////// split STD_EX in STD and the transformation matrix
775
776       L=transpose(L);
777       STD=L[st+1];
778       LT=L[1];
779
780
781       for(i=2; i<=s+st; i++)
782        {
783         if (i<=st)
784          {
785           LT=LT,L[i];
786          }
787         if (i>st+1)
788          {
789          STD=STD,L[i];
790          }
791        }
792
793       STD=transpose(STD);
794       STDFIN=matrix(STD);
795       COMPARE=insert(COMPARE,STDFIN);
796       LT=transpose(LT);
797
798      ////////////////////// compute the transformation matrices
799
800       if (flag mod 2 ==1)
801        {
802         TrafoL=transpose(LT)*TrafoL;
803        }
804       else
805        {
806         TrafoR=TrafoR*involution(LT,inv);
807        }
808
809
810      ///////////////////////// check whether the alg termined /////////////////
811      if(size(COMPARE)>=3)
812       {
813        if(STD==COMPARE[3]){finish=1;}
814       }
815       ////////////////////////////////// change to the opposite module
816       TSTD=transpose(STD);
817       STD=involution(TSTD,inv);
818       flag++;
819       dbprint(ppl,"Finished one while cycle");
820    }
821
822if (flag mod 2 ==0){ STD = involution(STD,inv);} else { STD = transpose(STD);  }
823
824    list REVERSE=TrafoL,STD,TrafoR;
825    return(REVERSE);
826}
827
828static proc divideByContent(matrix M)
829{
830//find first entrie not equal to zero
831int i,k;
832k=1;
833vector CON;
834   for(i=1;i<=ncols(M);i++)
835   {
837   }
838poly con=content(CON);
839matrix TL=1/con*freemodule(nrows(M));
840return(TL);
841}
842
843
844/////interesting examples for smith////////////////
845
846static proc Ex_One_wheeled_bicycle()
847{
848ring RA=(0,m), D, lp;
849matrix bicycle[2][3]=(1+m)*D^2, D^2, 1, D^2, D^2-1, 0;
850list s=smith(RA,bicycle, 1,0);
851print(s[2]);
852print(s[1]*bicycle*s[3]-s[2]);
853}
854
855
856static proc Ex_RLC-circuit()
857{
858ring RA=(0,m,R1,R2,L,C), D, lp;
859matrix  circuit[2][3]=D+1/(R1*C), 0, -1/(R1*C), 0, D+R2/L, -1/L;
860list s=smith(RA,circuit, 1,0);
861print(s[2]);
862print(s[1]*circuit*s[3]-s[2]);
863}
864
865
866static proc Ex_two_pendula()
867{
868ring RA=(0,m,M,L1,L2,m1,m2,g), D, lp;
869
870matrix 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,
871m2*L2^2*D^2-m2*L2*g,m2*L2*D^2,0;
872list s=smith(RA,pendula, 1,0);
873print(s[2]);
874print(s[1]*pendula*s[3]-s[2]);
875}
876
877
878
879static proc Ex_linerized_satellite_in_a_circular_equatorial_orbit()
880{
881ring RA=(0,m,omega,r,a,b), D, lp;
882
883matrix satellite[4][6]=
884D,-1,0,0,0,0,
885-3*omega^2,D,0,-2*omega*r,-a/m,0,
8860,0,D,-1,0,0,
8870,2*omega/r,0,D,0,-b/(m*r);
888
889list s=smith(RA,satellite, 1,0);
890print(s[2]);
891print(s[1]*satellite*s[3]-s[2]);
892}
893
895{
896ring RA=(0,M11,M12,M13,M21,M22,M31,M33,K1,K2), D, lp;
897
898matrix 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;
899list s=smith(RA,robot, 1,0);
900print(s[2]);
901print(s[1]*robot*s[3]-s[2]);
902}
903
904
905
906/////interesting examples for jacobson////////////////
907
908static proc Ex_compare_output_with_maple_package_JanetOre()
909{   ring w = 0,(x,d),Dp;
910    def W=nc_algebra(1,1);
911    setring W;
912    matrix m[3][3]=[d2,d+1,0],[d+1,0,d3-x2*d],[2d+1, d3+d2, d2];
913    list J=jacobson(W,m,0);
914    print(J[1]*m*J[3]);
915    print(J[2]);
916    print(J[1]);
917    print(J[3]);
918    print(J[1]*m*J[3]-J[2]);
919}
920
921
922static proc Ex_cyclic()
923{   ring w = 0,(x,d),Dp;
924    def W=nc_algebra(1,1);
925    setring W;
926    matrix m[3][3]=d,0,0,x*d+1,d,0,0,x*d,d;
927    list J=jacobson(W,m,0);
928    print(J[1]*m*J[3]);
929    print(J[2]);
930    print(J[1]);
931    print(J[3]);
932    print(J[1]*m*J[3]-J[2]);
933}
934
Note: See TracBrowser for help on using the repository browser.