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

spielwiese
Last change on this file since cb8103a was f7d745, checked in by Hans Schoenemann <hannes@…>, 13 years ago
format git-svn-id: file:///usr/local/Singular/svn/trunk@14246 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100755
File size: 29.7 KB
Line 
1//////////////////////////////////////////////////////////////////////////////
2version="$Id$";
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@*           Viktor Levandovskyy,      levandov@math.rwth-aachen.de
8
9OVERVIEW:
10We work over a ring R, that is an Euclidean principal ideal domain.
11If R is commutative, we suppose R to be a polynomial ring in one variable.
12If R is non-commutative, we suppose R to have two variables, say x and d.
13We treat then the basering as the Ore localization of R
14with respect to the mult. closed set S = K[x] without 0.
15Thus, we treat basering as principal ideal ring with d a polynomial
16variable and x an invertible one.
17Note, that in computations no division by x will actually happen.
18
19Given a rectangular matrix M over R, one can compute unimodular (that is
20invertible) square matrices U and V, such that U*M*V=D is a diagonal matrix.
21Depending on the ring, the diagonal entries of D have certain properties.
22
23We call a square matrix D as before 'a weak Jacobson normal form of M'.
24It is known, that over the first rational Weyl algebra K(x)<d>, D can be further
25transformed into a diagonal matrix (1,1,...,1,f,0,..,0), where f is in K(x)<d>.
26We call such a form of D the strong Jacobson normal form.
27The existence of strong form in not guaranteed if one works with algebra,
28which is not rational Weyl algebra.
29
30
31REFERENCES:
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
39PROCEDURES:
40smith(M[,eng1,eng2]);  compute the Smith Normal Form of M over commutative ring
41jacobson(M[,eng]); compute a weak Jacobson Normal Form of M over non-commutative ring
42divideUnits(L);     create ones out of units in the output of smith or jacobson
43
44SEE ALSO: control_lib
45
46KEYWORDS: Jacobson form; Jacobson normal form; Smith form; Smith normal form; matrix diagonalization
47
48";
49
50LIB "poly.lib";
51LIB "involut.lib"; // involution
52LIB "dmodapp.lib"; // engine
53LIB "qhmoduli.lib"; // Min
54LIB "latex.lib"; // tex output for printlevel=4
55
56proc tstjacobson()
57{
58  /* tests all procs for consistency */
59  example divideUnits;
60  example smith;
61  example jacobson;
62}
63
64proc divideUnits(list L)
65"USAGE: divideUnits(L); list L
66RETURN: matrix or list of matrices
67ASSUME: 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
71PURPOSE: divide out units from the diagonal and reflect this in transformation matrices
72EXAMPLE: 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}
128example
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
142proc smith(matrix MA, list #)
143"USAGE: smith(M[, eng1, eng2]);  M matrix, eng1 and eng2 are optional integers
144RETURN: matrix or list of matrices, depending on arguments
145ASSUME: Basering is a commutative polynomial ring in one variable
146PURPOSE: compute the Smith Normal Form of M with (optionally) transformation matrices
147THEORY: Groebner bases are used for the Smith form like in [2] and [3].
148NOTE: 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.
157DISPLAY: If @code{printlevel}=1, progress debug messages will be printed,
158@* if @code{printlevel}>=2, all the debug messages will be printed.
159EXAMPLE: example smith; shows examples
160SEE 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}
229example
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
240static proc diagonal_with_trafo( R, matrix MA, int B)
241{
242
243int 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
464static 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
532static 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
689proc jacobson(matrix MA, list #)
690 "USAGE:  jacobson(M, eng);  M matrix, eng an optional int
691RETURN: list
692ASSUME: Basering is a (non-commutative) ring in two variables.
693PURPOSE: compute a weak Jacobson normal form of M over the basering
694THEORY: Groebner bases and involutions are used, following [3]
695NOTE: 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.
701DISPLAY: If @code{printlevel}=1, progress debug messages will be printed,
702@* if @code{printlevel}>=2, all the debug messages will be printed.
703EXAMPLE: example jacobson; shows examples
704SEE 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
769list RUECK=CON*TL, CON*MAA, TRR;
770return(RUECK);
771}
772example
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
791proc 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
995if (flag mod 2 ==0){ STD = involution(STD,inv);} else { STD = transpose(STD);  }
996
997    list REVERSE=TrafoL,STD,TrafoR;
998    return(REVERSE);
999}
1000
1001static proc divideByContent(matrix M)
1002{
1003//find first entrie not equal to zero
1004int i,k;
1005k=1;
1006vector 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   }
1011poly con=content(CON);
1012matrix TL=1/con*freemodule(nrows(M));
1013return(TL);
1014}
1015
1016
1017/////interesting examples for smith////////////////
1018
1019/*
1020
1021//static proc Ex_One_wheeled_bicycle()
1022{
1023ring RA=(0,m), D, lp;
1024matrix bicycle[2][3]=(1+m)*D^2, D^2, 1, D^2, D^2-1, 0;
1025list s=smith(bicycle, 1,0);
1026print(s[2]);
1027print(s[1]*bicycle*s[3]-s[2]);
1028}
1029
1030
1031//static proc Ex_RLC-circuit()
1032{
1033ring RA=(0,m,R1,R2,L,C), D, lp;
1034matrix  circuit[2][3]=D+1/(R1*C), 0, -1/(R1*C), 0, D+R2/L, -1/L;
1035list s=smith(circuit, 1,0);
1036print(s[2]);
1037print(s[1]*circuit*s[3]-s[2]);
1038}
1039
1040
1041//static proc Ex_two_pendula()
1042{
1043ring RA=(0,m,M,L1,L2,m1,m2,g), D, lp;
1044matrix 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,
1045m2*L2^2*D^2-m2*L2*g,m2*L2*D^2,0;
1046list s=smith(pendula, 1,0);
1047print(s[2]);
1048print(s[1]*pendula*s[3]-s[2]);
1049}
1050
1051//static proc Ex_linerized_satellite_in_a_circular_equatorial_orbit()
1052{
1053ring RA=(0,m,omega,r,a,b), D, lp;
1054matrix satellite[4][6]=
1055D,-1,0,0,0,0,
1056-3*omega^2,D,0,-2*omega*r,-a/m,0,
10570,0,D,-1,0,0,
10580,2*omega/r,0,D,0,-b/(m*r);
1059list s=smith(satellite, 1,0);
1060print(s[2]);
1061print(s[1]*satellite*s[3]-s[2]);
1062}
1063
1064//static proc Ex_flexible_one_link_robot()
1065{
1066ring RA=(0,M11,M12,M13,M21,M22,M31,M33,K1,K2), D, lp;
1067matrix 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;
1068list s=smith(robot, 1,0);
1069print(s[2]);
1070print(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!
1154matrix m[2][2]=s,x,0,s; print(m);
1155
1156// more interesting:
1157matrix n[3][3]=s,x,0,0,s,x,s,0,x;
1158
1159// things blow up
1160matrix m[2][3] = s,x^2*s,x^3*s,s*x^2,s*x+1,(x+1)^3;
1161matrix m[2][3] = s,x^2*s,x^3*s,s*x^2,s*x+1,(x+s)^2; // variation
1162matrix 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
1163matrix 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
1166matrix 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*/
Note: See TracBrowser for help on using the repository browser.