source: git/Singular/LIB/multigrading.lib @ 78f84b8

spielwiese
Last change on this file since 78f84b8 was ea87a9, checked in by Hans Schoenemann <hannes@…>, 13 years ago
tabs and spaces git-svn-id: file:///usr/local/Singular/svn/trunk@13403 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 54.9 KB
Line 
1///////////////////////////////////////////////////////////////////
2version="$Id$";
3category="Combinatorial Commutative Algebra";
4info="
5LIBRARY:  multigrading.lib          Multigraded Rings
6
7AUTHORS:  Rene Birkner, rbirkner@math.fu-berlin.de
8@*        Lars Kastner, lkastner@math.fu-berlin.de
9@*        Oleksandr Motsak, U@D, where U={motsak}, D={mathematik.uni-kl.de}
10
11
12OVERVIEW: This library allows one to virtually add multigradings to Singular.
13For more see http://code.google.com/p/convex-singular/wiki/Multigrading
14For theoretical references see:
15E. Miller, B. Sturmfels: 'Combinatorial Commutative Algebra' and
16M. Kreuzer, L. Robbiano: 'Computational Commutative Algebra'.
17
18NOTE: 'mDegBasis' relies on 4ti2 for computing Hilbert Bases.
19
20PROCEDURES:
21setBaseMultigrading(M,T); attach multiweights/torsion matrices to the basering
22getVariableWeights([R]);  get matrix of multidegrees of vars attached to a ring
23getTorsion([R]);          get torsion matrix attached to a ring
24
25setModuleGrading(M,v);    attach multiweights of units to a module and return it
26getModuleGrading(M);      get multiweights of module units (attached to M)
27
28mDeg(A);                  compute the multidegree of A
29mDegBasis(d);             compute all monomials of multidegree d
30mDegPartition(p);         compute the multigraded-homogenous components of p
31
32isTorsionFree();          test whether the current multigrading is torsion free
33isTorsionElement(p);      test whether p has zero multidegree
34isHomogenous(a);          test whether 'a' is multigraded-homogenous
35
36equalMDeg(e1,e2[,V]);     test whether e1==e2 in the current multigrading
37
38mDegGroebner(M);          compute the multigraded GB/SB of M
39mDegSyzygy(M);            compute the multigraded syzygies of M
40mDegResolution(M,l[,m]);  compute the multigraded resolution of M
41
42defineHomogenous(p);      get a torsion matrix wrt which p becomes homogenous
43pushForward(f);           find the finest grading on the image ring, homogenizing f
44
45hermite(A);               compute the Hermite Normal Form of a matrix
46
47hilbertSeries(M);         compute the multigraded Hilbert Series of M
48
49           (parameters in square brackets [] are optional)
50
51KEYWORDS:  multigradeding, multidegree, multiweights, multigraded-homogenous
52";
53
54// finestMDeg(def r)
55// newMap(map F, intmat Q, list #)
56
57LIB "standard.lib"; // for groebner
58
59/******************************************************/
60proc setBaseMultigrading(intmat M, list #)
61"USAGE: setBaseMultigrading(M[, T]); M, T are integer matrices
62PURPOSE: attaches weights of variables and torsion to the basering.
63NOTE: M encodes the weights of variables column-wise.
64The torsion is given by the lattice spanned by the columns of the integer
65matrix T in Z^nrows(M) over Z.
66RETURN: nothing
67EXAMPLE: example setBaseMultigrading; shows an example
68"
69{
70  string attrMgrad   = "mgrad";
71  string attrTorsion = "torsion";
72  string attrTorsionHNF = "torsion_hermite";
73
74
75  attrib(basering, attrMgrad, M);
76
77  if( size(#) > 0 and typeof(#[1]) == "intmat" )
78  {
79    def T = #[1];
80  } else
81  {
82    intmat T[nrows(M)][1];
83  }
84
85  if( nrows(T) == nrows(M) )
86  {
87    attrib(basering, attrTorsion, T);
88//    def H;
89//    attrib(basering, attrTorsionHNF, H);
90  }
91  else
92  {
93    ERROR("Incompatible matrix sizes!");
94  }
95}
96example
97{
98  "EXAMPLE:"; echo=2;
99
100  ring R = 0, (x, y, z), dp;
101
102  // Weights of variables
103  intmat M[3][3] =
104    1, 0, 0,
105    0, 1, 0,
106    0, 0, 1;
107
108  // Torsion:
109  intmat L[3][2] =
110    1, 1,
111    1, 3,
112    1, 5;
113
114  // attaches M & L to R (==basering):
115  setBaseMultigrading(M, L); // Grading: Z^3/L
116
117  // Weights are accessible via "getVariableWeights()":
118  getVariableWeights();
119
120  // Test all possible usages:
121  (getVariableWeights() == M) && (getVariableWeights(R) == M) && (getVariableWeights(basering) == M);
122
123  // Torsion is accessible via "getTorsion()":
124  getTorsion();
125
126  // Test all possible usages:
127  (getTorsion() == L) && (getTorsion(R) == L) && (getTorsion(basering) == L);
128
129  // And its hermite NF via getTorsion("hermite"):
130  getTorsion("hermite");
131
132  // Test all possible usages:
133  intmat H = hermite(L);
134  (getTorsion("hermite") == H) && (getTorsion(R, "hermite") == H) && (getTorsion(basering, "hermite") == H);
135
136  kill L, M;
137
138  // ----------- isomorphic multigrading -------- //
139
140  // Weights of variables
141  intmat M[2][3] =
142    1, -2, 1,
143    1,  1, 0;
144
145  // Torsion:
146  intmat L[2][1] =
147    0,
148    2;
149
150  // attaches M & L to R (==basering):
151  setBaseMultigrading(M, L); // Grading: Z + (Z/2Z)
152
153  // Weights are accessible via "getVariableWeights()":
154  getVariableWeights() == M;
155
156  // Torsion is accessible via "getTorsion()":
157  getTorsion() == L;
158
159  kill L, M;
160  // ----------- extreme case ------------ //
161
162  // Weights of variables
163  intmat M[1][3] =
164    1,  -1, 10;
165
166  // Torsion:
167  intmat L[1][1] =
168    0;
169
170  // attaches M & L to R (==basering):
171  setBaseMultigrading(M); // Grading: Z^3
172
173  // Weights are accessible via "getVariableWeights()":
174  getVariableWeights() == M;
175
176  // Torsion is accessible via "getTorsion()":
177  getTorsion() == L;
178}
179
180
181/******************************************************/
182proc getVariableWeights(list #)
183"USAGE: getVariableWeights([R])
184PURPOSE: get associated multigrading matrix for the basering [or R]
185RETURN:  intmat, matrix of multidegrees of variables
186EXAMPLE: example getVariableWeights; shows an example
187"
188{
189  string attrMgrad = "mgrad";
190
191
192  if( size(#) > 0 )
193  {
194    if(( typeof(#[1]) == "ring" ) || ( typeof(#[1]) == "qring" ))
195    {
196      def R = #[1];
197    }
198    else
199    {
200      ERROR("Optional argument must be a ring!");
201    }
202  }
203  else
204  {
205    def R = basering;
206  }
207
208  def M = attrib(R, attrMgrad);
209  if( typeof(M) == "intmat"){ return (M); }
210  ERROR( "Sorry no multigrading matrix!" );
211}
212example
213{
214  "EXAMPLE:"; echo=2;
215
216  ring R = 0, (x, y, z), dp;
217
218  // Weights of variables
219  intmat M[3][3] =
220    1, 0, 0,
221    0, 1, 0,
222    0, 0, 1;
223
224  // Torsion:
225  intmat L[3][2] =
226    1, 1,
227    1, 3,
228    1, 5;
229
230  // attaches M & L to R (==basering):
231  setBaseMultigrading(M, L); // Grading: Z^3/L
232
233  // Weights are accessible via "getVariableWeights()":
234  getVariableWeights() == M;
235
236  kill L, M;
237
238  // ----------- isomorphic multigrading -------- //
239
240  // Weights of variables
241  intmat M[2][3] =
242    1, -2, 1,
243    1,  1, 0;
244
245  // Torsion:
246  intmat L[2][1] =
247    0,
248    2;
249
250  // attaches M & L to R (==basering):
251  setBaseMultigrading(M, L); // Grading: Z + (Z/2Z)
252
253  // Weights are accessible via "getVariableWeights()":
254  getVariableWeights() == M;
255
256  kill L, M;
257
258  // ----------- extreme case ------------ //
259
260  // Weights of variables
261  intmat M[1][3] =
262    1,  -1, 10;
263
264  // Torsion:
265  intmat L[1][1] =
266    0;
267
268  // attaches M & L to R (==basering):
269  setBaseMultigrading(M); // Grading: Z^3
270
271  // Weights are accessible via "getVariableWeights()":
272  getVariableWeights() == M;
273}
274
275/******************************************************/
276proc getTorsion(list #)
277"USAGE: getTorsion([R[,opt]])
278PURPOSE: get associated torsion matrix, i.e. generators (cols) of the Torsion group
279RETURN: intmat, the torsion matrix, or its hermite normal form
280        if an optional argument (\"hermite\") is given
281EXAMPLE: example getTorsion; shows an example
282"
283{
284  string attrTorsion    = "torsion";
285  string attrTorsionHNF = "torsion_hermite";
286
287  int i = 1;
288
289  if( size(#) >= i )
290  {
291    if( ( typeof(#[i]) == "ring" ) or ( typeof(#[i]) == "qring" ) )
292    {
293      def R = #[i];
294      i++;
295    }
296  }
297
298  if( !defined(R) )
299  {
300    def R = basering;
301  }
302
303  if( size(#) >= i )
304  {
305    if( #[i] == "hermite" )
306    {
307      if( typeof(attrib(R, attrTorsionHNF)) != "intmat" )
308      {
309        def M = getTorsion(R);
310        if( typeof(M) != "intmat")
311        {
312          ERROR( "Sorry no torsion matrix!" );
313        }
314        M = hermite(M);
315        attrib(R, attrTorsionHNF, M); // this might not work with R != basering...
316      }
317      return (attrib(R, attrTorsionHNF));
318    }
319  }
320
321  def M = attrib(R, attrTorsion);
322  if( typeof(M) != "intmat")
323  {
324    ERROR( "Sorry no torsion matrix!" );
325  }
326  return (M);
327}
328example
329{
330  "EXAMPLE:"; echo=2;
331
332  ring R = 0, (x, y, z), dp;
333
334  // Weights of variables
335  intmat M[3][3] =
336    1, 0, 0,
337    0, 1, 0,
338    0, 0, 1;
339
340  // Torsion:
341  intmat L[3][2] =
342    1, 1,
343    1, 3,
344    1, 5;
345
346  // attaches M & L to R (==basering):
347  setBaseMultigrading(M, L); // Grading: Z^3/L
348
349  // Torsion is accessible via "getTorsion()":
350  getTorsion() == L;
351
352  // its hermite NF:
353  print(getTorsion("hermite"));
354
355  kill L, M;
356
357  // ----------- isomorphic multigrading -------- //
358
359  // Weights of variables
360  intmat M[2][3] =
361    1, -2, 1,
362    1,  1, 0;
363
364  // Torsion:
365  intmat L[2][1] =
366    0,
367    2;
368
369  // attaches M & L to R (==basering):
370  setBaseMultigrading(M, L); // Grading: Z + (Z/2Z)
371
372  // Torsion is accessible via "getTorsion()":
373  getTorsion() == L;
374
375  // its hermite NF:
376  print(getTorsion("hermite"));
377
378  kill L, M;
379
380  // ----------- extreme case ------------ //
381
382  // Weights of variables
383  intmat M[1][3] =
384    1,  -1, 10;
385
386  // Torsion:
387  intmat L[1][1] =
388    0;
389
390  // attaches M & L to R (==basering):
391  setBaseMultigrading(M); // Grading: Z^3
392
393  // Torsion is accessible via "getTorsion()":
394  getTorsion() == L;
395
396  // its hermite NF:
397  print(getTorsion("hermite"));
398}
399
400/******************************************************/
401proc getModuleGrading(def m)
402"USAGE: getModuleGrading(m), 'm' module/vector
403RETURN: integer matrix of the multiweights of free module generators attached to 'm'
404EXAMPLE: example getModuleGrading; shows an example
405"
406{
407  string attrModuleGrading = "genWeights";
408
409  //  print(m);  typeof(m);  attrib(m);
410
411  def V = attrib(m, attrModuleGrading);
412
413  if( typeof(V) != "intmat" )
414  {
415    if( (typeof(m) == "ideal") or (typeof(m) == "poly") )
416    {
417      intmat M = getVariableWeights();
418      intmat VV[nrows(M)][1];
419      return (VV);
420    }
421
422    ERROR("Sorry: vector or module need module-grading-matrix! See 'getModuleGrading'.");
423  }
424
425  if( nrows(V) != nrows(getVariableWeights()) )
426  {
427    ERROR("Sorry wrong height of V: " + string(nrows(V)));
428  }
429
430  if( ncols(V) < nrows(m) )
431  {
432    ERROR("Sorry wrong width of V: " + string(ncols(V)));
433  }
434
435  return (V);
436}
437example
438{
439  "EXAMPLE:"; echo=2;
440
441   ring R = 0, (x,y), dp;
442   intmat M[2][2]=
443     1, 1,
444     0, 2;
445   intmat T[2][5]=
446     1,  2,  3,  4, 0,
447     0, 10, 20, 30, 1;
448
449   setBaseMultigrading(M, T);
450
451   ideal I = x, y, xy^5;
452   isHomogenous(I);
453
454   intmat V = mDeg(I); print(V);
455
456   module S = syz(I); print(S);
457
458   S = setModuleGrading(S, V);
459
460   getModuleGrading(S) == V;
461
462   vector v = setModuleGrading(S[1], V);
463   getModuleGrading(v) == V;
464   isHomogenous(v);
465   print( mDeg(v) );
466
467   isHomogenous(S);
468   print( mDeg(S) );
469}
470
471/******************************************************/
472proc setModuleGrading(def m, intmat G)
473"USAGE: setModuleGrading(m, G), m module/vector, G intmat
474PURPOSE: attaches the multiweights of free module generators to 'm'
475WARNING: The method does not verify whether the multigrading makes the
476         module/vector homogenous. One can do that using isHomogenous(m).
477EXAMPLE: example setModuleGrading; shows an example
478"
479{
480  string attrModuleGrading = "genWeights";
481
482  intmat R = getVariableWeights();
483
484  if(nrows(G) != nrows(R)){ ERROR("Incompatible gradings.");}
485  if(ncols(G) < nrows(m)){ ERROR("Multigrading does not fit to module.");}
486
487  attrib(m, attrModuleGrading, G);
488  return(m);
489}
490example
491{
492  "EXAMPLE:"; echo=2;
493
494   ring R = 0, (x,y), dp;
495   intmat M[2][2]=
496     1, 1,
497     0, 2;
498   intmat T[2][5]=
499     1,  2,  3,  4, 0,
500     0, 10, 20, 30, 1;
501
502   setBaseMultigrading(M, T);
503
504   ideal I = x, y, xy^5;
505   intmat V = mDeg(I);
506
507   // V == M; modulo T
508   print(V);
509
510   module S = syz(I);
511
512   S = setModuleGrading(S, V);
513   getModuleGrading(S) == V;
514
515   print(S);
516
517   vector v = S[1]; v = setModuleGrading(v, V);
518   getModuleGrading(v) == V;
519
520   print( mDeg(v) );
521
522   isHomogenous(S);
523
524   print( mDeg(S) );
525}
526
527
528
529
530/******************************************************/
531proc isTorsionFree()
532"USAGE: isTorsionFree()
533PURPOSE: Determines whether the multigrading attached to the current ring is torsion-free.
534RETURN: boolean, the result of the test
535EXAMPLE: example isTorsionFree; shows an example
536"
537{
538
539  intmat H = hermite(transpose(getTorsion("hermite"))); // TODO: ?cache it?  //******
540
541  int i, j;
542  int r = nrows(H);
543  int c = ncols(H);
544  int d = 1;
545  for( i = 1; (i <= c) && (i <= r); i++ )
546  {
547    for( j = i; (H[j, i] == 0)&&(j < r); j++ )
548    {
549    }
550
551    if(H[j, i]!=0)
552    {
553      d=d*H[j, i];
554    }
555  }
556
557  if( (d*d)==1 )
558  {
559    return(1==1);
560  }
561  return(0==1);
562}
563example
564{
565  "EXAMPLE:"; echo=2;
566
567   ring R = 0,(x,y),dp;
568   intmat M[2][2]=
569     1,0,
570     0,1;
571   intmat T[2][5]=
572     1, 2, 3, 4, 0,
573     0,10,20,30, 1;
574
575   setBaseMultigrading(M,T);
576
577   // Is the resulting group torsion free?
578   isTorsionFree();
579
580   kill R, M, T;
581   ///////////////////////////////////////////
582
583   ring R=0,(x,y,z),dp;
584   intmat A[3][3] =
585     1,0,0,
586     0,1,0,
587     0,0,1;
588   intmat B[3][4]=
589     3,3,3,3,
590     2,1,3,0,
591     1,2,0,3;
592   setBaseMultigrading(A,B);
593   // Is the resulting group torsion free?
594   isTorsionFree();
595
596   kill R, A, B;
597}
598
599
600
601/******************************************************/
602proc hermite(intmat A)
603"USAGE: hermite( A );
604PURPOSE: Computes the (lower triangular) Hermite Normal Form
605           of the matrix A by column operations.
606RETURN: intmat, the Hermite Normal Form of A
607EXAMPLE: example hermite; shows an example
608"
609{
610  intmat g;
611  int i,j,k, save, d, a1, a2, c, l;
612  c=0;
613  l=0;
614  for(i=1; ((i+l)<=nrows(A))&&((i+l)<=ncols(A)); i++)
615  {
616    //i;
617    if(A[i+l, i]==0)
618    {
619      for(j=i;j<=ncols(A); j++)
620      {
621        if(A[i+l, j]!=0)
622        {
623          for(k=1;k<=nrows(A);k++)
624          {
625            save=A[k, i];
626            A[k, i]=A[k, j];
627            A[k, j]=save;
628          }
629          break;
630        }
631        if(j==ncols(A)){ c=1; l=l+1; }
632      }
633    }
634
635    if((i+l)>nrows(A)){ break; }
636
637    if(A[i+l, i]<0)
638    {
639      for(k=1;k<=nrows(A);k++){ A[k, i]=A[k, i]*(-1); }
640    }
641
642    if(c==0)
643    {
644      for(j=i+1;j<=ncols(A);j++)
645      {
646        //A;
647        if(A[i+l, j]<0)
648        {
649          for(k=1;k<=nrows(A);k++){ A[k, j]=(-1)*A[k, j];}
650        }
651
652        if(A[i+l, i]==1)
653        {
654          d=A[i+l, j];
655          for(k=1;k<=nrows(A);k++)
656          {
657            A[k, j]=A[k, j]-d*A[k, i];
658          }
659        }
660        else
661        {
662          while((A[i+l, i] * A[i+l, j])!=0)
663          {
664            if(A[i+l, i]> A[i+l, j])
665            {
666
667              for(k=1;k<=nrows(A);k++)
668              {
669                save=A[k, i];
670                A[k, i]=A[k, j];
671                A[k, j]=save;
672              }
673            }
674            a1=A[i+l, j]%A[i+l,i];
675            a2=A[i+l, j]-a1;
676            d=a2/A[i+l, i];
677            for(k=1;k<=nrows(A);k++)
678            {
679              A[k, j]=A[k, j]- d*A[k, i];
680            }
681          }
682        }
683      }
684      for(j=1;j<i;j++)
685      {
686        a1=A[i+l, j]%A[i+l,i];
687        d=(A[i+l, j]-a1)/A[i+l, i];
688        for(k=1;k<=nrows(A);k++){ A[k, j]=A[k, j]-d*A[k, i];}
689      }
690    }
691    c=0;
692  }
693  return( A);
694}
695example
696{
697  "EXAMPLE:"; echo=2;
698
699   intmat M[2][5] =
700     1, 2, 3, 4, 0,
701     0,10,20,30, 1;
702
703   // Hermite Normal Form of M:
704   print(hermite(M));
705
706   intmat T[3][4] =
707     3,3,3,3,
708     2,1,3,0,
709     1,2,0,3;
710
711   // Hermite Normal Form of T:
712   print(hermite(T));
713
714   intmat A[4][5] =
715     1,2,3,2,2,
716     1,2,3,4,0,
717     0,5,4,2,1,
718     3,2,4,0,2;
719
720   // Hermite Normal Form of A:
721   print(hermite(A));
722}
723
724
725/******************************************************/
726proc isTorsionElement(intvec mdeg)
727"USAGE: isTorsionElement(intvec mdeg);
728PURPOSE: For a integer vector mdeg representing the multidegree of some polynomial
729or vector this method computes if the multidegree is contained in the torsion
730group, i.e. if it is zero in the multigrading.
731EXAMPLE: example isTorsionElement; shows an example
732"
733{
734  intmat H = getTorsion("hermite");
735  int x, k, i;
736
737  int r = nrows(H);
738  int c = ncols(H);
739
740  int rr = nrows(mdeg);
741
742  for( i = 1; (i<=r) && (i<=c); i++)
743  {
744    if(H[i, i]!=0)
745    {
746      x = mdeg[i]%H[i, i];
747
748      if(x!=0)
749      {
750        return(1==0);
751      }
752
753      x = mdeg[i] / H[i,i];
754
755      for( k=1; k <= rr; k++)
756      {
757        mdeg[k] = mdeg[k] - x*H[k,i];
758      }
759    }
760  }
761
762  return( mdeg == 0 );
763
764}
765example
766{
767 "EXAMPLE:"; echo=2;
768
769  ring r = 0,(x,y,z),dp;
770
771  intmat g[2][3]=
772    1,0,1,
773    0,1,1;
774  intmat t[2][1]=
775    -2,
776    1;
777
778  setBaseMultigrading(g,t);
779
780  poly a = x10yz;
781  poly b = x8y2z;
782  poly c = x4z2;
783  poly d = y5;
784  poly e = x2y2;
785  poly f = z2;
786
787  intvec v1 = mDeg(a) - mDeg(b);
788  v1;
789  isTorsionElement(v1);
790
791  intvec v2 = mDeg(a) - mDeg(c);
792  v2;
793  isTorsionElement(v2);
794
795  intvec v3 = mDeg(e) - mDeg(f);
796  v3;
797  isTorsionElement(v3);
798
799  intvec v4 = mDeg(c) - mDeg(d);
800  v4;
801  isTorsionElement(v4);
802}
803
804
805/******************************************************/
806proc defineHomogenous(poly f, list #)
807"USAGE: defineHomogenous(f[, G]); polynomial f, integer matrix G
808PURPOSE: Yields a matrix which has to be appended to the torsion matrix to make the
809polynomial f homogenous in the grading by grad.
810EXAMPLE: example defineHomogenous; shows an example
811"
812{
813  if( size(#) > 0 )
814  {
815    if( typeof(#[1]) == "intmat" )
816    {
817      intmat grad = #[1];
818    }
819  }
820
821  if( !defined(grad) )
822  {
823    intmat grad = getVariableWeights();
824  }
825
826  intmat newtor[nrows(grad)][size(f)-1];
827  int i,j;
828  intvec l = grad*leadexp(f);
829  intvec v;
830  for(i=2; i <= size(f); i++)
831  {
832    v = grad * leadexp(f[i]) - l;
833    for( j=1; j<=size(v); j++)
834    {
835      newtor[j,i-1] = v[j];
836    }
837  }
838  return(newtor);
839}
840example
841{
842  "EXAMPLE:"; echo=2;
843
844  ring r =0,(x,y,z),dp;
845  intmat grad[2][3] =
846    1,0,1,
847    0,1,1;
848
849  setBaseMultigrading(grad);
850
851  poly f = x2y3-z5+x-3zx;
852
853  intmat M = defineHomogenous(f);
854  M;
855  defineHomogenous(f, grad) == M;
856
857  isHomogenous(f);
858  setBaseMultigrading(grad, M);
859  isHomogenous(f);
860}
861
862/******************************************************/
863proc pushForward(map f)
864"USAGE: pushForward(f);
865PURPOSE: Computes the finest grading of the image ring which makes the map f
866a map of graded rings. The group map between the two grading groups is given
867by transpose( (Id, 0) ). Pay attention that the group spanned by the columns of
868the torsion matrix may not be a subgroup of the grading group. Still all columns
869are needed to find the correct image of the preimage gradings.
870EXAMPLE: example pushForward; shows an example
871"
872{
873
874  int k,i,j;
875//  f;
876
877//  listvar();
878  def pre = preimage(f);
879
880//  "pre: ";  pre;
881
882  intmat oldgrad=getVariableWeights(pre);
883  intmat oldtor=getTorsion(pre);
884
885  int n=nvars(pre);
886  int np=nvars(basering);
887  int p=nrows(oldgrad);
888  int pp=p+np;
889
890  intmat newgrad[pp][np];
891
892  for(i=1;i<=np;i++){ newgrad[p+i,i]=1;}
893
894  //newgrad;
895
896
897
898  list newtor;
899  intmat toadd;
900  int columns=0;
901
902  intmat toadd1[pp][n];
903  intvec v;
904  poly im;
905
906  for(i=1;i<=p;i++){
907    for(j=1;j<=n;j++){ toadd1[i,j]=oldgrad[i,j];}
908  }
909
910  for(i=1;i<=n;i++){
911    im=f[i];
912    //im;
913    toadd = defineHomogenous(im, newgrad);
914    newtor=insert(newtor,toadd);
915    columns=columns+ncols(toadd);
916
917    v=leadexp(f[i]);
918    for(j=p+1;j<=p+np;j++){ toadd1[j,i]=-v[j-p];}
919  }
920
921  newtor=insert(newtor,toadd1);
922  columns=columns+ncols(toadd1);
923
924
925  if(typeof(basering)=="qring"){
926    //"Entering qring";
927    ideal a=ideal(basering);
928    for(i=1;i<=size(a);i++){
929      toadd = defineHomogenous(a[i], newgrad);
930      //toadd;
931      columns=columns+ncols(toadd);
932      newtor=insert(newtor,toadd);
933    }
934  }
935
936  //newtor;
937  intmat imofoldtor[pp][ncols(oldtor)];
938  for(i=1; i<=nrows(oldtor);i++){
939    for(j=1; j<=ncols(oldtor); j++){
940      imofoldtor[i,j]=oldtor[i,j];
941    }
942  }
943
944  columns=columns+ncols(oldtor);
945  newtor=insert(newtor, imofoldtor);
946
947  intmat torsion[pp][columns];
948  columns=0;
949  for(k=1;k<=size(newtor);k++){
950    for(i=1;i<=pp;i++){
951      for(j=1;j<=ncols(newtor[k]);j++){torsion[i,j+columns]=newtor[k][i,j];}
952    }
953    columns=columns+ncols(newtor[k]);
954  }
955
956  torsion=hermite(torsion);
957  intmat result[pp][pp];
958  for(i=1;i<=pp;i++){
959    for(j=1;j<=pp;j++){result[i,j]=torsion[i,j];}
960  }
961
962  setBaseMultigrading(newgrad, result);
963
964}
965example
966{
967  "EXAMPLE:"; echo=2;
968
969  ring r = 0,(x,y,z),dp;
970
971
972
973  // Setting degrees for preimage ring.;
974  intmat grad[3][3] =
975    1,0,0,
976    0,1,0,
977    0,0,1;
978
979  setBaseMultigrading(grad);
980
981  // grading on r:
982  getVariableWeights();
983  getTorsion();
984
985  // only for the purpose of this example
986  if( voice > 1 ){ /*keepring(r);*/ export(r); }
987
988  ring R = 0,(a,b),dp;
989  ideal i = a2-b2+a6-b5+ab3,a7b+b15-ab6+a6b6;
990
991  // The quotient ring by this ideal will become our image ring.;
992  qring Q = std(i);
993
994  listvar();
995
996  map f = r,-a2b6+b5+a3b+a2+ab,-a2b7-3a2b5+b4+a,a6-b6-b3+a2; f;
997
998
999  // TODO: Unfortunately this is not a very spectacular example...:
1000  // Pushing forward f:
1001  pushForward(f);
1002
1003  // due to pushForward we have got new grading on Q
1004  getVariableWeights();
1005  getTorsion();
1006
1007
1008  // only for the purpose of this example
1009  if( voice > 1 ){ kill r; }
1010
1011}
1012
1013
1014/******************************************************/
1015proc equalMDeg(intvec exp1, intvec exp2, list #)
1016"USAGE: equalMDeg(exp1, exp2[, V]); intvec exp1, exp2, intmat V
1017PURPOSE: Tests if the exponent vectors of two monomials (given by exp1 and exp2)
1018represent the same multidegree.
1019NOTE: the integer matrix V encodes multidegrees of module components,
1020if module component is present in exp1 and exp2
1021EXAMPLE: example equalMDeg; shows an example
1022"
1023{
1024  if( size(exp1) != size(exp2) )
1025  {
1026    ERROR("Sorry: we cannot compare exponents comming from a polynomial and a vector yet!");
1027  }
1028
1029  if( exp1 == exp2)
1030  {
1031    return (1==1);
1032  }
1033
1034
1035
1036  intmat M = getVariableWeights();
1037
1038  if( nrows(exp1) > ncols(M) ) // vectors => last exponent is the module component!
1039  {
1040    if( (size(#) == 0) or (typeof(#[1])!="intmat") )
1041    {
1042      ERROR("Sorry: wrong or missing module-unit-weights-matrix V!");
1043    }
1044    intmat V = #[1];
1045
1046    // typeof(V); print(V);
1047
1048    int N = ncols(M);
1049    int r = nrows(M);
1050
1051    intvec d = intvec(exp1[1..N]) - intvec(exp2[1..N]);
1052    intvec dm = intvec(V[1..r, exp1[N+1]]) - intvec(V[1..r, exp2[N+1]]);
1053
1054    intvec difference = M * d + dm;
1055  }
1056  else
1057  {
1058    intvec d = (exp1 - exp2);
1059    intvec difference = M * d;
1060  }
1061
1062  if (isFreeRepresented()) // no torsion!?
1063  {
1064    return ( difference == 0);
1065  }
1066  return ( isTorsionElement( difference ) );
1067}
1068example
1069{
1070  "EXAMPLE:"; echo=2;
1071
1072  ring r = 0,(x,y,z),dp;
1073
1074  intmat g[2][3]=
1075    1,0,1,
1076    0,1,1;
1077
1078  intmat t[2][1]=
1079    -2,
1080    1;
1081
1082  setBaseMultigrading(g,t);
1083
1084  poly a = x10yz;
1085  poly b = x8y2z;
1086  poly c = x4z2;
1087  poly d = y5;
1088  poly e = x2y2;
1089  poly f = z2;
1090
1091
1092  equalMDeg(leadexp(a), leadexp(b));
1093  equalMDeg(leadexp(a), leadexp(c));
1094  equalMDeg(leadexp(a), leadexp(d));
1095  equalMDeg(leadexp(a), leadexp(e));
1096  equalMDeg(leadexp(a), leadexp(f));
1097
1098  equalMDeg(leadexp(b), leadexp(c));
1099  equalMDeg(leadexp(b), leadexp(d));
1100  equalMDeg(leadexp(b), leadexp(e));
1101  equalMDeg(leadexp(b), leadexp(f));
1102
1103  equalMDeg(leadexp(c), leadexp(d));
1104  equalMDeg(leadexp(c), leadexp(e));
1105  equalMDeg(leadexp(c), leadexp(f));
1106
1107  equalMDeg(leadexp(d), leadexp(e));
1108  equalMDeg(leadexp(d), leadexp(f));
1109
1110  equalMDeg(leadexp(e), leadexp(f));
1111
1112}
1113
1114
1115
1116/******************************************************/
1117static proc isFreeRepresented()
1118"check whether the base muligrading is torsion-free (it is zero).
1119"
1120{
1121  intmat T = getTorsion();
1122
1123  intmat Z[nrows(T)][ncols(T)];
1124
1125  return (T == Z); // no torsion!
1126}
1127
1128
1129/******************************************************/
1130proc isHomogenous(def a, list #)
1131"USAGE: isHomogenous(a[, f]); a polynomial/vector/ideal/module
1132RETURN: boolean, TRUE if a is (multi)homogenous, and FALSE otherwise
1133EXAMPLE: example isHomogenous; shows an example
1134"
1135{
1136  if( (typeof(a) == "poly") or (typeof(a) == "vector") )
1137  {
1138    return ( size(mDegPartition(a)) <= 1 )
1139  }
1140
1141  if( typeof(a) == "vector" or typeof(a) == "module" )
1142  {
1143    intmat V = getModuleGrading(a);
1144  }
1145
1146  if( (typeof(a) == "ideal") or (typeof(a) == "module") )
1147  {
1148    if(size(#) > 0)
1149    {
1150      if (#[1] == "checkGens")
1151      {
1152        def aa;
1153        for( int i = ncols(a); i > 0; i-- )
1154        {
1155          aa = a[i];
1156
1157          if(defined(V))
1158          {
1159            aa = setModuleGrading(aa, V);
1160          }
1161
1162          if(!isHomogenous(aa))
1163          {
1164            return(0==1);
1165          }
1166        }
1167        return(1==1);
1168      }
1169    }
1170
1171    def g = groebner(a); // !!!!
1172
1173    def b, aa; int j;
1174    for( int i = ncols(a); i > 0; i-- )
1175    {
1176      aa = a[i];
1177
1178      if(defined(V))
1179      {
1180        aa = setModuleGrading(aa, V);
1181      }
1182
1183      b = mDegPartition(aa);
1184      for( j = ncols(b); j > 0; j-- )
1185      {
1186        if(NF(b[j],g) != 0)
1187        {
1188          return(0==1);
1189        }
1190      }
1191    }
1192    return(1==1);
1193  }
1194}
1195example
1196{
1197  "EXAMPLE:"; echo=2;
1198
1199  ring r = 0,(x,y,z),dp;
1200
1201  //Grading and Torsion matrices:
1202  intmat M[3][3] =
1203    1,0,0,
1204    0,1,0,
1205    0,0,1;
1206
1207  intmat T[3][1] =
1208    1,2,3;
1209
1210  setBaseMultigrading(M,T);
1211
1212  attrib(r);
1213
1214  poly f = x-yz;
1215
1216  mDegPartition(f);
1217  print(mDeg(_));
1218
1219  isHomogenous(f);   // f: is not homogenous
1220
1221  poly g = 1-xy2z3;
1222  isHomogenous(g); // g: is homogenous
1223  mDegPartition(g);
1224
1225  kill T;
1226  /////////////////////////////////////////////////////////
1227  // new Torsion matrix:
1228  intmat T[3][4] =
1229    3,3,3,3,
1230    2,1,3,0,
1231    1,2,0,3;
1232
1233  setBaseMultigrading(M,T);
1234
1235  f;
1236  isHomogenous(f);
1237  mDegPartition(f);
1238
1239  // ---------------------
1240  g;
1241  isHomogenous(g);
1242  mDegPartition(g);
1243
1244  kill r, T, M;
1245
1246  ring R = 0, (x,y,z), dp;
1247
1248  intmat A[2][3] =
1249    0,0,1,
1250    3,2,1;
1251  intmat T[2][1] =
1252    -1,
1253     4;
1254  setBaseMultigrading(A, T);
1255
1256  isHomogenous(ideal(x2 - y3 -xy +z, x*y-z, x^3 - y^2*z + x^2 -y^3)); // 1
1257  isHomogenous(ideal(x2 - y3 -xy +z, x*y-z, x^3 - y^2*z + x^2 -y^3), "checkGens");
1258  isHomogenous(ideal(x+y, x2 - y2)); // 0
1259
1260  // Degree partition:
1261  mDegPartition(x2 - y3 -xy +z);
1262  mDegPartition(x3 -y2z + x2 -y3 + z + 1);
1263
1264
1265  module N = gen(1) + (x+y) * gen(2), z*gen(3);
1266
1267  intmat V[2][3] = 0; // 1, 2, 3,  4, 5, 6; //  column-wise weights of components!!??
1268
1269  vector v1, v2;
1270
1271  v1 = setModuleGrading(N[1], V); v1;
1272  mDegPartition(v1);
1273  print( mDeg(_) );
1274
1275  v2 = setModuleGrading(N[2], V); v2;
1276  mDegPartition(v2);
1277  print( mDeg(_) );
1278
1279  N = setModuleGrading(N, V);
1280  isHomogenous(N);
1281  print( mDeg(N) );
1282
1283  ///////////////////////////////////////
1284
1285  V =
1286    1, 2, 3,
1287    4, 5, 6;
1288
1289  v1 = setModuleGrading(N[1], V); v1;
1290  mDegPartition(v1);
1291  print( mDeg(_) );
1292
1293  v2 = setModuleGrading(N[2], V); v2;
1294  mDegPartition(v2);
1295  print( mDeg(_) );
1296
1297  N = setModuleGrading(N, V);
1298  isHomogenous(N);
1299  print( mDeg(N) );
1300
1301  ///////////////////////////////////////
1302
1303  V =
1304    0, 0, 0,
1305    4, 1, 0;
1306
1307  N = gen(1) + x * gen(2), z*gen(3);
1308  N = setModuleGrading(N, V); print(N);
1309  isHomogenous(N);
1310  print( mDeg(N) );
1311  v1 = setModuleGrading(N[1], V); print(v1);
1312  mDegPartition(v1);
1313  print( mDeg(_) );
1314  N = setModuleGrading(N, V); print(N);
1315  isHomogenous(N);
1316  print( mDeg(N) );
1317}
1318
1319/******************************************************/
1320proc mDeg(def A)
1321"USAGE: mDeg(A); polynomial/vector/ideal/module A
1322PURPOSE: compute multidegree
1323EXAMPLE: example mDeg; shows an example
1324"
1325{
1326  if( defined(attrib(A, "grad")) > 0 )
1327  {
1328    return (attrib(A, "grad"));
1329  }
1330
1331  intmat M = getVariableWeights();
1332  int N = nvars(basering);
1333
1334  if( ncols(M) != N )
1335  {
1336    ERROR("Sorry wrong mgrad-size of M: " + string(ncols(M)));
1337  }
1338
1339  int r = nrows(M);
1340
1341  if( A == 0 )
1342  {
1343    intvec v; v[r] = 0;
1344    return (v);
1345  }
1346
1347  if( (typeof(A) == "vector") or (typeof(A) == "module") )
1348  {
1349    intmat V = getModuleGrading(A);
1350
1351    if( nrows(V) != r )
1352    {
1353      ERROR("Sorry wrong mgrad-size of V: " + string(nrows(V)));
1354    }
1355  }
1356
1357  intvec m; m[r] = 0;
1358
1359  if( typeof(A) == "poly" )
1360  {
1361    intvec v = leadexp(A); //  v;
1362    m = M * v;
1363
1364    // We assume homogeneous input!
1365    return(m);
1366
1367    A = A - lead(A);
1368    while( size(A) > 0 )
1369    {
1370      v = leadexp(A); //  v;
1371      m = max( m, M * v, r ); // ????
1372      A = A - lead(A);
1373    }
1374
1375    return(m);
1376  }
1377
1378
1379  if( typeof(A) == "vector" )
1380  {
1381    intvec v;
1382    v = leadexp(A); //  v;
1383    m = intvec(M * intvec(v[1..N])) + intvec(V[1..r, v[N+1]]);
1384
1385    // We assume homogeneous input!
1386    return(m);
1387
1388    A = A - lead(A);
1389    while( size(A) > 0 )
1390    {
1391      v = leadexp(A); //  v;
1392
1393      // intvec(M * intvec(v[1..N])) + intvec(V[1..r, v[N+1]]);
1394
1395      m = max( m, intvec(M * intvec(v[1..N])) + intvec(V[1..r, v[N+1]]), r ); // ???
1396
1397      A = A - lead(A);
1398    }
1399
1400    return(m);
1401  }
1402
1403  int i, j; intvec d;
1404
1405  if( typeof(A) == "ideal" )
1406  {
1407    intmat G[ r ] [ ncols(A)];
1408    for( i = ncols(A); i > 0; i-- )
1409    {
1410      d = mDeg( A[i] );
1411
1412      for( j = 1; j <= r; j++ ) // see ticket: 253
1413      {
1414        G[j, i] = d[j];
1415      }
1416    }
1417    return(G);
1418  }
1419
1420  if( typeof(A) == "module" )
1421  {
1422    intmat G[ r ] [ ncols(A)];
1423    vector v;
1424
1425    for( i = ncols(A); i > 0; i-- )
1426    {
1427      v = setModuleGrading(A[i], V);
1428
1429      // G[1..r, i]
1430      d = mDeg(v);
1431
1432      for( j = 1; j <= r; j++ ) // see ticket: 253
1433      {
1434        G[j, i] = d[j];
1435      }
1436
1437    }
1438
1439    return(G);
1440  }
1441
1442}
1443example
1444{
1445  "EXAMPLE:"; echo=2;
1446
1447  ring r = 0,(x, y), dp;
1448
1449  intmat A[2][2] = 1, 0, 0, 1;
1450  print(A);
1451
1452  intmat Ta[2][1] = 0, 3;
1453  print(Ta);
1454
1455  //   attrib(A, "torsion", Ta); // to think about
1456
1457//  "poly:";
1458  setBaseMultigrading(A);
1459
1460
1461  mDeg( x*x, A );
1462  mDeg( y*y*y, A );
1463
1464  setBaseMultigrading(A, Ta);
1465
1466  mDeg( x*x*y );
1467
1468  mDeg( y*y*y*x );
1469
1470  mDeg( x*y + x + 1 );
1471
1472  mDegPartition(x*y + x + 1);
1473
1474  print ( mDeg(0) );
1475  poly zero = 0;
1476  print ( mDeg(zero) );
1477
1478//  "ideal:";
1479
1480  ideal I = y*x*x, x*y*y*y;
1481  print( mDeg(I) );
1482
1483  print ( mDeg(ideal(0)) );
1484  print ( mDeg(ideal(0,0,0)) );
1485
1486//  "vectors:";
1487
1488  intmat B[2][2] = 0, 1, 1, 0;
1489  print(B);
1490
1491  mDeg( setModuleGrading(y*y*y*gen(2), B ));
1492  mDeg( setModuleGrading(x*x*gen(1), B ));
1493
1494
1495  vector V = x*gen(1) + y*gen(2);
1496  V = setModuleGrading(V, B);
1497  mDeg( V );
1498
1499  vector v1 = setModuleGrading([0, 0, 0], B);
1500  print( mDeg( v1 ) );
1501
1502  vector v2 = setModuleGrading([0], B);
1503  print( mDeg( v2 ) );
1504
1505//  "module:";
1506
1507  module D = x*gen(1), y*gen(2);
1508  D;
1509  D = setModuleGrading(D, B);
1510  print( mDeg( D ) );
1511
1512
1513  module DD = [0, 0],[0, 0, 0];
1514  DD = setModuleGrading(DD, B);
1515  print( mDeg( DD ) );
1516
1517  module DDD = [0, 0];
1518  DDD = setModuleGrading(DDD, B);
1519  print( mDeg( DDD ) );
1520
1521};
1522
1523
1524
1525
1526
1527/******************************************************/
1528proc mDegPartition(def p)
1529"USAGE: mDegPartition(def p), p polynomial/vector
1530RETURNS: an ideal/module consisting of multigraded-homogeneous parts of p
1531EXAMPLE: example mDegPartition; shows an example
1532"
1533{
1534  if( typeof(p) == "poly" )
1535  {
1536    ideal I;
1537    poly mp, t, tt;
1538  }
1539  else
1540  {
1541    if(  typeof(p) == "vector" )
1542    {
1543      module I;
1544      vector mp, t, tt;
1545    }
1546    else
1547    {
1548      ERROR("Wrong ARGUMENT type!");
1549    }
1550  }
1551
1552  if( typeof(p) == "vector" )
1553  {
1554    intmat V = getModuleGrading(p);
1555  }
1556  else
1557  {
1558    intmat V;
1559  }
1560
1561  if( size(p) > 1)
1562  {
1563    intvec m;
1564
1565    while( p != 0 )
1566    {
1567      m = leadexp(p);
1568      mp = lead(p);
1569      p = p - lead(p);
1570      tt = p; t = 0;
1571
1572      while( size(tt) > 0 )
1573      {
1574        // TODO: we do not cache matrices (M,T,H,V), which remain the same :(
1575        // TODO: we need some low-level procedure with all these arguments...!
1576        if( equalMDeg( leadexp(tt), m, V  ) )
1577        {
1578          mp = mp + lead(tt); // "mp", mp;
1579        }
1580        else
1581        {
1582          t = t + lead(tt);  //  "t", t;
1583        }
1584
1585        tt = tt - lead(tt);
1586      }
1587
1588      I[size(I)+1] = mp;
1589
1590      p = t;
1591    }
1592  }
1593  else
1594  {
1595    I[1] = p; // single monom
1596  }
1597
1598  if( typeof(I) == "module" )
1599  {
1600    I = setModuleGrading(I, V);
1601  }
1602
1603  return (I);
1604}
1605example
1606{
1607  "EXAMPLE:"; echo=2;
1608
1609  ring r = 0,(x,y,z),dp;
1610
1611  intmat g[2][3]=
1612    1,0,1,
1613    0,1,1;
1614  intmat t[2][1]=
1615    -2,
1616    1;
1617
1618  setBaseMultigrading(g,t);
1619
1620  poly f = x10yz+x8y2z-x4z2+y5+x2y2-z2+x17z3-y6;
1621
1622  mDegPartition(f);
1623
1624  vector v = xy*gen(1)-x3y2*gen(2)+x4y*gen(3);
1625  intmat B[2][3]=1,-1,-2,0,0,1;
1626  v = setModuleGrading(v,B);
1627  getModuleGrading(v);
1628
1629  mDegPartition(v, B);
1630}
1631
1632
1633
1634/******************************************************/
1635static proc unitMatrix(int n)
1636{
1637  intmat A[n][n];
1638
1639  for( int i = n; i > 0; i-- )
1640  {
1641    A[i,i] = 1;
1642  }
1643
1644  return (A);
1645}
1646
1647
1648
1649/******************************************************/
1650static proc finestMDeg(def r)
1651"
1652USAGE: finestMDeg(r); ring r
1653RETURN: ring, r endowed with the finest multigrading
1654TODO: not yet...
1655"
1656{
1657  def save = basering;
1658  setring (r);
1659
1660  // in basering
1661      ideal I = ideal(basering);
1662
1663  int n = 0; int i; poly p;
1664  for( i = ncols(I); i > 0; i-- )
1665  {
1666    p = I[i];
1667    if( size(p) > 1 )
1668    {
1669      n = n + (size(p) - 1);
1670    }
1671    else
1672    {
1673      I[i] = 0;
1674    }
1675  }
1676
1677  int N = nvars(basering);
1678  intmat A = unitMatrix(N);
1679
1680
1681
1682  if( n > 0)
1683  {
1684
1685    intmat L[N][n];
1686    //  list L;
1687    int j = n;
1688
1689    for(  i = ncols(I); i > 0; i-- )
1690    {
1691      p = I[i];
1692
1693      if( size(p) > 1 )
1694      {
1695        intvec m0 = leadexp(p);
1696        p = p - lead(p);
1697
1698        while( size(p) > 0 )
1699        {
1700          L[ 1..N, j ] = leadexp(p) - m0;
1701          p = p - lead(p);
1702          j--;
1703        }
1704      }
1705    }
1706
1707    print(L);
1708    setBaseMultigrading(A, L);
1709  }
1710  else
1711  {
1712    setBaseMultigrading(A);
1713  }
1714
1715  //  ERROR("nope");
1716
1717  //  ring T = integer, (x), (C, dp);
1718
1719  setring(save);
1720  return (r);
1721}
1722example
1723{
1724  "EXAMPLE:"; echo=2;
1725
1726  ring r = 0,(x, y), dp;
1727  qring q  = std(x^2 - y);
1728
1729  finestMDeg(q);
1730
1731}
1732
1733
1734
1735
1736/******************************************************/
1737static proc newMap(map F, intmat Q, list #)
1738"
1739USAGE: newMap(F, Q[, P]); map F, intmat Q[, intmat P]
1740PURPOSE: endowe the map F with the integer matrices P [and Q]
1741"
1742{
1743  attrib(F, "Q", Q);
1744
1745  if( size(#) > 0 and typeof(#[1]) == "intmat" )
1746  {
1747    attrib(F, "P", #[1]);
1748  }
1749  return (F);
1750}
1751
1752/******************************************************/
1753static proc matrix2intmat( matrix M )
1754{
1755  execute( "intmat A[ "+ string(nrows(M)) + "]["+ string(ncols(M)) + "] = " + string(M) + ";" );
1756  return (A);
1757}
1758
1759
1760/******************************************************/
1761static proc leftKernelZ(intmat M)
1762"USAGE:  leftKernel(M);   M a matrix
1763RETURN:  module
1764PURPOSE: computes left kernel of matrix M (a module of all elements v such that vM=0)
1765EXAMPLE: example leftKernel; shows an example
1766"
1767{
1768  if( nameof(basering) != "basering" )
1769  {
1770    def save = basering;
1771  }
1772
1773  ring r = integer, (x), dp;
1774
1775
1776  //  basering;
1777  module N = matrix((M)); // transpose
1778  //  print(N);
1779
1780  def MM = modulo( N, std(0) ) ;
1781  //  print(MM);
1782
1783  intmat R = (  matrix2intmat( MM ) ); // transpose
1784
1785  if( defined(save) > 0 )
1786  {
1787    setring save;
1788  }
1789
1790  kill r;
1791  return( R );
1792}
1793example
1794{
1795  "EXAMPLE:"; echo=2;
1796
1797  ring r= 0,(x,y,z),dp;
1798  matrix M[3][1] = x,y,z;
1799  print(M);
1800  matrix L = leftKernel(M);
1801  print(L);
1802  // check:
1803  print(L*M);
1804};
1805
1806
1807
1808/******************************************************/
1809// the following is taken from "sing4ti2.lib" as we need 'hilbert' from 4ti2
1810
1811static proc hilbert4ti2intmat(intmat A, list #)
1812"USAGE:  hilbert4ti2(A[,i]);
1813@*       A=intmat
1814@*       i=int
1815ASSUME:  - A is a matrix with integer entries which describes the lattice
1816@*         as ker(A), if second argument is not present, and
1817@*         as the left image Im(A) = {zA : z \in ZZ^k}, if second argument is a positive integer
1818@*       - number of variables of basering equals number of columns of A
1819@*         (for ker(A)) resp. of rows of A (for Im(A))
1820CREATE:  temporary files sing4ti2.mat, sing4ti2.lat, sing4ti2.mar
1821@*       in the current directory (I/O files for communication with 4ti2)
1822NOTE:    input rules for 4ti2 also apply to input to this procedure
1823@*       hence ker(A)={x|Ax=0} and Im(A)={xA}
1824RETURN:  toric ideal specified by Hilbert basis thereof
1825EXAMPLE: example graver4ti2; shows an example
1826"
1827{
1828//--------------------------------------------------------------------------
1829// Initialization and Sanity Checks
1830//--------------------------------------------------------------------------
1831   int i,j;
1832   int nr=nrows(A);
1833   int nc=ncols(A);
1834   string fileending="mat";
1835   if (size(#)!=0)
1836   {
1837//--- default behaviour: use ker(A) as lattice
1838//--- if #[1]!=0 use Im(A) as lattice
1839      if(typeof(#[1])!="int")
1840      {
1841         ERROR("optional parameter needs to be integer value");
1842      }
1843      if(#[1]!=0)
1844      {
1845         fileending="lat";
1846      }
1847   }
1848//--- we should also be checking whether all entries are indeed integers
1849//--- or whether there are fractions, but in this case the error message
1850//--- of 4ti2 is printed directly
1851
1852//--------------------------------------------------------------------------
1853// preparing input file for 4ti2
1854//--------------------------------------------------------------------------
1855   link eing=":w sing4ti2."+fileending;
1856   string eingstring=string(nr)+" "+string(nc);
1857   write(eing,eingstring);
1858   for(i=1;i<=nr;i++)
1859   {
1860      kill eingstring;
1861      string eingstring;
1862      for(j=1;j<=nc;j++)
1863      {
1864        //          if(g(A[i,j])>0)||(char(basering)!=0)||(npars(basering)>0))
1865        //          {
1866        //             ERROR("Input to hilbert4ti2 needs to be a matrix with integer entries");
1867        //          }
1868        eingstring=eingstring+string(A[i,j])+" ";
1869      }
1870      write(eing, eingstring);
1871   }
1872   close(eing);
1873
1874//----------------------------------------------------------------------
1875// calling 4ti2 and converting output
1876// Singular's string is too clumsy for this, hence we first prepare
1877// using standard unix commands
1878//----------------------------------------------------------------------
1879   j=system("sh","hilbert -q -n sing4ti2"); ////////// be quiet + no loggin!!!
1880
1881   j=system("sh", "awk \'BEGIN{ORS=\",\";}{print $0;}\' sing4ti2.hil " +
1882                "| sed s/[\\\ \\\t\\\v\\\f]/,/g " +
1883                "| sed s/,+/,/g|sed s/,,/,/g " +
1884                "| sed s/,,/,/g " +
1885                "> sing4ti2.converted" );
1886   if( defined(keepfiles) <= 0)
1887   {
1888      j=system("sh",("rm -f sing4ti2.hil sing4ti2."+fileending));
1889   }
1890//----------------------------------------------------------------------
1891// reading output of 4ti2
1892//----------------------------------------------------------------------
1893   link ausg=":r sing4ti2.converted";
1894//--- last entry ideal(0) is used to tie the list to the basering
1895//--- it will not be processed any further
1896
1897   string s = read(ausg);
1898   string ergstr = "intvec erglist = " + s + "0;";
1899   execute(ergstr);
1900
1901   //   print(erglist);
1902
1903   int Rnc = erglist[1];
1904   int Rnr = erglist[2];
1905
1906   intmat R[Rnr][Rnc];
1907
1908   int k = 3;
1909
1910   for(i=1;i<=Rnc;i++)
1911   {
1912     for(j=1;j<=Rnr;j++)
1913     {
1914       //       "i: ", i, ", j: ", j, ", v: ", erglist[k];
1915       R[j, i] = erglist[k];
1916       k = k + 1;
1917     }
1918   }
1919
1920   return (R);
1921//--- get rid of leading entry 0;
1922//   toric=toric[2..ncols(toric)];
1923//   return(toric);
1924}
1925// A nice example here is the 3x3 Magic Squares
1926example
1927{
1928  "EXAMPLE:"; echo=2;
1929
1930   ring r=0,(x1,x2,x3,x4,x5,x6,x7,x8,x9),dp;
1931   intmat M[7][9]=
1932      1, 1, 1, -1, -1, -1, 0, 0, 0,
1933      1, 1, 1,  0,  0,  0,-1,-1,-1,
1934      0, 1, 1, -1,  0,  0,-1, 0, 0,
1935      1, 0, 1,  0, -1,  0, 0,-1, 0,
1936      1, 1, 0,  0,  0, -1, 0, 0,-1,
1937      0, 1, 1,  0, -1,  0, 0, 0,-1,
1938      1, 1, 0,  0, -1,  0,-1, 0, 0;
1939   hilbert4ti2intmat(M);
1940   hermite(M);
1941}
1942
1943/////////////////////////////////////////////////////////////////////////////
1944static proc getMonomByExponent(intvec exp)
1945{
1946  int n = nvars(basering);
1947
1948  if( nrows(exp) < n )
1949  {
1950    n = nrows(exp);
1951  }
1952
1953  poly m = 1; int e;
1954
1955  for( int i = 1; i <= n; i++ )
1956  {
1957    e = exp[i];
1958    if( e < 0 )
1959    {
1960      ERROR("Negative exponent!!!");
1961    }
1962
1963    m = m * (var(i)^e);
1964  }
1965
1966  return (m);
1967
1968}
1969
1970/******************************************************/
1971proc mDegBasis(intvec d)
1972"
1973USAGE: multidegree d
1974ASSUME: current ring is multigraded, monomial ordering is global
1975PURPOSE: compute all monomials of multidegree d
1976EXAMPLE: example mDegBasis; shows an example
1977"
1978{
1979  def R = basering;  // setring R;
1980
1981  intmat M = getVariableWeights(R);
1982
1983  //  print(M);
1984
1985  int nr = nrows(M);
1986  int nc = ncols(M);
1987
1988  intmat A[nr][nc+1];
1989  A[1..nr, 1..nc] = M[1..nr, 1..nc];
1990  //typeof(A[1..nr, nc+1]);
1991  if( nr==1)
1992  {
1993    A[1..nr, nc+1]=-d[1];
1994  }
1995  else
1996  {
1997    A[1..nr, nc+1] = -d;
1998  }
1999
2000  intmat T = getTorsion(R);
2001
2002  if( isFreeRepresented() )
2003  {
2004    intmat B = hilbert4ti2intmat(A);
2005
2006    //      matrix B = unitMatrix(nrows(T));
2007  }
2008  else
2009  {
2010    int n = ncols(T);
2011
2012    nc = ncols(A);
2013
2014    intmat AA[nr][nc + 2 * n];
2015    AA[1..nr, 1.. nc] = A[1..nr, 1.. nc];
2016    AA[1..nr, nc + (1.. n)] = T[1..nr, 1.. n];
2017    AA[1..nr, nc + n + (1.. n)] = -T[1..nr, 1.. n];
2018
2019
2020    //      print ( AA );
2021
2022    intmat K = leftKernelZ(( AA ) ); //
2023
2024    //      print(K);
2025
2026    intmat KK[nc][ncols(K)] = K[ 1.. nc, 1.. ncols(K) ];
2027
2028    //      print(KK);
2029    //      "!";
2030
2031    intmat B = hilbert4ti2intmat(transpose(KK), 1);
2032
2033    //      "!";      print(B);
2034
2035  }
2036
2037
2038  //  print(A);
2039
2040
2041
2042  int i;
2043  int nnr = nrows(B);
2044  int nnc = ncols(B);
2045  ideal I, J;
2046  if(nnc==0){
2047    I=0;
2048    return(I);
2049  }
2050  I[nnc] = 0;
2051  J[nnc] = 0;
2052
2053  for( i = 1; i <= nnc; i++ )
2054  {
2055    //      "i: ", i;    B[nnr, i];
2056
2057    if( B[nnr, i] == 1)
2058    {
2059      // intvec(B[1..nnr-1, i]);
2060      I[i] = getMonomByExponent(intvec(B[1..nnr-1, i]));
2061    }
2062    else
2063    {
2064      if( B[nnr, i] == 0)
2065      {
2066        // intvec(B[1..nnr-1, i]);
2067        J[i] = getMonomByExponent(intvec(B[1..nnr-1, i]));
2068      }
2069    }
2070    //      I[i];
2071  }
2072
2073  ideal Q = (ideal(basering));
2074
2075  if ( size(Q) > 0 )
2076  {
2077    I = NF( I, lead(Q) );
2078    J = NF( J, lead(Q) ); // Global ordering!!!
2079  }
2080
2081  I = simplify(I, 2); // d
2082  J = simplify(J, 2); // d
2083
2084  attrib(I, "ZeroPart", J);
2085
2086  return (I);
2087
2088  //  setring ;
2089}
2090example
2091{
2092  "EXAMPLE:"; echo=2;
2093
2094  ring R = 0, (x, y), dp;
2095
2096  intmat g1[2][2]=1,0,0,1;
2097  intmat t[2][1]=2,0;
2098  intmat g2[2][2]=1,1,1,1;
2099  intvec v1=4,0;
2100  intvec v2=4,4;
2101
2102  intmat g3[1][2]=1,1;
2103  setBaseMultigrading(g3);
2104  intvec v3=4:1;
2105  v3;
2106  mDegBasis(v3);
2107
2108  setBaseMultigrading(g1,t);
2109  mDegBasis(v1);
2110  setBaseMultigrading(g2);
2111  mDegBasis(v2);
2112
2113  intmat M[2][2] = 1, -1, -1, 1;
2114  intvec d = -2, 2;
2115
2116  setBaseMultigrading(M);
2117
2118  mDegBasis(d);
2119  attrib(_, "ZeroPart");
2120
2121  kill R;
2122  ring R = 0, (x, y, z), dp;
2123
2124  intmat M[2][3] = 1, -2, 1,     1, 1, 0;
2125
2126  intmat T[2][1] = 0, 2;
2127
2128  intvec d = 4, 1;
2129
2130  setBaseMultigrading(M, T);
2131
2132  mDegBasis(d);
2133  attrib(_, "ZeroPart");
2134
2135
2136  kill R;
2137
2138  ring R = 0, (x, y, z), dp;
2139  qring Q = std(ideal( y^6+ x*y^3*z-x^2*z^2 ));
2140
2141
2142  intmat M[2][3] = 1, 1, 2,     2, 1, 1;
2143  //  intmat T[2][1] = 0, 2;
2144
2145  setBaseMultigrading(M);
2146
2147  intvec d = 6, 6;
2148  mDegBasis(d);
2149  attrib(_, "ZeroPart");
2150
2151
2152
2153  kill R;
2154  ring R = 0, (x, y, z), dp;
2155  qring Q = std(ideal( x*z^3 - y *z^6, x*y*z  - x^4*y^2 ));
2156
2157
2158  intmat M[2][3] = 1, -2, 1,     1, 1, 0;
2159  intmat T[2][1] = 0, 2;
2160
2161  intvec d = 4, 1;
2162
2163  setBaseMultigrading(M, T);
2164
2165  mDegBasis(d);
2166  attrib(_, "ZeroPart");
2167}
2168
2169
2170proc mDegSyzygy(def I)
2171"USAGE: mDegSyzygy(I); I is a poly/vector/ideal/module
2172PURPOSE: computes the multigraded syzygy module of I
2173RETURNS: module, the syzygy module of I
2174NOTE: generators of I must be multigraded homogeneous
2175EXAMPLE: example mDegSyzygy; shows an example
2176"
2177{
2178  if( isHomogenous(I, "checkGens") == 0)
2179  {
2180    ERROR ("Sorry: inhomogenous input!");
2181  }
2182  module S = syz(I);
2183  S = setModuleGrading(S, mDeg(I));
2184  return (S);
2185}
2186example
2187{
2188  "EXAMPLE:"; echo=2;
2189
2190
2191  ring r = 0,(x,y,z,w),dp;
2192  intmat MM[2][4]=
2193    1,1,1,1,
2194    0,1,3,4;
2195  setBaseMultigrading(MM);
2196  module M = ideal(  xw-yz, x2z-y3, xz2-y2w, yw2-z3);
2197
2198
2199  intmat v[2][nrows(M)]=
2200    1,
2201    0;
2202
2203  M = setModuleGrading(M, v);
2204
2205  isHomogenous(M);
2206  "Multidegrees: "; print(mDeg(M));
2207  // Let's compute syzygies!
2208  def S = mDegSyzygy(M); S;
2209  "Module Units Multigrading: "; print( getModuleGrading(S) );
2210  "Multidegrees: "; print(mDeg(S));
2211
2212  isHomogenous(S);
2213}
2214
2215
2216proc mDegGroebner(def I)
2217"USAGE: mDegGroebner(I); I is a poly/vector/ideal/module
2218PURPOSE: computes the multigraded standard/groebner basis of I
2219NOTE: I must be multigraded homogeneous
2220RETURNS: ideal/module, the computed basis
2221EXAMPLE: example mDegGroebner; shows an example
2222"
2223{
2224  if( isHomogenous(I) == 0)
2225  {
2226    ERROR ("Sorry: inhomogenous input!");
2227  }
2228
2229  def S = groebner(I);
2230
2231  if( typeof(I) == "module" or typeof(I) == "vector" )
2232  {
2233    S = setModuleGrading(S, getModuleGrading(I));
2234  }
2235
2236  return(S);
2237}
2238example
2239{
2240  "EXAMPLE:"; echo=2;
2241
2242  ring r = 0,(x,y,z,w),dp;
2243
2244  intmat MM[2][4]=
2245    1,1,1,1,
2246    0,1,3,4;
2247
2248  setBaseMultigrading(MM);
2249
2250
2251  module M = ideal(  xw-yz, x2z-y3, xz2-y2w, yw2-z3);
2252
2253
2254  intmat v[2][nrows(M)]=
2255    1,
2256    0;
2257
2258  M = setModuleGrading(M, v);
2259
2260
2261  /////////////////////////////////////////////////////////////////////////////
2262  // GB:
2263  M = mDegGroebner(M); M;
2264  "Module Units Multigrading: "; print( getModuleGrading(M) );
2265  "Multidegrees: "; print(mDeg(M));
2266
2267  isHomogenous(M);
2268
2269  /////////////////////////////////////////////////////////////////////////////
2270  // Let's compute Syzygy!
2271  def S = mDegSyzygy(M); S;
2272  "Module Units Multigrading: "; print( getModuleGrading(S) );
2273  "Multidegrees: "; print(mDeg(S));
2274
2275  isHomogenous(S);
2276
2277  /////////////////////////////////////////////////////////////////////////////
2278  // GB:
2279  S = mDegGroebner(S); S;
2280  "Module Units Multigrading: "; print( getModuleGrading(S) );
2281  "Multidegrees: "; print(mDeg(S));
2282
2283  isHomogenous(S);
2284}
2285
2286
2287/******************************************************/
2288proc mDegResolution(def I, int ll, list #)
2289"USAGE: mDegResolution(I,l,[f]); I is poly/vector/ideal/module; l,f are integers
2290PURPOSE: computes the multigraded resolution of I of the length l,
2291or the whole resolution if l is zero. Returns minimal resolution if an optional
2292argument 1 is supplied
2293NOTE: input must have multigraded-homogeneous generators.
2294The returned list is truncated beginning with the first zero differential.
2295RETURNS: list, the computed resolution
2296EXAMPLE: example mDegResolution; shows an example
2297"
2298{
2299  if( isHomogenous(I, "checkGens") == 0)
2300  {
2301    ERROR ("Sorry: inhomogenous input!");
2302  }
2303
2304  def R = res(I, ll, #); list L = R; int l = size(L);
2305
2306  if( (typeof(I) == "module") or (typeof(I) == "vector") )
2307  {
2308    L[1] = setModuleGrading(L[1], getModuleGrading(I));
2309  }
2310
2311  int i;
2312  for( i = 2; i <= l; i++ )
2313  {
2314    if( size(L[i]) > 0 )
2315    {
2316      L[i] = setModuleGrading( L[i], mDeg(L[i-1]) );
2317    } else
2318    {
2319      return (L[1..(i-1)]);
2320    }
2321  }
2322
2323  return (L);
2324
2325
2326}
2327example
2328{
2329  "EXAMPLE:"; echo=2;
2330
2331  ring r = 0,(x,y,z,w),dp;
2332
2333  intmat M[2][4]=
2334    1,1,1,1,
2335    0,1,3,4;
2336
2337  setBaseMultigrading(M);
2338
2339
2340  module m= ideal(  xw-yz, x2z-y3, xz2-y2w, yw2-z3);
2341
2342  isHomogenous(ideal(  xw-yz, x2z-y3, xz2-y2w, yw2-z3), "checkGens");
2343
2344  ideal A = xw-yz, x2z-y3, xz2-y2w, yw2-z3;
2345
2346  int j;
2347
2348  for(j=1; j<=ncols(A); j++)
2349  {
2350    mDegPartition(A[j]);
2351  }
2352
2353  intmat v[2][1]=
2354    1,
2355    0;
2356
2357  m = setModuleGrading(m, v);
2358
2359  // Let's compute Syzygy!
2360  def S = mDegSyzygy(m); S;
2361  "Module Units Multigrading: "; print( getModuleGrading(S) );
2362  "Multidegrees: "; print(mDeg(S));
2363
2364  /////////////////////////////////////////////////////////////////////////////
2365
2366  S = mDegGroebner(S); S;
2367  "Module Units Multigrading: "; print( getModuleGrading(S) );
2368  "Multidegrees: "; print(mDeg(S));
2369
2370  /////////////////////////////////////////////////////////////////////////////
2371
2372  def L = mDegResolution(m, 0, 1);
2373
2374  for( j =1; j<=size(L); j++)
2375  {
2376    "----------------------------------- ", j, " -----------------------------";
2377    L[j];
2378    "Module Multigrading: "; print( getModuleGrading(L[j]) );
2379    "Multigrading: "; print(mDeg(L[j]));
2380  }
2381
2382  /////////////////////////////////////////////////////////////////////////////
2383
2384  L = mDegResolution(maxideal(1), 0, 1);
2385
2386  for( j =1; j<=size(L); j++)
2387  {
2388    "----------------------------------- ", j, " -----------------------------";
2389    L[j];
2390    "Module Multigrading: "; print( getModuleGrading(L[j]) );
2391    "Multigrading: "; print(mDeg(L[j]));
2392  }
2393
2394  kill v;
2395
2396
2397  def h = hilbertSeries(m);
2398  setring h;
2399
2400  numerator1;
2401  factorize(numerator1);
2402
2403  denominator1;
2404  factorize(denominator1);
2405
2406  numerator2;
2407  factorize(numerator2);
2408
2409  denominator2;
2410  factorize(denominator2);
2411}
2412
2413/******************************************************/
2414proc hilbertSeries(def I)
2415"USAGE: hilbertSeries(I); I is poly/vector/ideal/module
2416PURPOSE: computes the multigraded Hilbert Series of M
2417NOTE: input must have multigraded-homogeneous generators.
2418Multigrading should be positive.
2419RETURNS: a ring in variables t_(i), s_(i), with polynomials
2420numerator1 and denominator1 and muturally prime numerator2
2421and denominator2, quotients of which give the series.
2422EXAMPLE: example hilbertSeries; shows an example
2423"
2424{
2425
2426  if( !isFreeRepresented() )
2427  {
2428    ERROR("SORRY: ONLY TORSION-FREE CASE (POSITIVE GRADING)");
2429  }
2430
2431  int i, j, k, v;
2432
2433  intmat M = getVariableWeights();
2434
2435  int cc = ncols(M);
2436  int n = nrows(M);
2437
2438  if( n == 0 )
2439  {
2440    ERROR("Error: wrong Variable Weights?");
2441  }
2442
2443  list RES = mDegResolution(I,0,1);
2444
2445  int l = size(RES);
2446
2447  list L; L[l + 1] = 0;
2448
2449  if(typeof(I) == "ideal")
2450  {
2451    intmat zeros[n][1];
2452    L[1] = zeros;
2453  }
2454  else
2455  {
2456    L[1] = getModuleGrading(RES[1]);
2457  }
2458
2459  for( j = 1; j <= l; j++)
2460  {
2461    L[j + 1] = mDeg(RES[j]);
2462  }
2463
2464  l++;
2465
2466  ring R = 0,(t_(1..n),s_(1..n)),dp;
2467
2468  ideal units;
2469  for( i=n; i>=1; i--)
2470  {
2471    units[i] = (var(i) * var(n + i) - 1);
2472  }
2473
2474  qring Q = std(units);
2475
2476  // TODO: should not it be a quotient ring depending on Torsion???
2477  // I am not sure about what to do in the torsion case, but since
2478  // we want to evaluate the polynomial at certain points to get
2479  // a dimension we need uniqueness for this. I think we would lose
2480  // this uniqueness if switching to this torsion ring.
2481
2482  poly monom, summand, numerator;
2483  poly denominator = 1;
2484
2485  for( i = 1; i <= cc; i++)
2486  {
2487    monom = 1;
2488
2489    for( k = 1; k <= n; k++)
2490    {
2491      v = M[k,i];
2492
2493      if(v >= 0)
2494      {
2495        monom = monom * (var(k)^(v));
2496      }
2497      else
2498      {
2499        monom = monom * (var(n+k)^(-v));
2500      }
2501    }
2502
2503    if( monom == 1)
2504    {
2505      ERROR("Multigrading not positive.");
2506    }
2507
2508    denominator = denominator * (1 - monom);
2509  }
2510
2511  for( j = 1; j<= l; j++)
2512  {
2513    summand = 0;
2514    M = L[j];
2515
2516    for( i = 1; i <= ncols(M); i++)
2517    {
2518      monom = 1;
2519      for( k = 1; k <= n; k++)
2520      {
2521        v = M[k,i];
2522        if( v > 0 )
2523        {
2524          monom = monom * (var(k)^v);
2525        }
2526        else
2527        {
2528          monom = monom * (var(n+k)^(-v));
2529        }
2530      }
2531      summand = summand + monom;
2532    }
2533    numerator = numerator - (-1)^j * summand;
2534  }
2535
2536  if( denominator == 0 )
2537  {
2538    ERROR("Multigrading not positive.");
2539  }
2540
2541  poly denominator1 = denominator;
2542  poly numerator1 = numerator;
2543
2544  export denominator1;
2545  export numerator1;
2546
2547  if( numerator != 0 )
2548  {
2549    poly d = gcd(denominator, numerator);
2550
2551    poly denominator2 = denominator/d;
2552    poly numerator2 = numerator/d;
2553
2554    if( gcd(denominator2, numerator2) != 1 )
2555    {
2556      ERROR("Sorry: gcd should be 1 (after dividing out gcd)! Something went wrong!");
2557    }
2558  }
2559  else
2560  {
2561    poly denominator2 = denominator;
2562    poly numerator2 = numerator;
2563  }
2564
2565
2566  export denominator2;
2567  export numerator2;
2568
2569  " ------------ ";
2570  "This proc returns a ring with polynomials called 'numerator1/2' and 'denominator1/2'!";
2571  "They represent the first and the second Hilbert Series.";
2572  "The s_(i)-variables are defined to be the inverse of the t_(i)-variables.";
2573  " ------------ ";
2574
2575  return(Q);
2576}
2577example
2578{
2579  "EXAMPLE:"; echo=2;
2580
2581  ring r = 0,(x,y,z,w),dp;
2582  intmat g[2][4]=
2583    1,1,1,1,
2584    0,1,3,4;
2585  setBaseMultigrading(g);
2586
2587  module M = ideal(xw-yz, x2z-y3, xz2-y2w, yw2-z3);
2588  intmat V[2][1]=
2589    1,
2590    0;
2591
2592  M = setModuleGrading(M, V);
2593
2594  def h = hilbertSeries(M); setring h;
2595
2596  factorize(numerator2);
2597  factorize(denominator2);
2598
2599  kill g, h; setring r;
2600
2601  intmat g[2][4]=
2602    1,2,3,4,
2603    0,0,5,8;
2604
2605  setBaseMultigrading(g);
2606
2607  ideal I = x^2, y, z^3;
2608  I = std(I);
2609  def L = mDegResolution(I, 0, 1);
2610
2611  for( int j = 1; j<=size(L); j++)
2612  {
2613    "----------------------------------- ", j, " -----------------------------";
2614    L[j];
2615    "Module Multigrading: "; print( getModuleGrading(L[j]) );
2616    "Multigrading: "; print(mDeg(L[j]));
2617  }
2618
2619  mDeg(I);
2620  def h = hilbertSeries(I); setring h;
2621
2622  factorize(numerator2);
2623  factorize(denominator2);
2624
2625  kill r, h, g, V;
2626  ////////////////////////////////////////////////
2627  ring R = 0,(x,y,z),dp;
2628  intmat W[2][3] =
2629     1,1, 1,
2630     0,0,-1;
2631  setBaseMultigrading(W);
2632  ideal I = x3y,yz2,y2z,z4;
2633
2634  def h = hilbertSeries(I); setring h;
2635
2636  factorize(numerator2);
2637  factorize(denominator2);
2638
2639  kill R, W, h;
2640  ////////////////////////////////////////////////
2641  ring R = 0,(x,y,z,a,b,c),dp;
2642  intmat W[2][6] =
2643     1,1, 1,1,1,1,
2644     0,0,-1,0,0,0;
2645  setBaseMultigrading(W);
2646  ideal I = x3y,yz2,y2z,z4;
2647
2648  def h = hilbertSeries(I); setring h;
2649
2650  factorize(numerator2);
2651  factorize(denominator2);
2652
2653  kill R, W, h;
2654  ////////////////////////////////////////////////
2655  // This is example 5.3.9. from Robbianos book.
2656
2657  ring R = 0,(x,y,z,w),dp;
2658  intmat W[1][4] =
2659     1,1, 1,1;
2660  setBaseMultigrading(W);
2661  ideal I = z3,y3zw2,x2y4w2xyz2;
2662
2663  hilb(std(I));
2664
2665  def h = hilbertSeries(I); setring h;
2666
2667  numerator1;
2668  denominator1;
2669
2670  factorize(numerator2);
2671  factorize(denominator2);
2672
2673
2674  kill h;
2675  ////////////////////////////////////////////////
2676  setring R;
2677
2678  ideal I2 = x2,y2,z2; I2;
2679
2680  hilb(std(I2));
2681
2682  def h = hilbertSeries(I2); setring h;
2683
2684  numerator1;
2685  denominator1;
2686
2687
2688  kill h;
2689  ////////////////////////////////////////////////
2690  setring R;
2691
2692  W = 2,2,2,2;
2693
2694  setBaseMultigrading(W);
2695
2696  getVariableWeights();
2697
2698  intvec w = 2,2,2,2;
2699
2700  hilb(std(I2), 1, w);
2701
2702  kill w;
2703
2704
2705  def h = hilbertSeries(I2); setring h;
2706
2707
2708  numerator1; denominator1;
2709  kill h;
2710
2711
2712  kill R, W;
2713
2714  ////////////////////////////////////////////////
2715  ring R = 0,(x),dp;
2716  intmat W[1][1] =
2717     1;
2718  setBaseMultigrading(W);
2719
2720  ideal I;
2721
2722  I = 1; I;
2723
2724  hilb(std(I));
2725
2726  def h = hilbertSeries(I); setring h;
2727
2728  numerator1; denominator1;
2729
2730  kill h;
2731  ////////////////////////////////////////////////
2732  setring R;
2733
2734  I = x; I;
2735
2736  hilb(std(I));
2737
2738  def h = hilbertSeries(I); setring h;
2739
2740  numerator1; denominator1;
2741
2742  kill h;
2743  ////////////////////////////////////////////////
2744  setring R;
2745
2746  I = x^5; I;
2747
2748  hilb(std(I));
2749  hilb(std(I), 1);
2750
2751  def h = hilbertSeries(I); setring h;
2752
2753  numerator1; denominator1;
2754
2755
2756  kill h;
2757  ////////////////////////////////////////////////
2758  setring R;
2759
2760  I = x^10; I;
2761
2762  hilb(std(I));
2763
2764  def h = hilbertSeries(I); setring h;
2765
2766  numerator1; denominator1;
2767
2768  kill h;
2769  ////////////////////////////////////////////////
2770  setring R;
2771
2772  module M = 1;
2773
2774  M = setModuleGrading(M, W);
2775
2776
2777  hilb(std(M));
2778
2779  def h = hilbertSeries(M); setring h;
2780
2781  numerator1; denominator1;
2782
2783  kill h;
2784  ////////////////////////////////////////////////
2785  setring R;
2786
2787  kill M; module M = x^5*gen(1);
2788//  intmat V[1][3] = 0; // TODO: this would lead to a wrong result!!!?
2789  intmat V[1][1] = 0; // all gen(i) of degree 0!
2790
2791  M = setModuleGrading(M, V);
2792
2793  hilb(std(M));
2794
2795  def h = hilbertSeries(M); setring h;
2796
2797  numerator1; denominator1;
2798
2799  kill h;
2800  ////////////////////////////////////////////////
2801  setring R;
2802
2803  module N = x^5*gen(3);
2804
2805  kill V;
2806
2807  intmat V[1][3] = 0; // all gen(i) of degree 0!
2808
2809  N = setModuleGrading(N, V);
2810
2811  hilb(std(N));
2812
2813  def h = hilbertSeries(N); setring h;
2814
2815  numerator1; denominator1;
2816
2817  kill h;
2818  ////////////////////////////////////////////////
2819  setring R;
2820
2821
2822  module S = M + N;
2823
2824  S = setModuleGrading(S, V);
2825
2826  hilb(std(S));
2827
2828  def h = hilbertSeries(S); setring h;
2829
2830  numerator1; denominator1;
2831
2832  kill h;
2833
2834  kill V;
2835  kill R, W;
2836
2837}
2838
2839
2840
2841///////////////////////////////////////////////////////////////////////////////
2842// testing for consistency of the library:
2843proc testMultigradingLib ()
2844{
2845  echo = 2; printlevel = 3;
2846  example setBaseMultigrading;
2847  example setModuleGrading;
2848
2849  example getVariableWeights;
2850  example getTorsion;
2851  example getModuleGrading;
2852
2853
2854  example mDeg;
2855  example mDegPartition;
2856
2857
2858  example hermite;
2859  example isHomogenous;
2860  example isTorsionFree;
2861  example pushForward;
2862  example defineHomogenous;
2863
2864  example equalMDeg;
2865  example isTorsionElement;
2866
2867  example mDegResolution;
2868
2869  "// ******************* example hilbertSeries ************************//";
2870  example hilbertSeries;
2871
2872
2873// example mDegBasis; // needs 4ti2!
2874
2875  "The End!";
2876
2877}
Note: See TracBrowser for help on using the repository browser.