source: git/Singular/LIB/multigrading.lib @ 6f84e21

spielwiese
Last change on this file since 6f84e21 was 6f84e21, checked in by Oleksandr Motsak <motsak@…>, 12 years ago
fix: "? wrong type declaration. type 'help poly;'" was due to making 'denominator' & 'numerator' - reserved words (kernel functions)
  • Property mode set to 100644
File size: 109.7 KB
Line 
1// -*- Mode: C++; tab-width: 2; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2// Vi-modeline: vim: filetype=c:syntax:shiftwidth=2:tabstop=8:textwidth=0:expandtab
3///////////////////////////////////////////////////////////////////
4version="$Id$";
5category="Combinatorial Commutative Algebra";
6info="
7LIBRARY:  multigrading.lib          Multigraded Rings
8
9AUTHORS:  Benjamin Bechtold, benjamin.bechtold@googlemail.com
10@*        Rene Birkner, rbirkner@math.fu-berlin.de
11@*        Lars Kastner, lkastner@math.fu-berlin.de
12@*        Simon Keicher, keicher@mail.mathematik.uni-tuebingen.de
13@*        Oleksandr Motsak, U@D, where U={motsak}, D={mathematik.uni-kl.de}
14@*        Anna-Lena Winz, anna-lena.winz@math.fu-berlin.de
15
16OVERVIEW: This library allows one to virtually add multigradings to Singular:
17grade multivariate polynomial rings with arbitrary (fin. gen. Abelian) groups. 
18For more see http://code.google.com/p/convex-singular/wiki/Multigrading
19For theoretical references see:
20@* E. Miller, B. Sturmfels: 'Combinatorial Commutative Algebra'
21and
22@* M. Kreuzer, L. Robbiano: 'Computational Commutative Algebra'.
23
24NOTE: 'multiDegBasis' relies on 4ti2 for computing Hilbert Bases.
25All groups are finitely generated Abelian
26
27PROCEDURES:
28setBaseMultigrading(M,L); attach multiweights/grading group matrices to the basering
29getVariableWeights([R]);  get matrix of multidegrees of vars attached to a ring
30
31getGradingGroup([R]);     get grading group attached to a ring
32getLattice([R[,choice]]); get grading group' lattice attached to a ring (or its NF)
33
34createGroup(S,L);          create a group generated by S, with relations L
35createQuotientGroup(L);    create a group generated by the unit matrix whith relations L
36createTorsionFreeGroup(S); create a group generated by S which is torsionfree
37printGroup(G);             print a group
38
39isGroup(G);                   test whether G is a valid group
40isGroupHomomorphism(L1,L2,A); test wheter A defines a group homomrphism from L1 to L2
41
42isGradedRingHomomorphism(R,f,A);  test graded ring homomorph
43createGradedRingHomomorphism(R,f,A);  create a graded ring homomorph
44
45setModuleGrading(M,v);    attach multiweights of units to a module and return it
46getModuleGrading(M);      get multiweights of module units (attached to M)
47
48isSublattice(A,B);        test whether A is a sublattice of B
49imageLattice(P,L);        computes an integral basis for P(L)
50intRank(A);               computes the rank of the intmat A
51kernelLattice(P);         computes an integral basis for the kernel of the linear map P.
52latticeBasis(B);          computes an integral basis of the lattice B
53preimageLattice(P,L);     computes an integral basis for the preimage of the lattice L under the linear map P.
54projectLattice(B);        computes a linear map of lattices having the primitive span of B as its kernel.
55intersectLattices(A,B);   computes an integral basis for the intersection of the lattices A and B.
56isIntegralSurjective(P);  test whether the map P of lattices is surjective.
57isPrimitiveSublattice(A); test whether A generates a primitive sublattice.
58intInverse(A);            computes the integral inverse matrix of the intmat A
59integralSection(P);       for a given linear surjective map P of lattices this procedure returns an integral section of P.
60primitiveSpan(A);         computes a basis for the minimal primitive sublattice that contains the given vectors (by A).
61
62factorgroup(G,H);         create the group G mod H
63productgroup(G,H);        create the group G x H
64
65multiDeg(A);                  compute the multidegree of A
66multiDegBasis(d);             compute all monomials of multidegree d
67multiDegPartition(p);         compute the multigraded-homogeneous components of p
68
69isTorsionFree();          test whether the current multigrading is  free
70isPositive();             test whether the current multigrading is positive
71isZeroElement(p);         test whether p has zero multidegree
72areZeroElements(M);       test whether an integer matrix M considered as a collection of columns has zero multidegree
73isHomogeneous(a);         test whether 'a' is multigraded-homogeneous
74
75equalMultiDeg(e1,e2[,V]);     test whether e1==e2 in the current multigrading
76
77multiDegGroebner(M);          compute the multigraded GB/SB of M
78multiDegSyzygy(M);            compute the multigraded syzygies of M
79multiDegModulo(I,J);          compute the multigraded 'modulo' module of I and J
80multiDegResolution(M,l[,m]);  compute the multigraded resolution of M
81multiDegTensor(m,n);          compute the tensor product of multigraded modules m,n
82multiDegTor(i,m,n);           compute the Tor_i(m,n) for multigraded modules m,n
83
84defineHomogeneous(p);     get a grading group wrt which p becomes homogeneous
85pushForward(f);           find the finest grading on the image ring, homogenizing f
86gradiator(h);             coarsens grading of the ring until h becomes homogeneous
87
88hermiteNormalForm(A);     compute the Hermite Normal Form of a matrix
89smithNormalForm(A,#);     compute matrices D,P,Q with D=P*A*Q and D is the smith normal form of A
90
91hilbertSeries(M);         compute the multigraded Hilbert Series of M
92
93lll(A);                   applies LLL(.) of lll.lib which only works for lists on a matrix A
94
95           (parameters in square brackets [] are optional)
96
97KEYWORDS:  multigrading, multidegree, multiweights, multigraded-homogeneous, integral linear algebra
98";
99
100/// evalHilbertSeries(h,v);   evaluate hilberts series h by substituting v[i] for t_(i) (too experimentall)
101
102// finestMDeg(def r)
103// newMap(map F, intmat Q, list #)
104
105LIB "standard.lib"; // for groebner
106// LIB "lll.lib"; // for lll_matrix // no need now, right?
107LIB "matrix.lib"; // for multiDegTor
108
109/******************************************************/
110
111static proc concatintmat(intmat A, intmat B)
112{
113
114 if ( nrows(A) != nrows(B) )
115   {
116     ERROR("matrices A and B have different number of rows.");
117   }
118
119 intmat At = transpose(A);
120 intmat Bt = transpose(B);
121
122 intmat Ct[nrows(At) + nrows(Bt)][ncols(At)] = At, Bt;
123
124 return(transpose(Ct));
125}
126
127
128/******************************************************/
129proc createGradedRingHomomorphism(def src, ideal Im, def A)
130"USAGE: createGradedRingHomomorphism(R, f, A); ring R, ideal f, group homomorphism A
131PURPOSE: create a multigraded group ring homomorphism defined by
132a ring map from R to the current ring, given by generators images f
133and a group homomorphism A between grading groups
134RETURN: graded ring homorphism
135EXAMPLE: example createGradedRingHomomorphism; shows an example
136"
137{
138  string isGRH = "isGRH";
139
140  if( !isGradedRingHomomorphism(src, Im, A) )
141  {
142    ERROR("Input data is wrong");
143  }
144
145  list h;
146  h[3] = A;
147
148//  map f = src, Im;
149  h[2] = Im; // f?
150  h[1] = src;
151
152  attrib(h, isGRH, (1==1)); // mark it "a graded ring homomorphism"
153
154  return(h);
155}
156example
157{
158  "EXAMPLE:"; echo=2;
159
160  ring r = 0, (x, y, z), dp;
161  intmat S1[3][3] =
162                   1, 0, 0,
163                   0, 1, 0,
164                   0, 0, 1;
165  intmat L1[3][1] =
166                   0,
167                   0,
168                   0;
169
170  def G1 = createGroup(S1, L1); // (S1 + L1)/L1
171  printGroup(G1);
172
173  setBaseMultigrading(S1, L1); // to change...
174
175  ring R = 0, (a, b, c), dp;
176  intmat S2[2][3] =
177                   1, 0,
178                   0, 1;
179  intmat L2[2][1] =
180                   0,
181                   2;
182
183  def G2 = createGroup(S2, L2);
184  printGroup(G2);
185
186  setBaseMultigrading(S2, L2); // to change...
187
188
189  map F = r, a, b, c;
190  intmat A[nrows(L2)][nrows(L1)] =
191                                  1, 0, 0,
192                                  3, 2, -6;
193
194  // graded ring homomorphism is given by (compatible):
195  print(F);
196  print(A);
197
198  isGradedRingHomomorphism(r, ideal(F), A);
199  def h = createGradedRingHomomorphism(r, ideal(F), A);
200
201  print(h);
202}
203
204
205/******************************************************/
206proc isGradedRingHomomorphism(def src, ideal Im, def A)
207"USAGE: isGradedRingHomomorphism(R, f, A); ring R, ideal f, group homomorphism A
208PURPOSE: test a multigraded group ring homomorphism defined by
209a ring map from R to the current ring, given by generators images f
210and a group homomorphism A between grading groups
211RETURN: int, 1 for TRUE, 0 otherwise
212EXAMPLE: example isGradedRingHomomorphism; shows an example
213"
214{
215  def dst = basering;
216
217  intmat result_degs = multiDeg(Im);
218//  print(result_degs);
219
220  setring src;
221
222  intmat input_degs = multiDeg(maxideal(1));
223//  print(input_degs);
224
225  def image_degs = A * input_degs;
226//  print( image_degs );
227
228  def df = image_degs - result_degs;
229//  print(df);
230
231  setring dst;
232
233  return (areZeroElements( df ));
234}
235example
236{
237  "EXAMPLE:"; echo=2;
238
239  ring r = 0, (x, y, z), dp;
240  intmat S1[3][3] =
241    1, 0, 0,
242    0, 1, 0,
243    0, 0, 1;
244  intmat L1[3][1] =
245    0,
246    0,
247    0;
248
249  def G1 = createGroup(S1, L1); // (S1 + L1)/L1
250  printGroup(G1);
251
252  setBaseMultigrading(S1, L1); // to change...
253
254  ring R = 0, (a, b, c), dp;
255  intmat S2[2][3] =
256    1, 0,
257    0, 1;
258  intmat L2[2][1] =
259    0,
260    2;
261
262  def G2 = createGroup(S2, L2);
263  printGroup(G2);
264
265  setBaseMultigrading(S2, L2); // to change...
266
267
268  map F = r, a, b, c;
269  intmat A[nrows(L2)][nrows(L1)] =
270      1, 0, 0,
271      3, 2, -6;
272
273  // graded ring homomorphism is given by (compatible):
274  print(F);
275  print(A);
276
277  isGradedRingHomomorphism(r, ideal(F), A);
278  def h = createGradedRingHomomorphism(r, ideal(F), A);
279
280  print(h);
281
282  // not a homo..
283  intmat B[nrows(L2)][nrows(L1)] =
284     1, 1, 1,
285     0, 0, 0;
286  print(B);
287
288  isGradedRingHomomorphism(r, ideal(F), B); // FALSE: there is no such homomorphism!
289  // Therefore: the following command should return an error
290  // createGradedRingHomomorphism(r, ideal(F), B);
291
292}
293
294
295proc createQuotientGroup(intmat L)
296"USAGE: createGroup(L); L is an integer matrix
297PURPOSE: create the group of the form (I+L)/L,
298where I is the square identity matrix of size nrows(L) x nrows(L)
299NOTE: L specifies relations between free generators of Z^nrows(L)
300RETURN: group
301EXAMPLE: example createQuotientGroup; shows an example
302"
303{
304  int r = nrows(L); int i;
305  intmat S[r][r]; // SQUARE!!!
306  for(i = r; i > 0; i--){ S[i, i] = 1; }
307  return (createGroup(S,L));
308}
309example
310{
311  "EXAMPLE:"; echo=2;
312
313  intmat I[3][3] =
314                  1, 0, 0,
315                  0, 1, 0,
316                  0, 0, 1;
317
318  intmat L[3][2] =
319                  1, 1,
320                  1, 3,
321                  1, 5;
322
323
324  // The group Z^3 / L can be constructed as follows:
325
326  // shortcut:
327  def G = createQuotientGroup(L);
328  printGroup(G);
329
330  // the general way:
331  def GG = createGroup(I, L); // (I+L)/L
332  printGroup(GG);
333}
334
335proc createTorsionFreeGroup(intmat S)
336"USAGE: createTorsionFreeGroup(S); S is an integer matrix
337PURPOSE: create the free subgroup generated by S within the
338free Abelian group of rank nrows(S)
339RETURN: group
340EXAMPLE: example createTorsionFreeGroup; shows an example
341"
342{
343  int r = nrows(S); int i;
344  intmat L[r][1] = 0;
345  return (createGroup(S,L));
346}
347example
348{
349  "EXAMPLE:"; echo=2;
350
351  // ----------- extreme case ------------ //
352  intmat S[1][3] =
353                  1,  -1, 10;
354
355  // Torsion:
356  intmat L[1][1] =
357                  0;
358
359
360  // The free subgroup generated by elements of S within Z^1
361  // can be constructed as follows:
362
363  // shortcut:
364  def G = createTorsionFreeGroup(S);
365  printGroup(G);
366
367  // the general way:
368  def GG = createGroup(S, L); // (S+L)/L
369  printGroup(GG);
370}
371
372
373/******************************************************/
374proc createGroup(intmat S, intmat L)
375"USAGE: createGroup(S, L); S, L are integer matrices
376PURPOSE: create the group of the form (S+L)/L, i.e.
377S specifies generators, L specifies relations.
378RETURN: group
379EXAMPLE: example createGroup; shows an example
380"
381{
382  string isGroup = "isGroup";
383  string attrGroupHNF = "hermite";
384  string attrGroupSNF = "smith";
385
386
387/*
388  if( size(#) > 0 )
389  {
390    if( typeof(#[1]) == "intmat" )
391    {
392      intmat S = #[1];
393    } else { ERROR("Wrong optional argument: 1"); }
394
395    if( size(#) > 1 )
396    {
397      if( typeof(#[2]) == "intmat" )
398      {
399        intmat L = #[2];
400      } else { ERROR("Wrong optional argument: 2"); }
401    }
402  }
403*/
404
405  if( nrows(L) != nrows(S) )
406  {
407    ERROR("Incompatible matrices!");
408  }
409
410  def H = attrib(L, attrGroupHNF);
411  if( typeof(H) != "intmat")
412  {
413    attrib(L, attrGroupHNF, hermiteNormalForm(L));
414  } else { kill H; }
415
416  def HH = attrib(L, attrGroupSNF);
417  if( typeof(HH) != "intmat")
418  {
419    attrib(L, attrGroupSNF, smithNormalForm(L));
420  } else { kill HH; }
421
422  list G; // Please, note the order: Generators + Relations:
423  G[1] = S;
424  G[2] = L;
425  // And now a quick-and-dirty fix of Singular inability to handle attribs of attribs:
426  // For the use of a group as an attribute for multigraded rings
427  G[3] = attrib(L, attrGroupHNF); 
428  G[4] = attrib(L, attrGroupSNF);
429 
430
431  attrib(G, isGroup, (1==1)); // mark it "a group"
432
433  return (G);
434}
435example
436{
437  "EXAMPLE:"; echo=2;
438
439  intmat S[3][3] =
440    1, 0, 0,
441    0, 1, 0,
442    0, 0, 1;
443
444  intmat L[3][2] =
445    1, 1,
446    1, 3,
447    1, 5;
448
449  def G = createGroup(S, L); // (S+L)/L
450
451  printGroup(G);
452
453  kill S, L, G;
454
455  /////////////////////////////////////////////////
456  intmat S[2][3] =
457    1, -2, 1,
458    1,  1, 0;
459
460  intmat L[2][1] =
461    0,
462    2;
463
464  def G = createGroup(S, L); // (S+L)/L
465
466  printGroup(G);
467
468  kill S, L, G;
469
470  // ----------- extreme case ------------ //
471  intmat S[1][3] =
472    1,  -1, 10;
473
474  // Torsion:
475  intmat L[1][1] =
476    0;
477
478  def G = createGroup(S, L); // (S+L)/L
479
480  printGroup(G);
481}
482
483
484/******************************************************/
485proc printGroup(def G)
486"USAGE: printGroup(G); G is a group
487PURPOSE: prints the group G
488RETURN: nothing
489EXAMPLE: example printGroup; shows an example
490"
491{
492  "Generators: ";
493  print(G[1]);
494
495  "Relations: ";
496  print(G[2]);
497
498//  attrib(G[2]);
499}
500example
501{
502  "EXAMPLE:"; echo=2;
503
504  intmat S[3][3] =
505                  1, 0, 0,
506                  0, 1, 0,
507                  0, 0, 1;
508
509  intmat L[3][2] =
510                  1, 1,
511                  1, 3,
512                  1, 5;
513
514  def G = createGroup(S, L); // (S+L)/L
515  printGroup(G);
516
517}
518
519/******************************************************/
520static proc areIsomorphicGroups(def G, def H)
521"USAGE: areIsomorphicGroups(G, H); G and H are groups
522PURPOSE: Check whether G and H define isomorphic groups.
523RETURN: int, 1 for TRUE, 0 otherwise
524EXAMPLE: example areIsomorphicGroups; shows an example
525"
526{
527  ERROR("areIsomorphicGroups: Not yet implemented!");
528  return (1); // TRUE
529}
530example
531{
532  "EXAMPLE:"; echo=2;
533
534  // TODO!
535
536}
537
538/******************************************************/
539proc isGroup(def G)
540"USAGE: isGroup(G); G a list
541PURPOSE: checks whether G is a valid group
542NOTE: G should be created by createGroup
543(or createQuotientGroup, createTorsionFreeGroup)
544RETURN: int, 1 if G is a valid group and 0 otherwise
545EXAMPLE: example isGroup; shows an example
546"
547{
548  string isGroup = "isGroup";
549
550  // valid?
551  if( typeof(G) != "list" ){ return(0); }
552
553  def a = attrib(G, isGroup);
554
555///// TODO for Hans: fix attr^2 bug in Singular!
556
557//  if( !defined(a) ) { return(0); }
558//  if( typeof(a) != "int" ) { return(0); }
559  if( defined(a) ){ if(typeof(a) == "int") { return(a); } }
560
561
562  if( (size(G) != 2) && (size(G) != 4) ){ return(0); }
563  if( typeof(G[1]) != "intmat" ){ return(0); }
564  if( typeof(G[2]) != "intmat" ){ return(0); }
565  if( nrows(G[1]) != nrows(G[2]) ){ return(0); }
566
567  return(1);
568}
569example
570{
571  "EXAMPLE:"; echo=2;
572
573  intmat S[3][3] =
574                  1, 0, 0,
575                  0, 1, 0,
576                  0, 0, 1;
577
578  intmat L[3][2] =
579                  1, 1,
580                  1, 3,
581                  1, 5;
582
583  def G = createGroup(S, L); // (S+L)/L
584
585  isGroup(G);
586 
587  printGroup(G);
588
589}
590
591
592
593/******************************************************/
594proc setBaseMultigrading(intmat M, list #)
595"USAGE: setBaseMultigrading(M[, G]); M is an integer matrix, G is a group (or lattice)
596PURPOSE: attaches weights of variables and grading group to the basering.
597NOTE: M encodes the weights of variables column-wise.
598RETURN: nothing
599EXAMPLE: example setBaseMultigrading; shows an example
600"
601{
602  string attrMgrad   = "mgrad";
603  string attrGradingGroup = "gradingGroup";
604
605  int i = 1;
606  if( size(#) >= i )
607  {
608    def a = #[i];
609    if( typeof(a) == "intmat" )
610    {
611      def L = createGroup(M, a);
612      i++;
613    }
614
615    if( isGroup(a) )
616    {
617      def L = a;
618
619      if( !isSublattice(M, L[1]) )
620      {
621        ERROR("Multigrading is not contained in the grading group!");
622      }
623      i++;
624    }
625    if( i == 1 ){ ERROR("Wrong arguments: no group given?"); }
626    kill a;
627  }
628  else
629  {
630    def L = createTorsionFreeGroup(M);
631  }
632
633
634  attrib(basering, attrMgrad, M);
635  attrib(basering, attrGradingGroup, L);
636
637  ideal Q = ideal(basering); // quotient ideal is assumed to be a GB!
638  if( !isHomogeneous(Q, "checkGens") ) // easy now, but would be hard before setting ring attributes!
639  {
640    "Warning: your quotient ideal is not homogenous (multigrading was set anyway)!";
641  }
642
643}
644example
645{
646  "EXAMPLE:"; echo=2;
647
648  ring R = 0, (x, y, z), dp;
649
650  // Weights of variables
651  intmat M[3][3] =
652    1, 0, 0,
653    0, 1, 0,
654    0, 0, 1;
655
656  // GradingGroup:
657  intmat L[3][2] =
658    1, 1,
659    1, 3,
660    1, 5;
661
662  // attaches M & L to R (==basering):
663  setBaseMultigrading(M, L); // Grading: Z^3/L
664
665  // Weights are accessible via "getVariableWeights()":
666  getVariableWeights();
667
668  // Test all possible usages:
669  (getVariableWeights() == M) && (getVariableWeights(R) == M) && (getVariableWeights(basering) == M);
670
671  // Grading group is accessible via "getLattice()":
672  getLattice();
673
674  // Test all possible usages:
675  (getLattice() == L) && (getLattice(R) == L) && (getLattice(basering) == L);
676
677  // And its hermite NF via getLattice("hermite"):
678  getLattice("hermite");
679
680  // Test all possible usages:
681  intmat H = hermiteNormalForm(L);
682  (getLattice("hermite") == H) && (getLattice(R, "hermite") == H) && (getLattice(basering, "hermite") == H);
683
684  kill L, M;
685
686  // ----------- isomorphic multigrading -------- //
687
688  // Weights of variables
689  intmat M[2][3] =
690    1, -2, 1,
691    1,  1, 0;
692
693  // Torsion:
694  intmat L[2][1] =
695    0,
696    2;
697
698  // attaches M & L to R (==basering):
699  setBaseMultigrading(M, L); // Grading: Z + (Z/2Z)
700
701  // Weights are accessible via "getVariableWeights()":
702  getVariableWeights() == M;
703
704  // Torsion is accessible via "getLattice()":
705  getLattice() == L;
706
707  kill L, M;
708  // ----------- extreme case ------------ //
709
710  // Weights of variables
711  intmat M[1][3] =
712    1,  -1, 10;
713
714  // Torsion:
715  intmat L[1][1] =
716    0;
717
718  // attaches M & L to R (==basering):
719  setBaseMultigrading(M); // Grading: Z^3
720
721  // Weights are accessible via "getVariableWeights()":
722  getVariableWeights() == M;
723
724  // Torsion is accessible via "getLattice()":
725  getLattice() == L;
726}
727
728
729/******************************************************/
730proc getVariableWeights(list #)
731"USAGE: getVariableWeights([R])
732PURPOSE: get associated multigrading matrix for the basering [or R]
733RETURN:  intmat, matrix of multidegrees of variables
734EXAMPLE: example getVariableWeights; shows an example
735"
736{
737  string attrMgrad = "mgrad";
738
739
740  if( size(#) > 0 )
741  {
742    if(( typeof(#[1]) == "ring" ) || ( typeof(#[1]) == "qring" ))
743    {
744      def R = #[1];
745    }
746    else
747    {
748      ERROR("Optional argument must be a ring!");
749    }
750  }
751  else
752  {
753    def R = basering;
754  }
755
756  def M = attrib(R, attrMgrad);
757  if( typeof(M) == "intmat"){ return (M); }
758  ERROR( "Sorry no multigrading matrix!" );
759}
760example
761{
762  "EXAMPLE:"; echo=2;
763
764  ring R = 0, (x, y, z), dp;
765
766  // Weights of variables
767  intmat M[3][3] =
768    1, 0, 0,
769    0, 1, 0,
770    0, 0, 1;
771
772  // Grading group:
773  intmat L[3][2] =
774    1, 1,
775    1, 3,
776    1, 5;
777
778  // attaches M & L to R (==basering):
779  setBaseMultigrading(M, L); // Grading: Z^3/L
780
781  // Weights are accessible via "getVariableWeights()":
782  getVariableWeights() == M;
783
784  kill L, M;
785
786  // ----------- isomorphic multigrading -------- //
787
788  // Weights of variables
789  intmat M[2][3] =
790    1, -2, 1,
791    1,  1, 0;
792
793  // Grading group:
794  intmat L[2][1] =
795    0,
796    2;
797
798  // attaches M & L to R (==basering):
799  setBaseMultigrading(M, L); // Grading: Z + (Z/2Z)
800
801  // Weights are accessible via "getVariableWeights()":
802  getVariableWeights() == M;
803
804  kill L, M;
805
806  // ----------- extreme case ------------ //
807
808  // Weights of variables
809  intmat M[1][3] =
810    1,  -1, 10;
811
812  // Grading group:
813  intmat L[1][1] =
814    0;
815
816  // attaches M & L to R (==basering):
817  setBaseMultigrading(M); // Grading: Z^3
818
819  // Weights are accessible via "getVariableWeights()":
820  getVariableWeights() == M;
821}
822
823
824proc getGradingGroup(list #)
825"USAGE: getGradingGroup([R])
826PURPOSE: get associated grading group
827RETURN: group, the grading group
828EXAMPLE: example getGradingGroup; shows an example
829"
830{
831  string attrGradingGroup    = "gradingGroup";
832
833  int i = 1;
834
835  if( size(#) >= i )
836  {
837    if( ( typeof(#[i]) == "ring" ) or ( typeof(#[i]) == "qring" ) )
838    {
839      def R = #[i];
840      i++;
841    }
842  }
843
844  if( i == 1 )
845  {
846    def R = basering;
847  }
848
849  def G = attrib(R, attrGradingGroup);
850
851  if( !isGroup(G) )
852  {
853    ERROR("Sorry no grading group!");
854  }
855
856  return(G);
857}
858example
859{
860  "EXAMPLE:"; echo=2;
861
862  ring R = 0, (x, y, z), dp;
863
864  // Weights of variables
865  intmat M[3][3] =
866    1, 0, 0,
867    0, 1, 0,
868    0, 0, 1;
869
870  // Torsion:
871  intmat L[3][2] =
872    1, 1,
873    1, 3,
874    1, 5;
875
876  // attaches M & L to R (==basering):
877  setBaseMultigrading(M, L); // Grading: Z^3/L
878
879  def G = getGradingGroup();
880
881  printGroup( G );
882
883  G[1] == M; G[2] == L;
884
885  kill L, M, G;
886
887  // ----------- isomorphic multigrading -------- //
888
889  // Weights of variables
890  intmat M[2][3] =
891    1, -2, 1,
892    1,  1, 0;
893
894  // Torsion:
895  intmat L[2][1] =
896    0,
897    2;
898
899  // attaches M & L to R (==basering):
900  setBaseMultigrading(M, L); // Grading: Z + (Z/2Z)
901
902  def G = getGradingGroup();
903
904  printGroup( G );
905
906  G[1] == M; G[2] == L;
907
908  kill L, M, G;
909  // ----------- extreme case ------------ //
910
911  // Weights of variables
912  intmat M[1][3] =
913    1,  -1, 10;
914
915  // Torsion:
916  intmat L[1][1] =
917    0;
918
919  // attaches M & L to R (==basering):
920  setBaseMultigrading(M); // Grading: Z^3
921
922  def G = getGradingGroup();
923
924  printGroup( G );
925
926  G[1] == M; G[2] == L;
927
928  kill L, M, G;
929}
930
931
932/******************************************************/
933proc getLattice(list #)
934"USAGE: getLattice([R[,opt]])
935PURPOSE: get associated grading group matrix, i.e. generators (cols) of the grading group
936RETURN: intmat, the grading group matrix, or
937its hermite normal form if an optional argument (\"hermiteNormalForm\") is given or
938smith normal form if an optional argument (\"smith\") is given
939EXAMPLE: example getLattice; shows an example
940"
941{
942  int i = 1;
943  if( size(#) >= i )
944  {
945    def a = #[i];
946    if( ( typeof(a) == "ring" ) or ( typeof(a) == "qring" ) )
947    {
948      i++;
949    }
950    kill a;
951  }
952
953  string attrGradingGroupHNF = "hermite";
954  string attrGradingGroupSNF = "smith";
955
956  def G = getGradingGroup(#);
957
958//  printGroup(G);
959
960
961
962  def T = G[2]; 
963
964  if( size(#) >= i )
965  {
966    def a = #[i];
967
968    if( typeof(a) != "string" )
969    {
970      ERROR("Sorry wrong arguments!");
971    }
972   
973    if( a == "hermite" )
974    {
975      def M = attrib(T, attrGradingGroupHNF);
976      if( typeof(M) != "intmat" )
977      {
978        if( size(G) > 2 )
979        {
980          M = G[3];
981        } else
982        {
983          M = hermiteNormalForm(T);
984        }
985      }
986      return (M);
987    }
988
989    if( a == "smith" )
990    {
991      def M = attrib(T, attrGradingGroupSNF);
992      if( typeof(M) != "intmat" )
993      {
994        if( size(G) > 2 )
995        {
996          M = G[4];
997        } else
998        {
999          M = smithNormalForm(T);
1000        }
1001      }
1002      return (M);
1003    }
1004  }
1005
1006  return(T);
1007}
1008example
1009{
1010  "EXAMPLE:"; echo=2;
1011
1012  ring R = 0, (x, y, z), dp;
1013
1014  // Weights of variables
1015  intmat M[3][3] =
1016    1, 0, 0,
1017    0, 1, 0,
1018    0, 0, 1;
1019
1020  // Torsion:
1021  intmat L[3][2] =
1022    1, 1,
1023    1, 3,
1024    1, 5;
1025
1026  // attaches M & L to R (==basering):
1027  setBaseMultigrading(M, L); // Grading: Z^3/L
1028
1029  // Torsion is accessible via "getLattice()":
1030  getLattice() == L;
1031
1032  // its hermite NF:
1033  print(getLattice("hermite"));
1034
1035  kill L, M;
1036
1037  // ----------- isomorphic multigrading -------- //
1038
1039  // Weights of variables
1040  intmat M[2][3] =
1041    1, -2, 1,
1042    1,  1, 0;
1043
1044  // Torsion:
1045  intmat L[2][1] =
1046    0,
1047    2;
1048
1049  // attaches M & L to R (==basering):
1050  setBaseMultigrading(M, L); // Grading: Z + (Z/2Z)
1051
1052  // Torsion is accessible via "getLattice()":
1053  getLattice() == L;
1054
1055  // its hermite NF:
1056  print(getLattice("hermite"));
1057
1058  kill L, M;
1059
1060  // ----------- extreme case ------------ //
1061
1062  // Weights of variables
1063  intmat M[1][3] =
1064    1,  -1, 10;
1065
1066  // Torsion:
1067  intmat L[1][1] =
1068    0;
1069
1070  // attaches M & L to R (==basering):
1071  setBaseMultigrading(M); // Grading: Z^3
1072
1073  // Torsion is accessible via "getLattice()":
1074  getLattice() == L;
1075
1076  // its hermite NF:
1077  print(getLattice("hermite"));
1078}
1079
1080proc getGradedGenerator(def m, int i)
1081"USAGE: getGradedGenerator(M, i), 'M' module/ideal, 'i' int
1082RETURN: returns the i-th generator of M, endowed with the module grading from M
1083EXAMPLE: example getGradedGenerator; shows an example
1084"
1085{
1086  if( typeof(m) == "ideal" )
1087  {
1088    return (m[i]);
1089  }
1090
1091  if( typeof(m) == "module" )
1092  {
1093    def v = getModuleGrading(m);
1094
1095    return ( setModuleGrading(m[i],v) );
1096  }
1097
1098  ERROR("m is expected to be an ideal or a module");
1099}
1100example
1101{
1102  "EXAMPLE:"; echo=2;
1103
1104  ring r = 0,(x,y,z,w),dp;
1105  intmat MM[2][4]=
1106                  1,1,1,1,
1107                  0,1,3,4;
1108  setBaseMultigrading(MM);
1109 
1110  module M = ideal(  xw-yz, x2z-y3, xz2-y2w, yw2-z3);
1111
1112
1113  intmat v[2][nrows(M)]=
1114                        1,
1115                        0;
1116
1117  M = setModuleGrading(M, v);
1118
1119  isHomogeneous(M);
1120  "Multidegrees: "; print(multiDeg(M));
1121
1122  // Let's compute syzygies!
1123  def S = multiDegSyzygy(M); S;
1124  "Module Units Multigrading: "; print( getModuleGrading(S) );
1125  "Multidegrees: "; print(multiDeg(S));
1126
1127  isHomogeneous(S);
1128
1129  // same as S[1] together with the induced module weighting
1130  def v = getGradedGenerator(S, 1);
1131  print(v);
1132  print(setModuleGrading(v));
1133 
1134  isHomogeneous(v);
1135
1136  isHomogeneous(S[1]);
1137}
1138
1139/******************************************************/
1140proc getModuleGrading(def m)
1141"USAGE: getModuleGrading(m), 'm' module/vector
1142RETURN: integer matrix of the multiweights of free module generators attached to 'm'
1143EXAMPLE: example getModuleGrading; shows an example
1144"
1145{
1146  string attrModuleGrading = "genWeights";
1147
1148  //  print(m);  typeof(m);  attrib(m);
1149
1150  def V = attrib(m, attrModuleGrading);
1151
1152  if( typeof(V) != "intmat" )
1153  {
1154    if( (typeof(m) == "ideal") or (typeof(m) == "poly") )
1155    {
1156      intmat M = getVariableWeights();
1157      intmat VV[nrows(M)][1];
1158      return (VV);
1159    }
1160
1161    ERROR("Sorry: vector or module need module-grading-matrix! See 'getModuleGrading'.");
1162  }
1163
1164  if( nrows(V) != nrows(getVariableWeights()) )
1165  {
1166    ERROR("Sorry wrong height of V: " + string(nrows(V)));
1167  }
1168
1169  if( ncols(V) < nrows(m) )
1170  {
1171    ERROR("Sorry wrong width of V: " + string(ncols(V)));
1172  }
1173
1174  return (V);
1175}
1176example
1177{
1178  "EXAMPLE:"; echo=2;
1179
1180   ring R = 0, (x,y), dp;
1181   intmat M[2][2]=
1182     1, 1,
1183     0, 2;
1184   intmat T[2][5]=
1185     1,  2,  3,  4, 0,
1186     0, 10, 20, 30, 1;
1187
1188   setBaseMultigrading(M, T);
1189
1190   ideal I = x, y, xy^5;
1191   isHomogeneous(I);
1192
1193   intmat V = multiDeg(I); print(V);
1194
1195   module S = syz(I); print(S);
1196
1197   S = setModuleGrading(S, V);
1198
1199   getModuleGrading(S) == V;
1200
1201   vector v = getGradedGenerator(S, 1);
1202   getModuleGrading(v) == V;
1203   isHomogeneous(v);
1204   print( multiDeg(v) );
1205
1206   isHomogeneous(S);
1207   print( multiDeg(S) );
1208}
1209
1210/******************************************************/
1211proc setModuleGrading(def m, intmat G)
1212"USAGE: setModuleGrading(m, G), m module/vector, G intmat
1213PURPOSE: attaches the multiweights of free module generators to 'm'
1214WARNING: The method does not verify whether the multigrading makes the
1215         module/vector homogeneous. One can do that using isHomogeneous(m).
1216EXAMPLE: example setModuleGrading; shows an example
1217"
1218{
1219  string attrModuleGrading = "genWeights";
1220
1221  intmat R = getVariableWeights();
1222
1223  if(nrows(G) != nrows(R)){ ERROR("Incompatible gradings.");}
1224  if(ncols(G) < nrows(m)){ ERROR("Multigrading does not fit to module.");}
1225
1226  attrib(m, attrModuleGrading, G);
1227  return(m);
1228}
1229example
1230{
1231  "EXAMPLE:"; echo=2;
1232
1233   ring R = 0, (x,y), dp;
1234   intmat M[2][2]=
1235     1, 1,
1236     0, 2;
1237   intmat T[2][5]=
1238     1,  2,  3,  4, 0,
1239     0, 10, 20, 30, 1;
1240
1241   setBaseMultigrading(M, T);
1242
1243   ideal I = x, y, xy^5;
1244   intmat V = multiDeg(I);
1245
1246   // V == M; modulo T
1247   print(V);
1248
1249   module S = syz(I);
1250
1251   S = setModuleGrading(S, V);
1252   getModuleGrading(S) == V;
1253
1254   print(S);
1255
1256   vector v = getGradedGenerator(S, 1);
1257   getModuleGrading(v) == V;
1258
1259   print( multiDeg(v) );
1260
1261   isHomogeneous(S);
1262
1263   print( multiDeg(S) );
1264}
1265
1266
1267proc multiDegTensor(module m, module n)
1268"
1269USAGE: multiDegTensor(m, n), m,n  modules or matrices.
1270PURPOSE: Computes the multigraded tensor product of to multigraded modules.
1271RETURN: A module.
1272EXAMPLE: example multiDegTensor; shows an example
1273"
1274{
1275  matrix M = m;
1276  matrix N = n;
1277  intmat gm = getModuleGrading(m);
1278  intmat gn = getModuleGrading(n);
1279  int grows = nrows(gm);
1280  int mr = nrows(M);
1281  int mc = ncols(M);
1282  if(rank(M) == 0){ mc = 0;}
1283  int nr = nrows(N);
1284  int nc = ncols(N);
1285  if(rank(N) == 0){ nc = 0;}
1286  intmat gresult[nrows(gm)][mr*nr];
1287  matrix result[mr*nr][mr*nc+mc*nr];
1288  int i, j;
1289  int column = 1;
1290  for(i = 1; i<=mr; i++){
1291    for(j = 1; j<=nr; j++){
1292      gresult[1..grows,(i-1)*nr+j] = gm[1..grows,i]+gn[1..grows,j];
1293    }
1294  }
1295  //gresult;
1296  if( nc!=0 ){
1297    for(i = 1; i<=mr; i++)
1298    {
1299      result[((i-1)*nr+1)..(i*nr),((i-1)*nc+1)..(i*nc)] = N[1..nr,1..nc];
1300    }
1301  }
1302  list rownumbers, colnumbers;
1303  //print(result);
1304  if( mc!=0 ){
1305    for(j = 1; j<=nr; j++)
1306    {
1307      rownumbers = nr*(0..(mr-1))+j*(1:mr);
1308      colnumbers = ((mr*nc+j):mc)+nr*(0..(mc-1));
1309      result[rownumbers[1..mr],colnumbers[1..mc] ] = M[1..mr,1..mc];
1310    }
1311  }
1312  module res = result;
1313  res = setModuleGrading(res, gresult);
1314  //getModuleGrading(res);
1315  return(res);
1316}
1317example
1318{
1319"EXAMPLE: ";echo=2;
1320ring r = 0,(x),dp;
1321intmat g[2][1]=1,1;
1322setBaseMultigrading(g);
1323matrix m[5][3]=1,2,3,4,5,6,7,8,9,10,11,12,13,14,15;
1324matrix n[3][2]=x,x2,x3,x4,x5,x6;
1325module mm = m;
1326module nn = n;
1327intmat gm[2][5]=1,2,3,4,5,0,0,0,0,0;
1328intmat gn[2][3]=0,0,0,1,2,3;
1329mm = setModuleGrading(mm, gm);
1330nn = setModuleGrading(nn, gn);
1331module mmtnn = multiDegTensor(mm, nn);
1332print(mmtnn);
1333getModuleGrading(mmtnn);
1334LIB "homolog.lib";
1335module tt = tensorMod(mm,nn);
1336print(tt);
1337
1338kill m, mm, n, nn, gm, gn;
1339
1340matrix m[7][3] = x, x-1,x+2, 3x, 4x, x5, x6, x-7, x-8, 9, 10, 11x, 12 -x, 13x, 14x, x15, (x-4)^2, x17, 18x, 19x, 20x, 21x;
1341matrix n[2][4] = 1, 2, 3, 4, x, x2, x3, x4;
1342module mm = m;
1343module nn = n;
1344print(mm);
1345print(nn);
1346intmat gm[2][7] = 1, 2, 3, 4, 5, 6, 7, 0, 0, 0, 0, 0, 0, 0;
1347intmat gn[2][2] = 0, 0, 1, 2;
1348mm = setModuleGrading(mm, gm);
1349nn = setModuleGrading(nn, gn);
1350module mmtnn = multiDegTensor(mm, nn);
1351print(mmtnn);
1352getModuleGrading(mmtnn);
1353matrix a = mmtnn;
1354matrix b = tensorMod(mm, nn);
1355print(a-b);
1356
1357}
1358
1359proc multiDegTor(int i, module m, module n)
1360{
1361  def res = multiDegResolution(n, 0, 1);
1362  //print(res);
1363  list l = res;
1364  if(size(l)<i){ return(0);}
1365  else
1366  {
1367
1368    matrix fd[nrows(m)][0];
1369    matrix fd2[nrows(l[i+1])][0];
1370    matrix fd3[nrows(l[i])][0];
1371
1372    module freedim = fd;
1373    module freedim2 = fd2;
1374    module freedim3 = fd3;
1375
1376    freedim = setModuleGrading(freedim,getModuleGrading(m));
1377    freedim2 = setModuleGrading(freedim2,getModuleGrading(l[i+1]));
1378    freedim3 = setModuleGrading(freedim3, getModuleGrading(l[i]));
1379
1380    module mimag = multiDegTensor(freedim3, m);
1381    //"mimag ok.";
1382    module mf = multiDegTensor(l[i], freedim);
1383    //"mf ok.";
1384    module mim1 = multiDegTensor(freedim2 ,m);
1385    module mim2 = multiDegTensor(l[i+1],freedim);
1386    //"mim1+2 ok.";
1387    module mker = multiDegModulo(mf,mimag);
1388    //"mker ok.";
1389    module mim = mim1,mim2;
1390    mim = setModuleGrading(mim, getModuleGrading(mim1));
1391    //"mim: r: ",nrows(mim)," c: ",ncols(mim);
1392    //"mim1: r: ",nrows(mim1)," c: ",ncols(mim1);
1393    //"mim2: r: ",nrows(mim2)," c: ",ncols(mim2);
1394    //matrix mimmat = mim;
1395    //matrix mimmat1[16][4]=mimmat[1..16,25..28];
1396    //print(mimmat1-matrix(mim2));
1397    return(multiDegModulo(mker,mim));
1398    //return(0);
1399  }
1400  return(0);
1401}
1402example
1403{
1404"EXAMPLE: ";echo=2;
1405LIB "homolog.lib";
1406ring r = 0,(x_(1..4)),dp;
1407intmat g[2][4]=1,1,0,0,0,1,1,-1;
1408setBaseMultigrading(g);
1409ideal i = maxideal(1);
1410module m = multiDegSyzygy(i);
1411module rt = Tor(2,m,m);
1412module multiDegT = multiDegTor(2,m,m);
1413print(matrix(rt)-matrix(multiDegT));
1414/*
1415ring r = 0,(x),dp;
1416intmat g[2][1]=1,1;
1417setBaseMultigrading(g);
1418matrix m[5][3]=1,2,3,4,5,6,7,8,9,10,11,12,13,14,15;
1419matrix n[3][2]=x,x2,x3,x4,x5,x6;
1420module mm = m;
1421module nn = n;
1422intmat gm[2][5]=1,1,1,1,1,1,1,1,1,1,1;
1423intmat gn[2][3]=0,-2,-4,0,-2,-4;
1424mm = setModuleGrading(mm, gm);
1425nn = setModuleGrading(nn, gn);
1426isHomogeneous(mm,"checkGens");
1427isHomogeneous(nn,"checkGens");
1428multiDegTor(1,mm, nn);
1429
1430kill m, mm, n, nn, gm, gn;
1431
1432matrix m[7][3] = x, x-1,x+2, 3x, 4x, x5, x6, x-7, x-8, 9, 10, 11x, 12 -x, 13x, 14x, x15, (x-4)^2, x17, 18x, 19x, 20x, 21x;
1433matrix n[2][4] = 1, 2, 3, 4, x, x2, x3, x4;
1434module mm = m;
1435module nn = n;
1436print(mm);
1437print(nn);
1438intmat gm[2][7] = 1, 2, 3, 4, 5, 6, 7, 0, 0, 0, 0, 0, 0, 0;
1439intmat gn[2][2] = 0, 0, 1, 2;
1440mm = setModuleGrading(mm, gm);
1441nn = setModuleGrading(nn, gn);
1442module mmtnn = multiDegTensor(mm, nn);
1443*/
1444}
1445
1446
1447/******************************************************/
1448proc isGroupHomomorphism(def L1, def L2, intmat A)
1449"USAGE: isGoupHomomorphism(L1,L2,A); L1 and L2 are groups, A is an integer matrix
1450PURPOSE: checks whether A defines a group homomorphism phi: L1 --> L2
1451RETURN: int, 1 if A defines the homomorphism and 0 otherwise
1452EXAMPLE: example isGroupHomomorphism; shows an example
1453"
1454{
1455  // TODO: L1, L2
1456  if( (ncols(A) != nrows(L1)) or (nrows(A) != nrows(L2)) )
1457  {
1458    ERROR("Incompatible sizes!");
1459  }
1460
1461  intmat im = A * L1;
1462
1463  return  (areZeroElements(im, L2));
1464}
1465example
1466{
1467  "EXAMPLE:"; echo=2;
1468
1469   intmat L1[4][1]=
1470     0,
1471     0,
1472     0,
1473     2;
1474
1475   intmat L2[3][2]=
1476     0, 0,
1477     2, 0,
1478     0, 3;
1479
1480  intmat A[3][4] =
1481     1, 2, 3, 0,
1482     7, 0, 0, 0,
1483     1, 2, 0, 3;
1484  print( A );
1485
1486  isGroupHomomorphism(L1, L2, A);
1487
1488  intmat B[3][4] =
1489     1, 2, 3, 0,
1490     7, 0, 0, 0,
1491     1, 2, 0, 2;
1492  print( B );
1493
1494  isGroupHomomorphism(L1, L2, B); // Not a homomorphism!
1495}
1496
1497/******************************************************/
1498proc isTorsionFree()
1499"USAGE: isTorsionFree()
1500PURPOSE: Determines whether the multigrading attached to the current ring is free.
1501RETURN: boolean, the result of the test
1502EXAMPLE: example isTorsionFree; shows an example
1503"
1504{
1505  intmat H = smithNormalForm(getLattice()); // TODO: ?cache it?  //******
1506
1507  int i, j;
1508  int r = nrows(H);
1509  int c = ncols(H);
1510  int d = 1;
1511  for( i = 1; (i <= c) && (i <= r); i++ )
1512  {
1513    for( j = i; (H[j, i] == 0)&&(j < r); j++ )
1514    {
1515    }
1516
1517    if(H[j, i]!=0)
1518    {
1519      d=d*H[j, i];
1520    }
1521  }
1522
1523  if( (d*d)==1 )
1524  {
1525    return(1==1);
1526  }
1527  return(0==1);
1528}
1529example
1530{
1531  "EXAMPLE:"; echo=2;
1532
1533   ring R = 0,(x,y),dp;
1534   intmat M[2][2]=
1535     1,0,
1536     0,1;
1537   intmat T[2][5]=
1538     1, 2, 3, 4, 0,
1539     0,10,20,30, 1;
1540
1541   setBaseMultigrading(M,T);
1542
1543   // Is the resulting group  free?
1544   isTorsionFree();
1545
1546   kill R, M, T;
1547   ///////////////////////////////////////////
1548
1549   ring R=0,(x,y,z),dp;
1550   intmat A[3][3] =
1551     1,0,0,
1552     0,1,0,
1553     0,0,1;
1554   intmat B[3][4]=
1555     3,3,3,3,
1556     2,1,3,0,
1557     1,2,0,3;
1558   setBaseMultigrading(A,B);
1559   // Is the resulting group  free?
1560   isTorsionFree();
1561
1562   kill R, A, B;
1563}
1564
1565
1566static proc gcdcomb(int a, int b)
1567{
1568  // a;
1569  // b;
1570  intvec av = a,1,0;
1571  intvec bv = b,0,1;
1572  intvec save;
1573  while(av[1]*bv[1] != 0)
1574  {
1575    bv = bv - (bv[1] - bv[1]%av[1]) div av[1] * av;
1576    save = bv;
1577    bv = av;
1578    av = save;
1579  }
1580  if(bv[1] < 0)
1581  {
1582    bv = -bv;
1583  }
1584  return(bv);
1585}
1586
1587
1588proc lll(def A)
1589"
1590The lll algorithm of lll.lib only works for lists of vectors.
1591Maybe one should rescript it for matrices. This method will
1592convert a matrix to a list, plug it into lll and make the result
1593a matrix and return it.
1594"
1595{
1596  if(typeof(A) == "list")
1597  {
1598    int sizeA= size (A);
1599    if (sizeA == 0)
1600    {
1601      return (A);
1602    }
1603    if (typeof (A [1]) != "intvec")
1604    {
1605      ERROR("Unrecognized type.");
1606    }
1607    int columns= size (A [1]);
1608    int i;
1609    for (i= 2; i <= sizeA; i++)
1610    {
1611      if (typeof (A[i]) != "intvec")
1612      {
1613        ERROR("Unrecognized type.");
1614      }
1615      if (size (A [i]) != columns)
1616      {
1617        ERROR ("expected equal dimension");
1618      }
1619    }
1620    int j;
1621    intmat m [columns] [sizeA];
1622    for (i= 1; i <= sizeA; i++)
1623    {
1624      for (j= 1; j <= columns; j++)
1625      {
1626        m[i,j]= A[i] [j];
1627      }
1628    }
1629    m= system ("LLL", m);
1630    list result= list();
1631    intvec buf;
1632
1633    for (i= 1; i <= sizeA; i++)
1634    {
1635      buf = intvec (m[i , 1..columns]);
1636      result= result+ list (buf);
1637
1638    }
1639    return(result);
1640  }
1641  else
1642  {
1643    if(typeof(A) == "intmat")
1644    {
1645      A= system ("LLL", A);
1646      return(A);
1647    }
1648    else
1649    {
1650      ERROR("Unrecognized type.");
1651    }
1652  }
1653}
1654example
1655{
1656  "EXAMPLE:"; echo=2;
1657
1658  ring R = 0,x,dp;
1659  intmat m[5][5] =
1660                  13,25,37,83,294,
1661                  12,-33,9,0,64,
1662                  77,12,34,6,1,
1663                  43,2,88,91,100,
1664                  -46,32,37,42,15;
1665  lll(m);
1666 
1667  list l =
1668      intvec(13,25,37, 83, 294),
1669      intvec(12, -33, 9,0,64),
1670      intvec (77,12,34,6,1),
1671      intvec (43,2,88,91,100),
1672      intvec (-46,32,37,42,15);
1673  lll(l);
1674}
1675
1676
1677proc smithNormalForm(intmat A, list #)
1678"USAGE: smithNormalForm(A[,opt]); intmat A
1679PURPOSE: Computes the Smith Normal Form of A
1680RETURN: if no optional argument is given: intmat, the Smith Normal Form of A,
1681otherwise: a list of 3 integer matrices P, D Q, such that D == P*A*Q.
1682EXAMPLE: example smithNormalForm; shows an example
1683"
1684{
1685  list l1 = hermiteNormalForm(A, 5);
1686  // l1;
1687  intmat B = transpose(l1[1]);
1688  list l2 = hermiteNormalForm(B, 5);
1689  // l2;
1690  intmat P = transpose(l2[2]);
1691  intmat D = transpose(l2[1]);
1692  intmat Q = l1[2];
1693  int cc = ncols(D);
1694  int rr = nrows(D);
1695  intmat transform;
1696  int k = 1;
1697  int a, b, c;
1698  // D;
1699  intvec v;
1700  if((cc==1)||(rr==1)){
1701    if(size(#)==0)
1702    {
1703      return(D);
1704    } else
1705    {
1706      return(list(P,D,Q));
1707    }
1708  }
1709  while(D[k+1,k+1] !=0){
1710    if(D[k+1,k+1]%D[k,k]!=0){
1711      b = D[k, k]; c = D[k+1, k+1];
1712      v = gcdcomb(D[k,k],D[k+1,k+1]);
1713      transform = unitMatrix(cc);
1714      transform[k+1,k] = 1;
1715      a = -v[3]*D[k+1,k+1] div v[1];
1716      transform[k, k+1] = a;
1717      transform[k+1, k+1] = a+1;
1718      //det(transform);
1719      D = D*transform;
1720      Q = Q*transform;
1721      //D;
1722      transform = unitMatrix(rr);
1723      transform[k,k] = v[2];
1724      transform[k,k+1] = v[3];
1725      transform[k+1,k] = -c div v[1];
1726      transform[k+1,k+1] = b div v[1];
1727      D = transform * D;
1728      P = transform * P;
1729      //" ";
1730      //D;
1731      //"small transform: ", det(transform);
1732      //transform;
1733      k=0;
1734    }
1735    k++;
1736    if((k==rr) || (k==cc)){
1737      break;
1738    }
1739  }
1740  //"here is the size ",size(#);
1741  if(size(#) == 0){
1742    return(D);
1743  } else {
1744    return(list(P, D, Q));
1745  }
1746}
1747example
1748{
1749  "EXAMPLE:"; echo=2;
1750
1751
1752  intmat A[5][7] =
1753                  1,0,1,0,-2,9,-71,
1754                  0,-24,248,-32,-96,448,-3496,
1755                  0,4,-42,4,-8,30,-260,
1756                  0,0,0,18,-90,408,-3168,
1757                  0,0,0,-32,224,-1008,7872;
1758
1759  print( smithNormalForm(A) );
1760
1761  list l = smithNormalForm(A, 5);
1762
1763  l;
1764
1765  l[1]*A*l[3];
1766
1767  det(l[1]);
1768  det(l[3]);
1769}
1770
1771
1772/******************************************************/
1773proc hermiteNormalForm(intmat A, list #)
1774"USAGE: hermiteNormalForm( A );
1775PURPOSE: Computes the (lower triangular) Hermite Normal Form
1776           of the matrix A by column operations.
1777RETURN: intmat, the Hermite Normal Form of A
1778EXAMPLE: example hermiteNormalForm; shows an example
1779"
1780{
1781
1782  int row, column, i, j;
1783  int rr = nrows(A);
1784  int cc = ncols(A);
1785  intvec savev, gcdvec, v1, v2;
1786  intmat q = unitMatrix(cc);
1787  intmat transform;
1788  column = 1;
1789  for(row = 1; (row<=rr)&&(column<=cc); row++)
1790    {
1791      if(A[row,column]==0)
1792        {
1793          for(j = column; j<=cc; j++)
1794            {
1795              if(A[row, j]!=0)
1796                {
1797                  transform = unitMatrix(cc);
1798                  transform[j,j] = 0;
1799                  transform[column, column] = 0;
1800                  transform[column,j] = 1;
1801                  transform[j,column] = 1;
1802                  q = q*transform;
1803                  A = A*transform;
1804                  break;
1805                }
1806            }
1807        }
1808      if(A[row,column] == 0)
1809        {
1810          row++;
1811          continue;
1812        }
1813      for(j = column+1; j<=cc; j++)
1814        {
1815          if(A[row, j]!=0)
1816            {
1817              gcdvec = gcdcomb(A[row,column],A[row,j]);
1818              // gcdvec;
1819              // typeof(A[1..rr,column]);
1820              v1 = A[1..rr,column];
1821              v2 = A[1..rr,j];
1822              transform = unitMatrix(cc);
1823              transform[j,j] = v1[row] div gcdvec[1];
1824              transform[column, column] = gcdvec[2];
1825              transform[column,j] = -v2[row] div gcdvec[1];
1826              transform[j,column] = gcdvec[3];
1827              q = q*transform;
1828              A = A*transform;
1829              // A;
1830            }
1831        }
1832      if(A[row,column]<0)
1833        {
1834          transform = unitMatrix(cc);
1835          transform[column,column] = -1;
1836          q = q*transform;
1837          A = A*transform;
1838        }
1839      for( j=1; j<column; j++){
1840        if(A[row, j]!=0){
1841          transform = unitMatrix(cc);
1842          transform[column, j] = (-A[row,j]+A[row, j]%A[row, column]) div A[row, column];
1843          if(A[row,j]<0){
1844            transform[column,j]=transform[column,j]+1;}
1845          q = q*transform;
1846          A = A*transform;
1847        }
1848      }
1849      column++;
1850    }
1851  if(size(#) > 0){
1852    return(list(A, q));
1853  }
1854  return(A);
1855}
1856example
1857{
1858  "EXAMPLE:"; echo=2;
1859
1860   intmat M[2][5] =
1861     1, 2, 3, 4, 0,
1862     0,10,20,30, 1;
1863
1864   // Hermite Normal Form of M:
1865   print(hermiteNormalForm(M));
1866
1867   intmat T[3][4] =
1868     3,3,3,3,
1869     2,1,3,0,
1870     1,2,0,3;
1871
1872   // Hermite Normal Form of T:
1873   print(hermiteNormalForm(T));
1874
1875   intmat A[4][5] =
1876     1,2,3,2,2,
1877     1,2,3,4,0,
1878     0,5,4,2,1,
1879     3,2,4,0,2;
1880
1881   // Hermite Normal Form of A:
1882   print(hermiteNormalForm(A));
1883}
1884
1885proc areZeroElements(intmat m, list #)
1886"USAGE: areZeroElements(D, [T]); intmat D, group T
1887PURPOSE: For a integer matrix D, considered column-wise as a set of
1888   integer vecors representing the multidegree of some polynomial
1889   or vector this method checks whether all these multidegrees
1890   are contained in the grading group
1891   group (either set globally or given as an optional argument),
1892i.e. if they all are zero in the multigrading.
1893EXAMPLE: example areZeroElements; shows an example
1894"
1895{
1896  int r = nrows(m);
1897  int i = ncols(m);
1898
1899  intvec v;
1900
1901  for( ; i > 0; i-- )
1902  {
1903    v = m[1..r, i];
1904    if( !isZeroElement(v, #) )
1905    {
1906      return (0);
1907    }
1908  }
1909  return(1);
1910}
1911example
1912{
1913  "EXAMPLE:"; echo=2;
1914
1915  ring r = 0,(x,y,z),dp;
1916
1917  intmat S[2][3]=
1918    1,0,1,
1919    0,1,1;
1920
1921  intmat L[2][1]=
1922    2,
1923    2;
1924
1925  setBaseMultigrading(S,L);
1926
1927  poly a = 1;
1928  poly b = xyz;
1929
1930  ideal I = a, b;
1931  print(multiDeg(I));
1932 
1933  intmat m[5][2]=multiDeg(a),multiDeg(b); m=transpose(m);
1934
1935  print(multiDeg(a));
1936  print(multiDeg(b));
1937
1938  print(m);
1939 
1940  areZeroElements(m);
1941
1942  intmat LL[2][1]=
1943     1,
1944    -1;
1945
1946  areZeroElements(m,LL);
1947}
1948
1949
1950/******************************************************/
1951proc isZeroElement(intvec mdeg, list #)
1952"USAGE: isZeroElement(d, [T]); intvec d, group T
1953PURPOSE: For a integer vector 'd' representing the multidegree of some polynomial
1954or vector this method computes if the multidegree is contained in the grading group
1955group (either set globally or given as an optional argument), i.e. if it is zero in the multigrading.
1956EXAMPLE: example isZeroElement; shows an example
1957"
1958{
1959  int i = 1;
1960  if( size(#) >= i )
1961  {
1962    def a = #[1];
1963    if( typeof(a) == "intmat" )
1964    {
1965      intmat H = hermiteNormalForm(a);
1966      i++;
1967    }
1968    if( typeof(a) == "list" )
1969    {
1970      list L = a;
1971      intmat H = attrib(L, "hermite"); // todo
1972      i++;
1973    }
1974    kill a;
1975  }
1976 
1977  if( i == 1 )
1978  {
1979    intmat H = getLattice("hermite");
1980  }
1981
1982  int x, k, row;
1983
1984  int r = nrows(H);
1985  int c = ncols(H);
1986
1987  int rr = nrows(mdeg);
1988  row = 1;
1989  intvec v;
1990  for(i=1; (i<=r)&&(row<=r)&&(i<=c); i++)
1991  {
1992    while((H[row,i]==0)&&(row<=r))
1993    {
1994      row++;
1995      if(row == (r+1)){
1996        break;
1997      }
1998    }
1999    if(row<=r){
2000      if(H[row,i]!=0)
2001      {
2002        v = H[1..r,i];
2003        mdeg = mdeg-(mdeg[row]-mdeg[row]%v[row]) div v[row]*v;
2004      }
2005    }
2006  }
2007  return( mdeg == 0 );
2008
2009}
2010example
2011{
2012 "EXAMPLE:"; echo=2;
2013
2014  ring r = 0,(x,y,z),dp;
2015
2016  intmat g[2][3]=
2017    1,0,1,
2018    0,1,1;
2019  intmat t[2][1]=
2020    -2,
2021    1;
2022
2023  intmat tt[2][1]=
2024    1,
2025    -1;
2026
2027  setBaseMultigrading(g,t);
2028
2029  poly a = x10yz;
2030  poly b = x8y2z;
2031  poly c = x4z2;
2032  poly d = y5;
2033  poly e = x2y2;
2034  poly f = z2;
2035
2036  intvec v1 = multiDeg(a) - multiDeg(b);
2037  v1;
2038  isZeroElement(v1);
2039  isZeroElement(v1, tt);
2040
2041  intvec v2 = multiDeg(a) - multiDeg(c);
2042  v2;
2043  isZeroElement(v2);
2044  isZeroElement(v2, tt);
2045
2046  intvec v3 = multiDeg(e) - multiDeg(f);
2047  v3;
2048  isZeroElement(v3);
2049  isZeroElement(v3, tt);
2050
2051  intvec v4 = multiDeg(c) - multiDeg(d);
2052  v4;
2053  isZeroElement(v4);
2054  isZeroElement(v4, tt);
2055}
2056
2057
2058/******************************************************/
2059proc defineHomogeneous(poly f, list #)
2060"USAGE: defineHomogeneous(f[, G]); polynomial f, integer matrix G
2061PURPOSE: Yields a matrix which has to be appended to the grading group matrix to make the
2062polynomial f homogeneous in the grading by grad.
2063EXAMPLE: example defineHomogeneous; shows an example
2064"
2065{
2066  int i = 1;
2067  if( size(#) >= i )
2068  {
2069    def a = #[1];
2070    if( typeof(a) == "intmat" )
2071    {
2072      intmat grad = a;
2073      i++;
2074    }
2075    kill a;
2076  }
2077
2078  if( i == 1 )
2079  {
2080    intmat grad = getVariableWeights();
2081  }
2082
2083  intmat newgg[nrows(grad)][size(f)-1];
2084  int j;
2085  intvec l = grad*leadexp(f);
2086  intvec v;
2087  for(i=2; i <= size(f); i++)
2088  {
2089    v = grad * leadexp(f[i]) - l;
2090    for( j=1; j<=size(v); j++)
2091    {
2092      newgg[j,i-1] = v[j];
2093    }
2094  }
2095  return(newgg);
2096}
2097example
2098{
2099  "EXAMPLE:"; echo=2;
2100
2101  ring r =0,(x,y,z),dp;
2102  intmat grad[2][3] =
2103    1,0,1,
2104    0,1,1;
2105
2106  setBaseMultigrading(grad);
2107
2108  poly f = x2y3-z5+x-3zx;
2109
2110  intmat M = defineHomogeneous(f);
2111  M;
2112  defineHomogeneous(f, grad) == M;
2113
2114  isHomogeneous(f);
2115  setBaseMultigrading(grad, M);
2116  isHomogeneous(f);
2117}
2118
2119
2120proc gradiator(def h)
2121"PURPOSE: coarsens the grading of the basering until the polynom or ideal h becomes homogeneous."
2122{
2123  if(typeof(h)=="poly"){
2124    intmat W = getVariableWeights();
2125    intmat L = getLattice();
2126    intmat toadd = defineHomogeneous(h);
2127    //h;
2128    //toadd;
2129    if(ncols(toadd) == 0)
2130    {
2131      return(1==1);
2132    }
2133    int rr = nrows(W);
2134    intmat newL[rr][ncols(L)+ncols(toadd)];
2135    newL[1..rr,1..ncols(L)] = L[1..rr,1..ncols(L)];
2136    newL[1..rr,(ncols(L)+1)..(ncols(L)+ncols(toadd))] = toadd[1..rr,1..ncols(toadd)];
2137    setBaseMultigrading(W,newL);
2138    return(1==1);
2139  }
2140  if(typeof(h)=="ideal"){
2141    int i;
2142    def s = (1==1);
2143    for(i=1;i<=size(h);i++){
2144      s = s && gradiator(h[i]);
2145    }
2146    return(s);
2147  }
2148  return(1==0);
2149}
2150example
2151{
2152  "EXAMPLE:"; echo=2;
2153
2154ring r = 0,(x,y,z),dp;
2155intmat g[2][3] = 1,0,1,0,1,1;
2156intmat l[2][1] = 3,0;
2157
2158setBaseMultigrading(g,l);
2159
2160getLattice();
2161
2162ideal i = -y5+x4,
2163          y6+xz,
2164          x2y;
2165gradiator(i);
2166getLattice();
2167isHomogeneous(i);
2168}
2169
2170
2171proc pushForward(map f)
2172"USAGE: pushForward(f);
2173PURPOSE: Computes the finest grading of the image ring which makes the map f
2174a map of graded rings. The group map between the two grading groups is given
2175by transpose( (Id, 0) ). Pay attention that the group spanned by the columns of
2176the grading group matrix may not be a subgroup of the grading group. Still all columns
2177are needed to find the correct image of the preimage gradings.
2178EXAMPLE: example pushForward; shows an example
2179"
2180{
2181
2182  int k,i,j;
2183//  f;
2184
2185//  listvar();
2186  def pre = preimage(f);
2187
2188//  "pre: ";  pre;
2189
2190  intmat oldgrad=getVariableWeights(pre);
2191  intmat oldlat=getLattice(pre);
2192
2193  int n=nvars(pre);
2194  int np=nvars(basering);
2195  int p=nrows(oldgrad);
2196  int pp=p+np;
2197
2198  intmat newgrad[pp][np];
2199
2200  //This will set the finest grading on the image ring. We will proceed by coarsening this grading until f becomes homogeneous.
2201  for(i=1;i<=np;i++){ newgrad[p+i,i]=1;}
2202
2203  //newgrad;
2204
2205
2206
2207  list newlat;
2208  intmat toadd;
2209  int columns=0;
2210
2211  intmat toadd1[pp][n];
2212  intvec v;
2213  poly im;
2214
2215  for(i=1;i<=p;i++){
2216    for(j=1;j<=n;j++){ toadd1[i,j]=oldgrad[i,j];}
2217  }
2218 
2219  // This will make the images of homogeneous elements homogeneous, namely the variables of the preimage ring.
2220  for(i=1;i<=n;i++){
2221    im=f[i];
2222    //im;
2223    toadd = defineHomogeneous(im, newgrad);
2224    newlat=insert(newlat,toadd);
2225    columns=columns+ncols(toadd);
2226
2227    v=leadexp(f[i]);
2228    for(j=p+1;j<=p+np;j++){ toadd1[j,i]=-v[j-p];}
2229  }
2230
2231  newlat=insert(newlat,toadd1);
2232  columns=columns+ncols(toadd1);
2233
2234  //If the image ring is a quotient ring by some ideal, we have to coarsen the grading in order to make the ideal homogeneous.
2235  if(typeof(basering)=="qring"){
2236    //"Entering qring";
2237    ideal a=ideal(basering);
2238    for(i=1;i<=size(a);i++){
2239      toadd = defineHomogeneous(a[i], newgrad);
2240      //toadd;
2241      columns=columns+ncols(toadd);
2242      newlat=insert(newlat,toadd);
2243    }
2244  }
2245
2246  //The grading group of the preimage ring might not have been torsion free. We have to add this torsion to the grading group of the image ring.
2247  intmat imofoldlat[pp][ncols(oldlat)];
2248  for(i=1; i<=nrows(oldlat);i++){
2249    for(j=1; j<=ncols(oldlat); j++){
2250      imofoldlat[i,j]=oldlat[i,j];
2251    }
2252  }
2253
2254  columns=columns+ncols(oldlat);
2255  newlat=insert(newlat, imofoldlat);
2256
2257  intmat gragr[pp][columns];
2258  columns=0;
2259  for(k=1;k<=size(newlat);k++){
2260    for(i=1;i<=pp;i++){
2261      for(j=1;j<=ncols(newlat[k]);j++){gragr[i,j+columns]=newlat[k][i,j];}
2262    }
2263    columns=columns+ncols(newlat[k]);
2264  }
2265 
2266  //The following is just for reducing the size of the matrices.
2267  gragr=hermiteNormalForm(gragr);
2268  intmat result[pp][pp];
2269  for(i=1;i<=pp;i++){
2270    for(j=1;j<=pp;j++){result[i,j]=gragr[i,j];}
2271  }
2272
2273  setBaseMultigrading(newgrad, result);
2274
2275}
2276example
2277{
2278  "EXAMPLE:"; echo=2;
2279
2280  ring r = 0,(x,y,z),dp;
2281
2282
2283
2284  // Setting degrees for preimage ring.;
2285  intmat grad[3][3] =
2286    1,0,0,
2287    0,1,0,
2288    0,0,1;
2289
2290  setBaseMultigrading(grad);
2291
2292  // grading on r:
2293  getVariableWeights();
2294  getLattice();
2295
2296  // only for the purpose of this example
2297  if( voice > 1 ){ /*keepring(r);*/ export(r); }
2298
2299  ring R = 0,(a,b),dp;
2300  ideal i = a2-b2+a6-b5+ab3,a7b+b15-ab6+a6b6;
2301
2302  // The quotient ring by this ideal will become our image ring.;
2303  qring Q = std(i);
2304
2305  listvar();
2306
2307  map f = r,-a2b6+b5+a3b+a2+ab,-a2b7-3a2b5+b4+a,a6-b6-b3+a2; f;
2308
2309
2310  // TODO: Unfortunately this is not a very spectacular example...:
2311  // Pushing forward f:
2312  pushForward(f);
2313
2314  // due to pushForward we have got new grading on Q
2315  getVariableWeights();
2316  getLattice();
2317
2318
2319  // only for the purpose of this example
2320  if( voice > 1 ){ kill r; }
2321
2322}
2323
2324
2325/******************************************************/
2326proc equalMultiDeg(intvec exp1, intvec exp2, list #)
2327"USAGE: equalMultiDeg(exp1, exp2[, V]); intvec exp1, exp2, intmat V
2328PURPOSE: Tests if the exponent vectors of two monomials (given by exp1 and exp2)
2329represent the same multidegree.
2330NOTE: the integer matrix V encodes multidegrees of module components,
2331if module component is present in exp1 and exp2
2332EXAMPLE: example equalMultiDeg; shows an example
2333"
2334{
2335  if( size(exp1) != size(exp2) )
2336  {
2337    ERROR("Sorry: we cannot compare exponents comming from a polynomial and a vector yet!");
2338  }
2339
2340  if( exp1 == exp2)
2341  {
2342    return (1==1);
2343  }
2344
2345
2346
2347  intmat M = getVariableWeights();
2348
2349  if( nrows(exp1) > ncols(M) ) // vectors => last exponent is the module component!
2350  {
2351    if( (size(#) == 0) or (typeof(#[1])!="intmat") )
2352    {
2353      ERROR("Sorry: wrong or missing module-unit-weights-matrix V!");
2354    }
2355    intmat V = #[1];
2356
2357    // typeof(V); print(V);
2358
2359    int N = ncols(M);
2360    int r = nrows(M);
2361
2362    intvec d = intvec(exp1[1..N]) - intvec(exp2[1..N]);
2363    intvec dm = intvec(V[1..r, exp1[N+1]]) - intvec(V[1..r, exp2[N+1]]);
2364
2365    intvec difference = M * d + dm;
2366  }
2367  else
2368  {
2369    intvec d = (exp1 - exp2);
2370    intvec difference = M * d;
2371  }
2372
2373  if (isFreeRepresented()) // no grading group!?
2374  {
2375    return ( difference == 0);
2376  }
2377  return ( isZeroElement( difference ) );
2378}
2379example
2380{
2381  "EXAMPLE:"; echo=2;printlevel=3;
2382
2383  ring r = 0,(x,y,z),dp;
2384
2385  intmat g[2][3]=
2386    1,0,1,
2387    0,1,1;
2388
2389  intmat t[2][1]=
2390    -2,
2391    1;
2392
2393  setBaseMultigrading(g,t);
2394
2395  poly a = x10yz;
2396  poly b = x8y2z;
2397  poly c = x4z2;
2398  poly d = y5;
2399  poly e = x2y2;
2400  poly f = z2;
2401
2402
2403  equalMultiDeg(leadexp(a), leadexp(b));
2404  equalMultiDeg(leadexp(a), leadexp(c));
2405  equalMultiDeg(leadexp(a), leadexp(d));
2406  equalMultiDeg(leadexp(a), leadexp(e));
2407  equalMultiDeg(leadexp(a), leadexp(f));
2408
2409  equalMultiDeg(leadexp(b), leadexp(c));
2410  equalMultiDeg(leadexp(b), leadexp(d));
2411  equalMultiDeg(leadexp(b), leadexp(e));
2412  equalMultiDeg(leadexp(b), leadexp(f));
2413
2414  equalMultiDeg(leadexp(c), leadexp(d));
2415  equalMultiDeg(leadexp(c), leadexp(e));
2416  equalMultiDeg(leadexp(c), leadexp(f));
2417
2418  equalMultiDeg(leadexp(d), leadexp(e));
2419  equalMultiDeg(leadexp(d), leadexp(f));
2420
2421  equalMultiDeg(leadexp(e), leadexp(f));
2422
2423}
2424
2425
2426
2427/******************************************************/
2428static proc isFreeRepresented()
2429"check whether the base muligrading is free (it is zero).
2430"
2431{
2432  intmat T = getLattice();
2433
2434  intmat Z[nrows(T)][ncols(T)];
2435
2436  return (T == Z); // no grading group!
2437}
2438
2439
2440/******************************************************/
2441proc isHomogeneous(def a, list #)
2442"USAGE: isHomogeneous(a[, f]); a polynomial/vector/ideal/module
2443RETURN: boolean, TRUE if a is (multi)homogeneous, and FALSE otherwise
2444EXAMPLE: example isHomogeneous; shows an example
2445"
2446{
2447  if( (typeof(a) == "poly") or (typeof(a) == "vector") )
2448  {
2449    return ( size(multiDegPartition(a)) <= 1 )
2450  }
2451
2452  if( (typeof(a) == "ideal") or (typeof(a) == "module") )
2453  {
2454    if(size(#) > 0)
2455    {
2456      if (#[1] == "checkGens")
2457      {
2458        def aa;
2459        for( int i = ncols(a); i > 0; i-- )
2460        {
2461          aa = getGradedGenerator(a, i);
2462
2463          if(!isHomogeneous(aa))
2464          {
2465            return(0==1);
2466          }
2467        }
2468        return(1==1);
2469      }
2470    }
2471
2472    def g = groebner(a); // !!!!
2473
2474    def b, aa; int j;
2475    for( int i = ncols(a); i > 0; i-- )
2476    {
2477      aa = getGradedGenerator(a, i);
2478
2479      b = multiDegPartition(aa);
2480      for( j = ncols(b); j > 0; j-- )
2481      {
2482        if(NF(b[j],g) != 0)
2483        {
2484          return(0==1);
2485        }
2486      }
2487    }
2488    return(1==1);
2489  }
2490}
2491example
2492{
2493  "EXAMPLE:"; echo=2;
2494
2495  ring r = 0,(x,y,z),dp;
2496
2497  //Grading and Torsion matrices:
2498  intmat M[3][3] =
2499    1,0,0,
2500    0,1,0,
2501    0,0,1;
2502
2503  intmat T[3][1] =
2504    1,2,3;
2505
2506  setBaseMultigrading(M,T);
2507
2508  attrib(r);
2509
2510  poly f = x-yz;
2511
2512  multiDegPartition(f);
2513  print(multiDeg(_));
2514
2515  isHomogeneous(f);   // f: is not homogeneous
2516
2517  poly g = 1-xy2z3;
2518  isHomogeneous(g); // g: is homogeneous
2519  multiDegPartition(g);
2520
2521  kill T;
2522  /////////////////////////////////////////////////////////
2523  // new Torsion matrix:
2524  intmat T[3][4] =
2525    3,3,3,3,
2526    2,1,3,0,
2527    1,2,0,3;
2528
2529  setBaseMultigrading(M,T);
2530
2531  f;
2532  isHomogeneous(f);
2533  multiDegPartition(f);
2534
2535  // ---------------------
2536  g;
2537  isHomogeneous(g);
2538  multiDegPartition(g);
2539
2540  kill r, T, M;
2541
2542  ring R = 0, (x,y,z), dp;
2543
2544  intmat A[2][3] =
2545    0,0,1,
2546    3,2,1;
2547  intmat T[2][1] =
2548    -1,
2549     4;
2550  setBaseMultigrading(A, T);
2551
2552  isHomogeneous(ideal(x2 - y3 -xy +z, x*y-z, x^3 - y^2*z + x^2 -y^3)); // 1
2553  isHomogeneous(ideal(x2 - y3 -xy +z, x*y-z, x^3 - y^2*z + x^2 -y^3), "checkGens");
2554  isHomogeneous(ideal(x+y, x2 - y2)); // 0
2555
2556  // Degree partition:
2557  multiDegPartition(x2 - y3 -xy +z);
2558  multiDegPartition(x3 -y2z + x2 -y3 + z + 1);
2559
2560
2561  module N = gen(1) + (x+y) * gen(2), z*gen(3);
2562
2563  intmat V[2][3] = 0; // 1, 2, 3,  4, 5, 6; //  column-wise weights of components!!??
2564
2565  vector v1, v2;
2566
2567  v1 = setModuleGrading(N[1], V); v1;
2568  multiDegPartition(v1);
2569  print( multiDeg(_) );
2570
2571  v2 = setModuleGrading(N[2], V); v2;
2572  multiDegPartition(v2);
2573  print( multiDeg(_) );
2574
2575  N = setModuleGrading(N, V);
2576  isHomogeneous(N);
2577  print( multiDeg(N) );
2578
2579  ///////////////////////////////////////
2580
2581  V =
2582    1, 2, 3,
2583    4, 5, 6;
2584
2585  v1 = setModuleGrading(N[1], V); v1;
2586  multiDegPartition(v1);
2587  print( multiDeg(_) );
2588
2589  v2 = setModuleGrading(N[2], V); v2;
2590  multiDegPartition(v2);
2591  print( multiDeg(_) );
2592
2593  N = setModuleGrading(N, V);
2594  isHomogeneous(N);
2595  print( multiDeg(N) );
2596
2597  ///////////////////////////////////////
2598
2599  V =
2600    0, 0, 0,
2601    4, 1, 0;
2602
2603  N = gen(1) + x * gen(2), z*gen(3);
2604  N = setModuleGrading(N, V); print(N);
2605  isHomogeneous(N);
2606  print( multiDeg(N) );
2607  v1 = getGradedGenerator(N,1); print(v1);
2608  multiDegPartition(v1);
2609  print( multiDeg(_) );
2610  N = setModuleGrading(N, V); print(N);
2611  isHomogeneous(N);
2612  print( multiDeg(N) );
2613}
2614
2615/******************************************************/
2616proc multiDeg(def A)
2617"USAGE: multiDeg(A); polynomial/vector/ideal/module A
2618PURPOSE: compute multidegree
2619EXAMPLE: example multiDeg; shows an example
2620"
2621{
2622  def a = attrib(A, "grad");
2623  if( typeof(a) == "intvec" || typeof(a) == "intmat" )
2624  {
2625    return (a);
2626  }
2627
2628  intmat M = getVariableWeights();
2629  int N = nvars(basering);
2630
2631  if( ncols(M) != N )
2632  {
2633    ERROR("Sorry wrong mgrad-size of M: " + string(ncols(M)));
2634  }
2635
2636  int r = nrows(M);
2637
2638  if( (typeof(A) == "vector") or (typeof(A) == "module") )
2639  {
2640    intmat V = getModuleGrading(A);
2641
2642    if( nrows(V) != r )
2643    {
2644      ERROR("Sorry wrong mgrad-size of V: " + string(nrows(V)));
2645    }
2646  }
2647
2648  if( A == 0 )
2649  {
2650    intvec v; v[r] = 0;
2651    return (v);
2652  }
2653
2654  intvec m; m[r] = 0;
2655
2656  if( typeof(A) == "poly" )
2657  {
2658    intvec v = leadexp(A); //  v;
2659    m = M * v;
2660
2661    // We assume homogeneous input!
2662    return(m);
2663
2664    A = A - lead(A);
2665    while( size(A) > 0 )
2666    {
2667      v = leadexp(A); //  v;
2668      m = max( m, M * v, r ); // ????
2669      A = A - lead(A);
2670    }
2671
2672    return(m);
2673  }
2674
2675
2676  if( typeof(A) == "vector" )
2677  {
2678    intvec v;
2679    v = leadexp(A); //  v;
2680    m = intvec(M * intvec(v[1..N])) + intvec(V[1..r, v[N+1]]);
2681
2682    // We assume homogeneous input!
2683    return(m);
2684
2685    A = A - lead(A);
2686    while( size(A) > 0 )
2687    {
2688      v = leadexp(A); //  v;
2689
2690      // intvec(M * intvec(v[1..N])) + intvec(V[1..r, v[N+1]]);
2691
2692      m = max( m, intvec(M * intvec(v[1..N])) + intvec(V[1..r, v[N+1]]), r ); // ???
2693
2694      A = A - lead(A);
2695    }
2696
2697    return(m);
2698  }
2699
2700  int i, j; intvec d;
2701
2702  if( typeof(A) == "ideal" )
2703  {
2704    intmat G[ r ] [ ncols(A)];
2705    for( i = ncols(A); i > 0; i-- )
2706    {
2707      d = multiDeg( A[i] );
2708
2709      for( j = 1; j <= r; j++ ) // see ticket: 253
2710      {
2711        G[j, i] = d[j];
2712      }
2713    }
2714    return(G);
2715  }
2716
2717  if( typeof(A) == "module" )
2718  {
2719    intmat G[ r ] [ ncols(A)];
2720    vector v;
2721
2722    for( i = ncols(A); i > 0; i-- )
2723    {
2724      v = getGradedGenerator(A, i);
2725
2726      // G[1..r, i]
2727      d = multiDeg(v);
2728
2729      for( j = 1; j <= r; j++ ) // see ticket: 253
2730      {
2731        G[j, i] = d[j];
2732      }
2733
2734    }
2735
2736    return(G);
2737  }
2738
2739}
2740example
2741{
2742  "EXAMPLE:"; echo=2;
2743
2744  ring r = 0,(x, y), dp;
2745
2746  intmat A[2][2] = 1, 0, 0, 1;
2747  print(A);
2748
2749  intmat Ta[2][1] = 0, 3;
2750  print(Ta);
2751
2752  //   attrib(A, "gradingGroup", Ta); // to think about
2753
2754//  "poly:";
2755  setBaseMultigrading(A);
2756
2757
2758  multiDeg( x*x, A );
2759  multiDeg( y*y*y, A );
2760
2761  setBaseMultigrading(A, Ta);
2762
2763  multiDeg( x*x*y );
2764
2765  multiDeg( y*y*y*x );
2766
2767  multiDeg( x*y + x + 1 );
2768
2769  multiDegPartition(x*y + x + 1);
2770
2771  print ( multiDeg(0) );
2772  poly zero = 0;
2773  print ( multiDeg(zero) );
2774
2775//  "ideal:";
2776
2777  ideal I = y*x*x, x*y*y*y;
2778  print( multiDeg(I) );
2779
2780  print ( multiDeg(ideal(0)) );
2781  print ( multiDeg(ideal(0,0,0)) );
2782
2783//  "vectors:";
2784
2785  intmat B[2][2] = 0, 1, 1, 0;
2786  print(B);
2787
2788  multiDeg( setModuleGrading(y*y*y*gen(2), B ));
2789  multiDeg( setModuleGrading(x*x*gen(1), B ));
2790
2791
2792  vector V = x*gen(1) + y*gen(2);
2793  V = setModuleGrading(V, B);
2794  multiDeg( V );
2795
2796  vector v1 = setModuleGrading([0, 0, 0], B);
2797  print( multiDeg( v1 ) );
2798
2799  vector v2 = setModuleGrading([0], B);
2800  print( multiDeg( v2 ) );
2801
2802//  "module:";
2803
2804  module D = x*gen(1), y*gen(2);
2805  D;
2806  D = setModuleGrading(D, B);
2807  print( multiDeg( D ) );
2808
2809
2810  module DD = [0, 0],[0, 0, 0];
2811  DD = setModuleGrading(DD, B);
2812  print( multiDeg( DD ) );
2813
2814  module DDD = [0, 0];
2815  DDD = setModuleGrading(DDD, B);
2816  print( multiDeg( DDD ) );
2817
2818};
2819
2820
2821
2822
2823
2824/******************************************************/
2825proc multiDegPartition(def p)
2826"USAGE: multiDegPartition(def p), p polynomial/vector
2827RETURNS: an ideal/module consisting of multigraded-homogeneous parts of p
2828EXAMPLE: example multiDegPartition; shows an example
2829"
2830{ // TODO: What about an ideal or module???
2831
2832  if( typeof(p) == "poly" )
2833  {
2834    ideal I;
2835    poly mp, t, tt;
2836    intmat V;
2837  }
2838  else
2839  {
2840    if(  typeof(p) == "vector" )
2841    {
2842      module I;
2843      vector mp, t, tt;
2844      intmat V = getModuleGrading(p);
2845    }
2846    else
2847    {
2848      ERROR("Wrong ARGUMENT type!");
2849    }
2850  }
2851
2852  if( size(p) > 1)
2853  {
2854    intvec m;
2855
2856    while( p != 0 )
2857    {
2858      m = leadexp(p);
2859      mp = lead(p);
2860      p = p - lead(p);
2861      tt = p; t = 0;
2862
2863      while( size(tt) > 0 )
2864      {
2865        // TODO: we do not cache matrices (M,T,H,V), which remain the same :(
2866        // TODO: we need some low-level procedure with all these arguments...!
2867        if( equalMultiDeg( leadexp(tt), m, V  ) )
2868        {
2869          mp = mp + lead(tt); // "mp", mp;
2870        }
2871        else
2872        {
2873          t = t + lead(tt);  //  "t", t;
2874        }
2875
2876        tt = tt - lead(tt);
2877      }
2878
2879      I[size(I)+1] = mp;
2880
2881      p = t;
2882    }
2883  }
2884  else
2885  {
2886    I[1] = p; // single monom
2887  }
2888
2889  if( typeof(I) == "module" )
2890  {
2891    I = setModuleGrading(I, V);
2892  }
2893
2894  return (I);
2895}
2896example
2897{
2898  "EXAMPLE:"; echo=2;
2899
2900  ring r = 0,(x,y,z),dp;
2901
2902  intmat g[2][3]=
2903    1,0,1,
2904    0,1,1;
2905  intmat t[2][1]=
2906    -2,
2907    1;
2908
2909  setBaseMultigrading(g,t);
2910
2911  poly f = x10yz+x8y2z-x4z2+y5+x2y2-z2+x17z3-y6;
2912
2913  multiDegPartition(f);
2914
2915  vector v = xy*gen(1)-x3y2*gen(2)+x4y*gen(3);
2916  intmat B[2][3]=1,-1,-2,0,0,1;
2917  v = setModuleGrading(v,B);
2918  getModuleGrading(v);
2919
2920  multiDegPartition(v, B);
2921}
2922
2923
2924
2925/******************************************************/
2926static proc unitMatrix(int n)
2927{
2928  intmat A[n][n];
2929
2930  for( int i = n; i > 0; i-- )
2931  {
2932    A[i,i] = 1;
2933  }
2934
2935  return (A);
2936}
2937
2938
2939
2940/******************************************************/
2941static proc finestMDeg(def r)
2942"
2943USAGE: finestMDeg(r); ring r
2944RETURN: ring, r endowed with the finest multigrading
2945TODO: not yet...
2946"
2947{
2948  def save = basering;
2949  setring (r);
2950
2951  // in basering
2952      ideal I = ideal(basering);
2953
2954  int n = 0; int i; poly p;
2955  for( i = ncols(I); i > 0; i-- )
2956  {
2957    p = I[i];
2958    if( size(p) > 1 )
2959    {
2960      n = n + (size(p) - 1);
2961    }
2962    else
2963    {
2964      I[i] = 0;
2965    }
2966  }
2967
2968  int N = nvars(basering);
2969  intmat A = unitMatrix(N);
2970
2971
2972
2973  if( n > 0)
2974  {
2975
2976    intmat L[N][n];
2977    //  list L;
2978    int j = n;
2979
2980    for(  i = ncols(I); i > 0; i-- )
2981    {
2982      p = I[i];
2983
2984      if( size(p) > 1 )
2985      {
2986        intvec m0 = leadexp(p);
2987        p = p - lead(p);
2988
2989        while( size(p) > 0 )
2990        {
2991          L[ 1..N, j ] = leadexp(p) - m0;
2992          p = p - lead(p);
2993          j--;
2994        }
2995      }
2996    }
2997
2998    print(L);
2999    setBaseMultigrading(A, L);
3000  }
3001  else
3002  {
3003    setBaseMultigrading(A);
3004  }
3005
3006  //  ERROR("nope");
3007
3008  //  ring T = integer, (x), (C, dp);
3009
3010  setring(save);
3011  return (r);
3012}
3013example
3014{
3015  "EXAMPLE:"; echo=2;
3016
3017  ring r = 0,(x, y), dp;
3018  qring q  = std(x^2 - y);
3019
3020  finestMDeg(q);
3021
3022}
3023
3024
3025
3026
3027/******************************************************/
3028static proc newMap(map F, intmat Q, list #)
3029"
3030USAGE: newMap(F, Q[, P]); map F, intmat Q[, intmat P]
3031PURPOSE: endowe the map F with the integer matrices P [and Q]
3032"
3033{
3034  attrib(F, "Q", Q);
3035
3036  if( size(#) > 0 and typeof(#[1]) == "intmat" )
3037  {
3038    attrib(F, "P", #[1]);
3039  }
3040  return (F);
3041}
3042
3043/******************************************************/
3044static proc matrix2intmat( matrix M )
3045{
3046  execute( "intmat A[ "+ string(nrows(M)) + "]["+ string(ncols(M)) + "] = " + string(M) + ";" );
3047  return (A);
3048}
3049
3050
3051/******************************************************/
3052static proc leftKernelZ(intmat M)
3053"USAGE:  leftKernel(M);   M a matrix
3054RETURN:  module
3055PURPOSE: computes left kernel of matrix M (a module of all elements v such that vM=0)
3056EXAMPLE: example leftKernel; shows an example
3057"
3058{
3059  int @bf = 0;
3060  if( nameof(basering) != "basering" )
3061  {
3062    @bf = 1;
3063    def @save@ = basering;
3064  }
3065
3066  ring r = integer, (x), dp;
3067
3068
3069  //  basering;
3070  module N = matrix((M)); // transpose
3071  //  print(N);
3072
3073  def MM = modulo( N, std(0) ) ;
3074  //  print(MM);
3075
3076  intmat R = (  matrix2intmat( MM ) ); // transpose
3077
3078  if( @bf == 1 )
3079  {
3080    setring @save@;
3081  }
3082
3083  kill r;
3084  return( R );
3085}
3086example
3087{
3088  "EXAMPLE:"; echo=2;
3089
3090  ring r= 0,(x,y,z),dp;
3091  matrix M[3][1] = x,y,z;
3092  print(M);
3093  matrix L = leftKernel(M);
3094  print(L);
3095  // check:
3096  print(L*M);
3097};
3098
3099
3100
3101/******************************************************/
3102// the following is taken from "sing4ti2.lib" as we need 'hilbert' from 4ti2
3103
3104static proc hilbert4ti2intmat(intmat A, list #)
3105"USAGE:  hilbert4ti2(A[,i]);
3106@*       A=intmat
3107@*       i=int
3108ASSUME:  - A is a matrix with integer entries which describes the lattice
3109@*         as ker(A), if second argument is not present, and
3110@*         as the left image Im(A) = {zA : z \in ZZ^k}, if second argument is a positive integer
3111@*       - number of variables of basering equals number of columns of A
3112@*         (for ker(A)) resp. of rows of A (for Im(A))
3113CREATE:  temporary files sing4ti2.mat, sing4ti2.lat, sing4ti2.mar
3114@*       in the current directory (I/O files for communication with 4ti2)
3115NOTE:    input rules for 4ti2 also apply to input to this procedure
3116@*       hence ker(A)={x|Ax=0} and Im(A)={xA}
3117RETURN:  toric ideal specified by Hilbert basis thereof
3118EXAMPLE: example graver4ti2; shows an example
3119"
3120{
3121   if( system("sh","which hilbert 2> /dev/null 1> /dev/null") != 0 )
3122   {
3123     ERROR("Sorry: cannot find 'hilbert' command from 4ti2. Please install 4ti2!");
3124   }
3125
3126//--------------------------------------------------------------------------
3127// Initialization and Sanity Checks
3128//--------------------------------------------------------------------------
3129   int i,j;
3130   int nr=nrows(A);
3131   int nc=ncols(A);
3132   string fileending="mat";
3133   if (size(#)!=0)
3134   {
3135//--- default behaviour: use ker(A) as lattice
3136//--- if #[1]!=0 use Im(A) as lattice
3137      if(typeof(#[1])!="int")
3138      {
3139         ERROR("optional parameter needs to be integer value");
3140      }
3141      if(#[1]!=0)
3142      {
3143         fileending="lat";
3144      }
3145   }
3146//--- we should also be checking whether all entries are indeed integers
3147//--- or whether there are fractions, but in this case the error message
3148//--- of 4ti2 is printed directly
3149
3150//--------------------------------------------------------------------------
3151// preparing input file for 4ti2
3152//--------------------------------------------------------------------------
3153   link eing=":w sing4ti2."+fileending;
3154   string eingstring=string(nr)+" "+string(nc);
3155   write(eing,eingstring);
3156   for(i=1;i<=nr;i++)
3157   {
3158      kill eingstring;
3159      string eingstring;
3160      for(j=1;j<=nc;j++)
3161      {
3162        //          if(g(A[i,j])>0)||(char(basering)!=0)||(npars(basering)>0))
3163        //          {
3164        //             ERROR("Input to hilbert4ti2 needs to be a matrix with integer entries");
3165        //          }
3166        eingstring=eingstring+string(A[i,j])+" ";
3167      }
3168      write(eing, eingstring);
3169   }
3170   close(eing);
3171
3172//----------------------------------------------------------------------
3173// calling 4ti2 and converting output
3174// Singular's string is too clumsy for this, hence we first prepare
3175// using standard unix commands
3176//----------------------------------------------------------------------
3177
3178
3179   j=system("sh","hilbert -q -n sing4ti2"); ////////// be quiet + no loggin!!!
3180
3181   j=system("sh", "awk \'BEGIN{ORS=\",\";}{print $0;}\' sing4ti2.hil " +
3182                "| sed s/[\\\ \\\t\\\v\\\f]/,/g " +
3183                "| sed s/,+/,/g|sed s/,,/,/g " +
3184                "| sed s/,,/,/g " +
3185                "> sing4ti2.converted" );
3186
3187
3188//----------------------------------------------------------------------
3189// reading output of 4ti2
3190//----------------------------------------------------------------------
3191   link ausg=":r sing4ti2.converted";
3192//--- last entry ideal(0) is used to tie the list to the basering
3193//--- it will not be processed any further
3194
3195   string s = read(ausg);
3196
3197   if( defined(keepfiles) <= 0)
3198   {
3199      j=system("sh",("rm -f sing4ti2.hil sing4ti2.converted sing4ti2."+fileending));
3200   }
3201
3202   string ergstr = "intvec erglist = " + s + "0;";
3203   execute(ergstr);
3204
3205   //   print(erglist);
3206
3207   int Rnc = erglist[1];
3208   int Rnr = erglist[2];
3209
3210   intmat R[Rnr][Rnc];
3211
3212   int k = 3;
3213
3214   for(i=1;i<=Rnc;i++)
3215   {
3216     for(j=1;j<=Rnr;j++)
3217     {
3218       //       "i: ", i, ", j: ", j, ", v: ", erglist[k];
3219       R[j, i] = erglist[k];
3220       k = k + 1;
3221     }
3222   }
3223
3224
3225
3226   return (R);
3227//--- get rid of leading entry 0;
3228//   toric=toric[2..ncols(toric)];
3229//   return(toric);
3230}
3231// A nice example here is the 3x3 Magic Squares
3232example
3233{
3234  "EXAMPLE:"; echo=2;
3235
3236   ring r=0,(x1,x2,x3,x4,x5,x6,x7,x8,x9),dp;
3237   intmat M[7][9]=
3238      1, 1, 1, -1, -1, -1, 0, 0, 0,
3239      1, 1, 1,  0,  0,  0,-1,-1,-1,
3240      0, 1, 1, -1,  0,  0,-1, 0, 0,
3241      1, 0, 1,  0, -1,  0, 0,-1, 0,
3242      1, 1, 0,  0,  0, -1, 0, 0,-1,
3243      0, 1, 1,  0, -1,  0, 0, 0,-1,
3244      1, 1, 0,  0, -1,  0,-1, 0, 0;
3245   hilbert4ti2intmat(M);
3246   hermiteNormalForm(M);
3247}
3248
3249/////////////////////////////////////////////////////////////////////////////
3250static proc getMonomByExponent(intvec exp)
3251{
3252  int n = nvars(basering);
3253
3254  if( nrows(exp) < n )
3255  {
3256    n = nrows(exp);
3257  }
3258
3259  poly m = 1; int e;
3260
3261  for( int i = 1; i <= n; i++ )
3262  {
3263    e = exp[i];
3264    if( e < 0 )
3265    {
3266      ERROR("Negative exponent!!!");
3267    }
3268
3269    m = m * (var(i)^e);
3270  }
3271
3272  return (m);
3273
3274}
3275
3276/******************************************************/
3277proc multiDegBasis(intvec d)
3278"USAGE: multidegree d
3279ASSUME: current ring is multigraded, monomial ordering is global
3280PURPOSE: compute all monomials of multidegree d
3281EXAMPLE: example multiDegBasis; shows an example
3282"
3283{
3284  def R = basering;  // setring R;
3285
3286  intmat M = getVariableWeights(R);
3287
3288  //  print(M);
3289
3290  int nr = nrows(M);
3291  int nc = ncols(M);
3292
3293  intmat A[nr][nc+1];
3294  A[1..nr, 1..nc] = M[1..nr, 1..nc];
3295  //typeof(A[1..nr, nc+1]);
3296  if( nr==1)
3297  {
3298    A[1..nr, nc+1]=-d[1];
3299  }
3300  else
3301  {
3302    A[1..nr, nc+1] = -d;
3303  }
3304
3305  intmat T = getLattice(R);
3306
3307  if( isFreeRepresented() )
3308  {
3309    intmat B = hilbert4ti2intmat(A);
3310
3311    //      matrix B = unitMatrix(nrows(T));
3312  }
3313  else
3314  {
3315    int n = ncols(T);
3316
3317    nc = ncols(A);
3318
3319    intmat AA[nr][nc + 2 * n];
3320    AA[1..nr, 1.. nc] = A[1..nr, 1.. nc];
3321    AA[1..nr, nc + (1.. n)] = T[1..nr, 1.. n];
3322    AA[1..nr, nc + n + (1.. n)] = -T[1..nr, 1.. n];
3323
3324
3325    //      print ( AA );
3326
3327    intmat K = leftKernelZ(( AA ) ); //
3328
3329    //      print(K);
3330
3331    intmat KK[nc][ncols(K)] = K[ 1.. nc, 1.. ncols(K) ];
3332
3333    //      print(KK);
3334    //      "!";
3335
3336    intmat B = hilbert4ti2intmat(transpose(KK), 1);
3337
3338    //      "!";      print(B);
3339
3340  }
3341
3342
3343  //  print(A);
3344
3345
3346
3347  int i;
3348  int nnr = nrows(B);
3349  int nnc = ncols(B);
3350  ideal I, J;
3351  if(nnc==0){
3352    I=0;
3353    return(I);
3354  }
3355  I[nnc] = 0;
3356  J[nnc] = 0;
3357
3358  for( i = 1; i <= nnc; i++ )
3359  {
3360    //      "i: ", i;    B[nnr, i];
3361
3362    if( B[nnr, i] == 1)
3363    {
3364      // intvec(B[1..nnr-1, i]);
3365      I[i] = getMonomByExponent(intvec(B[1..nnr-1, i]));
3366    }
3367    else
3368    {
3369      if( B[nnr, i] == 0)
3370      {
3371        // intvec(B[1..nnr-1, i]);
3372        J[i] = getMonomByExponent(intvec(B[1..nnr-1, i]));
3373      }
3374    }
3375    //      I[i];
3376  }
3377
3378  ideal Q = (ideal(basering));
3379
3380  if ( size(Q) > 0 )
3381  {
3382    I = NF( I, lead(Q) );
3383    J = NF( J, lead(Q) ); // Global ordering!!!
3384  }
3385
3386  I = simplify(I, 2); // d
3387  J = simplify(J, 2); // d
3388
3389  attrib(I, "ZeroPart", J);
3390
3391  return (I);
3392
3393  //  setring ;
3394}
3395example
3396{
3397  "EXAMPLE:"; echo=2;
3398
3399  ring R = 0, (x, y), dp;
3400
3401  intmat g1[2][2]=1,0,0,1;
3402  intmat l[2][1]=2,0;
3403  intmat g2[2][2]=1,1,1,1;
3404  intvec v1=4,0;
3405  intvec v2=4,4;
3406
3407  intmat g3[1][2]=1,1;
3408  setBaseMultigrading(g3);
3409  intvec v3=4:1;
3410  v3;
3411  multiDegBasis(v3);
3412
3413  setBaseMultigrading(g1,l);
3414  multiDegBasis(v1);
3415  setBaseMultigrading(g2);
3416  multiDegBasis(v2);
3417
3418  intmat M[2][2] = 1, -1, -1, 1;
3419  intvec d = -2, 2;
3420
3421  setBaseMultigrading(M);
3422
3423  multiDegBasis(d);
3424  attrib(_, "ZeroPart");
3425
3426  kill R, M, d;
3427  ring R = 0, (x, y, z), dp;
3428
3429  intmat M[2][3] = 1, -2, 1,     1, 1, 0;
3430
3431  intmat L[2][1] = 0, 2;
3432
3433  intvec d = 4, 1;
3434
3435  setBaseMultigrading(M, L);
3436
3437  multiDegBasis(d);
3438  attrib(_, "ZeroPart");
3439
3440
3441  kill R, M, d;
3442
3443  ring R = 0, (x, y, z), dp;
3444  qring Q = std(ideal( y^6+ x*y^3*z-x^2*z^2 ));
3445
3446
3447  intmat M[2][3] = 1, 1, 2,     2, 1, 1;
3448  //  intmat T[2][1] = 0, 2;
3449
3450  setBaseMultigrading(M); // BUG????
3451
3452  intvec d = 6, 6;
3453  multiDegBasis(d);
3454  attrib(_, "ZeroPart");
3455
3456
3457
3458  kill R, Q, M, d;
3459  ring R = 0, (x, y, z), dp;
3460  qring Q = std(ideal( x*z^3 - y *z^6, x*y*- x^4*y^2 ));
3461
3462
3463  intmat M[2][3] = 1, -2, 1,     1, 1, 0;
3464  intmat T[2][1] = 0, 2;
3465
3466  intvec d = 4, 1;
3467
3468  setBaseMultigrading(M, T); // BUG????
3469
3470  multiDegBasis(d);
3471  attrib(_, "ZeroPart");
3472}
3473
3474
3475proc multiDegSyzygy(def I)
3476"USAGE: multiDegSyzygy(I); I is a ideal or a module
3477PURPOSE: computes the multigraded syzygy module of I
3478RETURNS: module, the syzygy module of I
3479NOTE: generators of I must be multigraded homogeneous
3480EXAMPLE: example multiDegSyzygy; shows an example
3481"
3482{
3483  if( isHomogeneous(I, "checkGens") == 0)
3484  {
3485    ERROR ("Sorry: inhomogeneous input!");
3486  }
3487  module S = syz(I);
3488  S = setModuleGrading(S, multiDeg(I));
3489  return (S);
3490}
3491example
3492{
3493  "EXAMPLE:"; echo=2;
3494
3495  ring r = 0,(x,y,z,w),dp;
3496  intmat MM[2][4]=
3497    1,1,1,1,
3498    0,1,3,4;
3499  setBaseMultigrading(MM);
3500  module M = ideal(  xw-yz, x2z-y3, xz2-y2w, yw2-z3);
3501
3502
3503  intmat v[2][nrows(M)]=
3504    1,
3505    0;
3506
3507  M = setModuleGrading(M, v);
3508
3509  isHomogeneous(M);
3510  "Multidegrees: "; print(multiDeg(M));
3511  // Let's compute syzygies!
3512  def S = multiDegSyzygy(M); S;
3513  "Module Units Multigrading: "; print( getModuleGrading(S) );
3514  "Multidegrees: "; print(multiDeg(S));
3515
3516  isHomogeneous(S);
3517}
3518
3519
3520
3521proc multiDegModulo(def I, def J)
3522"USAGE: multiDegModulo(I); I, J are ideals or modules
3523PURPOSE: computes the multigraded 'modulo' module of I and J
3524RETURNS: module, see 'modulo' command
3525NOTE: I and J should have the same multigrading, and their
3526generators must be multigraded homogeneous
3527EXAMPLE: example multiDegModulo; shows an example
3528"
3529{
3530  if( (isHomogeneous(I, "checkGens") == 0) or (isHomogeneous(J, "checkGens") == 0) )
3531  {
3532    ERROR ("Sorry: inhomogeneous input!");
3533  }
3534  module K = modulo(I, J);
3535  K = setModuleGrading(K, multiDeg(I));
3536  return (K);
3537}
3538example
3539{
3540  "EXAMPLE:"; echo=2;
3541
3542  ring r = 0,(x,y,z),dp;
3543  intmat MM[2][3]=
3544    -1,1,1,
3545     0,1,3;
3546  setBaseMultigrading(MM);
3547
3548  ideal h1 = x, y, z;
3549  ideal h2 = x;
3550
3551  "Multidegrees: "; print(multiDeg(h1));
3552
3553  // Let's compute modulo(h1, h2):
3554  def K = multiDegModulo(h1, h2); K;
3555
3556  "Module Units Multigrading: "; print( getModuleGrading(K) );
3557  "Multidegrees: "; print(multiDeg(K));
3558
3559  isHomogeneous(K);
3560}
3561
3562
3563proc multiDegGroebner(def I)
3564"USAGE: multiDegGroebner(I); I is a poly/vector/ideal/module
3565PURPOSE: computes the multigraded standard/groebner basis of I
3566NOTE: I must be multigraded homogeneous
3567RETURNS: ideal/module, the computed basis
3568EXAMPLE: example multiDegGroebner; shows an example
3569"
3570{
3571  if( isHomogeneous(I) == 0)
3572  {
3573    ERROR ("Sorry: inhomogeneous input!");
3574  }
3575
3576  def S = groebner(I);
3577
3578  if( typeof(I) == "module" or typeof(I) == "vector" )
3579  {
3580    S = setModuleGrading(S, getModuleGrading(I));
3581  }
3582
3583  return(S);
3584}
3585example
3586{
3587  "EXAMPLE:"; echo=2;
3588
3589  ring r = 0,(x,y,z,w),dp;
3590
3591  intmat MM[2][4]=
3592    1,1,1,1,
3593    0,1,3,4;
3594
3595  setBaseMultigrading(MM);
3596
3597
3598  module M = ideal(  xw-yz, x2z-y3, xz2-y2w, yw2-z3);
3599
3600
3601  intmat v[2][nrows(M)]=
3602    1,
3603    0;
3604
3605  M = setModuleGrading(M, v);
3606
3607
3608  /////////////////////////////////////////////////////////////////////////////
3609  // GB:
3610  M = multiDegGroebner(M); M;
3611  "Module Units Multigrading: "; print( getModuleGrading(M) );
3612  "Multidegrees: "; print(multiDeg(M));
3613
3614  isHomogeneous(M);
3615
3616  /////////////////////////////////////////////////////////////////////////////
3617  // Let's compute Syzygy!
3618  def S = multiDegSyzygy(M); S;
3619  "Module Units Multigrading: "; print( getModuleGrading(S) );
3620  "Multidegrees: "; print(multiDeg(S));
3621
3622  isHomogeneous(S);
3623
3624  /////////////////////////////////////////////////////////////////////////////
3625  // GB:
3626  S = multiDegGroebner(S); S;
3627  "Module Units Multigrading: "; print( getModuleGrading(S) );
3628  "Multidegrees: "; print(multiDeg(S));
3629
3630  isHomogeneous(S);
3631}
3632
3633
3634/******************************************************/
3635proc multiDegResolution(def I, int ll, list #)
3636"USAGE: multiDegResolution(I,l,[f]); I is poly/vector/ideal/module; l,f are integers
3637PURPOSE: computes the multigraded resolution of I of the length l,
3638or the whole resolution if l is zero. Returns minimal resolution if an optional
3639argument 1 is supplied
3640NOTE: input must have multigraded-homogeneous generators.
3641The returned list is truncated beginning with the first zero differential.
3642RETURNS: list, the computed resolution
3643EXAMPLE: example multiDegResolution; shows an example
3644"
3645{
3646  if( isHomogeneous(I, "checkGens") == 0)
3647  {
3648    ERROR ("Sorry: inhomogeneous input!");
3649  }
3650
3651  def R = res(I, ll, #); list L = R; int l = size(L);
3652  def V = getModuleGrading(I);
3653  if( (typeof(I) == "module") or (typeof(I) == "vector") )
3654  {
3655    L[1] = setModuleGrading(L[1], V);
3656  }
3657
3658  int i;
3659  for( i = 2; i <= l; i++ )
3660  {
3661    if( size(L[i]) > 0 )
3662    {
3663      L[i] = setModuleGrading( L[i], multiDeg(L[i-1]) );
3664    } else
3665    {
3666      return (L[1..(i-1)]);
3667    }
3668  }
3669
3670  return (L);
3671
3672
3673}
3674example
3675{
3676  "EXAMPLE:"; echo=2;
3677
3678  ring r = 0,(x,y,z,w),dp;
3679
3680  intmat M[2][4]=
3681    1,1,1,1,
3682    0,1,3,4;
3683
3684  setBaseMultigrading(M);
3685
3686
3687  module m= ideal(  xw-yz, x2z-y3, xz2-y2w, yw2-z3);
3688
3689  isHomogeneous(ideal(  xw-yz, x2z-y3, xz2-y2w, yw2-z3), "checkGens");
3690
3691  ideal A = xw-yz, x2z-y3, xz2-y2w, yw2-z3;
3692
3693  int j;
3694
3695  for(j=1; j<=ncols(A); j++)
3696  {
3697    multiDegPartition(A[j]);
3698  }
3699
3700  intmat v[2][1]=
3701    1,
3702    0;
3703
3704  m = setModuleGrading(m, v);
3705
3706  // Let's compute Syzygy!
3707  def S = multiDegSyzygy(m); S;
3708  "Module Units Multigrading: "; print( getModuleGrading(S) );
3709  "Multidegrees: "; print(multiDeg(S));
3710
3711  /////////////////////////////////////////////////////////////////////////////
3712
3713  S = multiDegGroebner(S); S;
3714  "Module Units Multigrading: "; print( getModuleGrading(S) );
3715  "Multidegrees: "; print(multiDeg(S));
3716
3717  /////////////////////////////////////////////////////////////////////////////
3718
3719  def L = multiDegResolution(m, 0, 1);
3720
3721  for( j =1; j<=size(L); j++)
3722  {
3723    "----------------------------------- ", j, " -----------------------------";
3724    L[j];
3725    "Module Multigrading: "; print( getModuleGrading(L[j]) );
3726    "Multigrading: "; print(multiDeg(L[j]));
3727  }
3728
3729  /////////////////////////////////////////////////////////////////////////////
3730
3731  L = multiDegResolution(maxideal(1), 0, 1);
3732
3733  for( j =1; j<=size(L); j++)
3734  {
3735    "----------------------------------- ", j, " -----------------------------";
3736    L[j];
3737    "Module Multigrading: "; print( getModuleGrading(L[j]) );
3738    "Multigrading: "; print(multiDeg(L[j]));
3739  }
3740
3741  kill v;
3742
3743
3744  def h = hilbertSeries(m);
3745  setring h;
3746
3747  numerator1;
3748  factorize(numerator1);
3749
3750  denominator1;
3751  factorize(denominator1);
3752
3753  numerator2;
3754  factorize(numerator2);
3755
3756  denominator2;
3757  factorize(denominator2);
3758}
3759
3760/******************************************************/
3761proc hilbertSeries(def I)
3762"USAGE: hilbertSeries(I); I is poly/vector/ideal/module
3763PURPOSE: computes the multigraded Hilbert Series of I
3764NOTE: input must have multigraded-homogeneous generators.
3765Multigrading should be positive.
3766RETURNS: a ring in variables t_(i), s_(i), with polynomials
3767numerator1 and denominator1 and mutually prime numerator2
3768and denominator2, quotients of which give the series.
3769EXAMPLE: example hilbertSeries; shows an example
3770"
3771{
3772
3773  if( !isFreeRepresented() )
3774  {
3775    "Things might happen, since we are not  free.";
3776    //ERROR("SORRY: ONLY TORSION-FREE CASE (POSITIVE GRADING)");
3777  }
3778
3779  int i, j, k, v;
3780
3781  intmat M = getVariableWeights();
3782
3783  int cc = ncols(M);
3784  int n = nrows(M);
3785
3786  if( n == 0 )
3787  {
3788    ERROR("Error: wrong Variable Weights?");
3789  }
3790
3791  list RES = multiDegResolution(I,0,1);
3792
3793  int l = size(RES);
3794
3795  list L; L[l + 1] = 0;
3796
3797  if(typeof(I) == "ideal")
3798  {
3799    intmat zeros[n][1];
3800    L[1] = zeros;
3801  }
3802  else
3803  {
3804    L[1] = getModuleGrading(RES[1]);
3805  }
3806
3807  for( j = 1; j <= l; j++)
3808  {
3809    L[j + 1] = multiDeg(RES[j]);
3810  }
3811
3812  l++;
3813
3814  ring R = 0,(t_(1..n),s_(1..n)),dp;
3815
3816  ideal units;
3817  for( i=n; i>=1; i--)
3818  {
3819    units[i] = (var(i) * var(n + i) - 1);
3820  }
3821
3822  qring Q = std(units);
3823
3824  // TODO: should not it be a quotient ring depending on Torsion???
3825  // I am not sure about what to do in the torsion case, but since
3826  // we want to evaluate the polynomial at certain points to get
3827  // a dimension we need uniqueness for this. I think we would lose
3828  // this uniqueness if switching to this torsion ring.
3829
3830  poly monom, summand, @numerator;
3831  poly @denominator = 1;
3832
3833  for( i = 1; i <= cc; i++)
3834  {
3835    monom = 1;
3836    for( k = 1; k <= n; k++)
3837    {
3838      v = M[k,i];
3839
3840      if(v >= 0)
3841      {
3842        monom = monom * (var(k)^(v));
3843      }
3844      else
3845      {
3846        monom = monom * (var(n+k)^(-v));
3847      }
3848    }
3849
3850    if( monom == 1)
3851    {
3852      ERROR("Multigrading not positive.");
3853    }
3854
3855    @denominator = @denominator * (1 - monom);
3856  }
3857
3858  for( j = 1; j<= l; j++)
3859  {
3860    summand = 0;
3861    M = L[j];
3862
3863    for( i = 1; i <= ncols(M); i++)
3864    {
3865      monom = 1;
3866      for( k = 1; k <= n; k++)
3867      {
3868        v = M[k,i];
3869        if( v > 0 )
3870        {
3871          monom = monom * (var(k)^v);
3872        }
3873        else
3874        {
3875          monom = monom * (var(n+k)^(-v));
3876        }
3877      }
3878      summand = summand + monom;
3879    }
3880    @numerator = @numerator - (-1)^j * summand;
3881  }
3882
3883  if( @denominator == 0 )
3884  {
3885    ERROR("Multigrading not positive.");
3886  }
3887
3888  poly denominator1 = @denominator;
3889  poly numerator1 = @numerator;
3890
3891  export denominator1;
3892  export numerator1;
3893
3894  if( numerator != 0 )
3895  {
3896    poly d = gcd(@denominator, @numerator);
3897
3898    poly denominator2 = @denominator/d;
3899    poly numerator2 = @numerator/d;
3900
3901    if( gcd(denominator2, numerator2) != 1 )
3902    {
3903      ERROR("Sorry: gcd should be 1 (after dividing out gcd)! Something went wrong!");
3904    }
3905  }
3906  else
3907  {
3908    poly denominator2 = @denominator;
3909    poly numerator2 = @numerator;
3910  }
3911
3912
3913  export denominator2;
3914  export numerator2;
3915
3916  " ------------ ";
3917  "This proc returns a ring with polynomials called 'numerator1/2' and 'denominator1/2'!";
3918  "They represent the first and the second Hilbert Series.";
3919  "The s_(i)-variables are defined to be the inverse of the t_(i)-variables.";
3920  " ------------ ";
3921
3922  return(Q);
3923}
3924example
3925{
3926  "EXAMPLE:"; echo=2;
3927
3928  ring r = 0,(x,y,z,w),dp;
3929  intmat g[2][4]=
3930    1,1,1,1,
3931    0,1,3,4;
3932  setBaseMultigrading(g);
3933
3934  module M = ideal(xw-yz, x2z-y3, xz2-y2w, yw2-z3);
3935  intmat V[2][1]=
3936    1,
3937    0;
3938
3939  M = setModuleGrading(M, V);
3940
3941  def h = hilbertSeries(M); setring h;
3942
3943  factorize(numerator2);
3944  factorize(denominator2);
3945
3946  kill g, h; setring r;
3947
3948  intmat g[2][4]=
3949    1,2,3,4,
3950    0,0,5,8;
3951
3952  setBaseMultigrading(g);
3953
3954  ideal I = x^2, y, z^3;
3955  I = std(I);
3956  def L = multiDegResolution(I, 0, 1);
3957
3958  for( int j = 1; j<=size(L); j++)
3959  {
3960    "----------------------------------- ", j, " -----------------------------";
3961    L[j];
3962    "Module Multigrading: "; print( getModuleGrading(L[j]) );
3963    "Multigrading: "; print(multiDeg(L[j]));
3964  }
3965
3966  multiDeg(I);
3967  def h = hilbertSeries(I); setring h;
3968
3969  factorize(numerator2);
3970  factorize(denominator2);
3971
3972  kill r, h, g, V;
3973  ////////////////////////////////////////////////
3974  ring R = 0,(x,y,z),dp;
3975  intmat W[2][3] =
3976     1,1, 1,
3977     0,0,-1;
3978  setBaseMultigrading(W);
3979  ideal I = x3y,yz2,y2z,z4;
3980
3981  def h = hilbertSeries(I); setring h;
3982
3983  factorize(numerator2);
3984  factorize(denominator2);
3985
3986  kill R, W, h;
3987  ////////////////////////////////////////////////
3988  ring R = 0,(x,y,z,a,b,c),dp;
3989  intmat W[2][6] =
3990     1,1, 1,1,1,1,
3991     0,0,-1,0,0,0;
3992  setBaseMultigrading(W);
3993  ideal I = x3y,yz2,y2z,z4;
3994
3995  def h = hilbertSeries(I); setring h;
3996
3997  factorize(numerator2);
3998  factorize(denominator2);
3999
4000  kill R, W, h;
4001  ////////////////////////////////////////////////
4002  // This is example 5.3.9. from Robbianos book.
4003
4004  ring R = 0,(x,y,z,w),dp;
4005  intmat W[1][4] =
4006     1,1, 1,1;
4007  setBaseMultigrading(W);
4008  ideal I = z3,y3zw2,x2y4w2xyz2;
4009
4010  hilb(std(I));
4011
4012  def h = hilbertSeries(I); setring h;
4013
4014  numerator1;
4015  denominator1;
4016
4017  factorize(numerator2);
4018  factorize(denominator2);
4019
4020
4021  kill h;
4022  ////////////////////////////////////////////////
4023  setring R;
4024
4025  ideal I2 = x2,y2,z2; I2;
4026
4027  hilb(std(I2));
4028
4029  def h = hilbertSeries(I2); setring h;
4030
4031  numerator1;
4032  denominator1;
4033
4034
4035  kill h;
4036  ////////////////////////////////////////////////
4037  setring R;
4038
4039  W = 2,2,2,2;
4040
4041  setBaseMultigrading(W);
4042
4043  getVariableWeights();
4044
4045  intvec w = 2,2,2,2;
4046
4047  hilb(std(I2), 1, w);
4048
4049  kill w;
4050
4051
4052  def h = hilbertSeries(I2); setring h;
4053
4054
4055  numerator1; denominator1;
4056  kill h;
4057
4058
4059  kill R, W;
4060
4061  ////////////////////////////////////////////////
4062  ring R = 0,(x),dp;
4063  intmat W[1][1] =
4064     1;
4065  setBaseMultigrading(W);
4066
4067  ideal I;
4068
4069  I = 1; I;
4070
4071  hilb(std(I));
4072
4073  def h = hilbertSeries(I); setring h;
4074
4075  numerator1; denominator1;
4076
4077  kill h;
4078  ////////////////////////////////////////////////
4079  setring R;
4080
4081  I = x; I;
4082
4083  hilb(std(I));
4084
4085  def h = hilbertSeries(I); setring h;
4086
4087  numerator1; denominator1;
4088
4089  kill h;
4090  ////////////////////////////////////////////////
4091  setring R;
4092
4093  I = x^5; I;
4094
4095  hilb(std(I));
4096  hilb(std(I), 1);
4097
4098  def h = hilbertSeries(I); setring h;
4099
4100  numerator1; denominator1;
4101
4102
4103  kill h;
4104  ////////////////////////////////////////////////
4105  setring R;
4106
4107  I = x^10; I;
4108
4109  hilb(std(I));
4110
4111  def h = hilbertSeries(I); setring h;
4112
4113  numerator1; denominator1;
4114
4115  kill h;
4116  ////////////////////////////////////////////////
4117  setring R;
4118
4119  module M = 1;
4120
4121  M = setModuleGrading(M, W);
4122
4123
4124  hilb(std(M));
4125
4126  def h = hilbertSeries(M); setring h;
4127
4128  numerator1; denominator1;
4129
4130  kill h;
4131  ////////////////////////////////////////////////
4132  setring R;
4133
4134  kill M; module M = x^5*gen(1);
4135//  intmat V[1][3] = 0; // TODO: this would lead to a wrong result!!!?
4136  intmat V[1][1] = 0; // all gen(i) of degree 0!
4137
4138  M = setModuleGrading(M, V);
4139
4140  hilb(std(M));
4141
4142  def h = hilbertSeries(M); setring h;
4143
4144  numerator1; denominator1;
4145
4146  kill h;
4147  ////////////////////////////////////////////////
4148  setring R;
4149
4150  module N = x^5*gen(3);
4151
4152  kill V;
4153
4154  intmat V[1][3] = 0; // all gen(i) of degree 0!
4155
4156  N = setModuleGrading(N, V);
4157
4158  hilb(std(N));
4159
4160  def h = hilbertSeries(N); setring h;
4161
4162  numerator1; denominator1;
4163
4164  kill h;
4165  ////////////////////////////////////////////////
4166  setring R;
4167
4168
4169  module S = M + N;
4170
4171  S = setModuleGrading(S, V);
4172
4173  hilb(std(S));
4174
4175  def h = hilbertSeries(S); setring h;
4176
4177  numerator1; denominator1;
4178
4179  kill h;
4180
4181  kill V;
4182  kill R, W;
4183
4184}
4185
4186static proc evalHilbertSeries(def h, intvec v)
4187"
4188   TODO
4189   evaluate hilbert series h by substibuting v[i] for t_(i) (1/v[i] for s_(i))
4190  return: int (h(v))
4191"
4192{
4193  if( 2*size(v) != nvars(h) )
4194  {
4195    ERROR("Wrong input/size!");
4196  }
4197
4198  setring h;
4199
4200  if( defined(numerator2) and defined(denominator2) )
4201  {
4202    poly n = numerator2; poly d = denominator2;
4203  } else
4204  {
4205    poly n = numerator1; poly d = denominator1;
4206  }
4207
4208  int N = size(v);
4209  int i; number k;
4210  ideal V;
4211
4212  for( i = N; i > 0; i -- )
4213  {
4214    k = v[i];
4215    V[i] = var(i) - k;
4216  }
4217
4218  V = groebner(V);
4219
4220  n = NF(n, V);
4221  d = NF(d, V);
4222
4223  n;
4224  d;
4225
4226  if( d == 0 )
4227  {
4228    ERROR("Sorry: denominator is zero!");
4229  }
4230
4231  if( n == 0 )
4232  {
4233    return (0);
4234  }
4235
4236  poly g = gcd(n, d);
4237
4238  if( g != leadcoef(g) )
4239  {
4240    n = n / g;
4241    d = d / g;
4242  }
4243
4244  n;
4245  d;
4246
4247
4248  for( i = N; i > 0; i -- )
4249  {
4250    "i: ", i;
4251    n;
4252    d;
4253
4254    k = v[i];
4255    k;
4256
4257    n = subst(n, var(i), k);
4258    d = subst(d, var(i), k);
4259
4260    if( k != 0 )
4261    {
4262      k = 1/k;
4263      n = subst(n, var(N+i), k);
4264      d = subst(d, var(N+i), k);
4265    }
4266  }
4267
4268  n;
4269  d;
4270
4271  if( d == 0 )
4272  {
4273    ERROR("Sorry: denominator is zero!");
4274  }
4275
4276  if( n == 0 )
4277  {
4278    return (0);
4279  }
4280
4281  poly g = gcd(n, d);
4282
4283  if( g != leadcoef(g) )
4284  {
4285    n = n / g;
4286    d = d / g;
4287  }
4288
4289  n;
4290  d;
4291
4292  if( n != leadcoef(n) || d != leadcoef(d) )
4293  {
4294    ERROR("Sorry cannot completely evaluate. Partial result: (" + string(n) + ")/(" + string(d) + ")");
4295  }
4296
4297  n;
4298  d;
4299
4300  return (leadcoef(n)/leadcoef(d));
4301}
4302example
4303{
4304  "EXAMPLE:"; echo=2;
4305
4306  // TODO!
4307
4308}
4309
4310
4311proc isPositive()
4312"USAGE: isPositive()
4313PURPOSE: Computes whether the multigrading of the ring is positive.
4314For computation theorem 8.6 of the Miller/Sturmfels book is used.
4315RETURNS: true if the multigrading is positive
4316EXAMPLE: example isPositive; shows an example
4317"
4318{
4319ideal I = multiDegBasis(0);
4320ideal J = attrib(I,"ZeroPart");
4321/*
4322I am not quite sure what this ZeroPart is anymore. I thought it
4323should contain all monomials of degree 0, but then apparently 1 should
4324be contained. It makes sense to exclude 1, but was this also the intention?
4325*/
4326return(J==0);
4327}
4328example
4329{
4330  echo = 2; printlevel = 3;
4331  ring r = 0,(x,y),dp;
4332  intmat A[1][2]=-1,1;
4333  setBaseMultigrading(A);
4334  isPositive();
4335
4336  intmat B[1][2]=1,1;
4337  setBaseMultigrading(B);
4338  isPositive(B);
4339}
4340
4341///////////////////////////////////////////////////////////////////////////////
4342// testing for consistency of the library:
4343proc testMultigradingLib ()
4344{
4345  example setBaseMultigrading;
4346  example setModuleGrading;
4347
4348  example getVariableWeights;
4349  example getLattice;
4350  example getGradingGroup;
4351  example getModuleGrading;
4352
4353
4354  example multiDeg;
4355  example multiDegPartition;
4356
4357
4358  example hermiteNormalForm;
4359  example isHomogeneous;
4360  example isTorsionFree;
4361  example pushForward;
4362  example defineHomogeneous;
4363
4364  example equalMultiDeg;
4365  example isZeroElement;
4366
4367  example multiDegResolution;
4368
4369  "// ******************* example hilbertSeries ************************//";
4370  example hilbertSeries;
4371
4372
4373// example multiDegBasis; // needs 4ti2!
4374
4375  "The End!";
4376}
4377
4378
4379static proc multiDegTruncate(def M, intvec md)
4380{
4381  "d: ";
4382  print(md);
4383
4384  "M: ";
4385  module LL = M; // + L for d+1
4386  LL;
4387  print(multiDeg(LL));
4388
4389
4390  intmat V = getModuleGrading(M);
4391  intvec vi;
4392  int s = nrows(M);
4393  int r = nrows(V);
4394  int i;
4395  module L; def B;
4396  for (i=s; i>0; i--)
4397  {
4398    "comp: ", i;
4399    vi = V[1..r, i];
4400    "v[i]: "; vi;
4401
4402    B = multiDegBasis(md - vi); // ZeroPart is always the same...
4403    "B: "; B;
4404
4405    L = L, B*gen(i);
4406  }
4407  L = simplify(L, 2);
4408  L = setModuleGrading(L,V);
4409
4410  "L: "; L;
4411  print(multiDeg(L));
4412
4413  L = multiDegModulo(L, LL);
4414  L = multiDegGroebner(L);
4415//  L = minbase(prune(L));
4416
4417  "??????????";
4418  print(L);
4419  print(multiDeg(L));
4420
4421  V = getModuleGrading(L);
4422
4423  // take out other degrees
4424  for(i = ncols(L); i > 0; i-- )
4425  {
4426    if( !equalMultiDeg( multiDeg(getGradedGenerator(L, i)), md ) )
4427    {
4428      L[i] = 0;
4429    }
4430  }
4431
4432  L = simplify(L, 2);
4433  L = setModuleGrading(L, V);
4434  print(L);
4435  print(multiDeg(L));
4436
4437  return(L);
4438}
4439example
4440{
4441  "EXAMPLE:"; echo=2;
4442
4443  // TODO!
4444  ring r = 32003, (x,y), dp;
4445
4446  intmat M[2][2] =
4447    1, 0,
4448    0, 1;
4449
4450  setBaseMultigrading(M);
4451
4452  intmat V[2][1] =
4453    0,
4454    0;
4455
4456  "X:";
4457  module h1 = x;
4458  h1 = setModuleGrading(h1, V);
4459  multiDegTruncate(h1, multiDeg(x));
4460  multiDegTruncate(h1, multiDeg(y));
4461
4462  "XY:";
4463  module h2 = ideal(x, y);
4464  h2 = setModuleGrading(h2, V);
4465  multiDegTruncate(h2, multiDeg(x));
4466  multiDegTruncate(h2, multiDeg(y));
4467  multiDegTruncate(h2, multiDeg(xy));
4468}
4469
4470
4471/******************************************************/
4472/* Some functions on lattices.
4473TODO Tuebingen: - add functionality (see wiki) and
4474- adjust them to work for groups as well.*/
4475/******************************************************/
4476
4477
4478
4479/******************************************************/
4480proc imageLattice(intmat Q, intmat L)
4481"USAGE: imageLattice(Q,L); Q and L are of type intmat
4482PURPOSE: compute an integral basis for the image of the
4483lattice L under the homomorphism of lattices Q.
4484RETURN: intmat
4485EXAMPLE: example imageLattice; shows an example
4486"
4487{
4488  intmat Mul = Q*L;
4489  intmat LL = latticeBasis(Mul);
4490
4491  return(LL);
4492}
4493example
4494{
4495  "EXAMPLE:"; echo=2;
4496
4497  intmat Q[2][3] =
4498    1,2,3,
4499    3,2,1;
4500
4501  intmat L[3][2] =
4502    1,4,
4503    2,5,
4504    3,6;
4505
4506  // should be a 2x2 matrix with columns
4507  // [2,-14], [0,36]
4508  imageLattice(Q,L);
4509
4510}
4511
4512/******************************************************/
4513proc intRank(intmat A)
4514"USAGE: intRank(A); intmat A
4515PURPOSE: compute the rank of the integral matrix A
4516by computing a hermite normalform.
4517RETURNS: int
4518EXAMPLE: example intRank; shows an example
4519"
4520{
4521  intmat B = hermiteNormalForm(A);
4522
4523  // get number of zero columns
4524  int nzerocols = 0;
4525  int j;
4526  int i;
4527  int iszero;
4528  for ( j = 1; j <= ncols(B); j++ )
4529  {
4530    iszero = 1;
4531
4532    for ( i = 1; i <= nrows(B); i++ )
4533    {
4534      if ( B[i,j] != 0 )
4535      {
4536        iszero = 0;
4537        break;
4538      }
4539    }
4540
4541    if ( iszero == 1 )
4542    {
4543      nzerocols++;
4544    }
4545  }
4546
4547  // get number of zero rows
4548  int nzerorows = 0;
4549
4550  for ( i = 1; i <= nrows(B); i++ )
4551  {
4552    iszero = 1;
4553
4554    for ( j = 1; j <= ncols(B); j++ )
4555    {
4556      if ( B[i,j] != 0 )
4557      {
4558        iszero = 0;
4559        break;
4560      }
4561    }
4562
4563    if ( iszero == 1 )
4564    {
4565      nzerorows++;
4566    }
4567  }
4568
4569  int r = nrows(B) - nzerorows;
4570
4571  if ( ncols(B) - nzerocols < r )
4572  {
4573    r = ncols(B) - nzerocols;
4574  }
4575
4576  return(r);
4577}
4578example
4579{
4580  "EXAMPLE:"; echo=2;
4581
4582  intmat A[3][4] =
4583    1,0,1,0,
4584    1,2,0,0,
4585    0,0,0,0;
4586
4587  int r = intRank(A);
4588
4589  print(A);
4590  print(r); // Should be 2
4591
4592  // another example
4593  intmat B[2][2] =
4594    1,2,
4595    1,2;
4596
4597  int d = intRank(B);
4598
4599  print(B);
4600  print(d); // Should be 1
4601
4602  kill A, B, r, d;
4603
4604}
4605
4606/*****************************************************/
4607
4608proc isSublattice(intmat L, intmat S)
4609"USAGE: isSublattice(L, S); L, S are of tpye intmat
4610PURPOSE: checks whether the lattice created by L is a
4611sublattice of the lattice created by S.
4612The procedure checks whether each generator of L is
4613contained in S.
4614RETURN: integer, 0 if false, 1 if true
4615EXAMPLE: example isSublattice; shows an example
4616"
4617{
4618  int a,b,g,i,j,k;
4619  intmat Ker;
4620
4621  // check whether each column v of L is contained in
4622  // the lattice generated by S
4623  for ( i = 1; i <= ncols(L); i++ )
4624  {
4625
4626    // v is the i-th column of L
4627    intvec v;
4628     for ( j = 1; j <= nrows(L); j++ )
4629    {
4630      v[j] = L[j,i];
4631    }
4632
4633    // concatenate B = [S,v]
4634    intmat B[nrows(L)][ncols(S) + 1];
4635
4636    for ( a = 1; a <= nrows(S); a++ )
4637    {
4638      for ( b = 1; b <= ncols(S); b++ )
4639      {
4640        B[a,b] = S[a,b];
4641      }
4642    }
4643
4644    for ( a = 1; a <= size(v); a++ )
4645    {
4646      B[a,ncols(B)] = v[a];
4647    }
4648
4649
4650    // check gcd
4651    Ker = kernelLattice(B);
4652    k = nrows(Ker);
4653    list R; // R is the last row
4654
4655    for ( j = 1; j <= ncols(Ker); j++ )
4656    {
4657      R[j] = Ker[k,j];
4658    }
4659
4660    g = R[1];
4661
4662    for ( j = 2; j <= size(R); j++ )
4663    {
4664      g = gcd(g,R[j]);
4665    }
4666
4667    if ( g != 1 and g != -1 )
4668    {
4669      return(0);
4670    }
4671
4672    kill B, v, R;
4673
4674  }
4675
4676  return(1);
4677}
4678example
4679{
4680  "EXAMPLE:"; echo=2;
4681
4682  //ring R = 0,(x,y),dp;
4683  intmat S2[3][3]=
4684    0, 2, 3,
4685    0, 1, 1,
4686    3, 0, 2;
4687
4688  intmat S1[3][2]=
4689    0, 6,
4690    0, 2,
4691    3, 4;
4692
4693  isSublattice(S1,S2); // Yes!
4694
4695  intmat S3[3][1] =
4696    0,
4697    0,
4698    1;
4699
4700  not(isSublattice(S3,S2)); // Yes!
4701
4702}
4703
4704/******************************************************/
4705
4706proc latticeBasis(intmat B)
4707"USAGE: latticeBasis(B); intmat B
4708PURPOSE: compute an integral basis for the lattice defined by
4709the columns of B.
4710RETURNS: intmat
4711EXAMPLE: example latticeBasis; shows an example
4712"
4713{
4714  int n = ncols(B);
4715  int r = intRank(B);
4716
4717  if ( r == 0 )
4718  {
4719    intmat H[nrows(B)][1];
4720    int j;
4721
4722    for ( j = 1; j <= nrows(B); j++ )
4723    {
4724      H[j,1] = 0;
4725    }
4726  }
4727  else
4728  {
4729    intmat H = hermiteNormalForm(B);;
4730
4731    if (r < n)
4732    {
4733      // delete columns r+1 to n
4734      // should be identical with the function
4735      // H = submat(H,1..nrows(H),1..r);
4736      // for matrices
4737      intmat Hdel[nrows(H)][r];
4738      int k;
4739      int m;
4740
4741      for ( k = 1; k <= nrows(H); k++ )
4742      {
4743        for ( m = 1; m <= r; m++ )
4744        {
4745          Hdel[k,m] = H[k,m];
4746        }
4747      }
4748
4749      H = Hdel;
4750    }
4751  }
4752
4753  return(H);
4754}
4755example
4756{
4757  "EXAMPLE:"; echo=2;
4758
4759  intmat L[3][3] =
4760    1,4,8,
4761    2,5,10,
4762    3,6,12;
4763
4764  intmat B = latticeBasis(L); 
4765  print(B); // should result in a matrix whose columns generate the same lattice as  [1,2,3] and [0,3,6]:
4766
4767   // another example
4768  intmat C[2][4] =
4769    1,1,2,0,
4770    2,3,4,0;
4771
4772  // should result in a matrix whose
4773  // colums create the same  lattice as
4774  // [0,1],[1,0]
4775  intmat D = latticeBasis(C);
4776
4777  print(D);
4778
4779  kill B,L;
4780}
4781
4782/******************************************************/
4783
4784proc kernelLattice(def P)
4785"USAGE: kernelLattice(P); intmat P
4786PURPOSE: compute a integral basis for the kernel of the
4787homomorphism of lattices defined by the intmat P.
4788RETURNS: intmat
4789EXAMPLE: example kernelLattice; shows an example
4790"
4791{
4792  int n = ncols(P);
4793  int r = intRank(P);
4794
4795  if ( r == 0 )
4796  {
4797    intmat U = unitMatrix(n);
4798  }
4799  else
4800  {
4801    if ( r == n )
4802    {
4803      intmat U[n][1];  // now all entries are zero.
4804    }
4805    else
4806    {
4807      list L = hermiteNormalForm(P, "transform"); //hermite(P, "transform");  // now, Hermite = L[1] = A*L[2]
4808      intmat U = L[2];
4809
4810      // delete columns 1 to r
4811      // should be identical with the function
4812      // U = submat(U,1..nrows(U),r+1..);
4813      // for matrices
4814      intmat Udel[nrows(U)][ncols(U) - r];
4815      int k;
4816      int m;
4817
4818      for ( k = 1; k <= nrows(U); k++ )
4819      {
4820        for ( m = r + 1; m <= ncols(U); m++ )
4821        {
4822          Udel[k,m - r] = U[k,m];
4823        }
4824      }
4825
4826      U = Udel;
4827
4828    }
4829  }
4830
4831  return(U);
4832}
4833example
4834{
4835  "EXAMPLE"; echo = 2;
4836
4837  intmat LL[3][4] =
4838    1,4,7,10,
4839    2,5,8,11,
4840    3,6,9,12;
4841
4842  // should be a 4x2 matrix whose columns
4843  // generate the same lattice as [-1,2,-1,0],[2,-3,0,1]
4844  intmat B = kernelLattice(LL);
4845 
4846  print(B);
4847
4848  // another example
4849  intmat C[2][4] =
4850    1,0,2,0,
4851    0,1,2,0;
4852
4853  // should result in a matrix whose
4854  // colums create the same  lattice as
4855  // [-2,-2,1,0], [0,0,0,1]
4856  intmat D = kernelLattice(C);
4857
4858  print(D);
4859
4860  kill B;
4861
4862}
4863
4864/*****************************************************/
4865
4866proc preimageLattice(def P, def B)
4867"
4868USAGE: preimageLattice(P, B); intmat P, intmat B
4869PURPOSE: compute an integral basis for the preimage of B under
4870 the homomorphism of lattices defined by the intmat P.
4871RETURNS: intmat
4872EXAMPLE: example preimageLattice; shows an example
4873"
4874{
4875  // concatenate matrices: Con = [P,-B]
4876  intmat Con[nrows(P)][ncols(P) + ncols(B)];
4877  int i;
4878  int j;
4879
4880  for ( i = 1; i <= nrows(Con); i++ )
4881  {
4882    for ( j = 1; j <= ncols(P); j++ ) // P first
4883    {
4884      Con[i,j] = P[i,j];
4885    }
4886  }
4887
4888  for ( i = 1; i <= nrows(Con); i++ )
4889  {
4890    for ( j = 1; j <= ncols(B); j++ ) // now -B
4891    {
4892      Con[i,ncols(P) + j] = - B[i,j];
4893    }
4894  }
4895
4896  intmat L = kernelLattice(Con);
4897
4898  // delete rows ncols(P)+1 to nrows(L) out of L
4899  intmat Del[ncols(P)][ncols(L)];
4900  int k;
4901  int m;
4902
4903  for ( k = 1; k <= nrows(Del); k++ )
4904  {
4905    for ( m = 1; m <= ncols(Del); m++ )
4906    {
4907      Del[k,m] = L[k,m];
4908    }
4909  }
4910
4911  L = latticeBasis(Del);
4912
4913  return(L);
4914
4915}
4916example
4917{
4918  "EXAMPLE"; echo = 2;
4919
4920  intmat P[2][3] =
4921    2,6,10,
4922    4,8,12;
4923
4924  intmat B[2][1] =
4925    1,
4926    0;
4927
4928  // should be a (3x2)-matrix with columns e.g. [1,1,-1] and [0,3,-2] (the generated lattice should be identical)
4929  print(preimageLattice(P,B));
4930
4931  // another example
4932  intmat W[3][3] =
4933    1,0,0,
4934    0,1,1,
4935    0,2,0;
4936
4937  intmat Z[3][2] =
4938    1,0,
4939    0,1,
4940    0,0;
4941
4942  // should be a (3x2)-matrix with columns e.g. [1,0,0] and [0,0,-1] (the generated lattice should be identical)
4943  print(preimageLattice(W,Z));
4944
4945}
4946
4947/******************************************************/
4948proc isPrimitiveSublattice(intmat A);
4949"USAGE: isPrimitiveSublattice(A); intmat A
4950PURPOSE: check whether the given set of integral vectors in ZZ^m,
4951i.e. the columns of A, generate a primitive sublattice in ZZ^m
4952(a direct summand of ZZ^m).
4953RETURNS: int, where 0 is false and 1 is true.
4954EXAMPLE: example isPrimitiveSublattice; shows an example
4955"
4956{
4957  intmat B = smithNormalForm(A);
4958  int r = intRank(B);
4959
4960  if ( r == 0 )
4961  {
4962    return(1);
4963  }
4964
4965  if ( 1 < B[r,r] )
4966  {
4967    return(0);
4968  }
4969
4970  return(1);
4971}
4972example
4973{
4974  "EXAMPLE"; echo = 2;
4975
4976  intmat A[3][2] =
4977    1,4,
4978    2,5,
4979    3,6;
4980
4981  // should be 0
4982  int b = isPrimitiveSublattice(A);
4983  print(b);
4984
4985  // another example
4986
4987  intmat B[2][2] =
4988    1,0,
4989    0,1;
4990
4991  // should be 1
4992  int c = isPrimitiveSublattice(B);
4993  print(c);
4994
4995  kill A, b, B, c;
4996}
4997
4998/******************************************************/
4999proc isIntegralSurjective(intmat P);
5000"USAGE: isIntegralSurjective(P); intmat P
5001PURPOSE: test whether the given linear map P of lattices is
5002surjective.
5003RETURNS: int, where 0 is false and 1 is true.
5004EXAMPLE: example isIntegralSurjective; shows an example
5005"
5006{
5007  int r = intRank(P);
5008
5009  if ( r < nrows(P) )
5010  {
5011    return(0);
5012  }
5013
5014
5015  if ( isPrimitiveSublattice(P) == 1 )
5016  {
5017    return(1);
5018  }
5019
5020  return(0);
5021}
5022example
5023{
5024  "EXAMPLE"; echo = 2;
5025
5026  intmat A[2][3] =
5027    1,3,5,
5028    2,4,6;
5029
5030  // should be 0
5031  int b = isIntegralSurjective(A);
5032  print(b);
5033
5034  // another example
5035  intmat B[2][3] =
5036    1,1,5,
5037    2,3,6;
5038
5039  // should be 1
5040  int c = isIntegralSurjective(B);
5041  print(c);
5042
5043  kill A, b, B, c;
5044}
5045
5046/******************************************************/
5047proc projectLattice(intmat B)
5048"USAGE: projectLattice(B); intmat B
5049PURPOSE: A set of vectors in ZZ^m is given as the columns of B.
5050Then this function provides a linear map ZZ^m --> ZZ^n
5051having the primitive span of B its kernel.
5052RETURNS: intmat
5053EXAMPLE: example projectLattice; shows an example
5054"
5055{
5056  int n = nrows(B);
5057  int r = intRank(B);
5058
5059  if ( r == 0 )
5060  {
5061    intmat U = unitMatrix(n);
5062  }
5063  else
5064  {
5065    if ( r == n )
5066    {
5067      intmat U[1][n]; // U now is the n-dim zero-vector
5068    }
5069    else
5070    {
5071      // we want a matrix with column operations so we transpose
5072      intmat BB = transpose(B);
5073      list L = hermiteNormalForm(BB, "transform"); 
5074      intmat U = transpose(L[2]);
5075
5076 
5077      // delete rows 1 to r
5078      intmat Udel[nrows(U) - r][ncols(U)];
5079      int k;
5080      int m;
5081
5082      for ( k = 1; k <= nrows(U) - r ; k++ )
5083      {
5084        for ( m = 1; m <= ncols(U); m++ )
5085        {
5086          Udel[k,m] = U[k + r,m];
5087        }
5088      }
5089
5090      U = Udel;
5091
5092    }
5093  }
5094
5095  return(U);
5096}
5097example
5098{
5099  "EXAMPLE"; echo = 2;
5100
5101  intmat B[4][2] =
5102    1,5,
5103    2,6,
5104    3,7,
5105    4,8;
5106
5107  // should result in a (2x4)-matrix such that the corresponding lattice is created by
5108  // [-1, 2], [-2, 3], [-1, 0] and [0, 1]
5109  print(projectLattice(B));
5110
5111  // another example
5112
5113  intmat BB[4][2] =
5114    1,0,
5115    0,1,
5116    0,0,
5117    0,0;
5118
5119  // should result in a (2x4)-matrix such that the corresponding lattice is created by
5120  // [0,0],[0,0],[1,0],[0,1]
5121  print(projectLattice(BB));
5122
5123  // another example
5124
5125  intmat BBB[3][4] =
5126    1,0,1,2,
5127    1,1,0,0,
5128    3,0,0,3;
5129
5130  // should result in the (1x3)-matrix that consists of just zeros
5131  print(projectLattice(BBB));
5132
5133}
5134
5135/******************************************************/
5136proc intersectLattices(intmat A, intmat B)
5137"USAGE: intersectLattices(A, B); intmat A, intmat B
5138PURPOSE: compute an integral basis for the intersection of the
5139lattices A and B.
5140RETURNS: intmat
5141EXAMPLE: example intersectLattices; shows an example
5142"
5143{
5144  // concatenate matrices: Con = [A,-B]
5145  intmat Con[nrows(A)][ncols(A) + ncols(B)];
5146  int i;
5147  int j;
5148
5149  for ( i = 1; i <= nrows(Con); i++ )
5150  {
5151    for ( j = 1; j <= ncols(A); j++ ) // A first
5152    {
5153      Con[i,j] = A[i,j];
5154    }
5155  }
5156
5157  for ( i = 1; i <= nrows(Con); i++ )
5158  {
5159    for ( j = 1; j <= ncols(B); j++ ) // now -B
5160    {
5161      Con[i,ncols(A) + j] = - B[i,j];
5162    }
5163  }
5164
5165  intmat K = kernelLattice(Con);
5166
5167  // delete all rows in K from ncols(A)+1 onwards
5168  intmat Bas[ncols(A)][ncols(K)];
5169
5170  for ( i = 1; i <= nrows(Bas); i++ )
5171  {
5172    for ( j = 1; j <= ncols(Bas); j++ )
5173    {
5174      Bas[i,j] = K[i,j];
5175    }
5176  }
5177
5178  // take product in order to obtain the intersection
5179  intmat S = A * Bas;
5180  intmat Cut = hermiteNormalForm(S); //hermite(S);
5181  int r = intRank(Cut);
5182
5183  if ( r == 0 )
5184  {
5185    intmat Cutdel[nrows(Cut)][1]; // is now the zero-vector
5186
5187    Cut = Cutdel;
5188  }
5189  else
5190  {
5191    // delte columns from r+1 onwards
5192    intmat Cutdel[nrows(Cut)][r];
5193
5194    for ( i = 1; i <= nrows(Cutdel); i++ )
5195    {
5196      for ( j = 1; j <= r; j++ )
5197      {
5198        Cutdel[i,j] = Cut[i,j];
5199      }
5200    }
5201
5202    Cut = Cutdel;
5203  }
5204
5205  return(Cut);
5206}
5207example
5208{
5209  "EXAMPLE"; echo = 2;
5210
5211  intmat A[3][2] =
5212    1,4,
5213    2,5,
5214    3,6;
5215
5216  intmat B[3][2] =
5217    6,9,
5218    7,10,
5219    8,11;
5220
5221  // should result in a (3x2)-matrix with columns
5222  //  e.g. [0, 3, 6], [-3, 0, 3] (the lattice should be the same)
5223  print(intersectLattices(A,B));
5224
5225  // another example
5226  intmat C[2][3] =
5227    1,0,0,
5228    3,2,5;
5229
5230  intmat D[2][3] =
5231    4,5,0,
5232    0,5,0;
5233
5234  // should result in a (3x2)-matrix whose columns generate the
5235  // same lattice as [1,5], [0, 20]
5236  print(intersectLattices(C,D));
5237}
5238
5239////////////////////////////////////
5240
5241proc intInverse(intmat A);
5242"USAGE: intInverse(A); intmat A
5243PURPOSE: compute the integral inverse of the intmat A.
5244If det(A) is neither 1 nor -1 an error is returned.
5245RETURNS: intmat
5246EXAMPLE: example intInverse; shows an example
5247"
5248{
5249  int d = det(A);
5250
5251  if ( d * d != 1 ) // is d = 1 or -1? Else: error
5252  {
5253    ERROR("determinant of the given intmat has to be 1 or -1.");
5254  }
5255
5256  int c;
5257  int i,j;
5258  intmat C[nrows(A)][ncols(A)];
5259  intmat Ad;
5260  int s;
5261
5262  for ( i = 1; i <= nrows(C); i++ )
5263  {
5264    for ( j = 1; j <= ncols(C); j++ )
5265    {
5266      Ad = intAdjoint(A,i,j);
5267      s = 1;
5268
5269      if ( ((i + j) % 2) > 0 )
5270      {
5271        s = -1;
5272      }
5273
5274      C[i,j] = d * s * det(Ad); // mult by d is equal to div by det
5275    }
5276  }
5277
5278  C = transpose(C);
5279
5280  return(C);
5281}
5282example
5283{
5284  "EXAMPLE"; echo = 2;
5285
5286  intmat A[3][3] =
5287    1,1,3,
5288    3,2,0,
5289    0,0,1;
5290
5291  intmat B = intInverse(A);
5292
5293  // should be the unit matrix
5294  print(A * B);
5295
5296  // another example
5297  intmat C[2][2] =
5298    2,1,
5299    3,2;
5300 
5301  intmat D  = intInverse(C);
5302
5303  // should be the unit matrix
5304  print(C * D);
5305
5306  kill A, B, C, D;
5307}
5308
5309
5310/******************************************************/
5311static proc intAdjoint(intmat A, int indrow, int indcol)
5312"USAGE: intAdjoint(A); intmat A
5313PURPOSE: return the matrix where the given row and column are deleted.
5314RETURNS: intmat
5315EXAMPLE: example intAdjoint; shows an example
5316"
5317{
5318  int n = nrows(A);
5319  int m = ncols(A);
5320  int i, j;
5321  intmat B[n - 1][m - 1];
5322  int a, b;
5323
5324  for ( i = 1; i < indrow; i++ )
5325  {
5326    for ( j = 1; j < indcol; j++ )
5327    {
5328      B[i,j] = A[i,j];
5329    }
5330    for ( j = indcol + 1; j <= ncols(A); j++ )
5331    {
5332      B[i,j - 1] = A[i,j];
5333    }
5334  }
5335
5336  for ( i = indrow + 1; i <= nrows(A); i++ )
5337  {
5338    for ( j = 1; j < indcol; j++ )
5339    {
5340      B[i - 1,j] = A[i,j];
5341    }
5342    for ( j = indcol+1; j <= ncols(A); j++ )
5343    {
5344      B[i - 1,j - 1] = A[i,j];
5345    }
5346  }
5347
5348  return(B);
5349}
5350example
5351{
5352  "EXAMPLE"; echo = 2;
5353
5354  intmat A[2][3] =
5355    1,3,5,
5356    2,4,6;
5357
5358  intmat B = intAdjoint(A,2,2);
5359  print(B);
5360
5361  kill A,B;
5362}
5363
5364/******************************************************/
5365proc integralSection(intmat P);
5366"USAGE: integralSection(P); intmat P
5367PURPOSE: for a given linear surjective map P of lattices
5368 this procedure returns an integral section of P.
5369RETURNS: intmat
5370EXAMPLE: example integralSection; shows an example
5371"
5372{
5373  int m = nrows(P);
5374  int n = ncols(P);
5375
5376  if ( m == n )
5377  {
5378    intmat U = intInverse(P);
5379  }
5380  else
5381  {
5382    intmat U = (hermiteNormalForm(P, "transform"))[2];
5383
5384    // delete columns m+1 to n
5385    intmat Udel[nrows(U)][ncols(U) - (n - m)];
5386    int k;
5387    int z;
5388
5389    for ( k = 1; k <= nrows(U); k++ )
5390    {
5391      for ( z = 1; z <= m; z++ )
5392      {
5393        Udel[k,z] = U[k,z];
5394      }
5395    }
5396
5397    U = Udel;
5398  }
5399
5400  return(U);
5401}
5402example
5403{
5404  "EXAMPLE"; echo = 2;
5405
5406  intmat P[2][4] =
5407    1,3,4,6,
5408    2,4,5,7;
5409
5410  // should be a matrix with two columns
5411  // for example: [-2, 1, 0, 0], [3, -3, 0, 1]
5412  intmat U = integralSection(P);
5413
5414  print(U);
5415  print(P * U);
5416
5417  kill U;
5418}
5419
5420
5421
5422/******************************************************/
5423proc factorgroup(G,H)
5424"USAGE: factorgroup(G,H); list G, list H
5425PURPOSE: returns a representation of the factor group G mod H using the first isomorphism thm
5426RETURNS: list
5427EXAMPLE: example factorgroup(G,H); shows an example
5428"
5429{
5430  intmat S1 = G[1];
5431  intmat L1 = G[2];
5432  intmat S2 = H[1];
5433  intmat L2 = H[2];
5434
5435  // check whether G,H are subgroups of a common group, i.e. whether L1 and L2 span the same lattice
5436  if ( !isSublattice(L1,L2) || !isSublattice(L2,L1))
5437  {
5438    ERROR("G and H are not subgroups of a common group.");
5439  }
5440
5441  // check whether H is a subgroup of G, i.e. whether S2 is a sublattice of S1+L1
5442  intmat B = concatintmat(S1,L1); // check whether this gives the concatinated matrix
5443  if ( !isSublattice(S2,B) )
5444  {
5445    ERROR("H is not a subgroup of G");
5446  }
5447  // use first isomorphism thm to get the factor group
5448  intmat L = concatintmat(L1,S2); // check whether this gives the concatinated matrix
5449  list GmodH;
5450  GmodH[1]=S1;
5451  GmodH[2]=L;
5452  return(GmodH);
5453}
5454example
5455{
5456  "EXAMPLE"; echo = 2;
5457
5458  intmat S1[2][2] =
5459    1,0,
5460    0,1;
5461  intmat L1[2][1] =
5462    2,
5463    0;
5464
5465  intmat S2[2][1] =
5466    1,
5467    0;
5468  intmat L2[2][1] =
5469    2,
5470    0;
5471
5472  list G = createGroup(S1,L1);
5473  list H = createGroup(S2,L2);
5474
5475  list N = factorgroup(G,H);
5476  print(N);
5477
5478  kill G,H,N,S1,L1,S2,L2;
5479
5480}
5481
5482/******************************************************/
5483proc productgroup(G,H)
5484"USAGE: productgroup(G,H); list G, list H
5485PURPOSE: Returns a representation of the group G x H
5486RETURNS: list
5487EXAMPLE: example productgroup(G,H); shows an example
5488"
5489{
5490  intmat S1 = G[1];
5491  intmat L1 = G[2];
5492  intmat S2 = H[1];
5493  intmat L2 = H[2];
5494  intmat OS1[nrows(S1)][ncols(S2)];
5495  intmat OS2[nrows(S2)][ncols(S1)];
5496  intmat OL1[nrows(L1)][ncols(L2)];
5497  intmat OL2[nrows(L2)][ncols(L1)];
5498
5499  // concatinate matrices to get S
5500  intmat A = concatintmat(S1,OS1);
5501  intmat B = concatintmat(OS2,S2);
5502  intmat At = transpose(A);
5503  intmat Bt = transpose(B);
5504  intmat St = concatintmat(At,Bt);
5505  intmat S = transpose(St);
5506
5507  // concatinate matrices to get L
5508  intmat C = concatintmat(L1,OL1);
5509  intmat D = concatintmat(OL2,L2);
5510  intmat Ct = transpose(C);
5511  intmat Dt = transpose(D);
5512  intmat Lt = concatintmat(Ct,Dt);
5513  intmat L = transpose(Lt);
5514
5515  list GxH;
5516  GxH[1]=S;
5517  GxH[2]=L;
5518  return(GxH);
5519}
5520example
5521{
5522  "EXAMPLE"; echo = 2;
5523
5524  intmat S1[2][2] =
5525    1,0,
5526    0,1;
5527  intmat L1[2][1] =
5528    2,
5529    0;
5530
5531  intmat S2[2][2] =
5532    1,0,
5533    0,2;
5534  intmat L2[2][1] =
5535    0,
5536    3;
5537
5538  list G = createGroup(S1,L1);
5539  list H = createGroup(S2,L2);
5540
5541  list N = productgroup(G,H);
5542  print(N);
5543
5544  kill G,H,N,S1,L1,S2,L2;
5545
5546}
5547
5548/******************************************************/
5549proc primitiveSpan(intmat V);
5550"USAGE: primitiveSpan(V); intmat V
5551PURPOSE: compute an integral basis for the minimal primitive
5552sublattice that contains the given vectors, i.e. the columns of V.
5553RETURNS: int, where 0 is false and 1 is true.
5554EXAMPLE: example primitiveSpan; shows an example
5555"
5556{
5557  int n = ncols(V);
5558  int m = nrows(V);
5559  int r = intRank(V);
5560
5561
5562  if ( r == 0 )
5563  {
5564    intmat P[m][1]; //  this is the m-zero-vector now
5565  }
5566  else
5567  {
5568    list L = smithNormalForm(V, "transform"); // L = [A,S,B] where S is the smith-NF and S = A*S*B
5569    intmat P = intInverse(L[1]);
5570
5571//    print(L);
5572
5573    if ( r < m )
5574    {
5575      // delete columns r+1 to m in P:
5576      intmat Pdel[nrows(P)][r];
5577      int i,j;
5578
5579      for ( i = 1; i <= nrows(Pdel); i++ )
5580      {
5581        for ( j = 1; j <= ncols(Pdel); j++ )
5582        {
5583          Pdel[i,j] = P[i,j];
5584        }
5585      }
5586
5587      P = Pdel;
5588    }
5589  }
5590
5591  return(P);
5592}
5593example
5594{
5595  "EXAMPLE"; echo = 2;
5596
5597  intmat V[3][2] =
5598    1,4,
5599    2,5,
5600    3,6;
5601
5602  // should return a (3x2)-matrix whose columns
5603  // generate the same lattice as [1, 2, 3] and [0, 1, 2]
5604  intmat R = primitiveSpan(V);
5605  print(R);
5606
5607  // another example
5608  intmat W[2][2] =
5609    1,0,
5610    0,1;
5611
5612  // should return a (2x2)-matrix whose columns
5613  // generate the same lattice as [1, 0] and [0, 1]
5614  intmat S = primitiveSpan(W);
5615  print(S);
5616
5617  kill V, R, S, W;
5618}
5619
5620/***********************************************************/
Note: See TracBrowser for help on using the repository browser.