# source:git/Singular/LIB/jacobson.lib@36e346

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