source: git/Singular/LIB/multigrading.lib @ b840b1

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