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

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