Changeset b6ae8c in git


Ignore:
Timestamp:
Mar 15, 2011, 7:17:52 PM (13 years ago)
Author:
Oleksandr Motsak <motsak@…>
Branches:
(u'fieker-DuVal', '117eb8c30fc9e991c4decca4832b1d19036c4c65')(u'spielwiese', '3720ae8bfcff4a4649ee98a15552089151d2d59b')
Children:
199cd68d23f96331d3327f7a0b11a46accefa1db
Parents:
4f7d761d43192f9c176ef0f44b0f61e891545eb8
Message:
ADD: newest version of multigrading.lib

From: Oleksandr Motsak <motsak@mathematik.uni-kl.de>

git-svn-id: file:///usr/local/Singular/svn/trunk@13984 2c84dea3-7e68-4137-9b89-c4e89433aadc
File:
1 edited

Legend:

Unmodified
Added
Removed
  • Singular/LIB/multigrading.lib

    r4f7d76 rb6ae8c  
     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
    13///////////////////////////////////////////////////////////////////
    24version="$Id$";
     
    57LIBRARY:  multigrading.lib          Multigraded Rings
    68
    7 AUTHORS:  Rene Birkner, rbirkner@math.fu-berlin.de
     9AUTHORS:  Benjamin Bechtold, benjamin.bechtold@googlemail.com
     10@*        Rene Birkner, rbirkner@math.fu-berlin.de
    811@*        Lars Kastner, lkastner@math.fu-berlin.de
     12@*        Simon Keicher, keicher@mail.mathematik.uni-tuebingen.de
    913@*        Oleksandr Motsak, U@D, where U={motsak}, D={mathematik.uni-kl.de}
    1014
     
    1216OVERVIEW: This library allows one to virtually add multigradings to Singular.
    1317For more see http://code.google.com/p/convex-singular/wiki/Multigrading
    14 For theoretical references see:
     18For theoretical references see: 
    1519E. Miller, B. Sturmfels: 'Combinatorial Commutative Algebra' and
    1620M. Kreuzer, L. Robbiano: 'Computational Commutative Algebra'.
    1721
    1822NOTE: 'mDegBasis' relies on 4ti2 for computing Hilbert Bases.
     23All groups are finitely generated Abelian
    1924
    2025PROCEDURES:
    21 setBaseMultigrading(M,T); attach multiweights/torsion matrices to the basering
     26setBaseMultigrading(M,L); attach multiweights/grading group matrices to the basering
    2227getVariableWeights([R]);  get matrix of multidegrees of vars attached to a ring
    23 getTorsion([R]);          get torsion matrix attached to a ring
     28
     29getGradingGroup([R]);     get grading group attached to a ring
     30getLattice([R[,choice]]); get grading group' lattice attached to a ring (or its NF)
     31
     32createGroup(S,L);          create a group generated by S, with relations L
     33createQuotientGroup(L);    create a group generated by the unit matrix whith relations L
     34createTorsionFreeGroup(S); create a group generated by S which is torsionfree
     35printGroup(G);             print a group
     36
     37areIsomorphicGroups(G,H);     test wheter G an H are isomorphic groups (TODO Tuebingen)
     38isGroup(G);                   test whether G is a valid group
     39isGroupHomomorphism(L1,L2,A); test wheter A defines a group homomrphism from L1 to L2
     40
     41isGradedRingHomomorphism(R,f,A);  test graded ring homomorph
     42createGradedRingHomomorphism(R,f,A);  create a graded ring homomorph
    2443
    2544setModuleGrading(M,v);    attach multiweights of units to a module and return it
    2645getModuleGrading(M);      get multiweights of module units (attached to M)
    2746
     47isSublattice(A,B);        test whether A is a sublattice of B
     48imageLattice(P,L);        computes an integral basis for the image of the lattice L under the linear map P.
     49intRank(A);               computes the rank of the intmat A
     50kernelLattice(P);         computes an integral basis for the kernel of the linear map P.
     51latticeBasis(B);          compute an integral basis of the lattice B
     52preimageLattice(P,L);     computes an integral basis for the preimage of the lattice L under the linear map P.
     53projectLattices(B);       compute a linear map of lattices having the primitive span of B as its kernel.
     54intersectLattices(A,B);   compute an integral basis for the intersection of the lattices A and B.
     55isIntegralSurjective(P);  test whether the map P of lattices is surjective.
     56isPrimitiveSublattice(A); test whether A generates a primitive sublattice.
     57intInverse(A);            compute the integral inverse matrix of the intmat A
     58intAdjoint(A,i,j);        delete row i and column j of the intmat A.
     59integralSection(P);       for a given linear surjective map P of lattices this procedure returns an integral section of P.
     60primitiveSpan(A);         compute a basis for the minimal primitive sublattice that contains the given vectors (by A).
     61
     62factorgroup(G,H);         create the group G mod H
     63productgroup(G,H);        create the group G x H
     64
    2865mDeg(A);                  compute the multidegree of A
    2966mDegBasis(d);             compute all monomials of multidegree d
    30 mDegPartition(p);         compute the multigraded-homogenous components of p
    31 
    32 isTorsionFree();          test whether the current multigrading is torsion free
    33 isTorsionElement(p);      test whether p has zero multidegree
    34 isHomogenous(a);          test whether 'a' is multigraded-homogenous
     67mDegPartition(p);         compute the multigraded-homogeneous components of p
     68
     69isTorsionFree();          test whether the current multigrading is  free
     70isPositive();             test whether the current multigrading is positive
     71isZeroElement(p);         test whether p has zero multidegree
     72areZeroElements(M);       test whether an integer matrix M considered as a collection of columns has zero multidegree   
     73isHomogeneous(a);         test whether 'a' is multigraded-homogeneous
    3574
    3675equalMDeg(e1,e2[,V]);     test whether e1==e2 in the current multigrading
     
    3877mDegGroebner(M);          compute the multigraded GB/SB of M
    3978mDegSyzygy(M);            compute the multigraded syzygies of M
     79mDegModulo(I,J);          compute the multigraded 'modulo' module of I and J
    4080mDegResolution(M,l[,m]);  compute the multigraded resolution of M
    41 
    42 defineHomogenous(p);      get a torsion matrix wrt which p becomes homogenous
     81mDegTensor(m,n);          compute the tensor product of multigraded modules m,n
     82mDegTor(i,m,n);           compute the Tor_i(m,n) for multigraded modules m,n
     83
     84defineHomogeneous(p);     get a grading group wrt which p becomes homogeneous
    4385pushForward(f);           find the finest grading on the image ring, homogenizing f
    44 
    45 hermite(A);               compute the Hermite Normal Form of a matrix
     86gradiator(h);             coarsens grading of the ring until h becomes homogeneous
     87
     88hermiteNormalForm(A);     compute the Hermite Normal Form of a matrix
     89smithNormalForm(A,#);     compute matrices D,P,Q with D=P*A*Q and D is the smith normal form of A
    4690
    4791hilbertSeries(M);         compute the multigraded Hilbert Series of M
     92evalHilbertSeries(h,v);   evaluate hilberts series h by substituting v[i] for t_(i)
     93
     94lll(A);                   applies LLL(.) of lll.lib which only works for lists on a matrix A
    4895
    4996           (parameters in square brackets [] are optional)
    5097
    51 KEYWORDS:  multigradeding, multidegree, multiweights, multigraded-homogenous
     98KEYWORDS:  multigrading, multidegree, multiweights, multigraded-homogeneous, integral linear algebra
    5299";
    53100
     
    56103
    57104LIB "standard.lib"; // for groebner
     105LIB "lll.lib"; // for lll_matrix
     106LIB "matrix.lib"; // for mDegTor
     107
     108/******************************************************/
     109
     110static proc concatintmat(intmat A, intmat B)
     111{
     112
     113 if ( nrows(A) != nrows(B) )
     114   {
     115     ERROR("matrices A and B have different number of rows.");
     116   }
     117
     118 intmat At = transpose(A);
     119 intmat Bt = transpose(B);
     120
     121 intmat Ct[nrows(At) + nrows(Bt)][ncols(At)] = At, Bt;
     122
     123 return(transpose(Ct));
     124}
     125
     126
     127/******************************************************/
     128proc createGradedRingHomomorphism(def src, ideal Im, def A)
     129"USAGE: createGradedRingHomomorphism(R, f, A); ring R, ideal f, group homomorphism A
     130PURPOSE: create a multigraded group ring homomorphism defined by
     131a ring map from R to the current ring, given by generators images f
     132and a group homomorphism A between grading groups
     133RETURN: graded ring homorphism
     134EXAMPLE: example createGradedRingHomomorphism; shows an example
     135"
     136{
     137  string isGRH = "isGRH";
     138
     139  if( !isGradedRingHomomorphism(def src, ideal Im, def A) )
     140  {
     141    ERROR("Input data is wrong");
     142  }
     143 
     144  list h;
     145  h[3] = A;
     146 
     147//  map f = src, Im;
     148  h[2] = Im; // f?
     149  h[1] = src;
     150 
     151  attrib(h, isGRH, (1==1)); // mark it "a graded ring homomorphism"
     152
     153  return(h);
     154}
     155example
     156{
     157  "EXAMPLE:"; echo=2;
     158
     159  // TODO!
     160
     161}
     162
     163
     164/******************************************************/
     165proc isGradedRingHomomorphism(def src, ideal Im, def A)
     166"USAGE: isGradedRingHomomorphism(R, f, A); ring R, ideal f, group homomorphism A
     167PURPOSE: test a multigraded group ring homomorphism defined by
     168a ring map from R to the current ring, given by generators images f
     169and a group homomorphism A between grading groups
     170RETURN: int, 1 for TRUE, 0 otherwise
     171EXAMPLE: example isGradedRingHomomorphism; shows an example
     172"
     173{
     174  def dst = basering;
     175
     176  intmat result_degs = mDeg(Im);
     177  print(result_degs);
     178   
     179  setring src;
     180   
     181  intmat input_degs = mDeg(maxideal(1));
     182  print(input_degs);
     183
     184  def image_degs = A * input_degs;
     185  print( image_degs );
     186 
     187  def df = image_degs - result_degs;
     188  print(df);
     189 
     190  setring dst;
     191   
     192  return (areZeroElements( df ));
     193}
     194example
     195{
     196  "EXAMPLE:"; echo=2;
     197 
     198  ring r = 0, (x, y, z), dp;
     199  intmat S1[3][3] =
     200    1, 0, 0,
     201    0, 1, 0,
     202    0, 0, 1;
     203  intmat L1[3][1] =
     204    0,
     205    0,
     206    0;
     207   
     208  def G1 = createGroup(S1, L1); // (S1 + L1)/L1
     209  printGroup(G1);
     210   
     211  setBaseMultigrading(S1, L1); // to change...
     212
     213  ring R = 0, (a, b, c), dp;
     214  intmat S2[2][3] =
     215    1, 0,
     216    0, 1;
     217  intmat L2[2][1] =
     218    0,
     219    2;
     220
     221  def G2 = createGroup(S2, L2);
     222  printGroup(G2);
     223   
     224  setBaseMultigrading(S2, L2); // to change...
     225
     226
     227  map F = r, a, b, c;
     228  intmat A[nrows(L2)][nrows(L1)] =
     229      1, 0, 0,
     230      3, 2, -6;
     231   
     232  // graded ring homomorphism is given by (compatible):
     233  print(F);
     234  print(A);
     235
     236  isGradedRingHomomorphism(r, ideal(F), A);
     237  def h = createGradedRingHomomorphism(r, ideal(F), A);
     238
     239  print(h);
     240 
     241  // not a homo..
     242  intmat B[nrows(L2)][nrows(L1)] =
     243     1, 1, 1,
     244     0, 0, 0;
     245  print(B);
     246
     247  isGradedRingHomomorphism(r, ideal(F), B);
     248  def hh = createGradedRingHomomorphism(r, ideal(F), B);
     249 
     250  if( defined(hh) ) { ERROR("That input was not valid"); }
     251}
     252
     253
     254proc createQuotientGroup(intmat L)
     255"
     256L - relations
     257TODO: bad name
     258"
     259{
     260  int r = nrows(L); int i;
     261  intmat S[r][r]; // SQUARE!!!
     262  for(i = r; i > 0; i--){ S[i, i] = 1; }
     263  return (createGroup(S,L));
     264}
     265
     266proc createTorsionFreeGroup(intmat S)
     267"
     268S - generators
     269TODO: bad name
     270"
     271{
     272  int r = nrows(S); int i;
     273  intmat L[r][1] = 0;
     274  return (createGroup(S,L));
     275}
     276
     277
     278/******************************************************/
     279proc createGroup(intmat S, intmat L)
     280"USAGE: createGroup(S, L); S, L are integer matrices
     281PURPOSE: create the group of the form (S+L)/L, i.e.
     282S specifies generators, L specifies relations.
     283RETURN: group
     284EXAMPLE: example createGroup; shows an example
     285"
     286{
     287  string isGroup = "isGroup";
     288  string attrGroupHNF = "hermite";
     289  string attrGroupSNF = "smith";
     290
     291   
     292/*
     293  if( size(#) > 0 )
     294  {
     295    if( typeof(#[1]) == "intmat" )
     296    {
     297      intmat S = #[1];
     298    } else { ERROR("Wrong optional argument: 1"); }
     299
     300    if( size(#) > 1 )
     301    {
     302      if( typeof(#[2]) == "intmat" )
     303      {
     304        intmat L = #[2];
     305      } else { ERROR("Wrong optional argument: 2"); }
     306    }
     307  }
     308  if( !defined(S) )
     309  {}
     310*/   
     311
     312  if( nrows(L) != nrows(S) )
     313  {
     314    ERROR("Incompatible matrices!");
     315  }
     316 
     317  def H = attrib(L, attrGroupHNF);
     318  if( !defined(H) || typeof(H) != "intmat")
     319  {
     320    attrib(L, attrGroupHNF, hermiteNormalForm(L));
     321  } else { kill H; }
     322 
     323  def HH = attrib(L, attrGroupSNF);
     324  if( !defined(HH) || typeof(HH) != "intmat")
     325  {
     326    attrib(L, attrGroupSNF, smithNormalForm(L));
     327  } else { kill HH; }
     328
     329  list G; // Please, note the order: Generators + Relations:
     330  G[1] = S; G[2] = L;
     331
     332  attrib(G, isGroup, (1==1)); // mark it "a group"
     333
     334  return (G);
     335}
     336example
     337{
     338  "EXAMPLE:"; echo=2;
     339   
     340  intmat S[3][3] =
     341    1, 0, 0,
     342    0, 1, 0,
     343    0, 0, 1;
     344
     345  intmat L[3][2] =
     346    1, 1,
     347    1, 3,
     348    1, 5;
     349   
     350  def G = createGroup(S, L); // (S+L)/L
     351
     352  printGroup(G);
     353   
     354  kill S, L, G;
     355
     356  /////////////////////////////////////////////////
     357  intmat S[2][3] =
     358    1, -2, 1,
     359    1,  1, 0;
     360
     361  intmat L[2][1] =
     362    0,
     363    2;
     364
     365  def G = createGroup(S, L); // (S+L)/L
     366
     367  printGroup(G);
     368   
     369  kill S, L, G;
     370   
     371  // ----------- extreme case ------------ //
     372  intmat S[1][3] =
     373    1,  -1, 10;
     374
     375  // Torsion:
     376  intmat L[1][1] =
     377    0;
     378
     379  def G = createGroup(S, L); // (S+L)/L
     380
     381  printGroup(G);
     382}
     383
     384
     385/******************************************************/
     386proc printGroup(def G)
     387"USAGE: printGroup(G); G is a group
     388PURPOSE: prints the group G
     389RETURN: nothing
     390EXAMPLE: example printGroup; shows an example
     391"
     392{
     393  "Generators: ";
     394  print(G[1]);
     395
     396  "Relations: ";
     397  print(G[2]);
     398 
     399//  attrib(G[2]);
     400}
     401example
     402{
     403  "EXAMPLE:"; echo=2;
     404   
     405}
     406
     407/******************************************************/
     408proc areIsomorphicGroups(def G, def H)
     409"USAGE: areIsomorphicGroups(G, H); G and H are groups
     410PURPOSE: ?
     411RETURN: int, 1 for TRUE, 0 otherwise
     412EXAMPLE: example areIsomorphicGroups; shows an example
     413"
     414{
     415  return (1); // TRUE
     416}
     417example
     418{
     419  "EXAMPLE:"; echo=2;
     420   
     421}
     422
     423
     424proc isGroup(def G)
     425"test whether G is a valid group"
     426{
     427  string isGroup = "isGroup";
     428 
     429  // valid?
     430  if( typeof(G) != "list" ){ return(0); }
     431
     432  def a = attrib(G, isGroup);
     433
     434///// TODO for Hans: fix attr^2 bug in Singular!
     435
     436//  if( !defined(a) ) { return(0); }
     437//  if( typeof(a) != "int" ) { return(0); }
     438//  if( !a ){ return(0); }
     439
     440   
     441  if( size(G) != 2 ){ return(0); }
     442  if( typeof(G[1]) != "intmat" ){ return(0); }
     443  if( typeof(G[2]) != "intmat" ){ return(0); }
     444  if( nrows(G[1]) != nrows(G[2]) ){ return(0); }
     445   
     446  return(1==1);
     447}
     448
     449
    58450
    59451/******************************************************/
    60452proc setBaseMultigrading(intmat M, list #)
    61 "USAGE: setBaseMultigrading(M[, T]); M, T are integer matrices
    62 PURPOSE: attaches weights of variables and torsion to the basering.
     453"USAGE: setBaseMultigrading(M[, G]); M is an intege matrix, G is a group (or lattice)
     454PURPOSE: attaches weights of variables and grading group to the basering.
    63455NOTE: M encodes the weights of variables column-wise.
    64 The torsion is given by the lattice spanned by the columns of the integer
    65 matrix T in Z^nrows(M) over Z.
    66456RETURN: nothing
    67457EXAMPLE: example setBaseMultigrading; shows an example
     
    69459{
    70460  string attrMgrad   = "mgrad";
    71   string attrTorsion = "torsion";
    72   string attrTorsionHNF = "torsion_hermite";
    73 
    74 
    75   attrib(basering, attrMgrad, M);
    76 
    77   if( size(#) > 0 and typeof(#[1]) == "intmat" )
    78   {
    79     def T = #[1];
    80   } else
    81   {
    82     intmat T[nrows(M)][1];
    83   }
    84 
    85   if( nrows(T) == nrows(M) )
    86   {
    87     attrib(basering, attrTorsion, T);
    88 //    def H;
    89 //    attrib(basering, attrTorsionHNF, H);
    90   }
     461  string attrGradingGroup = "gradingGroup";
     462 
     463  if( size(#) > 0 )
     464  {
     465    if( typeof(#[1]) == "intmat" )
     466    {
     467      def L = createGroup(M, #[1]);
     468    }
     469
     470    if( isGroup(#[1]) )
     471    {
     472      def L = #[1];
     473
     474      if( !isSublattice(M, L[1]) )
     475      {
     476        ERROR("Multigrading is not contained in the grading group!");
     477      }
     478    }
     479  }
    91480  else
    92481  {
    93     ERROR("Incompatible matrix sizes!");
    94   }
     482    def L = createTorsionFreeGroup(M);
     483  }
     484
     485  if( !defined(L) ){ ERROR("Wrong arguments: no group given?"); }
     486   
     487
     488 
     489  attrib(basering, attrMgrad, M);
     490  attrib(basering, attrGradingGroup, L);
     491
     492  ideal Q = ideal(basering);
     493  if( !isHomogeneous(Q) ) // easy now, but would be hard before setting ring attributes!
     494  {
     495    "Warning: your quotient ideal is not homogenous (multigrading was set anyway)!";
     496  }
     497 
     498   
     499
     500
    95501}
    96502example
    97503{
    98504  "EXAMPLE:"; echo=2;
    99 
     505   
    100506  ring R = 0, (x, y, z), dp;
    101507
     
    106512    0, 0, 1;
    107513
    108   // Torsion:
     514  // GradingGroup:
    109515  intmat L[3][2] =
    110516    1, 1,
    111     1, 3,
     517    1, 3, 
    112518    1, 5;
    113 
     519   
    114520  // attaches M & L to R (==basering):
    115521  setBaseMultigrading(M, L); // Grading: Z^3/L
     
    117523  // Weights are accessible via "getVariableWeights()":
    118524  getVariableWeights();
    119 
    120   // Test all possible usages:
     525 
     526  // Test all possible usages: 
    121527  (getVariableWeights() == M) && (getVariableWeights(R) == M) && (getVariableWeights(basering) == M);
    122528
    123   // Torsion is accessible via "getTorsion()":
    124   getTorsion();
    125 
    126   // Test all possible usages:
    127   (getTorsion() == L) && (getTorsion(R) == L) && (getTorsion(basering) == L);
    128 
    129   // And its hermite NF via getTorsion("hermite"):
    130   getTorsion("hermite");
    131 
    132   // Test all possible usages:
    133   intmat H = hermite(L);
    134   (getTorsion("hermite") == H) && (getTorsion(R, "hermite") == H) && (getTorsion(basering, "hermite") == H);
     529  // Grading group is accessible via "getLattice()":
     530  getLattice();
     531 
     532  // Test all possible usages: 
     533  (getLattice() == L) && (getLattice(R) == L) && (getLattice(basering) == L);
     534
     535  // And its hermite NF via getLattice("hermite"):
     536  getLattice("hermite");
     537
     538  // Test all possible usages: 
     539  intmat H = hermiteNormalForm(L);
     540  (getLattice("hermite") == H) && (getLattice(R, "hermite") == H) && (getLattice(basering, "hermite") == H);
    135541
    136542  kill L, M;
     
    147553    0,
    148554    2;
    149 
     555   
    150556  // attaches M & L to R (==basering):
    151557  setBaseMultigrading(M, L); // Grading: Z + (Z/2Z)
     
    154560  getVariableWeights() == M;
    155561
    156   // Torsion is accessible via "getTorsion()":
    157   getTorsion() == L;
     562  // Torsion is accessible via "getLattice()":
     563  getLattice() == L;
    158564
    159565  kill L, M;
     
    167573  intmat L[1][1] =
    168574    0;
    169 
     575   
    170576  // attaches M & L to R (==basering):
    171577  setBaseMultigrading(M); // Grading: Z^3
     
    174580  getVariableWeights() == M;
    175581
    176   // Torsion is accessible via "getTorsion()":
    177   getTorsion() == L;
     582  // Torsion is accessible via "getLattice()":
     583  getLattice() == L;
    178584}
    179585
     
    208614  def M = attrib(R, attrMgrad);
    209615  if( typeof(M) == "intmat"){ return (M); }
    210   ERROR( "Sorry no multigrading matrix!" );
     616  ERROR( "Sorry no multigrading matrix!" ); 
    211617}
    212618example
    213619{
    214620  "EXAMPLE:"; echo=2;
    215 
     621   
    216622  ring R = 0, (x, y, z), dp;
    217623
     
    222628    0, 0, 1;
    223629
    224   // Torsion:
     630  // Grading group:
    225631  intmat L[3][2] =
    226632    1, 1,
    227     1, 3,
     633    1, 3, 
    228634    1, 5;
    229 
     635   
    230636  // attaches M & L to R (==basering):
    231637  setBaseMultigrading(M, L); // Grading: Z^3/L
     
    243649    1,  1, 0;
    244650
    245   // Torsion:
     651  // Grading group:
    246652  intmat L[2][1] =
    247653    0,
    248654    2;
    249 
     655   
    250656  // attaches M & L to R (==basering):
    251657  setBaseMultigrading(M, L); // Grading: Z + (Z/2Z)
     
    262668    1,  -1, 10;
    263669
    264   // Torsion:
     670  // Grading group:
    265671  intmat L[1][1] =
    266672    0;
    267 
     673   
    268674  // attaches M & L to R (==basering):
    269675  setBaseMultigrading(M); // Grading: Z^3
     
    273679}
    274680
    275 /******************************************************/
    276 proc getTorsion(list #)
    277 "USAGE: getTorsion([R[,opt]])
    278 PURPOSE: get associated torsion matrix, i.e. generators (cols) of the Torsion group
    279 RETURN: intmat, the torsion matrix, or its hermite normal form
    280         if an optional argument (\"hermite\") is given
    281 EXAMPLE: example getTorsion; shows an example
    282 "
    283 {
    284   string attrTorsion    = "torsion";
    285   string attrTorsionHNF = "torsion_hermite";
     681
     682proc getGradingGroup(list #)
     683"USAGE: getGradingGroup([R])
     684PURPOSE: get associated grading group
     685RETURN: group, the grading group
     686EXAMPLE: example getGradingGroup; shows an example
     687"
     688{
     689  string attrGradingGroup    = "gradingGroup";
    286690
    287691  int i = 1;
     
    293697      def R = #[i];
    294698      i++;
    295     }
    296   }
     699    } 
     700  } 
    297701
    298702  if( !defined(R) )
     
    301705  }
    302706
    303   if( size(#) >= i )
    304   {
    305     if( #[i] == "hermite" )
    306     {
    307       if( typeof(attrib(R, attrTorsionHNF)) != "intmat" )
    308       {
    309         def M = getTorsion(R);
    310         if( typeof(M) != "intmat")
    311         {
    312           ERROR( "Sorry no torsion matrix!" );
    313         }
    314         M = hermite(M);
    315         attrib(R, attrTorsionHNF, M); // this might not work with R != basering...
    316       }
    317       return (attrib(R, attrTorsionHNF));
    318     }
    319   }
    320 
    321   def M = attrib(R, attrTorsion);
    322   if( typeof(M) != "intmat")
    323   {
    324     ERROR( "Sorry no torsion matrix!" );
    325   }
    326   return (M);
     707  def G = attrib(R, attrGradingGroup);
     708 
     709  if( !isGroup(G) )
     710  {
     711    ERROR("Sorry no grading group!");
     712  }   
     713
     714  return(G);
    327715}
    328716example
    329717{
    330718  "EXAMPLE:"; echo=2;
    331 
     719   
    332720  ring R = 0, (x, y, z), dp;
    333721
     
    341729  intmat L[3][2] =
    342730    1, 1,
    343     1, 3,
     731    1, 3, 
    344732    1, 5;
    345 
     733   
    346734  // attaches M & L to R (==basering):
    347735  setBaseMultigrading(M, L); // Grading: Z^3/L
    348736
    349   // Torsion is accessible via "getTorsion()":
    350   getTorsion() == L;
    351 
    352   // its hermite NF:
    353   print(getTorsion("hermite"));
    354 
    355   kill L, M;
     737  def G = getGradingGroup();
     738
     739  printGroup( G );
     740 
     741  G[1] == M; G[2] == L;
     742
     743  kill L, M, G;
    356744
    357745  // ----------- isomorphic multigrading -------- //
     
    366754    0,
    367755    2;
    368 
     756   
    369757  // attaches M & L to R (==basering):
    370758  setBaseMultigrading(M, L); // Grading: Z + (Z/2Z)
    371759
    372   // Torsion is accessible via "getTorsion()":
    373   getTorsion() == L;
    374 
    375   // its hermite NF:
    376   print(getTorsion("hermite"));
    377 
    378   kill L, M;
    379 
     760  def G = getGradingGroup();
     761
     762  printGroup( G );
     763 
     764  G[1] == M; G[2] == L;
     765
     766  kill L, M, G;
    380767  // ----------- extreme case ------------ //
    381768
     
    387774  intmat L[1][1] =
    388775    0;
    389 
     776   
    390777  // attaches M & L to R (==basering):
    391778  setBaseMultigrading(M); // Grading: Z^3
    392779
    393   // Torsion is accessible via "getTorsion()":
    394   getTorsion() == L;
     780  def G = getGradingGroup();
     781
     782  printGroup( G );
     783 
     784  G[1] == M; G[2] == L;
     785
     786  kill L, M, G;
     787}
     788
     789
     790/******************************************************/
     791proc getLattice(list #)
     792"USAGE: getLattice([R[,opt]])
     793PURPOSE: get associated grading group matrix, i.e. generators (cols) of the grading group
     794RETURN: intmat, the grading group matrix, or
     795its hermite normal form if an optional argument (\"hermiteNormalForm\") is given or
     796smith normal form if an optional argument (\"smith\") is given
     797EXAMPLE: example getLattice; shows an example
     798"
     799{
     800  int i = 1;
     801  if( size(#) >= i )
     802  {
     803    if( ( typeof(#[i]) == "ring" ) or ( typeof(#[i]) == "qring" ) )
     804    {
     805      i++;
     806    }
     807  }
     808
     809  string attrGradingGroupHNF = "hermite";
     810  string attrGradingGroupSNF = "smith";
     811
     812  def G = getGradingGroup(#);
     813 
     814//  printGroup(G);
     815
     816
     817
     818  def T = G[2];
     819
     820  if( size(#) >= i )
     821  {
     822    if( #[i] == "hermite" )
     823    {
     824      def M = attrib(T, attrGradingGroupHNF);
     825      if( (!defined(M)) or (typeof(M) != "intmat") )
     826      {
     827        M = hermiteNormalForm(T);
     828      }
     829      return (M);
     830    }
     831
     832    if( #[i] == "smith" )
     833    {
     834      def M = attrib(T, attrGradingGroupSNF);
     835      if( (!defined(M)) or (typeof(M) != "intmat") )
     836      {
     837        M = smithNormalForm(T);
     838      }
     839      return (M);
     840    }
     841  }
     842
     843  return(T);
     844}
     845example
     846{
     847  "EXAMPLE:"; echo=2;
     848   
     849  ring R = 0, (x, y, z), dp;
     850
     851  // Weights of variables
     852  intmat M[3][3] =
     853    1, 0, 0,
     854    0, 1, 0,
     855    0, 0, 1;
     856
     857  // Torsion:
     858  intmat L[3][2] =
     859    1, 1,
     860    1, 3,
     861    1, 5;
     862   
     863  // attaches M & L to R (==basering):
     864  setBaseMultigrading(M, L); // Grading: Z^3/L
     865
     866  // Torsion is accessible via "getLattice()":
     867  getLattice() == L;
    395868
    396869  // its hermite NF:
    397   print(getTorsion("hermite"));
    398 }
     870  print(getLattice("hermite"));
     871
     872  kill L, M;
     873
     874  // ----------- isomorphic multigrading -------- //
     875
     876  // Weights of variables
     877  intmat M[2][3] =
     878    1, -2, 1,
     879    1,  1, 0;
     880
     881  // Torsion:
     882  intmat L[2][1] =
     883    0,
     884    2;
     885   
     886  // attaches M & L to R (==basering):
     887  setBaseMultigrading(M, L); // Grading: Z + (Z/2Z)
     888
     889  // Torsion is accessible via "getLattice()":
     890  getLattice() == L;
     891
     892  // its hermite NF:
     893  print(getLattice("hermite"));
     894
     895  kill L, M;
     896
     897  // ----------- extreme case ------------ //
     898
     899  // Weights of variables
     900  intmat M[1][3] =
     901    1,  -1, 10;
     902
     903  // Torsion:
     904  intmat L[1][1] =
     905    0;
     906   
     907  // attaches M & L to R (==basering):
     908  setBaseMultigrading(M); // Grading: Z^3
     909
     910  // Torsion is accessible via "getLattice()":
     911  getLattice() == L;
     912
     913  // its hermite NF:
     914  print(getLattice("hermite"));
     915}
     916
     917proc getGradedGenerator(def m, int i)
     918"
     919returns m[i], but with grading
     920"
     921{
     922  if( typeof(m) == "ideal" )
     923  {
     924    return (m[i]);
     925  }
     926
     927  if( typeof(m) == "module" )
     928  {
     929    def v = getModuleGrading(m);
     930   
     931    return ( setModuleGrading(m[i],v) );
     932  }
     933
     934  ERROR("m is expected to be an ideal or a module");
     935}
     936
    399937
    400938/******************************************************/
     
    410948
    411949  def V = attrib(m, attrModuleGrading);
    412 
     950 
    413951  if( typeof(V) != "intmat" )
    414952  {
     
    419957      return (VV);
    420958    }
    421 
     959     
    422960    ERROR("Sorry: vector or module need module-grading-matrix! See 'getModuleGrading'.");
    423961  }
     
    432970    ERROR("Sorry wrong width of V: " + string(ncols(V)));
    433971  }
    434 
     972   
    435973  return (V);
    436974}
     
    438976{
    439977  "EXAMPLE:"; echo=2;
    440 
     978   
    441979   ring R = 0, (x,y), dp;
    442980   intmat M[2][2]=
     
    446984     1,  2,  3,  4, 0,
    447985     0, 10, 20, 30, 1;
    448 
     986   
    449987   setBaseMultigrading(M, T);
    450 
     988   
    451989   ideal I = x, y, xy^5;
    452    isHomogenous(I);
    453 
     990   isHomogeneous(I);
     991 
    454992   intmat V = mDeg(I); print(V);
    455993
    456994   module S = syz(I); print(S);
    457 
     995   
    458996   S = setModuleGrading(S, V);
    459997
    460998   getModuleGrading(S) == V;
    461 
    462    vector v = setModuleGrading(S[1], V);
     999   
     1000   vector v = getGradedGenerator(S, 1);
    4631001   getModuleGrading(v) == V;
    464    isHomogenous(v);
    465    print( mDeg(v) );
    466 
    467    isHomogenous(S);
     1002   isHomogeneous(v);
     1003   print( mDeg(v) );   
     1004   
     1005   isHomogeneous(S);
    4681006   print( mDeg(S) );
    4691007}
     
    4741012PURPOSE: attaches the multiweights of free module generators to 'm'
    4751013WARNING: The method does not verify whether the multigrading makes the
    476          module/vector homogenous. One can do that using isHomogenous(m).
     1014         module/vector homogeneous. One can do that using isHomogeneous(m).
    4771015EXAMPLE: example setModuleGrading; shows an example
    4781016"
     
    4841022  if(nrows(G) != nrows(R)){ ERROR("Incompatible gradings.");}
    4851023  if(ncols(G) < nrows(m)){ ERROR("Multigrading does not fit to module.");}
    486 
     1024 
    4871025  attrib(m, attrModuleGrading, G);
    4881026  return(m);
     
    4911029{
    4921030  "EXAMPLE:"; echo=2;
    493 
     1031   
    4941032   ring R = 0, (x,y), dp;
    4951033   intmat M[2][2]=
     
    4991037     1,  2,  3,  4, 0,
    5001038     0, 10, 20, 30, 1;
    501 
     1039   
    5021040   setBaseMultigrading(M, T);
    503 
     1041   
    5041042   ideal I = x, y, xy^5;
    5051043   intmat V = mDeg(I);
    506 
     1044   
    5071045   // V == M; modulo T
    5081046   print(V);
    5091047
    5101048   module S = syz(I);
    511 
     1049   
    5121050   S = setModuleGrading(S, V);
    5131051   getModuleGrading(S) == V;
    5141052
    5151053   print(S);
    516 
    517    vector v = S[1]; v = setModuleGrading(v, V);
     1054   
     1055   vector v = getGradedGenerator(S, 1);
    5181056   getModuleGrading(v) == V;
    5191057
    520    print( mDeg(v) );
    521 
    522    isHomogenous(S);
     1058   print( mDeg(v) );   
     1059
     1060   isHomogeneous(S);
    5231061
    5241062   print( mDeg(S) );
     
    5261064
    5271065
    528 
     1066proc mDegTensor(module m, module n){
     1067  matrix M = m;
     1068  matrix N = n;
     1069  intmat gm = getModuleGrading(m);
     1070  intmat gn = getModuleGrading(n);
     1071  int grows = nrows(gm);
     1072  int mr = nrows(M);
     1073  int mc = ncols(M);
     1074  if(rank(M) == 0){ mc = 0;}
     1075  int nr = nrows(N);
     1076  int nc = ncols(N);
     1077  if(rank(N) == 0){ nc = 0;}
     1078  intmat gresult[nrows(gm)][mr*nr];
     1079  matrix result[mr*nr][mr*nc+mc*nr];
     1080  int i, j;
     1081  int column = 1;
     1082  for(i = 1; i<=mr; i++){
     1083    for(j = 1; j<=nr; j++){
     1084      gresult[1..grows,(i-1)*nr+j] = gm[1..grows,i]+gn[1..grows,j];
     1085    }
     1086  }
     1087  //gresult;
     1088  if( nc!=0 ){
     1089    for(i = 1; i<=mr; i++)
     1090    {
     1091      result[((i-1)*nr+1)..(i*nr),((i-1)*nc+1)..(i*nc)] = N[1..nr,1..nc];
     1092    }
     1093  }
     1094  list rownumbers, colnumbers;
     1095  //print(result);
     1096  if( mc!=0 ){
     1097    for(j = 1; j<=nr; j++)
     1098    {
     1099      rownumbers = nr*(0..(mr-1))+j*(1:mr);
     1100      colnumbers = ((mr*nc+j):mc)+nr*(0..(mc-1));
     1101      result[rownumbers[1..mr],colnumbers[1..mc] ] = M[1..mr,1..mc];
     1102    }
     1103  }
     1104  module res = result;
     1105  res = setModuleGrading(res, gresult);
     1106  //getModuleGrading(res);
     1107  return(res);
     1108}
     1109example
     1110{
     1111"EXAMPLE: ";echo=2;
     1112ring r = 0,(x),dp;
     1113intmat g[2][1]=1,1;
     1114setBaseMultigrading(g);
     1115matrix m[5][3]=1,2,3,4,5,6,7,8,9,10,11,12,13,14,15;
     1116matrix n[3][2]=x,x2,x3,x4,x5,x6;
     1117module mm = m;
     1118module nn = n;
     1119intmat gm[2][5]=1,2,3,4,5,0,0,0,0,0;
     1120intmat gn[2][3]=0,0,0,1,2,3;
     1121mm = setModuleGrading(mm, gm);
     1122nn = setModuleGrading(nn, gn);
     1123module mmtnn = mDegTensor(mm, nn);
     1124print(mmtnn);
     1125getModuleGrading(mmtnn);
     1126LIB "homolog.lib";
     1127module tt = tensorMod(mm,nn);
     1128print(tt);
     1129
     1130kill m, mm, n, nn, gm, gn;
     1131
     1132matrix 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;
     1133matrix n[2][4] = 1, 2, 3, 4, x, x2, x3, x4;
     1134module mm = m;
     1135module nn = n;
     1136print(mm);
     1137print(nn);
     1138intmat gm[2][7] = 1, 2, 3, 4, 5, 6, 7, 0, 0, 0, 0, 0, 0, 0;
     1139intmat gn[2][2] = 0, 0, 1, 2;
     1140mm = setModuleGrading(mm, gm);
     1141nn = setModuleGrading(nn, gn);
     1142module mmtnn = mDegTensor(mm, nn);
     1143print(mmtnn);
     1144getModuleGrading(mmtnn);
     1145matrix a = mmtnn;
     1146matrix b = tensorMod(mm, nn);
     1147print(a-b);
     1148
     1149}
     1150
     1151proc mDegTor(int i, module m, module n)
     1152{
     1153  def res = mDegResolution(n, 0, 1);
     1154  //print(res);
     1155  list l = res;
     1156  if(size(l)<i){ return(0);}
     1157  else
     1158  {
     1159   
     1160    matrix fd[nrows(m)][0];
     1161    matrix fd2[nrows(l[i+1])][0];
     1162    matrix fd3[nrows(l[i])][0];
     1163
     1164    module freedim = fd;
     1165    module freedim2 = fd2;
     1166    module freedim3 = fd3;
     1167
     1168    freedim = setModuleGrading(freedim,getModuleGrading(m));
     1169    freedim2 = setModuleGrading(freedim2,getModuleGrading(l[i+1]));
     1170    freedim3 = setModuleGrading(freedim3, getModuleGrading(l[i]));
     1171   
     1172    module mimag = mDegTensor(freedim3, m);
     1173    //"mimag ok.";
     1174    module mf = mDegTensor(l[i], freedim);
     1175    //"mf ok.";
     1176    module mim1 = mDegTensor(freedim2 ,m);
     1177    module mim2 = mDegTensor(l[i+1],freedim);
     1178    //"mim1+2 ok.";
     1179    module mker = mDegModulo(mf,mimag);
     1180    //"mker ok.";
     1181    module mim = mim1,mim2;
     1182    mim = setModuleGrading(mim, getModuleGrading(mim1));
     1183    //"mim: r: ",nrows(mim)," c: ",ncols(mim);
     1184    //"mim1: r: ",nrows(mim1)," c: ",ncols(mim1);
     1185    //"mim2: r: ",nrows(mim2)," c: ",ncols(mim2);
     1186    //matrix mimmat = mim;
     1187    //matrix mimmat1[16][4]=mimmat[1..16,25..28];
     1188    //print(mimmat1-matrix(mim2));
     1189    return(mDegModulo(mker,mim));
     1190    //return(0);
     1191  }
     1192  return(0);
     1193}
     1194example
     1195{
     1196"EXAMPLE: ";echo=2;
     1197LIB "homolog.lib";
     1198ring r = 0,(x_(1..4)),dp;
     1199intmat g[2][4]=1,1,0,0,0,1,1,-1;
     1200setBaseMultigrading(g);
     1201ideal i = maxideal(1);
     1202module m = mDegSyzygy(i);
     1203module rt = Tor(2,m,m);
     1204module mDegT = mDegTor(2,m,m);
     1205print(matrix(rt)-matrix(mDegT));
     1206/*
     1207ring r = 0,(x),dp;
     1208intmat g[2][1]=1,1;
     1209setBaseMultigrading(g);
     1210matrix m[5][3]=1,2,3,4,5,6,7,8,9,10,11,12,13,14,15;
     1211matrix n[3][2]=x,x2,x3,x4,x5,x6;
     1212module mm = m;
     1213module nn = n;
     1214intmat gm[2][5]=1,1,1,1,1,1,1,1,1,1,1;
     1215intmat gn[2][3]=0,-2,-4,0,-2,-4;
     1216mm = setModuleGrading(mm, gm);
     1217nn = setModuleGrading(nn, gn);
     1218isHomogeneous(mm,"checkGens");
     1219isHomogeneous(nn,"checkGens");
     1220mDegTor(1,mm, nn);
     1221
     1222kill m, mm, n, nn, gm, gn;
     1223
     1224matrix 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;
     1225matrix n[2][4] = 1, 2, 3, 4, x, x2, x3, x4;
     1226module mm = m;
     1227module nn = n;
     1228print(mm);
     1229print(nn);
     1230intmat gm[2][7] = 1, 2, 3, 4, 5, 6, 7, 0, 0, 0, 0, 0, 0, 0;
     1231intmat gn[2][2] = 0, 0, 1, 2;
     1232mm = setModuleGrading(mm, gm);
     1233nn = setModuleGrading(nn, gn);
     1234module mmtnn = mDegTensor(mm, nn);
     1235*/
     1236}
     1237
     1238
     1239/******************************************************/
     1240proc isGroupHomomorphism(def L1, def L2, intmat A)
     1241"USAGE: gisGoupHomomorphism(L1,L2,A); L1 and L2 are groups, A is an integer matrix
     1242PURPOSE: checks whether A defines a group homomorphism phi: L1 --> L2
     1243RETURN: int, 1 if A defines the homomorphism and 0 otherwise
     1244EXAMPLE: example isGroupHomomorphism; shows an example
     1245"
     1246{
     1247  // TODO: L1, L2
     1248  if( (ncols(A) != nrows(L1)) or (nrows(A) != nrows(L2)) )
     1249  {
     1250    ERROR("Incompatible sizes!");
     1251  }
     1252
     1253  intmat im = A * L1;
     1254 
     1255  return  (areZeroElements(im, L2));
     1256}
     1257example
     1258{
     1259  "EXAMPLE:"; echo=2;
     1260
     1261   intmat L1[4][1]=
     1262     0,
     1263     0,
     1264     0,
     1265     2;
     1266     
     1267   intmat L2[3][2]=
     1268     0, 0,
     1269     2, 0,
     1270     0, 3;
     1271
     1272  intmat A[3][4] =
     1273     1, 2, 3, 0,
     1274     7, 0, 0, 0,
     1275     1, 2, 0, 3;
     1276  print( A );
     1277
     1278  isGroupHomomorphism(L1, L2, A);
     1279
     1280  intmat B[3][4] =
     1281     1, 2, 3, 0,
     1282     7, 0, 0, 0,
     1283     1, 2, 0, 2;
     1284  print( B );
     1285
     1286  isGroupHomomorphism(L1, L2, B); // Not a homomorphism!
     1287}
    5291288
    5301289/******************************************************/
    5311290proc isTorsionFree()
    5321291"USAGE: isTorsionFree()
    533 PURPOSE: Determines whether the multigrading attached to the current ring is torsion-free.
     1292PURPOSE: Determines whether the multigrading attached to the current ring is free.
    5341293RETURN: boolean, the result of the test
    5351294EXAMPLE: example isTorsionFree; shows an example
    5361295"
    5371296{
    538 
    539   intmat H = hermite(transpose(getTorsion("hermite"))); // TODO: ?cache it?  //******
     1297  intmat H = smithNormalForm(getLattice()); // TODO: ?cache it?  //******
    5401298
    5411299  int i, j;
     
    5501308
    5511309    if(H[j, i]!=0)
    552     {
     1310    { 
    5531311      d=d*H[j, i];
    5541312    }
     
    5561314
    5571315  if( (d*d)==1 )
    558   {
     1316  { 
    5591317    return(1==1);
    5601318  }
     
    5721330     1, 2, 3, 4, 0,
    5731331     0,10,20,30, 1;
    574 
     1332   
    5751333   setBaseMultigrading(M,T);
    576 
    577    // Is the resulting group torsion free?
     1334   
     1335   // Is the resulting group free?
    5781336   isTorsionFree();
    5791337
     
    5821340
    5831341   ring R=0,(x,y,z),dp;
    584    intmat A[3][3] =
     1342   intmat A[3][3] = 
    5851343     1,0,0,
    5861344     0,1,0,
     
    5911349     1,2,0,3;
    5921350   setBaseMultigrading(A,B);
    593    // Is the resulting group torsion free?
     1351   // Is the resulting group free?
    5941352   isTorsionFree();
    5951353
     
    5981356
    5991357
     1358static proc gcdcomb(int a, int b)
     1359{
     1360  // a;
     1361  // b;
     1362  intvec av = a,1,0;
     1363  intvec bv = b,0,1;
     1364  intvec save;
     1365  while(av[1]*bv[1] != 0)
     1366  {
     1367    bv = bv - (bv[1] - bv[1]%av[1])/av[1] * av;
     1368    save = bv;
     1369    bv = av;
     1370    av = save;
     1371  }
     1372  if(bv[1] < 0)
     1373  {
     1374    bv = -bv;
     1375  }
     1376  return(bv);
     1377}
     1378
     1379
     1380proc lll(def A)
     1381"
     1382The lll algorithm of lll.lib only works for lists of vectors.
     1383Maybe one should rescript it for matrices. This method will
     1384convert a matrix to a list, plug it into lll and make the result
     1385a matrix and return it.
     1386"
     1387{
     1388  if(typeof(A) == "list")
     1389  {
     1390    int sizeA= size (A);
     1391    if (sizeA == 0)
     1392    {
     1393      return (A);
     1394    }
     1395    if (typeof (A [1]) != "intvec")
     1396    {
     1397      ERROR("Unrecognized type.");
     1398    }
     1399    int columns= size (A [1]);
     1400    int i;
     1401    for (i= 2; i <= sizeA; i++)
     1402    {
     1403      if (typeof (A[i]) != "intvec")
     1404      {
     1405        ERROR("Unrecognized type.");
     1406      }
     1407      if (size (A [i]) != columns)
     1408      {
     1409        ERROR ("expected equal dimension");
     1410      }
     1411    }
     1412    int j;
     1413    intmat m [columns] [sizeA];
     1414    for (i= 1; i <= sizeA; i++)
     1415    {
     1416      for (j= 1; j <= columns; j++)
     1417      {
     1418        m[i,j]= A[i] [j];
     1419      }
     1420    }
     1421    m= system ("LLL", m);
     1422    list result= list();
     1423
     1424    for (i= 1; i <= sizeA; i++)
     1425    {
     1426      intvec buf= intvec (m[i , 1]);
     1427      for (j= 2; j <= columns; j++)
     1428      {
     1429        buf= buf,intvec (m [i, j]);
     1430      }
     1431      result= result+ list (buf);
     1432     
     1433    }
     1434    return(result);
     1435  }
     1436  else
     1437  {
     1438    if(typeof(A) == "intmat")
     1439    {
     1440      A= system ("LLL", A);
     1441      return(A);
     1442    }
     1443    else
     1444    {
     1445      ERROR("Unrecognized type.");
     1446    }
     1447  }
     1448}
     1449
     1450example
     1451{
     1452
     1453"EXAMPLE:";
     1454
     1455ring R = 0,x,dp;
     1456intmat 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;
     1457lll(m);
     1458list 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);
     1459lll(l);
     1460}
     1461
     1462
     1463proc smithNormalForm(intmat A, list #)
     1464"
     1465This method returns 3 Matrices P, D and Q such that D = P*A*Q.
     1466WARNING: This might not be what you expect.
     1467"
     1468{
     1469  list l1 = hermiteNormalForm(A, 5);
     1470  // l1;
     1471  intmat B = transpose(l1[1]);
     1472  list l2 = hermiteNormalForm(B, 5);
     1473  // l2;
     1474  intmat P = transpose(l2[2]);
     1475  intmat D = transpose(l2[1]);
     1476  intmat Q = l1[2];
     1477  int cc = ncols(D);
     1478  int rr = nrows(D);
     1479  intmat transform;
     1480  int k = 1;
     1481  int a, b, c;
     1482  // D;
     1483  intvec v;
     1484  if((cc==1)||(rr==1)){
     1485    if(size(#)==0)
     1486    {
     1487      return(D);
     1488    } else
     1489    {
     1490      return(list(P,D,Q));
     1491    }
     1492  }
     1493  while(D[k+1,k+1] !=0){
     1494    if(D[k+1,k+1]%D[k,k]!=0){
     1495      b = D[k, k]; c = D[k+1, k+1];
     1496      v = gcdcomb(D[k,k],D[k+1,k+1]);
     1497      transform = unitMatrix(cc);
     1498      transform[k+1,k] = 1;
     1499      a = -v[3]*D[k+1,k+1]/v[1];
     1500      transform[k, k+1] = a;
     1501      transform[k+1, k+1] = a+1;
     1502      //det(transform);
     1503      D = D*transform;
     1504      Q = Q*transform;
     1505      //D;
     1506      transform = unitMatrix(rr);
     1507      transform[k,k] = v[2];
     1508      transform[k,k+1] = v[3];
     1509      transform[k+1,k] = -c/v[1];
     1510      transform[k+1,k+1] = b/v[1];
     1511      D = transform * D;
     1512      P = transform * P;
     1513      //" ";
     1514      //D;
     1515      //"small transform: ", det(transform);
     1516      //transform;
     1517      k=0;
     1518    }
     1519    k++;
     1520    if((k==rr) || (k==cc)){
     1521      break;
     1522    }
     1523  }
     1524  //"here is the size ",size(#);
     1525  if(size(#) == 0){
     1526    return(D);
     1527  } else {
     1528    return(list(P, D, Q));
     1529  }
     1530}
     1531example
     1532{
     1533"EXAMPLE: "; echo=2;
     1534
     1535intmat A[5][7] =
     15361,0,1,0,-2,9,-71,
     15370,-24,248,-32,-96,448,-3496,
     15380,4,-42,4,-8,30,-260,
     15390,0,0,18,-90,408,-3168,
     15400,0,0,-32,224,-1008,7872;
     1541
     1542list l = smithNormalForm(A, 5);
     1543
     1544l;
     1545l[1]*A*l[3];
     1546det(l[1]);
     1547det(l[3]);
     1548}
     1549
    6001550
    6011551/******************************************************/
    602 proc hermite(intmat A)
    603 "USAGE: hermite( A );
    604 PURPOSE: Computes the (lower triangular) Hermite Normal Form
     1552proc hermiteNormalForm(intmat A, list #)
     1553"USAGE: hermiteNormalForm( A );
     1554PURPOSE: Computes the (lower triangular) Hermite Normal Form 
    6051555           of the matrix A by column operations.
    6061556RETURN: intmat, the Hermite Normal Form of A
    607 EXAMPLE: example hermite; shows an example
    608 "
    609 {
    610   intmat g;
    611   int i,j,k, save, d, a1, a2, c, l;
    612   c=0;
    613   l=0;
    614   for(i=1; ((i+l)<=nrows(A))&&((i+l)<=ncols(A)); i++)
    615   {
    616     //i;
    617     if(A[i+l, i]==0)
    618     {
    619       for(j=i;j<=ncols(A); j++)
    620       {
    621         if(A[i+l, j]!=0)
     1557EXAMPLE: example hermiteNormalForm; shows an example
     1558"
     1559{
     1560 
     1561  int row, column, i, j;
     1562  int rr = nrows(A);
     1563  int cc = ncols(A);
     1564  intvec savev, gcdvec, v1, v2;
     1565  intmat q = unitMatrix(cc);
     1566  intmat transform;
     1567  column = 1;
     1568  for(row = 1; (row<=rr)&&(column<=cc); row++)
     1569    {
     1570      if(A[row,column]==0)
    6221571        {
    623           for(k=1;k<=nrows(A);k++)
    624           {
    625             save=A[k, i];
    626             A[k, i]=A[k, j];
    627             A[k, j]=save;
    628           }
    629           break;
     1572          for(j = column; j<=cc; j++)
     1573            {
     1574              if(A[row, j]!=0)
     1575                {
     1576                  transform = unitMatrix(cc);
     1577                  transform[j,j] = 0;
     1578                  transform[column, column] = 0;
     1579                  transform[column,j] = 1;
     1580                  transform[j,column] = 1;
     1581                  q = q*transform;
     1582                  A = A*transform;
     1583                  break;
     1584                }
     1585            }
    6301586        }
    631         if(j==ncols(A)){ c=1; l=l+1; }
    632       }
    633     }
    634 
    635     if((i+l)>nrows(A)){ break; }
    636 
    637     if(A[i+l, i]<0)
    638     {
    639       for(k=1;k<=nrows(A);k++){ A[k, i]=A[k, i]*(-1); }
    640     }
    641 
    642     if(c==0)
    643     {
    644       for(j=i+1;j<=ncols(A);j++)
    645       {
    646         //A;
    647         if(A[i+l, j]<0)
     1587      if(A[row,column] == 0)
    6481588        {
    649           for(k=1;k<=nrows(A);k++){ A[k, j]=(-1)*A[k, j];}
     1589          row++;
     1590          continue;
    6501591        }
    651 
    652         if(A[i+l, i]==1)
     1592      for(j = column+1; j<=cc; j++)
    6531593        {
    654           d=A[i+l, j];
    655           for(k=1;k<=nrows(A);k++)
    656           {
    657             A[k, j]=A[k, j]-d*A[k, i];
    658           }
     1594          if(A[row, j]!=0)
     1595            {
     1596              gcdvec = gcdcomb(A[row,column],A[row,j]);
     1597              // gcdvec;
     1598              // typeof(A[1..rr,column]);
     1599              v1 = A[1..rr,column];
     1600              v2 = A[1..rr,j];
     1601              transform = unitMatrix(cc);
     1602              transform[j,j] = v1[row]/gcdvec[1];
     1603              transform[column, column] = gcdvec[2];
     1604              transform[column,j] = -v2[row]/gcdvec[1];
     1605              transform[j,column] = gcdvec[3];
     1606              q = q*transform;
     1607              A = A*transform;
     1608              // A;
     1609            }
    6591610        }
    660         else
     1611      if(A[row,column]<0)
    6611612        {
    662           while((A[i+l, i] * A[i+l, j])!=0)
    663           {
    664             if(A[i+l, i]> A[i+l, j])
    665             {
    666 
    667               for(k=1;k<=nrows(A);k++)
    668               {
    669                 save=A[k, i];
    670                 A[k, i]=A[k, j];
    671                 A[k, j]=save;
    672               }
    673             }
    674             a1=A[i+l, j]%A[i+l,i];
    675             a2=A[i+l, j]-a1;
    676             d=a2/A[i+l, i];
    677             for(k=1;k<=nrows(A);k++)
    678             {
    679               A[k, j]=A[k, j]- d*A[k, i];
    680             }
    681           }
     1613          transform = unitMatrix(cc);
     1614          transform[column,column] = -1;
     1615          q = q*transform;
     1616          A = A*transform;
     1617        }
     1618      for( j=1; j<column; j++){
     1619        if(A[row, j]!=0){
     1620          transform = unitMatrix(cc);
     1621          transform[column, j] = (-A[row,j]+A[row, j]%A[row, column])/A[row, column];
     1622          if(A[row,j]<0){
     1623            transform[column,j]=transform[column,j]+1;}
     1624          q = q*transform;
     1625          A = A*transform;
    6821626        }
    6831627      }
    684       for(j=1;j<i;j++)
    685       {
    686         a1=A[i+l, j]%A[i+l,i];
    687         d=(A[i+l, j]-a1)/A[i+l, i];
    688         for(k=1;k<=nrows(A);k++){ A[k, j]=A[k, j]-d*A[k, i];}
    689       }
    690     }
    691     c=0;
    692   }
    693   return( A);
     1628      column++;
     1629    }
     1630  if(size(#) > 0){
     1631    return(list(A, q));
     1632  }
     1633  return(A);
    6941634}
    6951635example
     
    6971637  "EXAMPLE:"; echo=2;
    6981638
    699    intmat M[2][5] =
     1639   intmat M[2][5] = 
    7001640     1, 2, 3, 4, 0,
    7011641     0,10,20,30, 1;
    7021642
    7031643   // Hermite Normal Form of M:
    704    print(hermite(M));
    705 
    706    intmat T[3][4] =
     1644   print(hermiteNormalForm(M));
     1645
     1646   intmat T[3][4] = 
    7071647     3,3,3,3,
    7081648     2,1,3,0,
     
    7101650
    7111651   // Hermite Normal Form of T:
    712    print(hermite(T));
    713 
    714    intmat A[4][5] =
     1652   print(hermiteNormalForm(T));
     1653
     1654   intmat A[4][5] = 
    7151655     1,2,3,2,2,
    7161656     1,2,3,4,0,
     
    7191659
    7201660   // Hermite Normal Form of A:
    721    print(hermite(A));
    722 }
    723 
    724 
    725 /******************************************************/
    726 proc isTorsionElement(intvec mdeg)
    727 "USAGE: isTorsionElement(intvec mdeg);
    728 PURPOSE: For a integer vector mdeg representing the multidegree of some polynomial
    729 or vector this method computes if the multidegree is contained in the torsion
    730 group, i.e. if it is zero in the multigrading.
    731 EXAMPLE: example isTorsionElement; shows an example
    732 "
    733 {
    734   intmat H = getTorsion("hermite");
    735   int x, k, i;
    736 
    737   int r = nrows(H);
    738   int c = ncols(H);
    739 
    740   int rr = nrows(mdeg);
    741 
    742   for( i = 1; (i<=r) && (i<=c); i++)
    743   {
    744     if(H[i, i]!=0)
    745     {
    746       x = mdeg[i]%H[i, i];
    747 
    748       if(x!=0)
    749       {
    750         return(1==0);
    751       }
    752 
    753       x = mdeg[i] / H[i,i];
    754 
    755       for( k=1; k <= rr; k++)
    756       {
    757         mdeg[k] = mdeg[k] - x*H[k,i];
    758       }
    759     }
    760   }
    761 
    762   return( mdeg == 0 );
    763 
    764 }
     1661   print(hermiteNormalForm(A));
     1662}
     1663
     1664proc areZeroElements(intmat m, list #)
     1665"same as isZeroElement but for an integer matrix considered as a collection of columns"
     1666{
     1667  int r = nrows(m);
     1668  int i = ncols(m);
     1669
     1670  intvec v;
     1671 
     1672  for( ; i > 0; i-- )
     1673  {
     1674    v = m[1..r, i];
     1675    if( !isZeroElement(v, #) )
     1676    {
     1677      return (0);
     1678    }
     1679  }
     1680  return(1);
     1681}
     1682
    7651683example
    7661684{
     
    7761694    1;
    7771695
     1696  intmat tt[2][1]=
     1697    1,
     1698    -1;
     1699
    7781700  setBaseMultigrading(g,t);
    7791701
     
    7851707  poly f = z2;
    7861708
    787   intvec v1 = mDeg(a) - mDeg(b);
    788   v1;
    789   isTorsionElement(v1);
    790 
    791   intvec v2 = mDeg(a) - mDeg(c);
    792   v2;
    793   isTorsionElement(v2);
    794 
    795   intvec v3 = mDeg(e) - mDeg(f);
    796   v3;
    797   isTorsionElement(v3);
    798 
    799   intvec v4 = mDeg(c) - mDeg(d);
    800   v4;
    801   isTorsionElement(v4);
     1709 intmat m[5][2]=mDeg(a)-mDeg(b),mDeg(b)-mDeg(c),mDeg(c)-mDeg(d),mDeg(d)-mDeg(e),mDeg(e)-mDeg(f);
     1710 m=transpose(m);
     1711 areZeroElementes(m);
     1712 areZeroElementes(m,tt);
    8021713}
    8031714
    8041715
    8051716/******************************************************/
    806 proc defineHomogenous(poly f, list #)
    807 "USAGE: defineHomogenous(f[, G]); polynomial f, integer matrix G
    808 PURPOSE: Yields a matrix which has to be appended to the torsion matrix to make the
    809 polynomial f homogenous in the grading by grad.
    810 EXAMPLE: example defineHomogenous; shows an example
     1717proc isZeroElement(intvec mdeg, list #)
     1718"USAGE: isZeroElement(d, [T]); intvec d, group T
     1719PURPOSE: For a integer vector mdeg representing the multidegree of some polynomial
     1720or vector this method computes if the multidegree is contained in the grading group
     1721group (either set globally or given as an optional argument), i.e. if it is zero in the multigrading.
     1722EXAMPLE: example isZeroElement; shows an example
    8111723"
    8121724{
     
    8151727    if( typeof(#[1]) == "intmat" )
    8161728    {
    817       intmat grad = #[1];
    818     }
    819   }
    820 
    821   if( !defined(grad) )
    822   {
    823     intmat grad = getVariableWeights();
    824   }
    825 
    826   intmat newtor[nrows(grad)][size(f)-1];
    827   int i,j;
    828   intvec l = grad*leadexp(f);
     1729      intmat H = hermiteNormalForm(#[1]);
     1730    } else
     1731    {
     1732      if( typeof(#[1]) == "list" )
     1733      {
     1734        list L = #[1];
     1735        intmat H = attrib(L, "hermite"); // todo
     1736      }
     1737    }
     1738   
     1739  }
     1740  if( !defined(H) )
     1741  {
     1742    intmat H = getLattice("hermite");
     1743  }
     1744
     1745  int x, k, i, row;
     1746
     1747  int r = nrows(H);
     1748  int c = ncols(H);
     1749
     1750  int rr = nrows(mdeg);
     1751  row = 1;
    8291752  intvec v;
    830   for(i=2; i <= size(f); i++)
    831   {
    832     v = grad * leadexp(f[i]) - l;
    833     for( j=1; j<=size(v); j++)
    834     {
    835       newtor[j,i-1] = v[j];
    836     }
    837   }
    838   return(newtor);
    839 }
    840 example
    841 {
    842   "EXAMPLE:"; echo=2;
    843 
    844   ring r =0,(x,y,z),dp;
    845   intmat grad[2][3] =
    846     1,0,1,
    847     0,1,1;
    848 
    849   setBaseMultigrading(grad);
    850 
    851   poly f = x2y3-z5+x-3zx;
    852 
    853   intmat M = defineHomogenous(f);
    854   M;
    855   defineHomogenous(f, grad) == M;
    856 
    857   isHomogenous(f);
    858   setBaseMultigrading(grad, M);
    859   isHomogenous(f);
    860 }
    861 
    862 /******************************************************/
    863 proc pushForward(map f)
    864 "USAGE: pushForward(f);
    865 PURPOSE: Computes the finest grading of the image ring which makes the map f
    866 a map of graded rings. The group map between the two grading groups is given
    867 by transpose( (Id, 0) ). Pay attention that the group spanned by the columns of
    868 the torsion matrix may not be a subgroup of the grading group. Still all columns
    869 are needed to find the correct image of the preimage gradings.
    870 EXAMPLE: example pushForward; shows an example
    871 "
    872 {
    873 
    874   int k,i,j;
    875 //  f;
    876 
    877 //  listvar();
    878   def pre = preimage(f);
    879 
    880 //  "pre: ";  pre;
    881 
    882   intmat oldgrad=getVariableWeights(pre);
    883   intmat oldtor=getTorsion(pre);
    884 
    885   int n=nvars(pre);
    886   int np=nvars(basering);
    887   int p=nrows(oldgrad);
    888   int pp=p+np;
    889 
    890   intmat newgrad[pp][np];
    891 
    892   for(i=1;i<=np;i++){ newgrad[p+i,i]=1;}
    893 
    894   //newgrad;
    895 
    896 
    897 
    898   list newtor;
    899   intmat toadd;
    900   int columns=0;
    901 
    902   intmat toadd1[pp][n];
    903   intvec v;
    904   poly im;
    905 
    906   for(i=1;i<=p;i++){
    907     for(j=1;j<=n;j++){ toadd1[i,j]=oldgrad[i,j];}
    908   }
    909 
    910   for(i=1;i<=n;i++){
    911     im=f[i];
    912     //im;
    913     toadd = defineHomogenous(im, newgrad);
    914     newtor=insert(newtor,toadd);
    915     columns=columns+ncols(toadd);
    916 
    917     v=leadexp(f[i]);
    918     for(j=p+1;j<=p+np;j++){ toadd1[j,i]=-v[j-p];}
    919   }
    920 
    921   newtor=insert(newtor,toadd1);
    922   columns=columns+ncols(toadd1);
    923 
    924 
    925   if(typeof(basering)=="qring"){
    926     //"Entering qring";
    927     ideal a=ideal(basering);
    928     for(i=1;i<=size(a);i++){
    929       toadd = defineHomogenous(a[i], newgrad);
    930       //toadd;
    931       columns=columns+ncols(toadd);
    932       newtor=insert(newtor,toadd);
    933     }
    934   }
    935 
    936   //newtor;
    937   intmat imofoldtor[pp][ncols(oldtor)];
    938   for(i=1; i<=nrows(oldtor);i++){
    939     for(j=1; j<=ncols(oldtor); j++){
    940       imofoldtor[i,j]=oldtor[i,j];
    941     }
    942   }
    943 
    944   columns=columns+ncols(oldtor);
    945   newtor=insert(newtor, imofoldtor);
    946 
    947   intmat torsion[pp][columns];
    948   columns=0;
    949   for(k=1;k<=size(newtor);k++){
    950     for(i=1;i<=pp;i++){
    951       for(j=1;j<=ncols(newtor[k]);j++){torsion[i,j+columns]=newtor[k][i,j];}
    952     }
    953     columns=columns+ncols(newtor[k]);
    954   }
    955 
    956   torsion=hermite(torsion);
    957   intmat result[pp][pp];
    958   for(i=1;i<=pp;i++){
    959     for(j=1;j<=pp;j++){result[i,j]=torsion[i,j];}
    960   }
    961 
    962   setBaseMultigrading(newgrad, result);
    963 
    964 }
    965 example
    966 {
    967   "EXAMPLE:"; echo=2;
    968 
    969   ring r = 0,(x,y,z),dp;
    970 
    971 
    972 
    973   // Setting degrees for preimage ring.;
    974   intmat grad[3][3] =
    975     1,0,0,
    976     0,1,0,
    977     0,0,1;
    978 
    979   setBaseMultigrading(grad);
    980 
    981   // grading on r:
    982   getVariableWeights();
    983   getTorsion();
    984 
    985   // only for the purpose of this example
    986   if( voice > 1 ){ /*keepring(r);*/ export(r); }
    987 
    988   ring R = 0,(a,b),dp;
    989   ideal i = a2-b2+a6-b5+ab3,a7b+b15-ab6+a6b6;
    990 
    991   // The quotient ring by this ideal will become our image ring.;
    992   qring Q = std(i);
    993 
    994   listvar();
    995 
    996   map f = r,-a2b6+b5+a3b+a2+ab,-a2b7-3a2b5+b4+a,a6-b6-b3+a2; f;
    997 
    998 
    999   // TODO: Unfortunately this is not a very spectacular example...:
    1000   // Pushing forward f:
    1001   pushForward(f);
    1002 
    1003   // due to pushForward we have got new grading on Q
    1004   getVariableWeights();
    1005   getTorsion();
    1006 
    1007 
    1008   // only for the purpose of this example
    1009   if( voice > 1 ){ kill r; }
    1010 
    1011 }
    1012 
    1013 
    1014 /******************************************************/
    1015 proc equalMDeg(intvec exp1, intvec exp2, list #)
    1016 "USAGE: equalMDeg(exp1, exp2[, V]); intvec exp1, exp2, intmat V
    1017 PURPOSE: Tests if the exponent vectors of two monomials (given by exp1 and exp2)
    1018 represent the same multidegree.
    1019 NOTE: the integer matrix V encodes multidegrees of module components,
    1020 if module component is present in exp1 and exp2
    1021 EXAMPLE: example equalMDeg; shows an example
    1022 "
    1023 {
    1024   if( size(exp1) != size(exp2) )
    1025   {
    1026     ERROR("Sorry: we cannot compare exponents comming from a polynomial and a vector yet!");
    1027   }
    1028 
    1029   if( exp1 == exp2)
    1030   {
    1031     return (1==1);
    1032   }
    1033 
    1034 
    1035 
    1036   intmat M = getVariableWeights();
    1037 
    1038   if( nrows(exp1) > ncols(M) ) // vectors => last exponent is the module component!
    1039   {
    1040     if( (size(#) == 0) or (typeof(#[1])!="intmat") )
    1041     {
    1042       ERROR("Sorry: wrong or missing module-unit-weights-matrix V!");
    1043     }
    1044     intmat V = #[1];
    1045 
    1046     // typeof(V); print(V);
    1047 
    1048     int N = ncols(M);
    1049     int r = nrows(M);
    1050 
    1051     intvec d = intvec(exp1[1..N]) - intvec(exp2[1..N]);
    1052     intvec dm = intvec(V[1..r, exp1[N+1]]) - intvec(V[1..r, exp2[N+1]]);
    1053 
    1054     intvec difference = M * d + dm;
    1055   }
    1056   else
    1057   {
    1058     intvec d = (exp1 - exp2);
    1059     intvec difference = M * d;
    1060   }
    1061 
    1062   if (isFreeRepresented()) // no torsion!?
    1063   {
    1064     return ( difference == 0);
    1065   }
    1066   return ( isTorsionElement( difference ) );
    1067 }
    1068 example
    1069 {
    1070   "EXAMPLE:"; echo=2;
     1753  for(i=1; (i<=r)&&(row<=r)&&(i<=c); i++)
     1754  {
     1755    while((H[row,i]==0)&&(row<=r))
     1756    {
     1757      row++;
     1758      if(row == (r+1)){
     1759        break;
     1760      }
     1761    }
     1762    if(row<=r){
     1763      if(H[row,i]!=0)
     1764      {
     1765        v = H[1..r,i];
     1766        mdeg = mdeg-(mdeg[row]-mdeg[row]%v[row])/v[row]*v;
     1767      }
     1768    }
     1769  }
     1770  return( mdeg == 0 );
     1771
     1772}
     1773example
     1774{
     1775 "EXAMPLE:"; echo=2;
    10711776
    10721777  ring r = 0,(x,y,z),dp;
     
    10751780    1,0,1,
    10761781    0,1,1;
    1077 
    10781782  intmat t[2][1]=
    10791783    -2,
    10801784    1;
     1785
     1786  intmat tt[2][1]=
     1787    1,
     1788    -1;
    10811789
    10821790  setBaseMultigrading(g,t);
     
    10891797  poly f = z2;
    10901798
     1799  intvec v1 = mDeg(a) - mDeg(b);
     1800  v1;
     1801  isZeroElement(v1);
     1802  isZeroElement(v1, tt);
     1803 
     1804  intvec v2 = mDeg(a) - mDeg(c);
     1805  v2;
     1806  isZeroElement(v2);
     1807  isZeroElement(v2, tt);
     1808 
     1809  intvec v3 = mDeg(e) - mDeg(f);
     1810  v3;
     1811  isZeroElement(v3);
     1812  isZeroElement(v3, tt);
     1813 
     1814  intvec v4 = mDeg(c) - mDeg(d);
     1815  v4;
     1816  isZeroElement(v4);
     1817  isZeroElement(v4, tt);
     1818}
     1819
     1820
     1821/******************************************************/
     1822proc defineHomogeneous(poly f, list #)
     1823"USAGE: defineHomogeneous(f[, G]); polynomial f, integer matrix G
     1824PURPOSE: Yields a matrix which has to be appended to the grading group matrix to make the
     1825polynomial f homogeneous in the grading by grad.
     1826EXAMPLE: example defineHomogeneous; shows an example
     1827"
     1828{
     1829  if( size(#) > 0 )
     1830  {
     1831    if( typeof(#[1]) == "intmat" )
     1832    {
     1833      intmat grad = #[1];
     1834    }
     1835  }
     1836
     1837  if( !defined(grad) )
     1838  {
     1839    intmat grad = getVariableWeights();
     1840  }
     1841
     1842  intmat newgg[nrows(grad)][size(f)-1];
     1843  int i,j;
     1844  intvec l = grad*leadexp(f);
     1845  intvec v;
     1846  for(i=2; i <= size(f); i++)
     1847  {
     1848    v = grad * leadexp(f[i]) - l;
     1849    for( j=1; j<=size(v); j++)
     1850    {
     1851      newgg[j,i-1] = v[j];
     1852    }
     1853  }
     1854  return(newgg);
     1855}
     1856example
     1857{
     1858  "EXAMPLE:"; echo=2;
     1859
     1860  ring r =0,(x,y,z),dp;
     1861  intmat grad[2][3] =
     1862    1,0,1,
     1863    0,1,1;
     1864
     1865  setBaseMultigrading(grad);
     1866
     1867  poly f = x2y3-z5+x-3zx;
     1868
     1869  intmat M = defineHomogeneous(f);
     1870  M;
     1871  defineHomogeneous(f, grad) == M;
     1872 
     1873  isHomogeneous(f);
     1874  setBaseMultigrading(grad, M);
     1875  isHomogeneous(f);
     1876}
     1877
     1878
     1879proc gradiator(def h)
     1880PURPOSE: coarsens the grading of the basering until the polynom or ideal h becomes homogeneous.
     1881
     1882{
     1883  if(typeof(h)=="poly"){
     1884    intmat W = getVariableWeights();
     1885    intmat L = getLattice();
     1886    intmat toadd = defineHomogeneous(h);
     1887    //h;
     1888    //toadd;
     1889    if(ncols(toadd) == 0)
     1890    {
     1891      return(1==1);
     1892    }
     1893    int rr = nrows(W);
     1894    intmat newL[rr][ncols(L)+ncols(toadd)];
     1895    newL[1..rr,1..ncols(L)] = L[1..rr,1..ncols(L)];
     1896    newL[1..rr,(ncols(L)+1)..(ncols(L)+ncols(toadd))] = toadd[1..rr,1..ncols(toadd)];
     1897    setBaseMultigrading(W,newL);
     1898    return(1==1);
     1899  }
     1900  if(typeof(h)=="ideal"){
     1901    int i;
     1902    def s = (1==1);
     1903    for(i=1;i<=size(h);i++){
     1904      s = s && gradiator(h[i]);
     1905    }
     1906    return(s);
     1907  }
     1908  return(1==0);
     1909}
     1910example
     1911{
     1912"EXAMPLE:"; echo=2;
     1913ring r = 0,(x,y,z),dp;
     1914intmat g[2][3] = 1,0,1,0,1,1;
     1915intmat l[2][1] = 3,0;
     1916
     1917setBaseMultigrading(g,l);
     1918
     1919getLattice();
     1920
     1921ideal i = -y5+x4,
     1922          y6+xz,
     1923          x2y;
     1924gradiator(i);
     1925getLattice();
     1926isHomogeneous(i);
     1927}
     1928
     1929
     1930proc pushForward(map f)
     1931"USAGE: pushForward(f);
     1932PURPOSE: Computes the finest grading of the image ring which makes the map f
     1933a map of graded rings. The group map between the two grading groups is given
     1934by transpose( (Id, 0) ). Pay attention that the group spanned by the columns of
     1935the grading group matrix may not be a subgroup of the grading group. Still all columns
     1936are needed to find the correct image of the preimage gradings.
     1937EXAMPLE: example pushForward; shows an example
     1938"
     1939{
     1940
     1941  int k,i,j;
     1942//  f;
     1943
     1944//  listvar();
     1945  def pre = preimage(f);
     1946 
     1947//  "pre: ";  pre;
     1948
     1949  intmat oldgrad=getVariableWeights(pre);
     1950  intmat oldtor=getLattice(pre);
     1951
     1952  int n=nvars(pre);
     1953  int np=nvars(basering);
     1954  int p=nrows(oldgrad);
     1955  int pp=p+np;
     1956
     1957  intmat newgrad[pp][np];
     1958
     1959  for(i=1;i<=np;i++){ newgrad[p+i,i]=1;}
     1960
     1961  //newgrad;
     1962
     1963
     1964
     1965  list newtor;
     1966  intmat toadd;
     1967  int columns=0;
     1968
     1969  intmat toadd1[pp][n];
     1970  intvec v;
     1971  poly im;
     1972
     1973  for(i=1;i<=p;i++){
     1974    for(j=1;j<=n;j++){ toadd1[i,j]=oldgrad[i,j];}
     1975  }
     1976
     1977  for(i=1;i<=n;i++){
     1978    im=f[i];
     1979    //im;
     1980    toadd = defineHomogeneous(im, newgrad);
     1981    newtor=insert(newtor,toadd);
     1982    columns=columns+ncols(toadd);
     1983
     1984    v=leadexp(f[i]);
     1985    for(j=p+1;j<=p+np;j++){ toadd1[j,i]=-v[j-p];}
     1986  }
     1987
     1988  newtor=insert(newtor,toadd1);
     1989  columns=columns+ncols(toadd1);
     1990
     1991
     1992  if(typeof(basering)=="qring"){
     1993    //"Entering qring";
     1994    ideal a=ideal(basering);
     1995    for(i=1;i<=size(a);i++){
     1996      toadd = defineHomogeneous(a[i], newgrad);
     1997      //toadd;
     1998      columns=columns+ncols(toadd);
     1999      newtor=insert(newtor,toadd);
     2000    }
     2001  }
     2002
     2003  //newtor;
     2004  intmat imofoldtor[pp][ncols(oldtor)];
     2005  for(i=1; i<=nrows(oldtor);i++){
     2006    for(j=1; j<=ncols(oldtor); j++){
     2007      imofoldtor[i,j]=oldtor[i,j];
     2008    }
     2009  }
     2010
     2011  columns=columns+ncols(oldtor);
     2012  newtor=insert(newtor, imofoldtor);
     2013
     2014  intmat gragr[pp][columns];
     2015  columns=0;
     2016  for(k=1;k<=size(newtor);k++){
     2017    for(i=1;i<=pp;i++){
     2018      for(j=1;j<=ncols(newtor[k]);j++){gragr[i,j+columns]=newtor[k][i,j];}
     2019    }
     2020    columns=columns+ncols(newtor[k]);
     2021  }
     2022
     2023  gragr=hermiteNormalForm(gragr);
     2024  intmat result[pp][pp];
     2025  for(i=1;i<=pp;i++){
     2026    for(j=1;j<=pp;j++){result[i,j]=gragr[i,j];}
     2027  }
     2028
     2029  setBaseMultigrading(newgrad, result);
     2030
     2031}
     2032example
     2033{
     2034  "EXAMPLE:"; echo=2;
     2035
     2036  ring r = 0,(x,y,z),dp;
     2037 
     2038 
     2039
     2040  // Setting degrees for preimage ring.;
     2041  intmat grad[3][3] =
     2042    1,0,0,
     2043    0,1,0,
     2044    0,0,1;
     2045
     2046  setBaseMultigrading(grad);
     2047 
     2048  // grading on r:
     2049  getVariableWeights();
     2050  getLattice();
     2051
     2052  // only for the purpose of this example
     2053  if( voice > 1 ){ /*keepring(r);*/ export(r); }
     2054
     2055  ring R = 0,(a,b),dp;
     2056  ideal i = a2-b2+a6-b5+ab3,a7b+b15-ab6+a6b6;
     2057
     2058  // The quotient ring by this ideal will become our image ring.;
     2059  qring Q = std(i);
     2060
     2061  listvar();
     2062 
     2063  map f = r,-a2b6+b5+a3b+a2+ab,-a2b7-3a2b5+b4+a,a6-b6-b3+a2; f;
     2064
     2065 
     2066  // TODO: Unfortunately this is not a very spectacular example...:
     2067  // Pushing forward f:
     2068  pushForward(f);
     2069
     2070  // due to pushForward we have got new grading on Q
     2071  getVariableWeights();
     2072  getLattice();
     2073 
     2074
     2075  // only for the purpose of this example
     2076  if( voice > 1 ){ kill r; }
     2077
     2078}
     2079
     2080
     2081/******************************************************/
     2082proc equalMDeg(intvec exp1, intvec exp2, list #)
     2083"USAGE: equalMDeg(exp1, exp2[, V]); intvec exp1, exp2, intmat V
     2084PURPOSE: Tests if the exponent vectors of two monomials (given by exp1 and exp2)
     2085represent the same multidegree.
     2086NOTE: the integer matrix V encodes multidegrees of module components,
     2087if module component is present in exp1 and exp2
     2088EXAMPLE: example equalMDeg; shows an example
     2089"
     2090{
     2091  if( size(exp1) != size(exp2) )
     2092  {
     2093    ERROR("Sorry: we cannot compare exponents comming from a polynomial and a vector yet!");
     2094  }
     2095
     2096  if( exp1 == exp2)
     2097  {
     2098    return (1==1);
     2099  }
     2100
     2101
     2102
     2103  intmat M = getVariableWeights();
     2104
     2105  if( nrows(exp1) > ncols(M) ) // vectors => last exponent is the module component!
     2106  {
     2107    if( (size(#) == 0) or (typeof(#[1])!="intmat") )
     2108    {
     2109      ERROR("Sorry: wrong or missing module-unit-weights-matrix V!");
     2110    }
     2111    intmat V = #[1];
     2112
     2113    // typeof(V); print(V);
     2114
     2115    int N = ncols(M);
     2116    int r = nrows(M);
     2117
     2118    intvec d = intvec(exp1[1..N]) - intvec(exp2[1..N]);
     2119    intvec dm = intvec(V[1..r, exp1[N+1]]) - intvec(V[1..r, exp2[N+1]]);
     2120
     2121    intvec difference = M * d + dm;
     2122  }
     2123  else
     2124  {
     2125    intvec d = (exp1 - exp2);
     2126    intvec difference = M * d;
     2127  }
     2128
     2129  if (isFreeRepresented()) // no grading group!?
     2130  {
     2131    return ( difference == 0);
     2132  }
     2133  return ( isZeroElement( difference ) );
     2134}
     2135example
     2136{
     2137  "EXAMPLE:"; echo=2;printlevel=3;
     2138
     2139  ring r = 0,(x,y,z),dp;
     2140
     2141  intmat g[2][3]=
     2142    1,0,1,
     2143    0,1,1;
     2144
     2145  intmat t[2][1]=
     2146    -2,
     2147    1;
     2148
     2149  setBaseMultigrading(g,t);
     2150
     2151  poly a = x10yz;
     2152  poly b = x8y2z;
     2153  poly c = x4z2;
     2154  poly d = y5;
     2155  poly e = x2y2;
     2156  poly f = z2;
     2157
    10912158
    10922159  equalMDeg(leadexp(a), leadexp(b));
     
    11162183/******************************************************/
    11172184static proc isFreeRepresented()
    1118 "check whether the base muligrading is torsion-free (it is zero).
    1119 "
    1120 {
    1121   intmat T = getTorsion();
     2185"check whether the base muligrading is free (it is zero).
     2186"
     2187{
     2188  intmat T = getLattice();
    11222189
    11232190  intmat Z[nrows(T)][ncols(T)];
    11242191
    1125   return (T == Z); // no torsion!
     2192  return (T == Z); // no grading group!
    11262193}
    11272194
    11282195
    11292196/******************************************************/
    1130 proc isHomogenous(def a, list #)
    1131 "USAGE: isHomogenous(a[, f]); a polynomial/vector/ideal/module
    1132 RETURN: boolean, TRUE if a is (multi)homogenous, and FALSE otherwise
    1133 EXAMPLE: example isHomogenous; shows an example
     2197proc isHomogeneous(def a, list #)
     2198"USAGE: isHomogeneous(a[, f]); a polynomial/vector/ideal/module
     2199RETURN: boolean, TRUE if a is (multi)homogeneous, and FALSE otherwise
     2200EXAMPLE: example isHomogeneous; shows an example
    11342201"
    11352202{
     
    11372204  {
    11382205    return ( size(mDegPartition(a)) <= 1 )
    1139   }
    1140 
    1141   if( typeof(a) == "vector" or typeof(a) == "module" )
    1142   {
    1143     intmat V = getModuleGrading(a);
    11442206  }
    11452207
     
    11532215        for( int i = ncols(a); i > 0; i-- )
    11542216        {
    1155           aa = a[i];
    1156 
    1157           if(defined(V))
    1158           {
    1159             aa = setModuleGrading(aa, V);
    1160           }
    1161 
    1162           if(!isHomogenous(aa))
     2217          aa = getGradedGenerator(a, i);
     2218
     2219          if(!isHomogeneous(aa))
    11632220          {
    11642221            return(0==1);
     
    11712228    def g = groebner(a); // !!!!
    11722229
    1173     def b, aa; int j;
     2230    def b, aa; int j; 
    11742231    for( int i = ncols(a); i > 0; i-- )
    11752232    {
    1176       aa = a[i];
    1177 
    1178       if(defined(V))
    1179       {
    1180         aa = setModuleGrading(aa, V);
    1181       }
     2233      aa = getGradedGenerator(a, i);
    11822234
    11832235      b = mDegPartition(aa);
     
    11912243    }
    11922244    return(1==1);
    1193   }
     2245  } 
    11942246}
    11952247example
     
    12002252
    12012253  //Grading and Torsion matrices:
    1202   intmat M[3][3] =
     2254  intmat M[3][3] = 
    12032255    1,0,0,
    12042256    0,1,0,
     
    12172269  print(mDeg(_));
    12182270
    1219   isHomogenous(f);   // f: is not homogenous
     2271  isHomogeneous(f);   // f: is not homogeneous
    12202272
    12212273  poly g = 1-xy2z3;
    1222   isHomogenous(g); // g: is homogenous
     2274  isHomogeneous(g); // g: is homogeneous
    12232275  mDegPartition(g);
    12242276
     
    12262278  /////////////////////////////////////////////////////////
    12272279  // new Torsion matrix:
    1228   intmat T[3][4] =
     2280  intmat T[3][4] = 
    12292281    3,3,3,3,
    12302282    2,1,3,0,
    12312283    1,2,0,3;
    1232 
     2284 
    12332285  setBaseMultigrading(M,T);
    12342286
    12352287  f;
    1236   isHomogenous(f);
     2288  isHomogeneous(f);
    12372289  mDegPartition(f);
    12382290
    1239   // ---------------------
     2291  // --------------------- 
    12402292  g;
    1241   isHomogenous(g);
     2293  isHomogeneous(g);
    12422294  mDegPartition(g);
    12432295
     
    12462298  ring R = 0, (x,y,z), dp;
    12472299
    1248   intmat A[2][3] =
     2300  intmat A[2][3] = 
    12492301    0,0,1,
    12502302    3,2,1;
    1251   intmat T[2][1] =
    1252     -1,
     2303  intmat T[2][1] = 
     2304    -1, 
    12532305     4;
    12542306  setBaseMultigrading(A, T);
    12552307
    1256   isHomogenous(ideal(x2 - y3 -xy +z, x*y-z, x^3 - y^2*z + x^2 -y^3)); // 1
    1257   isHomogenous(ideal(x2 - y3 -xy +z, x*y-z, x^3 - y^2*z + x^2 -y^3), "checkGens");
    1258   isHomogenous(ideal(x+y, x2 - y2)); // 0
     2308  isHomogeneous(ideal(x2 - y3 -xy +z, x*y-z, x^3 - y^2*z + x^2 -y^3)); // 1
     2309  isHomogeneous(ideal(x2 - y3 -xy +z, x*y-z, x^3 - y^2*z + x^2 -y^3), "checkGens");
     2310  isHomogeneous(ideal(x+y, x2 - y2)); // 0
    12592311
    12602312  // Degree partition:
     
    12622314  mDegPartition(x3 -y2z + x2 -y3 + z + 1);
    12632315
    1264 
     2316 
    12652317  module N = gen(1) + (x+y) * gen(2), z*gen(3);
    12662318
    12672319  intmat V[2][3] = 0; // 1, 2, 3,  4, 5, 6; //  column-wise weights of components!!??
    1268 
     2320 
    12692321  vector v1, v2;
    1270 
     2322 
    12712323  v1 = setModuleGrading(N[1], V); v1;
    12722324  mDegPartition(v1);
     
    12782330
    12792331  N = setModuleGrading(N, V);
    1280   isHomogenous(N);
     2332  isHomogeneous(N);
    12812333  print( mDeg(N) );
    12822334
    1283   ///////////////////////////////////////
    1284 
    1285   V =
    1286     1, 2, 3,
     2335  /////////////////////////////////////// 
     2336
     2337  V = 
     2338    1, 2, 3, 
    12872339    4, 5, 6;
    12882340
     
    12962348
    12972349  N = setModuleGrading(N, V);
    1298   isHomogenous(N);
     2350  isHomogeneous(N);
    12992351  print( mDeg(N) );
    13002352
    1301   ///////////////////////////////////////
    1302 
    1303   V =
    1304     0, 0, 0,
     2353  /////////////////////////////////////// 
     2354
     2355  V = 
     2356    0, 0, 0, 
    13052357    4, 1, 0;
    13062358
    13072359  N = gen(1) + x * gen(2), z*gen(3);
    13082360  N = setModuleGrading(N, V); print(N);
    1309   isHomogenous(N);
     2361  isHomogeneous(N);
    13102362  print( mDeg(N) );
    1311   v1 = setModuleGrading(N[1], V); print(v1);
     2363  v1 = getGradedGenerator(N,1); print(v1);
    13122364  mDegPartition(v1);
    13132365  print( mDeg(_) );
    13142366  N = setModuleGrading(N, V); print(N);
    1315   isHomogenous(N);
     2367  isHomogeneous(N);
    13162368  print( mDeg(N) );
    13172369}
     
    13392391  int r = nrows(M);
    13402392
     2393  if( (typeof(A) == "vector") or (typeof(A) == "module") )
     2394  {
     2395    intmat V = getModuleGrading(A);
     2396 
     2397    if( nrows(V) != r )
     2398    {
     2399      ERROR("Sorry wrong mgrad-size of V: " + string(nrows(V)));
     2400    }
     2401  }
     2402                           
    13412403  if( A == 0 )
    13422404  {
     
    13452407  }
    13462408
    1347   if( (typeof(A) == "vector") or (typeof(A) == "module") )
    1348   {
    1349     intmat V = getModuleGrading(A);
    1350 
    1351     if( nrows(V) != r )
    1352     {
    1353       ERROR("Sorry wrong mgrad-size of V: " + string(nrows(V)));
    1354     }
    1355   }
    1356 
    13572409  intvec m; m[r] = 0;
    13582410
     
    13672419    A = A - lead(A);
    13682420    while( size(A) > 0 )
    1369     {
     2421    { 
    13702422      v = leadexp(A); //  v;
    13712423      m = max( m, M * v, r ); // ????
     
    13882440    A = A - lead(A);
    13892441    while( size(A) > 0 )
    1390     {
     2442    { 
    13912443      v = leadexp(A); //  v;
    13922444
     
    14132465      {
    14142466        G[j, i] = d[j];
    1415       }
     2467      }     
    14162468    }
    14172469    return(G);
     
    14252477    for( i = ncols(A); i > 0; i-- )
    14262478    {
    1427       v = setModuleGrading(A[i], V);
    1428 
    1429       // G[1..r, i]
     2479      v = getGradedGenerator(A, i);
     2480
     2481      // G[1..r, i] 
    14302482      d = mDeg(v);
    14312483
     
    14332485      {
    14342486        G[j, i] = d[j];
    1435       }
     2487      }     
    14362488
    14372489    }
     
    14392491    return(G);
    14402492  }
    1441 
     2493 
    14422494}
    14432495example
     
    14532505  print(Ta);
    14542506
    1455   //   attrib(A, "torsion", Ta); // to think about
     2507  //   attrib(A, "gradingGroup", Ta); // to think about
    14562508
    14572509//  "poly:";
     
    14632515
    14642516  setBaseMultigrading(A, Ta);
    1465 
     2517 
    14662518  mDeg( x*x*y );
    1467 
     2519 
    14682520  mDeg( y*y*y*x );
    1469 
     2521 
    14702522  mDeg( x*y + x + 1 );
    14712523
     
    14772529
    14782530//  "ideal:";
    1479 
     2531 
    14802532  ideal I = y*x*x, x*y*y*y;
    14812533  print( mDeg(I) );
     
    14852537
    14862538//  "vectors:";
    1487 
     2539 
    14882540  intmat B[2][2] = 0, 1, 1, 0;
    14892541  print(B);
    1490 
     2542 
    14912543  mDeg( setModuleGrading(y*y*y*gen(2), B ));
    14922544  mDeg( setModuleGrading(x*x*gen(1), B ));
    14932545
    1494 
     2546 
    14952547  vector V = x*gen(1) + y*gen(2);
    14962548  V = setModuleGrading(V, B);
     
    14992551  vector v1 = setModuleGrading([0, 0, 0], B);
    15002552  print( mDeg( v1 ) );
    1501 
     2553 
    15022554  vector v2 = setModuleGrading([0], B);
    15032555  print( mDeg( v2 ) );
    15042556
    15052557//  "module:";
    1506 
     2558 
    15072559  module D = x*gen(1), y*gen(2);
    15082560  D;
    15092561  D = setModuleGrading(D, B);
    15102562  print( mDeg( D ) );
    1511 
     2563 
    15122564
    15132565  module DD = [0, 0],[0, 0, 0];
     
    15182570  DDD = setModuleGrading(DDD, B);
    15192571  print( mDeg( DDD ) );
    1520 }
     2572
     2573};
     2574
     2575
     2576
     2577
    15212578
    15222579/******************************************************/
     
    15262583EXAMPLE: example mDegPartition; shows an example
    15272584"
    1528 {
    1529   if( typeof(p) == "poly" )
    1530   {
    1531     ideal I;
     2585{ // TODO: What about an ideal or module???
     2586
     2587  if( typeof(p) == "poly" )
     2588  {
     2589    ideal I;     
    15322590    poly mp, t, tt;
    15332591  }
     
    15362594    if(  typeof(p) == "vector" )
    15372595    {
    1538       module I;
     2596      module I;     
    15392597      vector mp, t, tt;
    15402598    }
     
    15472605  if( typeof(p) == "vector" )
    15482606  {
    1549     intmat V = getModuleGrading(p);
     2607    intmat V = getModuleGrading(p); 
    15502608  }
    15512609  else
     
    15542612  }
    15552613
    1556   if( size(p) > 1)
     2614  if( size(p) > 1) 
    15572615  {
    15582616    intvec m;
     
    15612619    {
    15622620      m = leadexp(p);
    1563       mp = lead(p);
     2621      mp = lead(p); 
    15642622      p = p - lead(p);
    15652623      tt = p; t = 0;
    15662624
    15672625      while( size(tt) > 0 )
    1568       {
     2626      { 
    15692627        // TODO: we do not cache matrices (M,T,H,V), which remain the same :(
    15702628        // TODO: we need some low-level procedure with all these arguments...!
    1571         if( equalMDeg( leadexp(tt), m, V  ) )
     2629        if( equalMDeg( leadexp(tt), m, V  ) ) 
    15722630        {
    15732631          mp = mp + lead(tt); // "mp", mp;
     
    16162674
    16172675  mDegPartition(f);
    1618 
     2676 
    16192677  vector v = xy*gen(1)-x3y2*gen(2)+x4y*gen(3);
    16202678  intmat B[2][3]=1,-1,-2,0,0,1;
    16212679  v = setModuleGrading(v,B);
    16222680  getModuleGrading(v);
    1623 
     2681 
    16242682  mDegPartition(v, B);
    16252683}
     
    16312689{
    16322690  intmat A[n][n];
    1633 
     2691 
    16342692  for( int i = n; i > 0; i-- )
    16352693  {
     
    16462704"
    16472705USAGE: finestMDeg(r); ring r
    1648 RETURN: ring, r endowed with the finest multigrading
     2706RETURN: ring, r endowed with the finest multigrading 
    16492707TODO: not yet...
    16502708"
     
    16752733
    16762734
    1677   if( n > 0)
    1678   {
    1679 
    1680     intmat L[N][n];
     2735  if( n > 0) 
     2736  {
     2737
     2738    intmat L[N][n]; 
    16812739    //  list L;
    16822740    int j = n;
     
    16862744      p = I[i];
    16872745
    1688       if( size(p) > 1 )
     2746      if( size(p) > 1 ) 
    16892747      {
    16902748        intvec m0 = leadexp(p);
     
    17012759
    17022760    print(L);
    1703     setBaseMultigrading(A, L);
    1704   }
     2761    setBaseMultigrading(A, L);     
     2762  } 
    17052763  else
    17062764  {
     
    17392797
    17402798  if( size(#) > 0 and typeof(#[1]) == "intmat" )
    1741   {
     2799  { 
    17422800    attrib(F, "P", #[1]);
    17432801  }
     
    17972855  // check:
    17982856  print(L*M);
    1799 }
     2857};
     2858
     2859
    18002860
    18012861/******************************************************/
     
    18192879"
    18202880{
     2881   if( system("sh","which hilbert 2> /dev/null 1> /dev/null") != 0 )
     2882   {
     2883     ERROR("Sorry: cannot find 'hilbert' command from 4ti2. Please install 4ti2!");
     2884   }
     2885   
    18212886//--------------------------------------------------------------------------
    18222887// Initialization and Sanity Checks
     
    18702935// using standard unix commands
    18712936//----------------------------------------------------------------------
     2937
     2938
    18722939   j=system("sh","hilbert -q -n sing4ti2"); ////////// be quiet + no loggin!!!
    18732940
    1874    j=system("sh", "awk \'BEGIN{ORS=\",\";}{print $0;}\' sing4ti2.hil " +
    1875                 "| sed s/[\\\ \\\t\\\v\\\f]/,/g " +
    1876                 "| sed s/,+/,/g|sed s/,,/,/g " +
    1877                 "| sed s/,,/,/g " +
     2941   j=system("sh", "awk \'BEGIN{ORS=\",\";}{print $0;}\' sing4ti2.hil " + 
     2942                "| sed s/[\\\ \\\t\\\v\\\f]/,/g " + 
     2943                "| sed s/,+/,/g|sed s/,,/,/g " + 
     2944                "| sed s/,,/,/g " + 
    18782945                "> sing4ti2.converted" );
    1879    if( defined(keepfiles) <= 0)
    1880    {
    1881       j=system("sh",("rm -f sing4ti2.hil sing4ti2."+fileending));
    1882    }
     2946
     2947
    18832948//----------------------------------------------------------------------
    18842949// reading output of 4ti2
     
    18892954
    18902955   string s = read(ausg);
     2956
     2957   if( defined(keepfiles) <= 0)
     2958   {
     2959      j=system("sh",("rm -f sing4ti2.hil sing4ti2.converted sing4ti2."+fileending));
     2960   }
     2961
    18912962   string ergstr = "intvec erglist = " + s + "0;";
    18922963   execute(ergstr);
    1893 
     2964 
    18942965   //   print(erglist);
    1895 
     2966 
    18962967   int Rnc = erglist[1];
    18972968   int Rnr = erglist[2];
    1898 
     2969   
    18992970   intmat R[Rnr][Rnc];
    19002971
     
    19102981     }
    19112982   }
     2983
     2984
    19122985
    19132986   return (R);
     
    19313004      1, 1, 0,  0, -1,  0,-1, 0, 0;
    19323005   hilbert4ti2intmat(M);
    1933    hermite(M);
     3006   hermiteNormalForm(M);
    19343007}
    19353008
     
    19623035
    19633036/******************************************************/
    1964 proc mDegBasis(intvec d)
     3037proc mDegBasis(intvec d) 
    19653038"
    19663039USAGE: multidegree d
     
    19913064  }
    19923065
    1993   intmat T = getTorsion(R);
     3066  intmat T = getLattice(R);
    19943067
    19953068  if( isFreeRepresented() )
     
    20063079
    20073080    intmat AA[nr][nc + 2 * n];
    2008     AA[1..nr, 1.. nc] = A[1..nr, 1.. nc];
    2009     AA[1..nr, nc + (1.. n)] = T[1..nr, 1.. n];
    2010     AA[1..nr, nc + n + (1.. n)] = -T[1..nr, 1.. n];
     3081    AA[1..nr, 1.. nc] = A[1..nr, 1.. nc]; 
     3082    AA[1..nr, nc + (1.. n)] = T[1..nr, 1.. n]; 
     3083    AA[1..nr, nc + n + (1.. n)] = -T[1..nr, 1.. n]; 
    20113084
    20123085
    20133086    //      print ( AA );
    20143087
    2015     intmat K = leftKernelZ(( AA ) ); //
     3088    intmat K = leftKernelZ(( AA ) ); // 
    20163089
    20173090    //      print(K);
     
    20223095    //      "!";
    20233096
    2024     intmat B = hilbert4ti2intmat(transpose(KK), 1);
     3097    intmat B = hilbert4ti2intmat(transpose(KK), 1); 
    20253098
    20263099    //      "!";      print(B);
     
    20333106
    20343107
    2035   int i;
     3108  int i; 
    20363109  int nnr = nrows(B);
    20373110  int nnc = ncols(B);
     
    20883161
    20893162  intmat g1[2][2]=1,0,0,1;
    2090   intmat t[2][1]=2,0;
     3163  intmat l[2][1]=2,0;
    20913164  intmat g2[2][2]=1,1,1,1;
    20923165  intvec v1=4,0;
    20933166  intvec v2=4,4;
    2094 
     3167 
    20953168  intmat g3[1][2]=1,1;
    20963169  setBaseMultigrading(g3);
     
    20983171  v3;
    20993172  mDegBasis(v3);
    2100 
    2101   setBaseMultigrading(g1,t);
     3173 
     3174  setBaseMultigrading(g1,l);
    21023175  mDegBasis(v1);
    21033176  setBaseMultigrading(g2);
    21043177  mDegBasis(v2);
    2105 
     3178 
    21063179  intmat M[2][2] = 1, -1, -1, 1;
    21073180  intvec d = -2, 2;
     
    21173190  intmat M[2][3] = 1, -2, 1,     1, 1, 0;
    21183191
    2119   intmat T[2][1] = 0, 2;
     3192  intmat L[2][1] = 0, 2;
    21203193
    21213194  intvec d = 4, 1;
    21223195
    2123   setBaseMultigrading(M, T);
     3196  setBaseMultigrading(M, L);
    21243197
    21253198  mDegBasis(d);
     
    21623235
    21633236proc mDegSyzygy(def I)
    2164 "USAGE: mDegSyzygy(I); I is a poly/vector/ideal/module
     3237"USAGE: mDegSyzygy(I); I is a ideal or a module
    21653238PURPOSE: computes the multigraded syzygy module of I
    21663239RETURNS: module, the syzygy module of I
     
    21693242"
    21703243{
    2171   if( isHomogenous(I, "checkGens") == 0)
    2172   {
    2173     ERROR ("Sorry: inhomogenous input!");
    2174   }
     3244  if( isHomogeneous(I, "checkGens") == 0)
     3245  { 
     3246    ERROR ("Sorry: inhomogeneous input!");
     3247  } 
    21753248  module S = syz(I);
    21763249  S = setModuleGrading(S, mDeg(I));
     
    21803253{
    21813254  "EXAMPLE:"; echo=2;
    2182 
     3255 
    21833256
    21843257  ring r = 0,(x,y,z,w),dp;
     
    21883261  setBaseMultigrading(MM);
    21893262  module M = ideal(  xw-yz, x2z-y3, xz2-y2w, yw2-z3);
    2190 
    2191 
     3263 
     3264 
    21923265  intmat v[2][nrows(M)]=
    21933266    1,
    21943267    0;
    2195 
     3268 
    21963269  M = setModuleGrading(M, v);
    21973270
    2198   isHomogenous(M);
     3271  isHomogeneous(M);
    21993272  "Multidegrees: "; print(mDeg(M));
    22003273  // Let's compute syzygies!
     
    22033276  "Multidegrees: "; print(mDeg(S));
    22043277
    2205   isHomogenous(S);
     3278  isHomogeneous(S);
     3279}
     3280
     3281
     3282
     3283proc mDegModulo(def I, def J)
     3284"USAGE: mDegModulo(I); I, J are ideals or modules
     3285PURPOSE: computes the multigraded 'modulo' module of I and J
     3286RETURNS: module, see 'modulo' command
     3287NOTE: I and J should have the same multigrading, and their
     3288generators must be multigraded homogeneous
     3289EXAMPLE: example mDegModulo; shows an example
     3290"
     3291{
     3292  if( (isHomogeneous(I, "checkGens") == 0) or (isHomogeneous(J, "checkGens") == 0) )
     3293  {
     3294    ERROR ("Sorry: inhomogeneous input!");
     3295  }
     3296  module K = modulo(I, J);
     3297  K = setModuleGrading(K, mDeg(I));
     3298  return (K);
     3299}
     3300example
     3301{
     3302  "EXAMPLE:"; echo=2;
     3303
     3304  ring r = 0,(x,y,z),dp;
     3305  intmat MM[2][3]=
     3306    -1,1,1,
     3307     0,1,3;
     3308  setBaseMultigrading(MM);
     3309
     3310  ideal h1 = x, y, z;
     3311  ideal h2 = x;
     3312
     3313  "Multidegrees: "; print(mDeg(h1));
     3314   
     3315  // Let's compute modulo(h1, h2):
     3316  def K = mDegModulo(h1, h2); K;
     3317
     3318  "Module Units Multigrading: "; print( getModuleGrading(K) );
     3319  "Multidegrees: "; print(mDeg(K));
     3320
     3321  isHomogeneous(K);
    22063322}
    22073323
     
    22103326"USAGE: mDegGroebner(I); I is a poly/vector/ideal/module
    22113327PURPOSE: computes the multigraded standard/groebner basis of I
    2212 NOTE: I must be multigraded homogeneous
     3328NOTE: I must be multigraded homogeneous 
    22133329RETURNS: ideal/module, the computed basis
    22143330EXAMPLE: example mDegGroebner; shows an example
    22153331"
    22163332{
    2217   if( isHomogenous(I) == 0)
    2218   {
    2219     ERROR ("Sorry: inhomogenous input!");
    2220   }
     3333  if( isHomogeneous(I) == 0)
     3334  { 
     3335    ERROR ("Sorry: inhomogeneous input!");
     3336  } 
    22213337
    22223338  def S = groebner(I);
    2223 
     3339 
    22243340  if( typeof(I) == "module" or typeof(I) == "vector" )
    22253341  {
    2226     S = setModuleGrading(S, getModuleGrading(I));
     3342    S = setModuleGrading(S, getModuleGrading(I));     
    22273343  }
    22283344
     
    22413357  setBaseMultigrading(MM);
    22423358
    2243 
     3359 
    22443360  module M = ideal(  xw-yz, x2z-y3, xz2-y2w, yw2-z3);
    2245 
    2246 
     3361 
     3362 
    22473363  intmat v[2][nrows(M)]=
    22483364    1,
    22493365    0;
    2250 
     3366 
    22513367  M = setModuleGrading(M, v);
    22523368
     
    22583374  "Multidegrees: "; print(mDeg(M));
    22593375
    2260   isHomogenous(M);
     3376  isHomogeneous(M);
    22613377
    22623378  /////////////////////////////////////////////////////////////////////////////
     
    22663382  "Multidegrees: "; print(mDeg(S));
    22673383
    2268   isHomogenous(S);
     3384  isHomogeneous(S);
    22693385
    22703386  /////////////////////////////////////////////////////////////////////////////
     
    22743390  "Multidegrees: "; print(mDeg(S));
    22753391
    2276   isHomogenous(S);
     3392  isHomogeneous(S);
    22773393}
    22783394
     
    22813397proc mDegResolution(def I, int ll, list #)
    22823398"USAGE: mDegResolution(I,l,[f]); I is poly/vector/ideal/module; l,f are integers
    2283 PURPOSE: computes the multigraded resolution of I of the length l,
    2284 or the whole resolution if l is zero. Returns minimal resolution if an optional
     3399PURPOSE: computes the multigraded resolution of I of the length l, 
     3400or the whole resolution if l is zero. Returns minimal resolution if an optional 
    22853401argument 1 is supplied
    22863402NOTE: input must have multigraded-homogeneous generators.
    2287 The returned list is truncated beginning with the first zero differential.
     3403The returned list is truncated beginning with the first zero differential.   
    22883404RETURNS: list, the computed resolution
    22893405EXAMPLE: example mDegResolution; shows an example
    22903406"
    22913407{
    2292   if( isHomogenous(I, "checkGens") == 0)
    2293   {
    2294     ERROR ("Sorry: inhomogenous input!");
    2295   }
     3408  if( isHomogeneous(I, "checkGens") == 0)
     3409  { 
     3410    ERROR ("Sorry: inhomogeneous input!");
     3411  } 
    22963412
    22973413  def R = res(I, ll, #); list L = R; int l = size(L);
    2298 
     3414  def V = getModuleGrading(I);
    22993415  if( (typeof(I) == "module") or (typeof(I) == "vector") )
    23003416  {
    2301     L[1] = setModuleGrading(L[1], getModuleGrading(I));
    2302   }
    2303 
    2304   int i;
     3417    L[1] = setModuleGrading(L[1], V);
     3418  }
     3419
     3420  int i; 
    23053421  for( i = 2; i <= l; i++ )
    23063422  {
     
    23133429    }
    23143430  }
    2315 
     3431 
    23163432  return (L);
    23173433
    2318 
     3434 
    23193435}
    23203436example
    23213437{
    23223438  "EXAMPLE:"; echo=2;
    2323 
     3439 
    23243440  ring r = 0,(x,y,z,w),dp;
    23253441
     
    23303446  setBaseMultigrading(M);
    23313447
    2332 
     3448 
    23333449  module m= ideal(  xw-yz, x2z-y3, xz2-y2w, yw2-z3);
    2334 
    2335   isHomogenous(ideal(  xw-yz, x2z-y3, xz2-y2w, yw2-z3), "checkGens");
    2336 
     3450 
     3451  isHomogeneous(ideal(  xw-yz, x2z-y3, xz2-y2w, yw2-z3), "checkGens");
     3452 
    23373453  ideal A = xw-yz, x2z-y3, xz2-y2w, yw2-z3;
    23383454
    23393455  int j;
    2340 
     3456 
    23413457  for(j=1; j<=ncols(A); j++)
    23423458  {
    23433459    mDegPartition(A[j]);
    23443460  }
    2345 
     3461 
    23463462  intmat v[2][1]=
    23473463    1,
    23483464    0;
    2349 
     3465 
    23503466  m = setModuleGrading(m, v);
    23513467
     
    23743490
    23753491  /////////////////////////////////////////////////////////////////////////////
    2376 
     3492 
    23773493  L = mDegResolution(maxideal(1), 0, 1);
    23783494
     
    23843500    "Multigrading: "; print(mDeg(L[j]));
    23853501  }
    2386 
     3502 
    23873503  kill v;
    2388 
     3504 
    23893505
    23903506  def h = hilbertSeries(m);
     
    23933509  numerator1;
    23943510  factorize(numerator1);
    2395 
     3511 
    23963512  denominator1;
    23973513  factorize(denominator1);
     
    24073523proc hilbertSeries(def I)
    24083524"USAGE: hilbertSeries(I); I is poly/vector/ideal/module
    2409 PURPOSE: computes the multigraded Hilbert Series of M
    2410 NOTE: input must have multigraded-homogeneous generators.
     3525PURPOSE: computes the multigraded Hilbert Series of M 
     3526NOTE: input must have multigraded-homogeneous generators. 
    24113527Multigrading should be positive.
    2412 RETURNS: a ring in variables t_(i), s_(i), with polynomials
    2413 numerator1 and denominator1 and muturally prime numerator2
     3528RETURNS: a ring in variables t_(i), s_(i), with polynomials 
     3529numerator1 and denominator1 and muturally prime numerator2 
    24143530and denominator2, quotients of which give the series.
    24153531EXAMPLE: example hilbertSeries; shows an example
    24163532"
    24173533{
    2418 
     3534   
    24193535  if( !isFreeRepresented() )
    24203536  {
    2421     ERROR("SORRY: ONLY TORSION-FREE CASE (POSITIVE GRADING)");
    2422   }
    2423 
     3537    "Things might happen, since we are not  free.";
     3538    //ERROR("SORRY: ONLY TORSION-FREE CASE (POSITIVE GRADING)");
     3539  }
     3540   
    24243541  int i, j, k, v;
    24253542
    24263543  intmat M = getVariableWeights();
    2427 
     3544 
    24283545  int cc = ncols(M);
    24293546  int n = nrows(M);
     
    24373554
    24383555  int l = size(RES);
    2439 
     3556 
    24403557  list L; L[l + 1] = 0;
    24413558
     
    24443561    intmat zeros[n][1];
    24453562    L[1] = zeros;
    2446   }
     3563  } 
    24473564  else
    24483565  {
     
    24543571    L[j + 1] = mDeg(RES[j]);
    24553572  }
    2456 
     3573 
    24573574  l++;
    24583575
    24593576  ring R = 0,(t_(1..n),s_(1..n)),dp;
    2460 
    2461   ideal units;
     3577 
     3578  ideal units; 
    24623579  for( i=n; i>=1; i--)
    24633580  {
    24643581    units[i] = (var(i) * var(n + i) - 1);
    24653582  }
    2466 
     3583 
    24673584  qring Q = std(units);
    2468 
     3585 
    24693586  // TODO: should not it be a quotient ring depending on Torsion???
    24703587  // I am not sure about what to do in the torsion case, but since
     
    24753592  poly monom, summand, numerator;
    24763593  poly denominator = 1;
    2477 
     3594 
    24783595  for( i = 1; i <= cc; i++)
    24793596  {
    24803597    monom = 1;
    2481 
    24823598    for( k = 1; k <= n; k++)
    24833599    {
     
    24873603      {
    24883604        monom = monom * (var(k)^(v));
    2489       }
     3605      } 
    24903606      else
    24913607      {
     
    24933609      }
    24943610    }
    2495 
     3611   
    24963612    if( monom == 1)
    24973613    {
     
    25013617    denominator = denominator * (1 - monom);
    25023618  }
    2503 
    2504   for( j = 1; j<= l; j++)
     3619 
     3620  for( j = 1; j<= l; j++) 
    25053621  {
    25063622    summand = 0;
     
    25163632        {
    25173633          monom = monom * (var(k)^v);
    2518         }
     3634        } 
    25193635        else
    25203636        {
     
    25263642    numerator = numerator - (-1)^j * summand;
    25273643  }
    2528 
     3644 
    25293645  if( denominator == 0 )
    25303646  {
    25313647    ERROR("Multigrading not positive.");
    2532   }
    2533 
     3648  } 
     3649 
    25343650  poly denominator1 = denominator;
    25353651  poly numerator1 = numerator;
     
    25653681  "The s_(i)-variables are defined to be the inverse of the t_(i)-variables.";
    25663682  " ------------ ";
    2567 
     3683 
    25683684  return(Q);
    25693685}
     
    25713687{
    25723688  "EXAMPLE:"; echo=2;
    2573 
     3689 
    25743690  ring r = 0,(x,y,z,w),dp;
    25753691  intmat g[2][4]=
     
    25773693    0,1,3,4;
    25783694  setBaseMultigrading(g);
    2579 
     3695 
    25803696  module M = ideal(xw-yz, x2z-y3, xz2-y2w, yw2-z3);
    25813697  intmat V[2][1]=
     
    25893705  factorize(numerator2);
    25903706  factorize(denominator2);
    2591 
     3707 
    25923708  kill g, h; setring r;
    25933709
     
    25953711    1,2,3,4,
    25963712    0,0,5,8;
    2597 
     3713 
    25983714  setBaseMultigrading(g);
    2599 
     3715 
    26003716  ideal I = x^2, y, z^3;
    26013717  I = std(I);
     
    26123728  mDeg(I);
    26133729  def h = hilbertSeries(I); setring h;
    2614 
     3730 
    26153731  factorize(numerator2);
    26163732  factorize(denominator2);
     
    26193735  ////////////////////////////////////////////////
    26203736  ring R = 0,(x,y,z),dp;
    2621   intmat W[2][3] =
     3737  intmat W[2][3] = 
    26223738     1,1, 1,
    26233739     0,0,-1;
    26243740  setBaseMultigrading(W);
    26253741  ideal I = x3y,yz2,y2z,z4;
    2626 
     3742 
    26273743  def h = hilbertSeries(I); setring h;
    2628 
     3744 
    26293745  factorize(numerator2);
    26303746  factorize(denominator2);
     
    26333749  ////////////////////////////////////////////////
    26343750  ring R = 0,(x,y,z,a,b,c),dp;
    2635   intmat W[2][6] =
     3751  intmat W[2][6] = 
    26363752     1,1, 1,1,1,1,
    26373753     0,0,-1,0,0,0;
    26383754  setBaseMultigrading(W);
    26393755  ideal I = x3y,yz2,y2z,z4;
    2640 
     3756 
    26413757  def h = hilbertSeries(I); setring h;
    2642 
     3758 
    26433759  factorize(numerator2);
    26443760  factorize(denominator2);
    2645 
     3761 
    26463762  kill R, W, h;
    26473763  ////////////////////////////////////////////////
    26483764  // This is example 5.3.9. from Robbianos book.
    2649 
     3765 
    26503766  ring R = 0,(x,y,z,w),dp;
    2651   intmat W[1][4] =
     3767  intmat W[1][4] = 
    26523768     1,1, 1,1;
    26533769  setBaseMultigrading(W);
     
    26553771
    26563772  hilb(std(I));
    2657 
     3773 
    26583774  def h = hilbertSeries(I); setring h;
    2659 
     3775 
    26603776  numerator1;
    26613777  denominator1;
     
    26633779  factorize(numerator2);
    26643780  factorize(denominator2);
    2665 
     3781 
    26663782
    26673783  kill h;
     
    26723788
    26733789  hilb(std(I2));
    2674 
     3790 
    26753791  def h = hilbertSeries(I2); setring h;
    26763792
     
    26823798  ////////////////////////////////////////////////
    26833799  setring R;
    2684 
     3800 
    26853801  W = 2,2,2,2;
    2686 
     3802 
    26873803  setBaseMultigrading(W);
    26883804
     
    26943810
    26953811  kill w;
    2696 
     3812 
    26973813
    26983814  def h = hilbertSeries(I2); setring h;
    26993815
    2700 
     3816 
    27013817  numerator1; denominator1;
    27023818  kill h;
    27033819
    2704 
     3820 
    27053821  kill R, W;
    27063822
     
    27163832
    27173833  hilb(std(I));
    2718 
     3834 
    27193835  def h = hilbertSeries(I); setring h;
    27203836
     
    27323848
    27333849  numerator1; denominator1;
     3850 
     3851  kill h; 
     3852  ////////////////////////////////////////////////
     3853  setring R;
     3854
     3855  I = x^5; I;
     3856
     3857  hilb(std(I));
     3858  hilb(std(I), 1);
     3859
     3860  def h = hilbertSeries(I); setring h;
     3861
     3862  numerator1; denominator1;
     3863 
     3864 
     3865  kill h; 
     3866  ////////////////////////////////////////////////
     3867  setring R;
     3868
     3869  I = x^10; I;
     3870
     3871  hilb(std(I));
     3872
     3873  def h = hilbertSeries(I); setring h;
     3874
     3875  numerator1; denominator1;
    27343876
    27353877  kill h;
     
    27373879  setring R;
    27383880
    2739   I = x^5; I;
    2740 
    2741   hilb(std(I));
    2742   hilb(std(I), 1);
    2743 
    2744   def h = hilbertSeries(I); setring h;
     3881  module M = 1;
     3882
     3883  M = setModuleGrading(M, W);
     3884
     3885 
     3886  hilb(std(M));
     3887
     3888  def h = hilbertSeries(M); setring h;
    27453889
    27463890  numerator1; denominator1;
    2747 
    27483891
    27493892  kill h;
     
    27513894  setring R;
    27523895
    2753   I = x^10; I;
    2754 
    2755   hilb(std(I));
    2756 
    2757   def h = hilbertSeries(I); setring h;
     3896  kill M; module M = x^5*gen(1);
     3897//  intmat V[1][3] = 0; // TODO: this would lead to a wrong result!!!?
     3898  intmat V[1][1] = 0; // all gen(i) of degree 0!
     3899
     3900  M = setModuleGrading(M, V);
     3901
     3902  hilb(std(M));
     3903
     3904  def h = hilbertSeries(M); setring h;
    27583905
    27593906  numerator1; denominator1;
     
    27613908  kill h;
    27623909  ////////////////////////////////////////////////
    2763   setring R;
    2764 
    2765   module M = 1;
    2766 
    2767   M = setModuleGrading(M, W);
    2768 
    2769 
    2770   hilb(std(M));
    2771 
    2772   def h = hilbertSeries(M); setring h;
     3910  setring R;
     3911
     3912  module N = x^5*gen(3);
     3913
     3914  kill V;
     3915 
     3916  intmat V[1][3] = 0; // all gen(i) of degree 0!
     3917
     3918  N = setModuleGrading(N, V);
     3919     
     3920  hilb(std(N));
     3921
     3922  def h = hilbertSeries(N); setring h;
    27733923
    27743924  numerator1; denominator1;
     
    27763926  kill h;
    27773927  ////////////////////////////////////////////////
    2778   setring R;
    2779 
    2780   kill M; module M = x^5*gen(1);
    2781 //  intmat V[1][3] = 0; // TODO: this would lead to a wrong result!!!?
    2782   intmat V[1][1] = 0; // all gen(i) of degree 0!
    2783 
    2784   M = setModuleGrading(M, V);
    2785 
    2786   hilb(std(M));
    2787 
    2788   def h = hilbertSeries(M); setring h;
    2789 
    2790   numerator1; denominator1;
    2791 
    2792   kill h;
    2793   ////////////////////////////////////////////////
    2794   setring R;
    2795 
    2796   module N = x^5*gen(3);
    2797 
    2798   kill V;
    2799 
    2800   intmat V[1][3] = 0; // all gen(i) of degree 0!
    2801 
    2802   N = setModuleGrading(N, V);
    2803 
    2804   hilb(std(N));
    2805 
    2806   def h = hilbertSeries(N); setring h;
    2807 
    2808   numerator1; denominator1;
    2809 
    2810   kill h;
    2811   ////////////////////////////////////////////////
    2812   setring R;
    2813 
    2814 
     3928  setring R;
     3929
     3930 
    28153931  module S = M + N;
    2816 
     3932 
    28173933  S = setModuleGrading(S, V);
    28183934
     
    28303946}
    28313947
    2832 
     3948proc evalHilbertSeries(def h, intvec v)
     3949"
     3950evaluate hilbert series h by substibuting v[i] for t_(i) (1/v[i] for s_(i))
     3951return: int (h(v))
     3952"
     3953{
     3954  if( 2*size(v) != nvars(h) )
     3955  {
     3956    ERROR("Wrong input/size!");
     3957  }
     3958
     3959  setring h;
     3960
     3961  if( defined(numerator2) and defined(denominator2) )
     3962  {
     3963    poly n = numerator2; poly d = denominator2;
     3964  } else
     3965  {
     3966    poly n = numerator1; poly d = denominator1;
     3967  }
     3968
     3969  int N = size(v);
     3970  int i; number k;
     3971  ideal V;
     3972
     3973  for( i = N; i > 0; i -- )
     3974  {
     3975    k = v[i];
     3976    V[i] = var(i) - k;
     3977  }
     3978   
     3979  V = groebner(V);
     3980   
     3981  n = NF(n, V);
     3982  d = NF(d, V);
     3983
     3984  n;
     3985  d;
     3986
     3987  if( d == 0 )
     3988  {
     3989    ERROR("Sorry: denominator is zero!");
     3990  }
     3991 
     3992  if( n == 0 )
     3993  {
     3994    return (0);
     3995  }
     3996
     3997  poly g = gcd(n, d);
     3998 
     3999  if( g != leadcoef(g) )
     4000  {
     4001    n = n / g;
     4002    d = d / g;
     4003  }
     4004
     4005  n;
     4006  d;
     4007   
     4008   
     4009  for( i = N; i > 0; i -- )
     4010  {
     4011    "i: ", i;
     4012    n;
     4013    d;
     4014     
     4015    k = v[i];
     4016    k;
     4017     
     4018    n = subst(n, var(i), k);
     4019    d = subst(d, var(i), k);
     4020   
     4021    if( k != 0 )
     4022    {
     4023      k = 1/k;
     4024      n = subst(n, var(N+i), k);
     4025      d = subst(d, var(N+i), k);
     4026    }
     4027  }
     4028
     4029  n;
     4030  d;
     4031
     4032  if( d == 0 )
     4033  {
     4034    ERROR("Sorry: denominator is zero!");
     4035  }
     4036 
     4037  if( n == 0 )
     4038  {
     4039    return (0);
     4040  }
     4041
     4042  poly g = gcd(n, d);
     4043 
     4044  if( g != leadcoef(g) )
     4045  {
     4046    n = n / g;
     4047    d = d / g;
     4048  }
     4049
     4050  n;
     4051  d;
     4052   
     4053  if( n != leadcoef(n) || d != leadcoef(d) )
     4054  {
     4055    ERROR("Sorry cannot completely evaluate. Partial result: (" + string(n) + ")/(" + string(d) + ")");
     4056  }
     4057
     4058  n;
     4059  d;
     4060
     4061  return (leadcoef(n)/leadcoef(d));
     4062}
     4063example
     4064{
     4065  "EXAMPLE:"; echo=2;
     4066
     4067  // TODO!
     4068
     4069}
     4070
     4071
     4072proc isPositive()
     4073"USAGE: isPositive()
     4074PURPOSE: Computes whether the multigrading of the ring is positive.
     4075For computation theorem 8.6 of the Miller/Sturmfels book is used.
     4076RETURNS: true if the multigrading is positive
     4077EXAMPLE: example isPositive; shows an example
     4078"
     4079{
     4080ideal I = mDegBasis(0);
     4081ideal J = attrib(I,"ZeroPart");
     4082/*
     4083I am not quite sure what this ZeroPart is anymore. I thought it
     4084should contain all monomials of degree 0, but then apparently 1 should
     4085be contained. It makes sense to exclude 1, but was this also the intention?
     4086*/
     4087return(J==0);
     4088}
     4089example
     4090{
     4091  echo = 2; printlevel = 3;
     4092  ring r = 0,(x,y),dp;
     4093  intmat A[1][2]=-1,1;
     4094  setBaseMultigrading(A);
     4095  isPositive();
     4096 
     4097  intmat A[1][2]=1,1;
     4098  setBaseMultigrading(A);
     4099  isPositive(A);
     4100}
    28334101
    28344102///////////////////////////////////////////////////////////////////////////////
     
    28364104proc testMultigradingLib ()
    28374105{
    2838   echo = 2; printlevel = 3;
    28394106  example setBaseMultigrading;
    28404107  example setModuleGrading;
    28414108
    28424109  example getVariableWeights;
    2843   example getTorsion;
     4110  example getLattice;
     4111  example getGradingGroup;
    28444112  example getModuleGrading;
    28454113
     
    28494117
    28504118
    2851   example hermite;
    2852   example isHomogenous;
     4119  example hermiteNormalForm;
     4120  example isHomogeneous;
    28534121  example isTorsionFree;
    28544122  example pushForward;
    2855   example defineHomogenous;
     4123  example defineHomogeneous;
    28564124
    28574125  example equalMDeg;
    2858   example isTorsionElement;
     4126  example isZeroElement;
    28594127
    28604128  example mDegResolution;
    2861 
    2862   "// ******************* example hilbertSeries ************************//";
     4129 
     4130  "// ******************* example hilbertSeries ************************//"; 
    28634131  example hilbertSeries;
    28644132
     
    28674135
    28684136  "The End!";
    2869 
    2870 }
     4137}
     4138
     4139
     4140proc mDegTruncate(def M, intvec md)
     4141{
     4142  "d: ";
     4143  print(md);
     4144 
     4145  "M: ";
     4146  module LL = M; // + L for d+1
     4147  LL;
     4148  print(mDeg(LL));
     4149
     4150
     4151  intmat V = getModuleGrading(M);
     4152  intvec vi;
     4153  int s = nrows(M);
     4154  int r = nrows(V);
     4155  int i;
     4156  module L; def B;
     4157  for (i=s; i>0; i--)
     4158  {
     4159    "comp: ", i;
     4160    vi = V[1..r, i];
     4161    "v[i]: "; vi;
     4162
     4163    B = mDegBasis(md - vi); // ZeroPart is always the same...
     4164    "B: "; B;
     4165
     4166    L = L, B*gen(i);
     4167  }
     4168  L = simplify(L, 2);
     4169  L = setModuleGrading(L,V);
     4170
     4171  "L: "; L;
     4172  print(mDeg(L));
     4173
     4174  L = mDegModulo(L, LL);
     4175  L = mDegGroebner(L);
     4176//  L = minbase(prune(L));
     4177
     4178  "??????????";
     4179  print(L);
     4180  print(mDeg(L));
     4181   
     4182  V = getModuleGrading(L);
     4183
     4184  // take out other degrees
     4185  for(i = ncols(L); i > 0; i-- )
     4186  {
     4187    if( !equalMDeg( mDeg(getGradedGenerator(L, i)), md ) )
     4188    {
     4189      L[i] = 0;
     4190    }
     4191  }
     4192 
     4193  L = simplify(L, 2);
     4194  L = setModuleGrading(L, V);
     4195  print(L);
     4196  print(mDeg(L));
     4197 
     4198  return(L);
     4199}
     4200example
     4201{
     4202  "EXAMPLE:"; echo=2;
     4203
     4204  // TODO!
     4205  ring r = 32003, (x,y), dp;
     4206 
     4207  intmat M[2][2] =
     4208    1, 0,
     4209    0, 1;
     4210
     4211  setBaseMultigrading(M);
     4212
     4213  intmat V[2][1] =
     4214    0,
     4215    0;
     4216 
     4217  "X:";
     4218  module h1 = x;
     4219  h1 = setModuleGrading(h1, V);
     4220  mDegTruncate(h1, mDeg(x));
     4221  mDegTruncate(h1, mDeg(y));
     4222
     4223  "XY:";
     4224  module h2 = ideal(x, y);
     4225  h2 = setModuleGrading(h2, V);
     4226  mDegTruncate(h2, mDeg(x));
     4227  mDegTruncate(h2, mDeg(y));
     4228  mDegTruncate(h2, mDeg(xy));
     4229}
     4230
     4231
     4232/******************************************************/
     4233/* Some functions on lattices.
     4234TODO Tuebingen: - add functionality (see wiki) and
     4235- adjust them to work for groups as well.*/
     4236/******************************************************/
     4237
     4238
     4239
     4240/******************************************************/
     4241proc imageLattice(intmat Q, intmat L)
     4242"USAGE: imageLattice(Q,L); Q and L are of type intmat
     4243PURPOSE: compute an integral basis for the image of the
     4244lattice L under the homomorphism of lattices Q.
     4245RETURN: intmat
     4246EXAMPLE: example imageLattice; shows an example
     4247"
     4248{
     4249  intmat Mul = Q*L;
     4250  intmat L = latticeBasis(Mul);
     4251
     4252  return(L);
     4253}
     4254example
     4255{
     4256  "EXAMPLE:"; echo=2;
     4257   
     4258  intmat Q[2][3] =
     4259    1,2,3,
     4260    3,2,1;
     4261
     4262  intmat L[3][2] =
     4263    1,4,
     4264    2,5,
     4265    3,6;
     4266
     4267  // should be a 2x2 matrix with columns
     4268  // [2,-14], [0,36]
     4269  imageLattice(Q,L);
     4270
     4271}
     4272
     4273/******************************************************/
     4274proc intRank(intmat A)
     4275"
     4276USAGE: intRank(A); intmat A
     4277PURPOSE: compute the rank of the integral matrix A
     4278by computing a hermite normalform.
     4279RETURNS: int
     4280EXAMPLE: example intRank; shows an example
     4281"
     4282{
     4283  intmat B = hermiteNormalForm(A);
     4284
     4285  // get number of zero columns
     4286  int nzerocols = 0;
     4287  int j;
     4288  int i;
     4289  int iszero;
     4290  for ( j = 1; j <= ncols(B); j++ )
     4291  {
     4292    iszero = 1;
     4293   
     4294    for ( i = 1; i <= nrows(B); i++ )
     4295    {
     4296      if ( B[i,j] != 0 )
     4297      {
     4298        iszero = 0;
     4299        break;
     4300      }
     4301    }
     4302   
     4303    if ( iszero == 1 )
     4304    {
     4305      nzerocols++;
     4306    }
     4307  }
     4308
     4309  // get number of zero rows
     4310  int nzerorows = 0;
     4311
     4312  for ( i = 1; i <= nrows(B); i++ )
     4313  {
     4314    iszero = 1;
     4315   
     4316    for ( j = 1; j <= ncols(B); j++ )
     4317    {
     4318      if ( B[i,j] != 0 )
     4319      {
     4320        iszero = 0;
     4321        break;
     4322      }
     4323    }
     4324   
     4325    if ( iszero == 1 )
     4326    {
     4327      nzerorows++;
     4328    }
     4329  }
     4330
     4331  int r = nrows(B) - nzerorows;
     4332
     4333  if ( ncols(B) - nzerocols < r )
     4334  {
     4335    r = ncols(B) - nzerocols;
     4336  }
     4337 
     4338  return(r);
     4339}
     4340example
     4341{
     4342
     4343  intmat A[3][4] =
     4344    1,0,1,0,
     4345    1,2,0,0,
     4346    0,0,0,0;
     4347
     4348  int r = intRank(A);
     4349
     4350  print(A);
     4351  print(r); // Should be 2
     4352
     4353  kill A;
     4354
     4355}
     4356
     4357/*****************************************************/
     4358
     4359proc isSublattice(intmat L, intmat S)
     4360"USAGE: isSublattice(L, S); L, S are of tpye intmat
     4361PURPOSE: checks whether the lattice created by L is a
     4362sublattice of the lattice created by S.
     4363The procedure checks whether each generator of L is
     4364contained in S.
     4365RETURN: 0 if false, 1 if true
     4366EXAMPLE: example isSublattice; shows an example
     4367"
     4368{
     4369  int a,b,g,i,j,k;
     4370  intmat Ker;
     4371   
     4372  // check whether each column v of L is contained in
     4373  // the lattice generated by S
     4374  for ( i = 1; i <= ncols(L); i++ )
     4375  {
     4376   
     4377    // v is the i-th column of L
     4378    intvec v;
     4379     for ( j = 1; j <= nrows(L); j++ )
     4380    {
     4381      v[j] = L[j,i];
     4382    }
     4383
     4384    // concatenate B = [S,v]
     4385    intmat B[nrows(L)][ncols(S) + 1];
     4386
     4387    for ( a = 1; a <= nrows(S); a++ )
     4388    {
     4389      for ( b = 1; b <= ncols(S); b++ )
     4390      {
     4391        B[a,b] = S[a,b];
     4392      }
     4393    }
     4394
     4395    for ( a = 1; a <= size(v); a++ )
     4396    {
     4397      B[a,ncols(B)] = v[a];
     4398    }
     4399
     4400   
     4401    // check gcd
     4402    Ker = kernelLattice(B);
     4403    k = nrows(Ker);
     4404    list R; // R is the last row
     4405
     4406    for ( j = 1; j <= ncols(Ker); j++ )
     4407    {
     4408      R[j] = Ker[k,j];
     4409    }
     4410
     4411    g = R[1];
     4412   
     4413    for ( j = 2; j <= size(R); j++ )
     4414    {
     4415      g = gcd(g,R[j]);
     4416    }
     4417
     4418    if ( g != 1 )
     4419    {
     4420      return(0);
     4421    }
     4422
     4423    kill B, v, R;
     4424
     4425  }
     4426
     4427  return(1);
     4428}
     4429example
     4430{
     4431  "EXAMPLE:"; echo=2;
     4432   
     4433  //ring R = 0,(x,y),dp;
     4434  intmat S2[3][2]=
     4435    0, 2, 3,
     4436    0, 1, 1,
     4437    3, 0, 2;
     4438
     4439  intmat S1[3][3]=
     4440    0, 6,
     4441    0, 2,
     4442    3, 4;
     4443
     4444  isSublattice(S1,S2); // Yes!
     4445
     4446  intmat S3[3][1] =
     4447    0,
     4448    0,
     4449    1;
     4450
     4451  not(isSublattice(S3,S2)); // Yes!
     4452
     4453}
     4454
     4455/******************************************************/
     4456
     4457proc latticeBasis(intmat B)
     4458"USAGE: latticeBasis(B); intmat B
     4459PURPOSE: compute an integral basis for the lattice defined by
     4460the columns of B.
     4461RETURNS: intmat
     4462EXAMPLE: example latticeBasis; shows an example
     4463"
     4464{
     4465  int n = ncols(B);
     4466  int r = intRank(B);
     4467 
     4468  if ( r == 0 )
     4469  {
     4470    intmat H[nrows(B)][1];
     4471    int j;
     4472
     4473    for ( j = 1; j <= nrows(B); j++ )
     4474    {
     4475      H[j,1] = 0;   
     4476    }
     4477  }
     4478  else
     4479  {
     4480    intmat H = hermiteNormalForm(B);;
     4481
     4482    if (r < n)
     4483    {
     4484      // delete columns r+1 to n
     4485      // should be identical with the function
     4486      // H = submat(H,1..nrows(H),1..r);
     4487      // for matrices
     4488      intmat Hdel[nrows(H)][r];
     4489      int k;
     4490      int m;
     4491     
     4492      for ( k = 1; k <= nrows(H); k++ )
     4493      {
     4494        for ( m = 1; m <= r; m++ )
     4495        {
     4496          Hdel[k,m] = H[k,m];
     4497        }
     4498      }
     4499
     4500      H = Hdel;     
     4501    }
     4502  }
     4503 
     4504  return(H); 
     4505}
     4506example
     4507{
     4508  "EXAMPLE:"; echo=2;
     4509 
     4510  intmat L[3][3] =
     4511    1,4,8,
     4512    2,5,10,
     4513    3,6,12;
     4514
     4515  intmat B = latticeBasis(B); // should be a matrix with columns [1,2,3] and [0,3,6]
     4516
     4517  kill B,L;
     4518}
     4519
     4520/******************************************************/
     4521
     4522proc kernelLattice(def P)
     4523"
     4524USAGE: kernelLattice(P); intmat P
     4525PURPOSE: compute a integral basis for the kernel of the
     4526homomorphism of lattices defined by the intmat P.
     4527RETURNS: intmat
     4528EXAMPLE: example kernelLattice; shows an example
     4529"
     4530{
     4531  int n = ncols(P);
     4532  int r = intRank(P);
     4533
     4534  if ( r == 0 )
     4535  {
     4536    intmat U = unitMatrix(n);
     4537  }
     4538  else
     4539  {
     4540    if ( r == n )
     4541    {
     4542      intmat U[n][1];  // now all entries are zero.
     4543    }
     4544    else
     4545    {
     4546      list L = hermiteNormalForm(P, "transform"); //hermite(P, "transform");  // now, Hermite = L[1] = A*L[2]
     4547      intmat U = L[2];
     4548
     4549      // delete columns 1 to r
     4550      // should be identical with the function
     4551      // U = submat(U,1..nrows(U),r+1..);
     4552      // for matrices
     4553      intmat Udel[nrows(U)][ncols(U) - r];
     4554      int k;
     4555      int m;
     4556     
     4557      for ( k = 1; k <= nrows(U); k++ )
     4558      {
     4559        for ( m = r + 1; m <= ncols(U); m++ )
     4560        {
     4561          Udel[k,m - r] = U[k,m];
     4562        }
     4563      }
     4564
     4565      U = Udel; 
     4566
     4567    }
     4568  }
     4569
     4570  return(U);
     4571}
     4572example
     4573{
     4574  "EXAMPLE"; echo = 2;
     4575
     4576  intmat LL[3][4] =
     4577    1,4,7,10,
     4578    2,5,8,11,
     4579    3,6,9,12;
     4580
     4581  // should be a 4x2 matrix with colums
     4582  // [-1,2,-1,0],[2,-3,0,1]
     4583  intmat B = kernelLattice(LL);
     4584
     4585  print(B);
     4586
     4587  kill B;
     4588
     4589}
     4590
     4591/*****************************************************/
     4592
     4593proc preimageLattice(def P, def B)
     4594"
     4595USAGE: preimageLattice(P, B); intmat P, intmat B
     4596PURPOSE: compute an integral basis for the preimage of B under
     4597 the homomorphism of lattices defined by the intmat P.
     4598RETURNS: intmat
     4599EXAMPLE: example preimageLattice; shows an example
     4600"
     4601{
     4602  // concatenate matrices: Con = [P,-B]
     4603  intmat Con[nrows(P)][ncols(P) + ncols(B)];
     4604  int i;
     4605  int j;
     4606
     4607  for ( i = 1; i <= nrows(Con); i++ )
     4608  {
     4609    for ( j = 1; j <= ncols(P); j++ ) // P first
     4610    {
     4611      Con[i,j] = P[i,j];
     4612    }
     4613  }
     4614
     4615  for ( i = 1; i <= nrows(Con); i++ )
     4616  {
     4617    for ( j = 1; j <= ncols(B); j++ ) // now -B
     4618    {
     4619      Con[i,ncols(P) + j] = - B[i,j];
     4620    }
     4621  }
     4622
     4623
     4624  print(Con);
     4625
     4626  intmat L = kernelLattice(Con);
     4627
     4628  print(L);
     4629  print(ncols(P));
     4630  print(ncols(L));
     4631
     4632  // delete rows ncols(P)+1 to nrows(L) out of L
     4633  intmat Del[ncols(P)][ncols(L)];
     4634  int k;
     4635  int m;
     4636     
     4637  for ( k = 1; k <= nrows(Del); k++ )
     4638  {
     4639    for ( m = 1; m <= ncols(Del); m++ )
     4640    {
     4641      Del[k,m] = L[k,m];
     4642    }
     4643  }
     4644 
     4645  L = latticeBasis(Del);
     4646
     4647  return(L);
     4648
     4649}
     4650example
     4651{
     4652  "EXAMPLE"; echo = 2;
     4653
     4654  intmat P[2][3] =
     4655    2,6,10,
     4656    4,8,12;
     4657
     4658  intmat B[2][1] =
     4659    1,
     4660    0;
     4661
     4662  // should be a 2x2 matrix with columns
     4663  // [1,1,-1] and [-2,1,0]
     4664  intmat L = preimageLattice(P,B);
     4665
     4666  kill B, P, L;
     4667
     4668}
     4669
     4670/******************************************************/
     4671proc isPrimitiveSublattice(intmat A);
     4672"USAGE: isPrimitiveSublattice(A); intmat A
     4673PURPOSE: check whether the given set of integral vectors in ZZ^m,
     4674i.e. the columns of A, generate a primitive sublattice in ZZ^m
     4675(a direct summand of ZZ^m).
     4676RETURNS: int, where 0 is false and 1 is true.
     4677EXAMPLE: example isPrimitiveSublattice; shows an example
     4678"
     4679{
     4680  intmat B = smithNormalForm(A);
     4681  int r = intRank(B);
     4682 
     4683  if ( r == 0 )
     4684  {
     4685    return(1);
     4686  }
     4687
     4688  if ( 1 < B[r,r] )
     4689  {
     4690    return(0);
     4691  }
     4692
     4693  return(1);
     4694}
     4695example
     4696{
     4697  "EXAMPLE"; echo = 2;
     4698 
     4699  intmat A[3][2] =
     4700    1,4,
     4701    2,5,
     4702    3,6;
     4703
     4704  // should be 0
     4705  int b = isPrimitiveSublattice(A);
     4706 
     4707  kill A,b;
     4708}
     4709
     4710/******************************************************/
     4711proc isIntegralSurjective(intmat P);
     4712"USAGE: isIntegralSurjective(P); intmat P
     4713PURPOSE: test whether the given linear map P of lattices is
     4714surjective.
     4715RETURNS: int, where 0 is false and 1 is true.
     4716EXAMPLE: example isIntegralSurjective; shows an example
     4717"
     4718{
     4719  int r = intRank(P);
     4720 
     4721  if ( r < nrows(P) )
     4722  {
     4723    return(0);
     4724  }
     4725
     4726  if ( isPrimitiveSublattice(A) == 1 )
     4727  {
     4728    return(1);
     4729  }
     4730
     4731  return(0);
     4732}
     4733example
     4734{
     4735  "EXAMPLE"; echo = 2;
     4736 
     4737  intmat A[3][2] =
     4738    1,3,5,
     4739    2,4,6;
     4740   
     4741  // should be 0
     4742  int b = isIntegralSurjective(A);
     4743  print(b);
     4744 
     4745  kill A,b;
     4746}
     4747
     4748/******************************************************/
     4749proc projectLattice(intmat B)
     4750"USAGE: projectLattice(B); intmat B
     4751PURPOSE: A set of vectors in ZZ^m is given as the columns of B.
     4752Then this function provides a linear map ZZ^m --> ZZ^n
     4753having the primitive span of B its kernel.
     4754RETURNS: intmat
     4755EXAMPLE: example projectLattice; shows an example
     4756"
     4757{
     4758  int n = nrows(B);
     4759  int r = intRank(B);
     4760
     4761  if ( r == 0 )
     4762  {
     4763    intmat U = unitMatrix(n);
     4764  }
     4765  else
     4766  {
     4767    if ( r == n )
     4768    {
     4769      intmat U[n][1]; // U now is the n-dim zero-vector
     4770    }
     4771    else
     4772    {
     4773      // we want a matrix with column operations so we transpose
     4774      list L = hermiteNormalForm(B, "transform"); //hermite(transpose(B), "transform");
     4775      intmat U = transpose(L[2]);
     4776
     4777      // delete rows 1 to r
     4778      intmat Udel[nrows(U) - r][ncols(U)];
     4779      int k;
     4780      int m;
     4781     
     4782      for ( k = 1; k <= nrows(U) - r ; k++ )
     4783      {
     4784        for ( m = 1; m <= ncols(U); m++ )
     4785        {
     4786          Udel[k,m] = U[k + r,m];
     4787        }
     4788      }
     4789
     4790      U = Udel; 
     4791     
     4792    }
     4793  }
     4794 
     4795  return(U);
     4796}
     4797example
     4798{
     4799  "EXAMPLE"; echo = 2;
     4800 
     4801  intmat B[4][2] =
     4802    1,5,
     4803    2,6,
     4804    3,7,
     4805    4,8;
     4806 
     4807  // should result in a (2x4)-matrix with columns
     4808  // [-1, 2], [2, −3], [-1, 0] and [0, 1].
     4809  intmat U = projectLattice(B);
     4810 
     4811}
     4812
     4813/******************************************************/
     4814proc intersectLattices(intmat A, intmat B)
     4815"USAGE: intersectLattices(A, B); intmat A, intmat B
     4816PURPOSE: compute an integral basis for the intersection of the
     4817lattices A and B.
     4818RETURNS: intmat
     4819EXAMPLE: example intersectLattices; shows an example
     4820"
     4821{
     4822  // concatenate matrices: Con = [A,-B]
     4823  intmat Con[nrows(A)][ncols(A) + ncols(B)];
     4824  int i;
     4825  int j;
     4826
     4827  for ( i = 1; i <= nrows(Con); i++ )
     4828  {
     4829    for ( j = 1; j <= ncols(A); j++ ) // A first
     4830    {
     4831      Con[i,j] = A[i,j];
     4832    }
     4833  }
     4834
     4835  for ( i = 1; i <= nrows(Con); i++ )
     4836  {
     4837    for ( j = 1; j <= ncols(B); j++ ) // now -B
     4838    {
     4839      Con[i,ncols(A) + j] = - B[i,j];
     4840    }
     4841  }
     4842
     4843  intmat K = kernelLattice(Con);
     4844
     4845  // delete all rows in K from ncols(A)+1 onwards
     4846  intmat Bas[ncols(A)][ncols(K)];
     4847 
     4848  for ( i = 1; i <= nrows(Bas); i++ )
     4849  {
     4850    for ( j = 1; j <= ncols(Bas); j++ )
     4851    {
     4852      Bas[i,j] = K[i,j];
     4853    }
     4854  }
     4855
     4856  // take product in order to obtain the intersection
     4857  intmat S = A * Bas;
     4858  intmat Cut = hermiteNormalForm(S); //hermite(S);
     4859  int r = intRank(Cut);
     4860
     4861  if ( r == 0 )
     4862  {
     4863    intmat Cutdel[nrows(Cut)][1]; // is now the zero-vector
     4864
     4865    Cut = Cutdel;
     4866  }
     4867  else
     4868  {
     4869    // delte columns from r+1 onwards
     4870    intmat Cutdel[nrows(Cut)][r];
     4871
     4872    for ( i = 1; i <= nrows(Cutdel); i++ )
     4873    {
     4874      for ( j = 1; j <= r; j++ )
     4875      {
     4876        Cutdel[i,j] = Cut[i,j];
     4877      }
     4878    }
     4879
     4880    Cut = Cutdel;
     4881  }
     4882
     4883  return(Cut);
     4884}
     4885example
     4886{
     4887  "EXAMPLE"; echo = 2;
     4888 
     4889  intmat A[3][2] =
     4890    1,4,
     4891    2,5,
     4892    3,6;
     4893
     4894  intmat B[3][2] =
     4895    6,9,
     4896    7,10,
     4897    8,11;
     4898 
     4899  // should result in a (2x4)-matrix with columns
     4900  // [3, 3, 3], [0, 3, 6]
     4901  intmat U = intersectLattices(A,B);
     4902 
     4903}
     4904
     4905proc intInverse(intmat A);
     4906"USAGE: intInverse(A); intmat A
     4907PURPOSE: compute the integral inverse of the intmat A.
     4908If det(A) is neither 1 nor -1 an error is returned.
     4909RETURNS: intmat
     4910EXAMPLE: example intInverse; shows an example
     4911"
     4912{
     4913  int d = det(A);
     4914 
     4915  if ( d * d != 1 ) // is d = 1 or -1? Else: error
     4916  {
     4917    ERROR("determinant of the given intmat has to be 1 or -1.");
     4918  }
     4919 
     4920  int c;
     4921  int i,j;
     4922  intmat C[nrows(A)][ncols(A)];
     4923  intmat Ad;
     4924  int s;
     4925
     4926  for ( i = 1; i <= nrows(C); i++ )
     4927  {
     4928    for ( j = 1; j <= ncols(C); j++ )
     4929    {
     4930      Ad = intAdjoint(A,i,j);
     4931      s = 1;
     4932   
     4933      if ( ((i + j) % 2) > 0 )
     4934      {
     4935        s = -1;
     4936      }
     4937
     4938      C[i,j] = d * s * intDet(Ad); // mult by d is equal to div by det
     4939    }
     4940  }
     4941
     4942  C = transpose(C);
     4943
     4944  return(C);
     4945}
     4946example
     4947{
     4948  "EXAMPLE"; echo = 2;
     4949
     4950  intmat A[3][3] =
     4951    1,1,3,
     4952    3,2,0,
     4953    0,0,1;
     4954
     4955  intmat B = intInverse(A);
     4956
     4957  // should be the unit matrix
     4958  print(A * B);
     4959
     4960  kill A,B;
     4961}
     4962
     4963
     4964/******************************************************/
     4965proc intAdjoint(intmat A, int indrow, int indcol)
     4966"USAGE: intAdjoint(A); intmat A
     4967PURPOSE: return the matrix where the given row and column are deleted.
     4968RETURNS: intmat
     4969EXAMPLE: example intAdjoint; shows an example
     4970"
     4971{
     4972  int n = nrows(A);
     4973  int m = ncols(A);
     4974  int i, j;
     4975  intmat B[n - 1][m - 1];
     4976  int a, b;
     4977
     4978  for ( i = 1; i < indrow; i++ )
     4979  {
     4980    for ( j = 1; j < indcol; j++ )
     4981    {
     4982      B[i,j] = A[i,j];
     4983    }
     4984    for ( j = indcol + 1; j <= ncols(A); j++ )
     4985    {
     4986      B[i,j - 1] = A[i,j];
     4987    }
     4988  }
     4989
     4990  for ( i = indrow + 1; i <= nrows(A); i++ )
     4991  {
     4992    for ( j = 1; j < indcol; j++ )
     4993    {
     4994      B[i - 1,j] = A[i,j];
     4995    }
     4996    for ( j = indcol+1; j <= ncols(A); j++ )
     4997    {
     4998      B[i - 1,j - 1] = A[i,j];
     4999    }
     5000  }
     5001
     5002  return(B);
     5003}
     5004example
     5005{
     5006  "EXAMPLE"; echo = 2;
     5007 
     5008  intmat A[2][3] =
     5009    1,3,5,
     5010    2,4,6;
     5011
     5012  intmat B = intAdjoint(A,2,2);
     5013  print(B);
     5014
     5015  kill A,B;
     5016}
     5017
     5018/******************************************************/
     5019proc integralSection(intmat P);
     5020"USAGE: integralSection(P); intmat P
     5021PURPOSE: for a given linear surjective map P of lattices
     5022 this procedure returns an integral section of P.
     5023RETURNS: intmat
     5024EXAMPLE: example integralSection; shows an example
     5025"
     5026{
     5027  int m = nrows(P);
     5028  int n = ncols(P);
     5029 
     5030  if ( m == n )
     5031  {
     5032    intmat U = intInverse(P);
     5033  }
     5034  else
     5035  {
     5036    intmat U = (hermiteNormalForm(P, "transform"))[2];
     5037   
     5038    // delete columns m+1 to n
     5039    intmat Udel[nrows(U)][ncols(U) - (n - m)];
     5040    int k;
     5041    int z;
     5042     
     5043    for ( k = 1; k <= nrows(U); k++ )
     5044    {
     5045      for ( z = 1; z <= m; z++ )
     5046      {
     5047        Udel[k,z] = U[k,z];
     5048      }
     5049    }
     5050   
     5051    U = Udel; 
     5052  }
     5053
     5054  return(U);
     5055}
     5056example
     5057{
     5058  "EXAMPLE"; echo = 2;
     5059
     5060  intmat P[2][4] =
     5061    1,3,4,6,
     5062    2,4,5,7;
     5063
     5064  // should be a matrix with two columns
     5065  // for example: [−2, 1, 0, 0], [3, −3, 0, 1]
     5066  intmat U = integralSection(P);
     5067
     5068  print(U);
     5069  print(P * U);
     5070
     5071  kill U; 
     5072}
     5073
     5074
     5075
     5076/******************************************************/
     5077proc factorgroup(G,H)
     5078"USAGE: factorgroup(G,H); list G, list H
     5079PURPOSE: returns a representation of the factor group G mod H using the first isomorphism thm
     5080RETURNS: list
     5081EXAMPLE: example factorgroup(G,H); shows an example
     5082"
     5083{
     5084  intmat S1 = G[1];
     5085  intmat L1 = G[2];
     5086  intmat S2 = H[1];
     5087  intmat L2 = H[2];
     5088
     5089  // check whether G,H are subgroups of a common group, i.e. whether L1 and L2 span the same lattice 
     5090  if ( !isSublattice(L1,L2) || !isSublattice(L2,L1))
     5091  {
     5092    ERROR("G and H are not subgroups of a common group.");
     5093  }
     5094
     5095  // check whether H is a subgroup of G, i.e. whether S2 is a sublattice of S1+L1
     5096  intmat B = concatintmat(S1,L1); // check whether this gives the concatinated matrix
     5097  if ( !isSublattice(S2,B) )
     5098  {
     5099    ERROR("H is not a subgroup of G");
     5100  }
     5101  // use first isomorphism thm to get the factor group
     5102  intmat L = concatintmat(L1,S2); // check whether this gives the concatinated matrix
     5103  list GmodH;
     5104  GmodH[1]=S1;
     5105  GmodH[2]=L;
     5106  return(GmodH);
     5107}
     5108example
     5109{
     5110  "EXAMPLE"; echo = 2;
     5111 
     5112  intmat S1[2][2] =
     5113    1,0,
     5114    0,1;
     5115  intmat L1[2][1] =
     5116    2,
     5117    0;
     5118
     5119  intmat S2[2][1] =
     5120    1,
     5121    0;
     5122  intmat L2[2][1] =
     5123    2,
     5124    0;
     5125
     5126  list G = createGroup(S1,L1);
     5127  list H = createGroup(S2,L2);
     5128
     5129  list N = factorgroup(G,H);
     5130  print(N);
     5131
     5132  kill G,H,N,S1,L1,S2,L2;
     5133   
     5134}
     5135
     5136/******************************************************/
     5137proc productgroup(G,H)
     5138"USAGE: productgroup(G,H); list G, list H
     5139PURPOSE: returns a representation of the G x H
     5140RETURNS: list
     5141EXAMPLE: example productgroup(G,H); shows an example
     5142"
     5143{
     5144  intmat S1 = G[1];
     5145  intmat L1 = G[2];
     5146  intmat S2 = H[1];
     5147  intmat L2 = H[2];
     5148  intmat OS1[nrows(S1)][ncols(S2)];
     5149  intmat OS2[nrows(S2)][ncols(S1)];
     5150  intmat OL1[nrows(L1)][ncols(L2)];
     5151  intmat OL2[nrows(L2)][ncols(L1)];
     5152
     5153  // concatinate matrices to get S
     5154  intmat A = concatintmat(S1,OS1);
     5155  intmat B = concatintmat(OS2,S2);
     5156  intmat At = transpose(A);
     5157  intmat Bt = transpose(B);
     5158  intmat St = concatintmat(At,Bt);
     5159  intmat S = transpose(St);
     5160
     5161  // concatinate matrices to get L
     5162  intmat C = concatintmat(L1,OL1);
     5163  intmat D = concatintmat(OL2,L2);
     5164  intmat Ct = transpose(C);
     5165  intmat Dt = transpose(D);
     5166  intmat Lt = concatintmat(Ct,Dt);
     5167  intmat L = transpose(Lt);
     5168
     5169  list GxH;
     5170  GxH[1]=S;
     5171  GxH[2]=L;
     5172  return(GxH);
     5173}
     5174example
     5175{
     5176  "EXAMPLE"; echo = 2;
     5177 
     5178  intmat S1[2][2] =
     5179    1,0,
     5180    0,1;
     5181  intmat L1[2][1] =
     5182    2,
     5183    0;
     5184
     5185  intmat S2[2][2] =
     5186    1,0,
     5187    0,2;
     5188  intmat L2[2][1] =
     5189    0,
     5190    3;
     5191
     5192  list G = createGroup(S1,L1);
     5193  list H = createGroup(S2,L2);
     5194
     5195  list N = productgroup(G,H);
     5196  print(N);
     5197
     5198  kill G,H,N,S1,L1,S2,L2;
     5199   
     5200}
     5201
     5202/******************************************************/
     5203proc primitiveSpan(intmat V);
     5204"USAGE: isIntegralSurjective(V); intmat V
     5205PURPOSE: compute an integral basis for the minimal primitive
     5206sublattice that contains the given vectors, i.e. the columns of V.
     5207RETURNS: int, where 0 is false and 1 is true.
     5208EXAMPLE: example isIntegralSurjective; shows an example
     5209"
     5210{
     5211  int n = ncols(V);
     5212  int m = nrows(V);
     5213  int r = intRank(V);
     5214
     5215 
     5216  if ( r == 0 )
     5217  {
     5218    intmat P[m][1]; //  this is the m-zero-vector now
     5219  }
     5220  else
     5221  {
     5222    list L = smithNormalForm(V, "transform"); // L = [A,S,B] where S is the smith-NF and S = A*S*B
     5223    intmat P = intInverse(L[1]);
     5224
     5225    print(L);
     5226   
     5227    if ( r < m )
     5228    {
     5229      // delete columns r+1 to m in P:
     5230      intmat Pdel[nrows(P)][r];
     5231      int i,j;
     5232
     5233      for ( i = 1; i <= nrows(Pdel); i++ )
     5234      {
     5235        for ( j = 1; j <= ncols(Pdel); j++ )
     5236        {
     5237          Pdel[i,j] = P[i,j];
     5238        }
     5239      }
     5240
     5241      P = Pdel;
     5242    }
     5243  }
     5244
     5245  return(P);
     5246}
     5247example
     5248{
     5249  "EXAMPLE"; echo = 2;
     5250 
     5251  intmat V[2][4] =
     5252    1,4,
     5253    2,5,
     5254    3,6;
     5255
     5256  // should return a (4x2)-matrix with columns
     5257  // [1, 2, 3] and [1, 1, 1] (or similar)
     5258  intmat R = primitiveSpan(V);
     5259  print(R);
     5260
     5261  kill V,R;
     5262}
     5263
     5264/***********************************************************/
Note: See TracChangeset for help on using the changeset viewer.