1#############################################################################
2##
3#W  affine-extra-4ti2gap.gi
4#W                          Manuel Delgado <mdelgado@fc.up.pt>
5#W                          Pedro Garcia-Sanchez <pedro@ugr.es>
6##
7#Y  Copyright 2015-- Centro de Matemática da Universidade do Porto, Portugal and Universidad de Granada, Spain
8#############################################################################
9InstallOtherMethod(DegreesOfPrimitiveElementsOfAffineSemigroup,
10        "Computes the set of primitive elements of an affine semigroup",
11        [IsAffineSemigroup],4,
12        function(a)
13    local  matrix, facs, mat, trunc, ls;
14
15    ls:=MinimalGenerators(a);
16
17    Info(InfoNumSgps,2,"Using 4ti2gap for Graver.");
18
19    mat:=TransposedMat(ls);
20    matrix := GraverBasis4ti2(["mat",mat]);
21
22    trunc:=function(ls)
23        return List(ls, y->Maximum(y,0));
24    end;
25
26    matrix:=Set(matrix,trunc);
27    return Union(Set(matrix, x->x*ls),ls);
28end);
29
30
31InstallOtherMethod(HilbertBasisOfSystemOfHomogeneousEquations,
32        "Computes a Hilbert basis of a systemd of linear Diophantine equations, some eventually in congruences.",
33        [IsRectangularTable,IsHomogeneousList],6,
34        function(ls,md)
35    local  homogeneous, withCongruences;
36
37    homogeneous:= function(l)
38        local  problem, matrix,mat,sign;
39
40        Info(InfoNumSgps,2,"Using 4ti2gap for Hilbert.");
41
42        #if not(IsRectangularTable(l)) then
43        #    Error("The argument must be a matrix.");
44        #fi;
45        if not(IsInt(l[1][1])) then
46            Error("The matrix must be of integers.");
47        fi;
48
49	mat:=l;
50        sign:=[List(l[1],_->1)];
51        problem:=["mat",mat, "sign", sign];
52
53        matrix := HilbertBasis4ti2(problem).zhom;
54        return matrix;
55
56    end;
57
58  withCongruences:=function(ls,md)
59      local l,n,m,diag,dim,d, hil, zero, leq;
60
61      leq:= function(v1,v2)
62          local v;
63          v:=v2-v1;
64          return (First(v,n->n<0)=fail);
65      end;
66
67      #if not(IsRectangularTable(ls)) then
68      #    Error("The first argument must be a matrix.");
69      #fi;
70
71      if not(IsListOfIntegersNS(md)) or ForAny(md, x->not(IsPosInt(x))) then
72          Error("The second argument must be a list of positive integers.");
73      fi;
74
75      n:=Length(ls);
76      dim:=Length(ls[1]);
77      m:=Length(md);
78      if m>n then
79          Error("There are more modulus than equations.");
80      fi;
81
82      diag:=Concatenation(md,List([1..n-m],_->0));
83      d:=DiagonalMat(diag);
84      l:=TransposedMat(Concatenation(TransposedMat(ls),d,-d));
85      zero:=List([1..dim],_->0);
86
87      hil:=Difference(List(homogeneous(l), x->x{[1..dim]}),[zero]);
88      return hil;
89
90      return Filtered(hil, y->Filtered(hil,x->leq(x,y))=[y]);
91  end;
92  ## end of local functions ...
93
94  #ls := arg[1][1];
95  #md := arg[1][2];
96  if md = [] then
97      return homogeneous(ls);
98  else
99      return withCongruences(ls,md);
100
101  fi;
102
103end);
104
105InstallOtherMethod(HilbertBasisOfSystemOfHomogeneousInequalities,
106        "Computes a Hilbert basis of l*x>=0, x>=0",
107        [IsRectangularTable],6,
108        function(l)
109    local  problem, matrix,mat,sign,rel;
110
111    Info(InfoNumSgps,2,"Using 4ti2gap for Hilbert.");
112
113    #if not(IsRectangularTable(l)) then
114    #    Error("The argument must be a matrix.");
115    #fi;
116    if not(IsInt(l[1][1])) then
117        Error("The matrix must be of integers.");
118    fi;
119
120    mat:=l;
121    sign:=[List(l[1],_->1)];
122    rel:=[List(l[1],_->">")];
123    problem:=["mat",mat,"rel",rel,"sign",sign];
124    matrix:=HilbertBasis4ti2(problem);
125    return matrix;
126
127end);
128
129
130InstallOtherMethod(FactorizationsVectorWRTList,
131        "Computes the factorizations of v in terms of the elments in ls",
132        [IsHomogeneousList,IsMatrix],6,
133        function(v,l)
134    local  matrix,mat,rhs,sign,problem, n;
135
136    Info(InfoNumSgps,2,"Using 4ti2gap for factorizations.");
137
138    if not(IsListOfIntegersNS(v)) then
139        Error("The first argument must be a list of integers.");
140    fi;
141
142    if not(IsInt(l[1][1])) then
143        Error("The matrix must be of integers.");
144    fi;
145
146    mat:=TransposedMat(Concatenation(l,[-v]));
147    if not(IsRectangularTable(mat)) then
148        Error("The list in the second argument must have the same length as all the lists in the first argument.");
149    fi;
150
151    sign:=[List(l,_->1)];
152    rhs:=[v];
153    problem:=["mat",TransposedMat(l),"sign",sign,"rhs",rhs];
154    matrix := ZSolve4ti2(problem);
155    return matrix.zinhom;
156
157end);
158
159
160InstallOtherMethod(GeneratorsOfKernelCongruence,
161        "Computes a set of generators of the kernel congruence of the monoid morphism associated to a matrix",
162        [IsRectangularTable],7,
163        function(m)
164    local positivenegative, gr;
165
166    positivenegative:=function(p)
167        local d1, d2;
168        d1:=List(p, i->Maximum(i,0));
169        d2:=List(p, i->-Minimum(0,i));
170        return [d1,d2];
171    end;
172
173    if not(ForAll(m, l->ForAll(l, x->(x=0) or IsPosInt(x)))) then
174        Error("The argument must be a matrix of nonnegative integer.");
175    fi;
176
177    gr:=GroebnerBasis4ti2(TransposedMat(m));
178    Info(InfoNumSgps,2,"4ti output:",gr);
179
180    return List(gr, x->positivenegative(x));
181end);
182
183############################################################
184# computes a canonical basis of the kernel congruence
185# of the monoid morphism associated to the matrix m with
186# nonnegative integer coefficients wrt the term ordering
187# the kernel is the pairs (x,y) such that xm=ym
188############################################################
189InstallMethod(CanonicalBasisOfKernelCongruence,
190"Computes a canonical basis for the congruence of of the monoid morphism associated to the matrix",
191	[IsRectangularTable, IsMonomialOrdering],7,
192  function(m,ord)
193    local positivenegative, gr, nord, to;
194
195  	positivenegative:=function(p)
196  		local d1, d2;
197  		d1:=List(p, i->Maximum(i,0));
198  		d2:=List(p, i->-Minimum(0,i));
199  		return [d1,d2];
200  	end;
201
202  	if not(ForAll(m, l->ForAll(l, x->(x=0) or IsPosInt(x)))) then
203  		Error("The argument must be a matrix of nonnegative integer.");
204  	fi;
205
206
207  	# trick taken from the package Singular
208  	nord := Name( ord );
209  	nord := nord{[ 1 .. Position( nord, '(' ) - 1 ]};
210  	if nord = "MonomialLexOrdering"  then
211  			to := "lex";
212  	elif nord = "MonomialGrevlexOrdering"  then
213  			to := "grevlex";
214  	elif nord = "MonomialGrlexOrdering"  then
215  			to := "grlex";
216  	else
217  			Error( "the ordering ", ord, " is not yet supported\n" );
218  	fi;
219
220  	gr:=GroebnerBasis4ti2(TransposedMat(m),to);
221  	Info(InfoNumSgps,2,"4ti output:",gr);
222
223  	return Set(gr, x->positivenegative(x));
224  end);
225
226############################################################
227# computes the Graver basis of matrix with integer entries
228############################################################
229InstallMethod(GraverBasis,
230        "Computes the Graver basis of the matrix",
231        [IsRectangularTable],8,
232  function(a)
233    #4ti2gap implementation
234    local gr;
235
236
237    if not(IsRectangularTable(a)) then
238      Error("The argument must be a matrix.");
239    fi;
240
241    if not(IsInt(a[1][1])) then
242      Error("The entries of the matrix must be integers.");
243    fi;
244
245    Info(InfoNumSgps,2,"Using 4ti2gap for Graver.");
246
247
248    gr := GraverBasis4ti2(["mat",a]);
249    return Union(gr,-gr);
250  end);
251
252
253
254InstallOtherMethod(MinimalPresentationOfAffineSemigroup,
255        "Computes a minimimal presentation of the affine semigroup",
256        [IsAffineSemigroup],6,
257        function(a)
258    local gens, positive, gr, candidates, pres, rclass,exps, c;
259
260    positive:=function(x)
261        local p,i;
262
263        p:=[];
264        for i in [1..Length(x)] do
265            p[i]:=Maximum(x[i],0);
266        od;
267
268        return p;
269    end;
270    if not(IsAffineSemigroup(a)) then
271        Error("The argument must be an affine semigroup.");
272    fi;
273
274    gens:=MinimalGenerators(a);
275
276    gr:=GroebnerBasis4ti2(TransposedMat(gens));
277    Info(InfoNumSgps,2,"4ti output:",gr);
278
279    candidates:=Set(gr,q->positive(q));
280    candidates:=Set(candidates,c->c*gens);
281    Info(InfoNumSgps,2, "Candidates to Betti elements",candidates);
282    pres:=[];
283    for c in candidates do
284        exps:=FactorizationsVectorWRTList(c,gens);
285        rclass:=RClassesOfSetOfFactorizations(exps);
286        if Length(rclass)>1 then
287            pres:=Concatenation(pres,List([2..Length(rclass)],
288                          i->[rclass[1][1],rclass[i][1]]));
289        fi;
290    od;
291    return pres;
292end);
293
294
295
296#####################################################################
297# Computes the omega-primality of v in the affine semigroup a
298#####################################################################
299InstallOtherMethod(OmegaPrimalityOfElementInAffineSemigroup,
300        "Computes the omega-primality of v in the affine semigroup a",
301        [IsHomogeneousList,IsAffineSemigroup],6,
302        function(v,a)
303    local  ls, n, mat,extfact,par,tot,le;
304
305    le:=function(a,b)  #ordinary partial order
306    	return ForAll(b-a,x-> x>=0);
307    end;
308
309    if not(IsAffineSemigroup(a)) then
310        Error("The second argument must be an affine semigroup");
311    fi;
312
313    if not(IsListOfIntegersNS(v)) then
314        Error("The first argument must be a list of integers.");
315    fi;
316
317    if not(ForAll(v, x-> x>=0)) then
318        Error("The first argument must be a list of on nonnegative integers.");
319    fi;
320
321    ls:=MinimalGenerators(a);
322    n:=Length(ls);
323    mat:=TransposedMat(Concatenation(ls,-ls,[-v]));
324
325    if not(IsRectangularTable(mat)) then
326        Error("The first argument has not the dimension of the second.");
327    fi;
328
329    extfact:=FactorizationsVectorWRTList(v,Concatenation(ls,-ls));
330
331    par:=Set(extfact, f->f{[1..n]});
332    tot:=Filtered(par, f-> Filtered(par, g-> le(g,f))=[f]);
333    Info(InfoNumSgps,2,"Minimals of v+ls =",tot);
334    if tot=[] then
335        return 0;
336    fi;
337    return Maximum(Set(tot, Sum));
338end);
339
340#####################################################################
341# Computes the omega-primality of v in the full affine semigroup a
342#####################################################################
343InstallOtherMethod(OmegaPrimalityOfElementInAffineSemigroup,
344        "Computes the omega-primality of v in the full affine semigroup a",
345        [IsHomogeneousList,IsAffineSemigroup and HasEquations],6,
346        function(v,a)
347    local  ls, n, mat,extfact,par,tot,le;
348
349    le:=function(a,b)  #ordinary partial order
350    	return ForAll(b-a,x-> x>=0);
351    end;
352
353    if not(IsAffineSemigroup(a)) then
354        Error("The second argument must be an affine semigroup");
355    fi;
356
357    if not(IsListOfIntegersNS(v)) then
358        Error("The first argument must be a list of integers.");
359    fi;
360
361    if not(ForAll(v, x-> x>=0)) then
362        Error("The first argument must be a list of on nonnegative integers.");
363    fi;
364
365    ls:=MinimalGenerators(a);
366    n:=Length(ls);
367    mat:=TransposedMat(Concatenation(ls,-ls,[-v]));
368
369    if not(IsRectangularTable(mat)) then
370        Error("The first argument has not the dimension of the second.");
371    fi;
372
373    Info(InfoNumSgps,2,"Using 4ti2gap with full affine semigroup");
374
375    extfact:=ZSolve4ti2(["mat",TransposedMat(ls),"rel",List(v,_->1),
376                     "sign",List([1..n],_->1),"rhs",v ]);
377
378    tot:=extfact.zinhom;
379    Info(InfoNumSgps,2,"Minimals of v+ls =",tot);
380    if tot=[] then
381        return 0;
382    fi;
383    return Maximum(Set(tot, Sum));
384end);
385
386#ZSolve4ti2(["mat",TransposedMat([[2,0],[0,2],[1,2],[2,1]]),"rel",[1,1],"sign",[1,1,1,1],"rhs",[[15,15]]]);
387
388########
389# Tame degree for full affine semigroups
390########
391InstallMethod(TameDegreeOfAffineSemigroup,
392        "Computes the tame degree of the full affine semigroup a",
393        [IsAffineSemigroup and HasEquations],2,
394        function(a)
395    local ls, min, tame, gen,m,n, facts, t, minfacts;
396
397    Info(InfoNumSgps,2,"Using 4ti2gap with full affine semigroup");
398
399    ls:=MinimalGenerators(a);
400    tame:=0;
401    n:=Length(ls);
402
403    for gen in ls do
404        minfacts:=ZSolve4ti2(["mat",TransposedMat(ls),"rel",List(gen,_->1),
405                     "sign",List([1..n],_->1),"rhs",gen ]).zinhom;
406        min:=List(minfacts, x->x*ls);
407        Info(InfoNumSgps,2,"Minimal elements of ",gen,"+a=",min);
408        for m in min do
409            facts:=FactorizationsVectorWRTList(m,ls);
410            t:=TameDegreeOfSetOfFactorizations(facts);
411            if t> tame then
412                tame:=t;
413                Info(InfoNumSgps,2,"Tame degree updated to ",tame);
414            fi;
415        od;
416    od;
417    return tame;
418
419end);
420