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