|
#############################################################################
##
#W polynomials.gi Manuel Delgado <mdelgado@fc.up.pt>
#W Pedro A. Garcia-Sanchez <pedro@ugr.es>
##
##
#Y Copyright 2015 by Manuel Delgado and Pedro Garcia-Sanchez
#Y We adopt the copyright regulations of GAP as detailed in the
#Y copyright notice in the GAP manual.
##
#############################################################################
#################################################
##
#F NumericalSemigroupPolynomial(s,x)
## s is a numerical semigroup, and x a variable (or a value)
## returns the polynomial (1-x)\sum_{s\in S}x^s
##
##################################################
InstallGlobalFunction(NumericalSemigroupPolynomial, function(s,x)
local gaps, p;
if not IsNumericalSemigroup(s) then
Error("The first argument must be a numerical semigroup.\n");
fi;
gaps:=GapsOfNumericalSemigroup(s);
p:=List(gaps, g-> x^g);
return 1+(x-1)*Sum(p);
end);
##############################################################
##
#F IsNumericalSemigroupPolynomial(f) detects
## if there exists S a numerical semigroup such that P_S(x)=f
## f is a polynomial
##############################################################
InstallGlobalFunction(IsNumericalSemigroupPolynomial, function(f)
local x, coef, gaps, small, s, p,c;
if not(IsUnivariatePolynomial(f)) then
Error("The argument must be a univariate polynomial");
fi;
if not(IsListOfIntegersNS(CoefficientsOfUnivariatePolynomial(f))) then
return false;
fi;
if not(IsOne(LeadingCoefficient(f))) then
return false;
fi;
x:=IndeterminateOfUnivariateRationalFunction(f);
p:=(f-1)/(x-1);
if not(IsUnivariatePolynomial(p)) then
return false;
fi;
coef:=CoefficientsOfUnivariatePolynomial(p);
if Set(coef)<>Set([0,1]) then
return false;
fi;
c:=Length(coef);
gaps:=Filtered([0..c-1], i->coef[i+1]<>0);
#Print("gaps ", gaps,"\n");
small:=Difference([0..c],gaps);
#Print("small ",small,"\n");
return RepresentsSmallElementsOfNumericalSemigroup(small);
end);
##############################################################
##
#F NumericalSemigroupFromNumericalSemigroupPolynomial(f) outputs
## a numerical semigroup S such that P_S(x)=f; error if no such
## S exists
## f is a polynomial
##############################################################
InstallGlobalFunction(NumericalSemigroupFromNumericalSemigroupPolynomial, function(f)
local x, coef, gaps,small, s, p,c;
if not(IsUnivariatePolynomial(f)) then
Error("The argument must be a univariate polynomial");
fi;
if not(IsListOfIntegersNS(CoefficientsOfUnivariatePolynomial(f))) then
Error("The argument is not a numerical semigroup polynomial");
fi;
if not(IsOne(LeadingCoefficient(f))) then
Error("The argument is not a numerical semigroup polynomial");
fi;
x:=IndeterminateOfUnivariateRationalFunction(f);
p:=(f-1)/(x-1);
if not(IsUnivariatePolynomial(p)) then
Error("The argument is not a numerical semigroup polynomial");
fi;
coef:=CoefficientsOfUnivariatePolynomial(p);
if Set(coef)<>Set([0,1]) then
Error("The argument is not a numerical semigroup polynomial");
fi;
c:=Length(coef);
gaps:=Filtered([0..c-1], i->coef[i+1]<>0);
small:=Difference([0..c],gaps);
if not(RepresentsSmallElementsOfNumericalSemigroup(small)) then
Error("The argument is not a numerical semigroup polynomial");
fi;
return NumericalSemigroupBySmallElements(small);
end);
###################################################
#F HilbertSeriesOfNumericalSemigroup(s,x)
## Computes the Hilber series of s in x : \sum_{s\in S}x^s
###################################################
InstallGlobalFunction(HilbertSeriesOfNumericalSemigroup,function(s,x)
local m,ap;
if not IsNumericalSemigroup(s) then
Error("The first argument must be a numerical semigroup.\n");
fi;
if HasAperyList(s) then #uses J.Ramirez-Alfonsin trick
m:=MultiplicityOfNumericalSemigroup(s);
ap:=AperyListOfNumericalSemigroup(s);
return 1/(1-x^m)*Sum(List(ap, w->x^w));
fi;
return NumericalSemigroupPolynomial(s,x)/(1-x);
end);
######################################################
##
#F Computes the Graeffe polynomial of p
## see for instance [BD-cyclotomic]
##
######################################################
InstallGlobalFunction(GraeffePolynomial,function(p)
local coef, h, g, i, f1,x;
if not(IsUnivariatePolynomial(p)) then
Error("The argument must be a univariate polynomial");
fi;
x:=IndeterminateOfUnivariateRationalFunction(p);
coef:=CoefficientsOfUnivariatePolynomial(p);
h:=[]; g:=[];
for i in [1..Length(coef)] do
if (i-1) mod 2 = 0 then
Add(g,coef[i]);
else
Add(h,coef[i]);
fi;
od;
#Print(g," ,",h,"\n");
g:=UnivariatePolynomial(Rationals,g);
h:=UnivariatePolynomial(Rationals,h);
f1:=Value(g,x)^2-x*Value(h,x)^2;
return f1/LeadingCoefficient(f1);
end);
#####################################################
##
#F IsCyclotomicPolynomial(f) detects
## if f is a cyclotomic polynomial using the method explained in
## BD-cyclotomic
#####################################################
InstallGlobalFunction(IsCyclotomicPolynomial,function(f)
local f1, x, f2, mf;
if not(IsUnivariatePolynomial(f)) then
Error("The argument must be a univariate polynomial");
fi;
if not(IsListOfIntegersNS(CoefficientsOfUnivariatePolynomial(f))) then
Error("The polynomial has not integer coefficients");
fi;
if not(IsOne(LeadingCoefficient(f))) then
return false;
fi;
if not(IsIrreducible(f)) then
return false;
fi;
x:=IndeterminateOfUnivariateRationalFunction(f);
f1:=GraeffePolynomial(f);
if (f=f1) then
return true;
fi;
mf:=Value(f,-x);# Print(f1, ", ",f, ",", mf,"\n");
if ((f1=mf) or (f1=-mf)) and IsCyclotomicPolynomial(mf/LeadingCoefficient(mf)) then
return true;
fi;
f2:=Set(Factors(f1))[1];
if f1=f2^2 and IsCyclotomicPolynomial(f2/LeadingCoefficient(f2)) then
return true;
fi;
return false;
end);
########################################################################
##
#F IsKroneckerPolynomial(f) decides if
## f is a Kronecker polynomial, that is, a monic polynomial with integer coefficients
## having all its roots in the unit circunference, equivalently, is a product of
## cyclotomic polynomials
## New implementation with A. Herrera-Poyatos that does not need factoring
#########################################################################
InstallGlobalFunction(IsKroneckerPolynomial,function(f)
local x, sf, f1, fs, fn, fp, factors_graeffe;
if not(IsUnivariatePolynomial(f)) then
Error("The argument must be a polynomial in one variable");
fi;
if IsZero(f) then
return false;
fi;
x:=IndeterminateOfLaurentPolynomial(f);
# Take the largest square free divisor.
sf:=f/Gcd(f,Derivative(f));
# Remove the factors x, x+1 and x-1.
if Value(sf, 0) = 0 then
sf := sf / x;
fi;
if Value(sf, 1) = 0 then
sf := sf/(x-1);
fi;
if Value(sf, -1) = 0 then
sf := sf/(x+1);
fi;
# Check if the polynomial is a constant.
if Degree(sf) = 0 then
return true;
fi;
# Check if the polynomial has even degree and is self reciprocal.
if Degree(sf) mod 2 <> 0 or not(IsSelfReciprocalUnivariatePolynomial(sf)) then
return false;
fi;
# Apply the graeffe method to sf. Note that if sf(+-x) = f1, then
# sf is Kronecker.
f1:= GraeffePolynomial(sf);
if sf = f1 then
return true;
fi;
if Value(sf,-x) = f1 then
return true;
fi;
# Assume that sf is Kronecker.
# Find the descomposition of sf = fs*fp*fn, where:
# - fs is the part of sf which verifies Graeffe(fs) = (g(x))^2
# - fp is the part of sf which verifies Graeffe(fs) = fs.
# - fn is the part of sf which vefifies Graeffe(fs)(x) = fs(-x).
fs := Value(Gcd(Derivative(f1), f1), x^2);
fp := Gcd(sf / fs, f1);
fn := sf / (fs * fp);
factors_graeffe := Difference([fs, fp, fn], [1+0*x]);
if Length(factors_graeffe) = 1 then
if IsOne(fs) then
# We must have f = fp or f = fs, but we obtained Graeffe(sf) != sf,
# Graeffe(sf) != sf(-x), a contradiction. Hence sf is not Kronecker.
return false;
else
# sf is Kronecker if and only if f1 / fs(\sqrt{x}) is Kronecker.
return IsKroneckerPolynomial(f1 / Gcd(f1,Derivative(f1)));
fi;
else
# sf is Kronecker if and only if fs, fp and fn are Kronecker.
return ForAll(factors_graeffe, IsKroneckerPolynomial);
fi;
end);
###########################################
##
#P IsCyclotomicNumericalSemigroup(s)
## Checks if the polynomial fo s is Kronecker
###########################################
InstallMethod(IsCyclotomicNumericalSemigroup,
"Detects if the polinomial semigroup of the semigroup is Kronecker",
[IsNumericalSemigroup],
function(s)
local x, f;
x:=X(Rationals,"x");
f:=NumericalSemigroupPolynomial(s,x);
return IsKroneckerPolynomial(f);# ForAll(Factors(f),IsCyclotomicPolynomial);
end);
InstallTrueMethod(IsCyclotomicNumericalSemigroup, IsACompleteIntersectionNumericalSemigroup);
#####################################################
##
#F IsSelfReciprocalUnivariatePolynomial(p)
## Checks if the univariate polynomial p is selfreciprocal
## New implementation by A. Herrera-Poyatos
#####################################################
InstallGlobalFunction(IsSelfReciprocalUnivariatePolynomial,function(p)
local cf;
if not(IsUnivariatePolynomial(p)) then
Error("The argument must be a polynomial\n");
fi;
# Check if the polynomial is self-reciprocal.
cf:=CoefficientsOfUnivariatePolynomial(p);
return cf=Reversed(cf);
end);
#################################################################
##
# F SemigroupOfValuesOfPlaneCurveWithSinglePlaceAtInfinity(f)
## Computes the semigroup of values {mult(f,g) | g curve} of a plane curve
## with one place at the infinity in the variables X(Rationals,1) and X(Rationals,2)
## f must be monic on X(Rationals(2))
## SemigroupOfValuesOfPlaneCurveWithSinglePlaceAtInfinity(f,"all")
## The same as above, but the output are the approximate roots and
## delta-sequence
##
#################################################################
InstallGlobalFunction(SemigroupOfValuesOfPlaneCurveWithSinglePlaceAtInfinity, function(arg)
local monomials, degree, degreepol, polexprc, polexprs, polexpr, app, tt1, sv;
#returns the set of monomials of a polynomial
monomials:=function(f)
local term, out,temp;
out:=[];
if IsZero(f) then
return ([0]);
fi;
temp:=f;
while not(IsZero(temp)) do
term:=LeadingTermOfPolynomial(temp,MonomialLexOrdering());
Add(out,term);
temp:=temp-term;
od;
return out;
end;
#returns the total weighted degree of a monomial wrt l
degree:=function(mon,l)
local n;
n:=Length(l);
if IsRat(mon) then
return 0;
fi;
return Sum(List([1..n], i->DegreeIndeterminate(mon,i)*l[i]));
end;
#returns the total weighted degree of a polynomial wrt l
degreepol:=function(pol,l)
return Maximum(List(monomials(pol), t->degree(t,l)));
end;
#finds an expression of f in terms of p, and outputs it as a polynomial in x, var is the list of variables in f and p
polexpr:=function(f,p,var,x)
local rem, quo, pr;
pr:=PolynomialDivisionAlgorithm(f,[p],MonomialLexOrdering(var));
rem:=[pr[1]];
quo:=pr[2][1];
while not(IsZero(quo)) do
pr:=PolynomialDivisionAlgorithm(quo,[p],MonomialLexOrdering(var));
Add(rem,pr[1]);
quo:=pr[2][1];
od;
return Sum(List([1..Length(rem)],i->rem[i]*x^(i-1)));
end;
#uses the previous function to apply recurrently the procedure for all polynomials in ps and variables in xs
polexprs:=function(f,ps,var,xs)
local i, len, out;
len:=Length(ps);
out:=f;
for i in [1..len] do
out:=polexpr(out,ps[i],var,xs[i]);
od;
return out;
end;
#expression of f in terms of p; both in variables x and y...
polexprc:=function(f,p)
local rem, quo, pr;
pr:=PolynomialDivisionAlgorithm(f,[p],MonomialLexOrdering([X(Rationals,2),X(Rationals,1)]));
rem:=[pr[1]];
quo:=Value(pr[2][1],[],[]);
while not(IsZero(quo)) do
pr:=PolynomialDivisionAlgorithm(quo,[p],MonomialLexOrdering([X(Rationals,2),X(Rationals,1)]));
Add(rem,pr[1]);
quo:=pr[2][1];
od;
return Reversed(rem);
end;
##TT1
tt1:=function(f,g,y)
local d;
d:=DegreeIndeterminate(f,y)/DegreeIndeterminate(g,y);
if d>0 then
return g+polexprc(f,g)[2]/d;
fi;
return f;
end;
#approximate root
app:=function(f,y,d)
local G,S,F,P,e,t,coef;
e:=DegreeIndeterminate(f,y)/d;
G:=tt1(f,y^e,y);
S:=polexprc(f,G);
P:=S[2];
#tschirnhausen transform
while not(IsZero(P)) do
F:=G+P/e;
G:=tt1(f,F,y);
S:=polexprc(f,G);
P:=S[2];
od;
return G;
end;
#semigroup of values
sv:=function(f)
local m, d,dd, n,g,it, max, lt, e, coef, deg, tsti,rd, var, id, R, aroots;
if not(IsPolynomial(f)) then
Error("The argument must be a polynomial\n");
fi;
R:=PolynomialRing(Rationals,[X(Rationals,1),X(Rationals,2)]);
if not(f in R) then
Error("The argument must be a polynomial in ", R,"\n");
fi;
coef:=PolynomialCoefficientsOfPolynomial(f,X(Rationals,1));
if not(IsConstantRationalFunction(coef[Length(coef)])) then
Error("The polynomial does not have a single place at infinity or the leading coefficient in ", X(Rationals,1)," is not a rational number\n");
fi;
coef:=PolynomialCoefficientsOfPolynomial(f,X(Rationals,2));
if not(IsConstantRationalFunction(coef[Length(coef)])) then
Error("The polynomial does not have a single place at infinity or or the leading coefficient in ", X(Rationals,2)," is not a rational number\n");
fi;
f:=f/coef[Length(coef)];
m:=[];
coef:=PolynomialCoefficientsOfPolynomial(f,X(Rationals,2));
Info(InfoNumSgps,2,"f ",f);
Info(InfoNumSgps,2,"f ",coef);
it:=coef[1];
if IsZero(it) then
Error("The polynomial is not irreducible\n");
fi;
m[1]:=DegreeIndeterminate(f,X(Rationals,2));
m[2]:=DegreeIndeterminate(PolynomialCoefficientsOfPolynomial(f,X(Rationals,2))[1],X(Rationals,1));
n:=2;
d:=Gcd(m);
e:=m[1]/d;
g:=f;
var:=[X(Rationals,2),X(Rationals,1)];
aroots:=[X(Rationals,2)];
while d>1 do
Add(var,X(Rationals,n+1));
rd:=app(f,X(Rationals,2),d);
Add(aroots,rd);
g:=polexprs(f,Reversed(aroots),var,Reversed(List([2..n+1],i->X(Rationals,i))));
coef:=PolynomialCoefficientsOfPolynomial(g, X(Rationals,n+1));
Info(InfoNumSgps,2,"We obtain coefficients: ",coef);
it:=coef[1];
Info(InfoNumSgps,2,"Independent term ", it);
max:=degreepol(it,m/Gcd(m));
if (max=0) then #or (d=Gcd(d,max)) or (max in m) then
Error("The polynomial is not irreducible or it has not a single place at infinity (deg 0)\n");
fi;
#irreducibility test
tsti:=First(coef{[2..Length(coef)]}, t->degreepol(t,m/d)>=max);
if tsti<>fail then
Info(InfoNumSgps,1,"Term ", tsti," produces error");
Error("The polynomial is not irreducible", "\n");
fi;
dd:=Gcd(d,max);
Info(InfoNumSgps,2,"New candidate: ",max);
if (d=Gcd(d,max)) or (max in m) then
Error("The polynomial is not irreducible or it has not a single place at infinity\n");
fi;
m[n+1]:=max;
if(m[n]*Gcd(m{[1..n-1]})<=max*d)then
Error("The polynomial is not irreducible (not subdescending sequence)", m,"\n");
fi;
e:=d/Gcd(d,max);
d:=dd;
n:=n+1;
od;
return [m,aroots];
end;
if Length(arg)=1 then
return NumericalSemigroup(sv(arg[1])[1]);
fi;
if Length(arg)=2 and arg[2] = "all" then
return sv(arg[1]);
fi;
Error("Wrong number of arguments\n");
end);
#########################################################################
##
#F IsDeltaSequence(l)
## tests whether or not l is a \delta-sequence (see for instancd [AGS14])
#########################################################################
InstallGlobalFunction(IsDeltaSequence,function(l)
local d,e,h;
if not(IsListOfIntegersNS(l)) then
Error("The argument must be a list of positive integers\n");
fi;
if Gcd(l)<>1 then
Info(InfoNumSgps,2,"The gcd of the list is not one");
return false;
fi;
h:=Length(l)-1;
d:=List([1..h+1], i->Gcd(l{[1..i]}));
e:=List([1..h], i->d[i]/d[i+1]);
if Length(l)=1 then
return true;
fi;
if (l[1]<=l[2]) or (d[2]=l[2]) then
Info(InfoNumSgps,2,"Either the first element is smaller than or equal to the second, or the gcd of the first to equals the second");
return false;
fi;
if Length(Set(d))<h+1 then
Info(InfoNumSgps,2,"The list of gcds is not strctitly decreasing");
return false;
fi;
return ForAll([1..h], i->l[i+1]/Gcd(l{[1..i+1]}) in NumericalSemigroup(l{[1..i]}/Gcd(l{[1..i]})))
and
ForAll([2..h], i->l[i]*Gcd(l{[1..i-1]}) >l[i+1]*Gcd(l{[1..i]}));
end);
#########################################################
##
#F DeltaSequencesWithFrobeniusNumber(f)
## Computes the list of delta-sequences with Frobenius number f
##
#########################################################
InstallGlobalFunction(DeltaSequencesWithFrobeniusNumber,function(f)
local abs, test;
if not(IsInt(f)) then
Error("The argument must be an integer\n");
fi;
if (f<-1) or IsEvenInt(f) then
return [];
fi;
if f=-1 then
return [[1]];
fi;
abs:=function(f)
local dkcand, dk, rk, fp, candr, bound, total;
if (f<-1) or (f mod 2=0) then
return [];
fi;
if (f=-1) then
return [[1]];
fi;
if (f=1) then
return [[2,3],[3,2]];
fi;
total:=[[f+2,2],[2,f+2]];
for rk in Reversed([2..f-1]) do
bound:=(Int((f+1)/(rk-1)+1));
dkcand:=Filtered([2..bound],d->(Gcd(d,rk)=1)and((f+rk) mod d=0));
for dk in dkcand do
fp:=(f+rk*(1-dk))/dk;
candr:=Filtered(abs(fp), l-> rk in NumericalSemigroup(l));
candr:=List(candr, l-> Concatenation(l*dk, [rk]));
candr:=Filtered(candr, test);
candr:=Filtered(candr,l->l[1]>l[2]);
total:=Union(total,candr);
od;
od;
return total;
end;
test:=function(l)
local len;
len:=Length(l);
return ForAll([3..len], k-> l[k-1]*Gcd(l{[1..k-2]})>l[k]*Gcd(l{[1..k-1]}));
end;
return Difference(abs(f),[[2,f+2]]);
end);
#####################################################
##
#F CurveAssociatedToDeltaSequence(l)
## computes the curve associated to a delta-sequence l in
## the variables X(Rationals,1) and X(Rationals,2)
## as explained in [AGS14]
##
#####################################################
InstallGlobalFunction(CurveAssociatedToDeltaSequence,function(l)
local n, d, f, fact,facts, g, e, k, ll, len, dd, x,y;
if not(IsDeltaSequence(l)) then
Error("The sequence is not valid\n");
fi;
x:=X(Rationals,1);
y:=X(Rationals,2);
len:=Length(l);
if len<2 then
Error("The sequence must have at least two elements\n");
fi;
if Gcd(l)<>1 then
Error("The sequence must have gcd equal to one\n");
fi;
if len=2 then
return y^(l[1])-x^(l[2]);
fi;
e:=List([1..len-1],k->Gcd(l{[1..k]})/Gcd(l{[1..k+1]}));
g:=[]; g[1]:=x; g[2]:=y;
for k in [2..len] do
d:=Gcd(l{[1..k-1]});
dd:=Gcd(l{[1..k]});
ll:=l{[1..k-1]}/d;
fact:=First(FactorizationsIntegerWRTList(l[k]/dd, ll),ff->ForAll([2..k-1], i->ff[i]<e[i-1]));
if fact=fail then
Error("The sequence is not valid\n");
fi;
g[k+1]:=g[k]^e[k-1]-Product(List([1..k-1], i->g[i]^fact[i]));
Info(InfoNumSgps,2,"Temporal curve: ",g[k+1]);
od;
return g[Length(g)];
end);
#################################################################
##
#F SemigroupOfValuesOfCurve_Local(arg)
## Computes the semigroup of values of R=K[pols],
## that is, the set of possible order of series in this ring
## pols is a set of polynomials in a single variable
## The semigroup of values is a numerical semigroup if l(K[[x]]/R) is finite
## If this length is not finite, the output is fail
## If the second argument "basis" is given, then the output is a basis B of
## R such that o(B) minimally generates o(R), and it is reduced
## If the second argument is an integer, then the output is a polynomial f in R
## with o(f) that value (if there is none, then the output is fail)
## Implementation based in [AGSM14]
###########################################################
InstallGlobalFunction(SemigroupOfValuesOfCurve_Local, function(arg)
local G, F, n, T, max, s, d, newG, c, kernel,var, tomoniclocal, order, subduction, t,tt, TT, narg,val,facts,pols, msg,reduce;
### the order of the series (polynomial)
order:=function(p)
local cl;
if IsInt(p) or IsZero(p) then
return 0;
fi;
cl:=CoefficientsOfUnivariatePolynomial(p);
return First([1..Length(cl)], i->cl[i]<>0)-1;
end;
#### transforms to monic
tomoniclocal:=function(p)
local cl;
if IsZero(p) then
return 0*X(Rationals,var);
fi;
if IsConstantRationalFunction(p) or IsRat(p) then
return 1;
fi;
cl:=CoefficientsOfUnivariatePolynomial(p);
return p/First(cl,i->i<>0);
end;
#### subduction: tries to express p as a polynomial in pols
#### pols are supposed to be monic
subduction:=function(p)
local initial, facts, ord;
ord:=order(p);
if d<>1 and (ord mod d)<>0 then
return p;
fi;
if order(p)=0 then
return 0*X(Rationals,var);
fi;
if order(p)>c then
if d=1 then
return 0*X(Rationals,var);
else
return p;
fi;
fi;
initial:=List(F,order);
facts:=FactorizationsIntegerWRTList(order(p), initial);
if facts=[] then
return p;
fi;
Info(InfoNumSgps,2,"Reducing ",p, " to ",
tomoniclocal(p-Product(List([1..Length(F)],i->F[i]^(facts[1][i])))));
return subduction(tomoniclocal(p-Product(List([1..Length(F)],i->F[i]^(facts[1][i])))));
end;
#### reduces the elements in the minimal basis F, p is monic
reduce:=function(p)
local ords, cfs, cf, fact, reduced, head, tail;
head:=X(Rationals,var)^order(p);
tail:=p-head;
ords:=List(F,order);
repeat
reduced:=false;
#Print(p," ");
cfs:=CoefficientsOfUnivariatePolynomial(tail);
cf:=First([1..Length(cfs)], i->(cfs[i]<>0) and ((i-1) in s));
if cf<>fail then
fact:=FactorizationsIntegerWRTList(cf-1,ords)[1];
tail:=(tail-cfs[cf]*Product(List([1..Length(F)],i->F[i]^(fact[i])))) mod X(Rationals,var)^c;
cfs:=CoefficientsOfUnivariatePolynomial(tail);
else
reduced:=true;
fi;
until reduced;
return head+tail;
end;
#### computes the relations among the polynomials in pols
kernel:=function(pols)
local p, msg, ed, mp, minp, eval, candidates, c, pres, rclass, exps;
eval:=function(pair)
local m1,m2;
m1:=tomoniclocal(Product(List([1..ed],i-> pols[i]^pair[1][i])));
m2:=tomoniclocal(Product(List([1..ed],i-> pols[i]^pair[2][i])));
return tomoniclocal(-m1+m2);
end;
msg:=List(pols,order);
ed:=Length(msg);
if ed=0 then
return [];
fi;
minp:=GeneratorsOfKernelCongruence(List(msg,x->[x]));
Info(InfoNumSgps,2,"The exponents of the binomials of the kernel are ",minp);
# Contrary to the global case, here it is faster not to compute a minimal presentation
# candidates:=Set(minp, x->x[1]*msg);
# pres:=[];
# for c in candidates do
# exps:=FactorizationsIntegerWRTList(c,msg);
# rclass:=RClassesOfSetOfFactorizations(exps);
# if Length(rclass)>1 then
# pres:=Concatenation(pres,List([2..Length(rclass)],
# i->[rclass[1][1],rclass[i][1]]));
# fi;
# od;
return List(minp,p->eval(p));
end; #end of kernel
narg:=Length(arg);
if narg=2 then
pols:=arg[1];
val:=arg[2];
if not(IsInt(val)) and not(val="basis") then
Error("The second argument must be an integer or the string 'basis'");
fi;
fi;
if narg=1 then
pols:=arg[1];
fi;
if narg>2 then
Error("Wrong number of arguments (two or one)");
fi;
if not(IsHomogeneousList(pols)) then
Error("The argument must be a list of polynomials.");
fi;
if not(ForAll(pols, IsUnivariatePolynomial)) then
Error("The argument must be a list of polynomials.");
fi;
if Length(Set(pols, IndeterminateOfLaurentPolynomial))<>1 then
Error("The arguments must be polynomials in the same variable; constants not allowed nor the empty list.");
fi;
var:=IndeterminateNumberOfLaurentPolynomial(pols[1]);
F:=ShallowCopy(pols);
Sort(F,function(a,b) return order(a)< order(b); end);
F:=List(F,tomoniclocal);
G:=List(F,order);
n:=0;
while true do
d:=Gcd(G);
s:=NumericalSemigroup(G/d);
c:=d*ConductorOfNumericalSemigroup(s);
T:=kernel(F);
T:=Filtered(T, x->not(IsZero(x)));
Info(InfoNumSgps,2,"The kernel evaluated in the polynomials is ",T);
T:=Set(T,subduction);
T:=Filtered(T, x->not(IsZero(x)));
Info(InfoNumSgps,2,"After subduction: ",T);
if Gcd(G) = 1 then
T:=Filtered(T, t->order(t)<c);
fi;
if T=[] or F=Union(F,T) then
d:=Gcd(G);
if d=1 then
s:=NumericalSemigroup(G);
if narg=1 then
return s;
fi;
msg:=MinimalGeneratingSystem(s);
F:=Filtered(F,f->order(f) in msg);
if val="basis" then
return List(F,reduce);
fi;
if IsInt(val) then
Info(InfoNumSgps,2,"The generators of the algebra are ",F);
facts:=FactorizationsIntegerWRTList(val,List(F,order));
if facts=[] then
return fail;
fi;
return reduce(Product(List([1..Length(facts[1])],i->F[i]^facts[1][i])));
fi;
fi;
Info(InfoNumSgps,2,"The monoid is not a numerical semigroup and it is generated by ",
G);
#d*MinimalGeneratingSystem(NumericalSemigroup(G/d)));
#return fail;
fi;
Info(InfoNumSgps,2,"Adding ",T," to my list of polynomials");
F:=Union(F,T);
newG:=Set(T,order);
G:=Union(G,newG);
Info(InfoNumSgps,2,"The set of possible values uptates to ",G);
n:=n+1;
d:=Gcd(G);
s:=NumericalSemigroup(G/d);
Info(InfoNumSgps,2,"Small elements ",d*SmallElements(s));
od;
return fail;
end);
#################################################################
##
#F SemigroupOfValuesOfCurve_Global(arg)
## Computes the semigroup of values of R=K[pols],
## that is, the set of possible degrees of polynomials in this ring
## pols is a set of polynomials in a single variable
## The semigroup of values is a numerical semigroup if l(K[x]/R) is finite
## If this length is not finite, the output is fail
## If the second argument "basis" is given, then the output is a basis B of
## R such that deg(B) minimally generates deg(R), and it is reduced
## If the second argument is an integer, then the output is a polynomial f in R
## with deg(f) that value (if there is none, then the output is fail)
## Implementation based in [AGSM14]
###########################################################
InstallGlobalFunction(SemigroupOfValuesOfCurve_Global, function(arg)
local G, F, n, T, max, s, d, newG, c, kernel,var, tomonicglobal, degree, subduction, pols, val, narg, degs, facts, reduce, msg;
### the degree of the polynomial
degree:=function(p)
return DegreeOfUnivariateLaurentPolynomial(p);
end;
#### transforms to monic
tomonicglobal:=function(p)
if IsZero(p) then
return 0*X(Rationals,1);
fi;
return p/LeadingCoefficient(p);
end;
#### subduction: tries to express p as a polynomial in pols
#### pols are supposed to be monic
subduction:=function(p)
local initial, facts;
if IsZero(p) then
return p;
fi;
if degree(p)=0 then
return 0*X(Rationals,1);
fi;
initial:=List(F,degree);
facts:=FactorizationsIntegerWRTList(degree(p), initial);
if facts=[] then
return p;
fi;
Info(InfoNumSgps,2,"Reducing ",p, " to ",
tomonicglobal(p-Product(List([1..Length(F)],i->F[i]^(facts[1][i])))));
return subduction(tomonicglobal(
p-Product(List([1..Length(F)],i->F[i]^(facts[1][i])))));
end;
#### reduces the elements in the minimal basis F, p is monic
reduce:=function(p)
local degs, cfs, cf, fact, reduced, head, tail;
head:=X(Rationals,var)^degree(p);
tail:=p-head;
degs:=List(F,degree);
repeat
reduced:=false;
#Print(p," ");
cfs:=CoefficientsOfUnivariatePolynomial(tail);
cf:=First([1..Length(cfs)], i->(cfs[i]<>0) and ((i-1) in s));
#Error("blabla");
#Print("cf",cf,"\n");
if cf<>fail then
fact:=FactorizationsIntegerWRTList(cf-1,degs)[1];
tail:=(tail-cfs[cf]*Product(List([1..Length(F)],i->F[i]^(fact[i]))));
cfs:=CoefficientsOfUnivariatePolynomial(tail);
else
reduced:=true;
fi;
until reduced;
return head+tail;
end;
#### computes the relations among the polynomials in pols
kernel:=function(pols)
local p, msg, ed, mp, minp, eval, candidates, c, pres, rclass, exps;
eval:=function(pair)
local m1,m2;
m1:=tomonicglobal(Product(List([1..ed],i-> pols[i]^pair[1][i])));
m2:=tomonicglobal(Product(List([1..ed],i-> pols[i]^pair[2][i])));
return tomonicglobal(-m1+m2);
end;
msg:=List(pols,degree);
ed:=Length(msg);
if ed=0 then
return [];
fi;
minp:=GeneratorsOfKernelCongruence(List(msg,x->[x]));
Info(InfoNumSgps,2,"The exponents of the binomials of the kernel are ",minp);
candidates:=Set(minp, x->x[1]*msg);
pres:=[];
for c in candidates do
exps:=FactorizationsIntegerWRTList(c,msg);
rclass:=RClassesOfSetOfFactorizations(exps);
if Length(rclass)>1 then
pres:=Concatenation(pres,List([2..Length(rclass)],
i->[rclass[1][1],rclass[i][1]]));
fi;
od;
return List(pres,p->tomonicglobal(eval(p)));
end; #end of kernel
narg:=Length(arg);
if narg>2 then
Error("Wrong number of arguments (two or one)");
fi;
if narg=1 then
pols:=arg[1];
fi;
if narg=2 then
pols:=arg[1];
val:=arg[2];
if not(IsInt(val)) and not(val="basis") then
Error("The second argument must be an integer or the string 'basis'");
fi;
fi;
if not(IsHomogeneousList(pols)) then
Error("The argument must be a list of polynomials.");
fi;
if not(ForAll(pols, IsUnivariatePolynomial)) then
Error("The argument must be a list of polynomials.");
fi;
if Length(Set(pols, IndeterminateOfLaurentPolynomial))<>1 then
Error("The arguments must be polynomials in the same variable; constants not allowed nor the empty list.");
fi;
var:=IndeterminateNumberOfLaurentPolynomial(pols[1]);
F:=ShallowCopy(pols);
Sort(F,function(a,b) return degree(a)< degree(b); end);
F:=List(F,tomonicglobal);
G:=List(F,degree);
n:=0;
while true do
T:=kernel(F);
T:=Filtered(T, x->not(IsZero(x)));
Info(InfoNumSgps,2,"The kernel evaluated in the polynomials is ",T);
T:=Set(T,subduction);
T:=Filtered(T, x->not(IsZero(x)));
Info(InfoNumSgps,2,"After subduction: ",T);
if Gcd(G) = 1 then
s:=NumericalSemigroup(G);
c:=ConductorOfNumericalSemigroup(s);
T:=Filtered(T, t->degree(t)<c);
fi;
if T=[] or F=Union(F,T) then
d:=Gcd(G);
if d=1 then
s:=NumericalSemigroup(G);
if narg=1 then
return s;
fi;
msg:=MinimalGeneratingSystem(s);
F:=Filtered(F,f->degree(f) in msg);
if val="basis" then
return List(F,reduce);
fi;
if IsInt(val) then
Info(InfoNumSgps,2,"The generators of the algebra are ",F);
facts:=FactorizationsIntegerWRTList(val,List(F,degree));
if facts=[] then
return fail;
fi;
return reduce(Product(List([1..Length(facts[1])],i->F[i]^facts[1][i])));
fi;
fi;
Info(InfoNumSgps,2,"The monoid is not a numerical semigroup and it is generated by ",
G);
#d*MinimalGeneratingSystem(NumericalSemigroup(G/d)));
return fail;
fi;
Info(InfoNumSgps,2,"Adding ",T," to my list of polynomials");
F:=Union(F,T);
newG:=Set(T,degree);
G:=Union(G,newG);
if narg=2 and val in G then
return First(F,f->degree(f)=val);
fi;
Info(InfoNumSgps,2,"The set of possible values uptates to ",G);
n:=n+1;
d:=Gcd(G);
s:=NumericalSemigroup(G/d);
Info(InfoNumSgps,2,"Small elements ",d*SmallElements(s));
od;
return fail;
end);
##################################################################
##
#F GeneratorsModule_Global(A,M)
##
## A and M are lists of polynomials in the same variable
## Computes a basis of the ideal MK[A], that is, a set F such that
## deg(F) generates the ideal deg(MK[A]) of deg(K[A]), where deg
## stands for degree
##################################################################
InstallGlobalFunction(GeneratorsModule_Global, function(Al,M)
local S, gens, gM, a, b, da, db, i, j, rs, rd, rel, fcta, fctb, C, pair, reduction, n, reduce, A, R, t;
R:=function(a,b,s)
local i, mg;
i:=IntersectionIdealsOfNumericalSemigroup(a+s,b+s);
mg:=MinimalGenerators(i);
return List(mg, m->[m-a,m-b]);
end;
# improve to reducing all terms, not only the leading term
reduce:=function(A,M,f)
local gens,geni,cand,d, fact, c, r, s,a, ds, cs, csa;
gens:=List(A, DegreeOfLaurentPolynomial);
s:=NumericalSemigroup(gens);
geni:=List(M,DegreeOfLaurentPolynomial);
if IsZero(f) then
return f;
fi;
r:=f;
cs:=CoefficientsOfUnivariatePolynomial(r);
ds:=Filtered([1..Length(cs)],i->not(IsZero(cs[i])))-1;
c:=First([1..Length(geni)], i->ForAny(ds, d->d-geni[i] in s));
if c<>fail then
d:=First(ds,x->x-geni[c] in s);
fi;
while c<>fail do
fact:=FactorizationsIntegerWRTList(d-geni[c],gens);
a:=M[c]*Product(List([1..Length(gens)],i->A[i]^fact[1][i]));
csa:=CoefficientsOfUnivariatePolynomial(a);
r:=csa[Length(csa)]*r-cs[d+1]*a;#r-cs[d+1]*a/csa[Length(csa)];
#Info(InfoNumSgps,2,"New reduction ",r," degree ",d," coeff ",cs[d+1]);
if IsZero(r) then
return r;
fi;
cs:=CoefficientsOfUnivariatePolynomial(r);
ds:=Filtered([1..Length(cs)],i->not(IsZero(cs[i])))-1;
c:=First([1..Length(geni)], i->ForAny(ds, d->d-geni[i] in s));
if c<>fail then
d:=First(ds,x->x-geni[c] in s);
fi;
od;
return r/cs[Length(cs)];
end;
if not(IsHomogeneousList(Al)) then
Error("The first argument must be a list of polynomials.");
fi;
if not(IsHomogeneousList(M)) then
Error("The second argument must be a list of polynomials.");
fi;
if not(ForAll(Union(Al,M), IsUnivariatePolynomial)) then
Error("The arguments must be lists of polynomials.");
fi;
if Length(Set(Union(Al,M), IndeterminateOfLaurentPolynomial))<>1 then
Error("The arguments must be lists of polynomials in the same variable; constants not allowed nor the empty list.");
fi;
t:=IndeterminateNumberOfLaurentPolynomial(Al[1]);
A:=SemigroupOfValuesOfCurve_Global(Al,"basis");#List(A, DegreeOfLaurentPolynomial);
gens:=List(A, DegreeOfLaurentPolynomial);
S:=NumericalSemigroup(gens);
n:=Length(A);
gM:=ShallowCopy(M);
C:=[];
for i in [1..Length(gM)] do
for j in [i+1..Length(gM)] do
Add(C,[gM[i],gM[j]]);
od;
od;
while C<>[] do
pair:=Remove(C,1);
a:=pair[1];
b:=pair[2];
da:=DegreeOfLaurentPolynomial(a);
db:=DegreeOfLaurentPolynomial(b);
rs:=R(da,db,S);
reduction:=true;
for rel in rs do
fcta:=FactorizationsIntegerWRTList(rel[1],gens)[1];
fctb:=FactorizationsIntegerWRTList(rel[2],gens)[1];
#rd:=reduce(A,gM,a*Product(List(fcta, e->t^e))-b*Product(List(fctb, e->t^e)));
rd:=reduce(A,gM,a*Product(List([1..n], i->A[i]^fcta[i]))-b*Product(List([1..n], i->A[i]^fctb[i])));
if not(IsZero(rd)) then
Info(InfoNumSgps,2,"new generator ",rd," of degreee ",DegreeIndeterminate(rd,t));
C:=Union(C,List(gM, x->[x,rd]));
Add(gM,rd);
reduction:=false;
fi;
od;
while not reduction do
#Info(InfoNumSgps,2,"Reducing...");
reduction:=true;
a:=First(gM, x->x<>reduce(A,Difference(gM,[x]),x));
if a<>fail then
rd:=reduce(A,Difference(gM,[a]),a);
if IsZero(rd) then
gM:=Difference(gM,[a]);
Info(InfoNumSgps,2,"Removing redundant generator: ",a);
else
gM:=Union(Difference(gM,[a]),[rd]);
Info(InfoNumSgps,2,"Replacing ",a," by ",rd);
fi;
reduction:=false;
fi;
od;
od;
Info(InfoNumSgps,2,"Reducing...");
reduction:=false;
while not reduction do
reduction:=true;
a:=First(gM, x->x<>reduce(A,Difference(gM,[x]),x));
if a<>fail then
rd:=reduce(A,Difference(gM,[a]),a);
if IsZero(rd) then
gM:=Difference(gM,[a]);
Info(InfoNumSgps,2,"Removing redundant:",a);
else
gM:=Union(Difference(gM,[a]),[rd]);
Info(InfoNumSgps,2,"Replacing ",a," by ",rd);
fi;
reduction:=false;
fi;
od;
return gM;
end);
##################################################################
##
#F GeneratorsKahlerDifferentials(A)
##
## A synonym for GeneratorsModule_Global(A,M), with M the set of
## derivatives of the elements in A
##################################################################
InstallGlobalFunction(GeneratorsKahlerDifferentials, function(A)
local M, t;
if not(IsHomogeneousList(A)) then
Error("The argument must be a list of polynomials.");
fi;
if not(ForAll(A, IsUnivariatePolynomial)) then
Error("The argument must be list of polynomials.");
fi;
if Length(Set(A, IndeterminateOfLaurentPolynomial))<>1 then
Error("The argument must be a list of polynomials in the same variable; constants not allowed nor the empty list.");
fi;
t:=IndeterminateNumberOfLaurentPolynomial(A[1]);
M:=List(A,p->Derivative(p,t));
return GeneratorsModule_Global(A,M);
end);
##################################################################
##
#O CyclotomicExponentSequence(s,k)
##
## s is a numerical semigroup and k a positive integer, the
## output is the list of the first k elements of the cyclotomic
## exponent sequence of s (see [C-GS-M]).
## The sequence will be truncated if the semigroup is cyclotomic
## and k is bigger than the last nonzero element in its sequence.
##################################################################
InstallMethod(CyclotomicExponentSequence,[IsNumericalSemigroup, IsPosInt],
function(s,k)
local invs, i, sec, e, n, left, right,p, x;
x:=Indeterminate(Rationals,"x");
p:=NumericalSemigroupPolynomial(s,x);
sec :=[];
left:=1;
right:=p;
for i in [1..k] do
invs:=p mod x^(i+2);
for n in [1..i-1] do
invs:=(invs*PowerMod((1-x^n),sec[n],x^(i+2))) mod x^(i+2);
od;
e:=CoefficientsOfUnivariatePolynomial(invs+x^(i+2))[i+1];
if e>0 then
right:=right*(1-x^i)^e;
else
left:=left*(1-x^i)^-e;
fi;
sec[i]:=e;
if left = right then
return -sec;
fi;
od;
return -sec;
end
);
##################################################################
##
#O WittCoefficients(p,k)
##
## p is a univariate polynomial with integer coefficientas and
## p(1)=1. Then p(x)=\prod_{n\ge 0}(1-x^n)^{e_n}.
## output is the list [e_1,..,e_k], and it is computed bu using
## [C-GS-HP-M]
##################################################################
InstallMethod(WittCoefficients,[IsUnivariatePolynomial, IsPosInt],
function(p,k)
local e, s, divisors, j, a;
a:=CoefficientsOfUnivariatePolynomial(p);
if a[1]<> 1 then
Error("The value of the first argument in 1 must be 1");
fi;
if not(IsSubset(Integers,a)) then
Error("The firs argument must be a polynomial with integer coefficients");
fi;
a:=Concatenation(a{[2..Length(a)]},List([Length(a)-1..k],_->0));
s:=[-a[1]];
for j in [2..k] do
s[j]:=-Reversed(a{[1..j-1]})*s{[1..j-1]}-j*a[j];
od;
e:=[];
for j in [1..k] do
divisors:=DivisorsInt(j);
e[j]:=1/j*s{divisors}*List(divisors, i->MoebiusMu(j/i));
od;
return e;
end);
##################################################################
##
#F LegendrianGenericNumericalSemigroup(n,m)
## n and m are coprime integers with m>=2n+1. The output is the
## semigroup of a generic element in the class of irreducible
## Legendrian singularities with equisingularity equial to the
## topological type of y^n=x^m.
################################################################
InstallGlobalFunction(LegendrianGenericNumericalSemigroup,
function(n,m)
local gamma, i, ni, wi, c, s, li, smg;
if not(IsPosInt(n) and IsPosInt(m) and Gcd(n,m)=1 and m>2*n) then
Error("The arguments are two positive coprime integers such that the second is larger than twice the first");
fi;
s:=NumericalSemigroup(n,m-n);
c:=Conductor(s);
gamma:=MultipleOfNumericalSemigroup(NumericalSemigroup(1),n,c);
i:=0;
while(i<>m-n) do
i:=Maximum(Difference(s,gamma));
smg:=SmallElements(gamma);
li:=Difference([i..Conductor(gamma)], smg);
ni:=Minimum(NrRestrictedPartitions(i,[n,m,m-n]),Length(li));
wi:=li[ni];
gamma:=NumericalSemigroupBySmallElements(Union(smg,[i..wi]));
od;
return gamma;
end);
[ Dauer der Verarbeitung: 0.48 Sekunden
(vorverarbeitet)
]
|
