|
#############################################################################
##
## This file is part of GAP, a system for computational discrete algebra.
## This file's authors include Alexander Hulpke, David Joyner.
##
## Copyright of GAP belongs to its developers, whose names are too numerous
## to list here. Please refer to the COPYRIGHT file for details.
##
## SPDX-License-Identifier: GPL-2.0-or-later
##
## This file contains the implementations for monomial orderings and Groebner
## bases.
# convenience function to catch cases
BindGlobal("GBasisFixOrderingArgument",function(a)
if IsFunction(a) then
a:=a();
fi;
return a;
end);
BindGlobal("MakeMonomialOrdering",function(name,ercomp)
local obj;
obj:=Objectify(NewType(MonomialOrderingsFamily,IsMonomialOrderingDefaultRep),
rec());
SetName(obj,name);
SetMonomialExtrepComparisonFun(obj,ercomp);
return obj;
end);
BindGlobal("InstallMonomialOrdering",function(ord,ordfun,idxord,
# for external interfaces. We use the names as in Singular
# this might change in future versions
type
)
InstallGlobalFunction(ord,function(arg)
local idx,nam,ordname,neword,ov;
if Length(arg)=1 and IsList(arg[1]) then
idx:=arg[1];
else
idx:=arg;
fi;
nam:=String(idx);
if Length(idx)>0 and not IsInt(idx[1]) then
idx:=List(idx,IndeterminateNumberOfUnivariateRationalFunction);
fi;
if Length(idx)>0 then
ov:=Set(idx);
else
ov:=true;
fi;
if IsSSortedList(idx) then
idx:=[];
fi;
if IsString(ord) then
ordname:=ord;
else
ordname:=NameFunction(ord);
fi;
if Length(idx)>0 then
# get an indexed ordering
nam:=Concatenation(ordname,"(",nam,")");
neword:=MakeMonomialOrdering(nam,
function(a,b)
# for each variable give its index position.
return idxord(a,b,List([1..Maximum(idx)],i->Position(idx,i)));
end);
neword!.idxarrangement:=Immutable(idx);
neword!.type:=Immutable(type);
SetOccuringVariableIndices(neword,ov);
return neword;
else
nam:=Concatenation(ordname,"()");
neword:=MakeMonomialOrdering(nam,ordfun);
neword!.idxarrangement:=Immutable(idx);
neword!.type:=Immutable(type);
SetOccuringVariableIndices(neword,ov);
return neword;
fi;
end);
end);
InstallMethod(MonomialComparisonFunction,"default: use extrep",true,
[IsMonomialOrderingDefaultRep],0,
function(ord)
local fun;
fun:=MonomialExtrepComparisonFun(ord);
return function(a,b)
local e, f, le, lf,fam;
fam:=FamilyObj(a);
e:=ExtRepPolynomialRatFun(a);
f:=ExtRepPolynomialRatFun(b);
repeat
if Length(e)=0 and Length(f)>0 then return true;
elif Length(f)=0 then return false;
fi;
if Length(e)=2 then
le:=1;
else
le:=LeadingMonomialPosExtRep(fam,e,fun);
fi;
if Length(f)=2 then
lf:=1;
else
lf:=LeadingMonomialPosExtRep(fam,f,fun);
fi;
if fun(e[le],f[lf]) then
return true;
elif e[le]<>f[lf] then
return false;
fi;
e:=e{Concatenation([1..le-1],[le+2..Length(e)])};
f:=f{Concatenation([1..lf-1],[lf+2..Length(f)])};
until false;
end;
end);
#############################################################################
##
#F MonomialLexOrdering()
#F MonomialLexOrdering(<vari>)
##
## This function creates a lexicographic ordering for monomials. If <vari>
## is given and is a list of variables or variable indices, this is used as
## the underlying order of variables.
InstallMonomialOrdering(MonomialLexOrdering,
#ordinary case
function(a,b)
local l,m,i;
l:=Length(a);
m:=Length(b);
i:=1;
while i<l and i<m do
# in this ordering x_1>x_2
if a[i]>b[i] then
return true;
elif a[i]<b[i] then
return false;
elif a[i+1]<b[i+1] then
return true;
elif a[i+1]>b[i+1] then
return false;
fi;
i:=i+2;
od;
if i<m then
return true;
fi;
# bigger or equal
return false;
end,
# indexed case
function(a,b,idx)
local l, m, min, i, am, ap, ai, has, bi,ret;
l:=Length(a);
m:=Length(b);
# as variables are not necessarily stored in the right order, we'll have
# to run through by variables one by one, The biggest index variables
# first
min:=0;
repeat
# get the smallest position (i.e. largest) remaining variable in a
i:=1;
am:=infinity;
ap:=0;
while i<l do
ai:=idx[a[i]];
if ai>min and ai<am then
# smaller pos (i.e. larger) variable found
am:=ai;
ap:=i;
fi;
i:=i+2;
od;
if am=infinity then
# no variable left in a. Test what is left in b
i:=1;
has:=false;
while i<m do
bi:=idx[b[i]];
if bi>min then
# b has variables left. So a is smaller
return true;
fi;
i:=i+2;
od;
# b also not. thus a must be equal to b
return false;
fi;
# now search for am in b
i:=1;
has:=false;
while i<m do
bi:=idx[b[i]];
if bi>min then
if bi<am then
# b has a larger (smaller pos) variable left. Thus b is bigger.
return true;
elif bi=am then
# the same variable occurs.
has:=true;
if a[ap+1]>b[i+1] then
# a exponent is bigger
ret:=false; # unless a smaller (larger pos) variable found
elif a[ap+1]<b[i+1] then
# b exponent is bigger
ret:=true; # unless a smaller (larger pos) variable found
else
# same exponents -- cannot decide on this variable
ret:=fail;
fi;
fi;
fi;
i:=i+2;
od;
if not has then
# b has no variable as large as am. thus a is bigger
return false;
elif ret<>fail then
# b has no larger variable than am buty has am. Use the comparison
return ret;
fi;
min:=am; # will increase until no variable in a left
until false;
end,
MakeImmutable("lp"));
#############################################################################
##
#F MonomialGrlexOrdering()
#F MonomialGrlexOrdering(<vari>)
##
## This function creates a degree/lexicographic ordering for monomials. If
## <vari> is given and is a list of variables or variable indices, this is
## used as the underlying order of variables.
InstallMonomialOrdering(MonomialGrlexOrdering,
function(a,b)
local s,t,i,l,m;
s:=0;
l:=Length(a);
m:=Length(b);
for i in [2,4..l] do
s:=s+a[i];
od;
t:=0;
for i in [2,4..m] do
t:=t+b[i];
od;
if s<t then
return true;
elif s>t then
return false;
fi;
# Lex
i:=1;
while i<l and i<m do
# in this ordering x_1>x_2
if a[i]>b[i] then
return true;
elif a[i]<b[i] then
return false;
elif a[i+1]<b[i+1] then
return true;
elif a[i+1]>b[i+1] then
return false;
fi;
i:=i+2;
od;
if i<m then
return true;
fi;
# bigger or equal
return false;
end,
# indexed case
function(a,b,idx)
local s, l, m, t, min, i, am, ap, ai, has, bi, ret;
s:=0;
l:=Length(a);
m:=Length(b);
for i in [2,4..l] do
s:=s+a[i];
od;
t:=0;
for i in [2,4..m] do
t:=t+b[i];
od;
if s<t then
return true;
elif s>t then
return false;
fi;
# Lex
# as variables are not necessarily stored in the right order, we'll have
# to run through by variables one by one, The biggest index variables
# first
min:=0;
repeat
# get the smallest position (i.e. largest) remaining variable in a
i:=1;
am:=infinity;
ap:=0;
while i<l do
ai:=idx[a[i]];
if ai>min and ai<am then
# smaller pos (i.e. larger) variable found
am:=ai;
ap:=i;
fi;
i:=i+2;
od;
if am=infinity then
# no variable left in a. Test what is left in b
i:=1;
has:=false;
while i<m do
bi:=idx[b[i]];
if bi>min then
# b has variables left. So a is smaller
return true;
fi;
i:=i+2;
od;
# b also not. thus a must be equal to b
return false;
fi;
# now search for am in b
i:=1;
has:=false;
while i<m do
bi:=idx[b[i]];
if bi>min then
if bi<am then
# b has a larger (smaller pos) variable left. Thus b is bigger.
return true;
elif bi=am then
# the same variable occurs.
has:=true;
if a[ap+1]>b[i+1] then
# a exponent is bigger
ret:=false; # unless a smaller (larger pos) variable found
elif a[ap+1]<b[i+1] then
# b exponent is bigger
ret:=true; # unless a smaller (larger pos) variable found
else
# same exponents -- cannot decide on this variable
ret:=fail;
fi;
fi;
fi;
i:=i+2;
od;
if not has then
# b has no variable as large as am. thus a is bigger
return false;
elif ret<>fail then
# b has no larger variable than am buty has am. Use the comparison
return ret;
fi;
min:=am; # will increase until no variable in a left
until false;
end,
MakeImmutable("Dp"));
#############################################################################
##
#F MonomialGrevlexOrdering()
#F MonomialGrevlexOrdering(<vari>)
##
## This function creates a degree/reverse lexicographic ordering for
## monomials. If <vari> is given and is a list of variables or variable
## indices, this is used as the underlying order of variables.
InstallMonomialOrdering(MonomialGrevlexOrdering,
#ordinary
function(a,b)
local s,t,i,j,l,m;
s:=0;
l:=Length(a);
m:=Length(b);
for i in [2,4..l] do
s:=s+a[i];
od;
t:=0;
for i in [2,4..m] do
t:=t+b[i];
od;
if s<t then
return true;
elif s>t then
return false;
fi;
# Revlex
i:=l-1;
j:=m-1;
while i>0 and j>0 do
# in this ordering x_1>x_2
if a[i]>b[j] then
# case x/0 -- a not bigger
return true;
elif a[i]<b[j] then
# case 0/x a bigger
return false;
elif a[i+1]<b[j+1] then
# a-b is negative
return false;
elif a[i+1]>b[j+1] then
# a-b is positive
return true;
fi;
i:=i-2;
j:=j-2;
od;
if i>0 then
#a-part left
return true;
fi;
# bigger or equal
return false;
end,
# indexed case
function(a,b,idx)
Error("indexed grevlex not yet implemented");
end,
MakeImmutable("dp"));
#############################################################################
##
#F EliminationOrdering(<elim>)
#F EliminationOrdering(<elim>,<rest>)
##
## This function creates an elimination ordering for eliminating the
## variables in <elim>. Two monomials are compared first by the exponent
## vectors for the variables listed in <elim> (a lexicographic comparison
## with respect to the ordering indicated in <elim>).
## If these submonomial are equal, the submonomials given by the other
## variables are compared by a graded lexicographic ordering (with respect
## to the variable order given in <rest>, if called with two parameters).
##
## Both <elim> and <rest> may be a list of variables of a list of variable
## indices.
InstallGlobalFunction(EliminationOrdering,function(arg)
local elimvar, nam, ord1, othvar, ord2,ov,neword;
elimvar:=arg[1];
nam:=ShallowCopy(String(elimvar));
if not IsInt(elimvar[1]) then
elimvar:=List(elimvar,IndeterminateNumberOfUnivariateRationalFunction);
fi;
ov:=Set(elimvar);
ord1:=MonomialExtrepComparisonFun(MonomialLexOrdering(elimvar));
if Length(arg)>1 then
othvar:=arg[2];
Add(nam,',');
nam:=Concatenation(nam,String(othvar));
if not IsInt(othvar[1]) then
othvar:=List(othvar,IndeterminateNumberOfUnivariateRationalFunction);
fi;
ov:=Union(ov,othvar);
ord2:=MonomialExtrepComparisonFun(MonomialGrlexOrdering(othvar));
else
ov:=true;
ord2:=MonomialExtrepComparisonFun(MonomialGrlexOrdering());
fi;
neword:=MakeMonomialOrdering(
Concatenation("EliminationOrdering(",nam,")"),
function(a,b)
local la, lb, asel, bsel, ah, bh, i;
la:=Length(a);
lb:=Length(b);
asel:=[];
for i in [1,3..la-1] do
if a[i] in elimvar then
Add(asel,i);
Add(asel,i+1);
fi;
od;
bsel:=[];
for i in [1,3..lb-1] do
if b[i] in elimvar then
Add(bsel,i);
Add(bsel,i+1);
fi;
od;
ah:=a{asel};
bh:=b{bsel};
if ah=bh then
ah:=a{Difference([1..la],asel)};
bh:=b{Difference([1..lb],bsel)};
return ord2(ah,bh);
else
return ord1(ah,bh);
fi;
end);
SetOccuringVariableIndices(neword,ov);
return neword;
end);
InstallMethod(MonomialExtrepComparisonFun,
"functions are themselves -- for compatibility",true,[IsFunction],0,f->f);
InstallOtherMethod(OccuringVariableIndices,"polynomial",true,
[IsPolynomial],0,
function(p)
local ov, l, i, j;
p:=ExtRepPolynomialRatFun(p);
ov:=[];
for i in [1,3..Length(p)-1] do
l:=p[i];
for j in [1,3..Length(l)-1] do
AddSet(ov,l[j]);
od;
od;
return ov;
end);
InstallMethod(LeadingTermOfPolynomial,"with ordering",true,
[IsPolynomialFunction,IsMonomialOrdering],0,
function(p,order)
local e,fam,a;
order:=MonomialExtrepComparisonFun(order);
e:=ExtRepPolynomialRatFun(p);
fam:=FamilyObj(p);
a:=LeadingMonomialPosExtRep(fam,e,order);
return PolynomialByExtRepNC(fam,e{[a,a+1]});
end);
InstallOtherMethod(LeadingMonomialOfPolynomial,"with ordering",true,
[IsPolynomialFunction,IsMonomialOrdering],0,
function(p,order)
local e,fam,a;
if not IsPolynomial(p) then
TryNextMethod();
fi;
order:=MonomialExtrepComparisonFun(order);
e:=ExtRepPolynomialRatFun(p);
fam:=FamilyObj(p);
a:=LeadingMonomialPosExtRep(fam,e,order);
return PolynomialByExtRepNC(fam,[e[a],One(CoefficientsFamily(fam))]);
end);
InstallOtherMethod(LeadingCoefficientOfPolynomial,"with ordering",true,
[IsPolynomialFunction,IsMonomialOrdering],0,
function(p,order)
local e,fam,a;
if not IsPolynomial(p) then
TryNextMethod();
fi;
order:=MonomialExtrepComparisonFun(order);
e:=ExtRepPolynomialRatFun(p);
fam:=FamilyObj(p);
a:=LeadingMonomialPosExtRep(fam,e,order);
return e[a+1];
end);
# compute S-polynomial and return `0' if the leading monomials are coprime
BindGlobal("SPolynomial",function(p,q,order)
local a,b,c,d,e,f,i,j,ae,be,di,fam;
order:=MonomialExtrepComparisonFun(order);
e:=ExtRepPolynomialRatFun(p);
f:=ExtRepPolynomialRatFun(q);
fam:=FamilyObj(p);
a:=LeadingMonomialPosExtRep(fam,e,order);
b:=LeadingMonomialPosExtRep(fam,f,order);
ae:=e[a];
be:=f[b];
# compute the cofactors
i:=1;
j:=1;
c:=[];
d:=[];
while i<Length(ae) and j<Length(be) do
if ae[i]<be[j] then
# a has variable not in b
Append(d,ae{[i,i+1]});
i:=i+2;
elif ae[i]>be[j] then
# b has variable not in a
Append(c,be{[j,j+1]});
j:=j+2;
else
# variable in both: One larger?
di:=be[j+1]-ae[i+1];
if di>0 then
Append(c,[ae[i],di]);
elif di<0 then
Append(d,[be[j],-di]);
fi;
i:=i+2;
j:=j+2;
fi;
od;
# trailing entries
if i<Length(ae) then
Append(d,ae{[i..Length(ae)]});
elif j<Length(be) then
Append(c,be{[j..Length(be)]});
fi;
return PolynomialByExtRepNC(fam,[c,f[b+1]])*p
-PolynomialByExtRepNC(fam,[d,e[a+1]])*q;
end);
#############################################################################
##
#F PolynomialReduction(poly,plist,order) reduces poly with the ideal
## generated by plist, according to order
## The routine returns a list [remainder,[quotients]]
##
InstallGlobalFunction( PolynomialReduction, function(poly,plist,order)
local fam,quot,elist,lmp,lmo,lmc,x,y,z,mon,mon2,qmon,noreduce,
ep,pos,di,rem,qmex;
if IsMonomialOrdering(order) then
order:=MonomialExtrepComparisonFun(order);
fi;
fam:=FamilyObj(poly);
quot:=List(plist,i->Zero(poly));
elist:=List(plist,ExtRepPolynomialRatFun);
lmp:=List(elist,i->LeadingMonomialPosExtRep(fam,i,order));
lmo:=List([1..Length(lmp)],i->elist[i][lmp[i]]);
lmc:=List([1..Length(lmp)],i->elist[i][lmp[i]+1]);
ep:=ExtRepPolynomialRatFun(poly);
rem:=Zero(poly);
while Length(ep)>0 do
# now try whether we can reduce at x
noreduce:=true;
x:=LeadingMonomialPosExtRep(fam,ep,order);
mon:=ep[x];
y:=1;
while y<=Length(plist) and noreduce do
mon2:=lmo[y];
#check whether the monomial at position x is a multiple of the
#y-th leading monomial
z:=1;
pos:=1;
qmon:=[]; # potential quotient
noreduce:=false;
while noreduce=false and z<=Length(mon2) and pos<=Length(mon) do
if mon[pos]>mon2[z] then
noreduce:=true; # indet in mon2 does not occur in mon -> does not
# divide
elif mon[pos]<mon2[z] then
Append(qmon,mon{[pos,pos+1]}); # indet only in mon
pos:=pos+2;
else
# the indets are the same
di:=mon[pos+1]-mon2[z+1];
if di>0 then
#divides and there is remainder
Append(qmon,[mon[pos],di]);
elif di<0 then
noreduce:=true; # exponent to small
fi;
pos:=pos+2;
z:=z+2;
fi;
od;
# if there is a tail in mon2 left, cannot divide
if z<=Length(mon2) then
noreduce:=true;
fi;
y:=y+1;
od;
if noreduce then
mon:=PolynomialByExtRepNC(fam,[ep[x],ep[x+1]]);
rem:=rem+mon;
#Print("noreduce:",PolynomialByExtRepNC(fam,ep)," ",mon," ",rem,"\n");
ep:=ep{Difference([1..Length(ep)],[x,x+1])};
#Print(ep,"\n");
else
y:=y-1; # re-correct incremented numbers
# is there a tail in mon? if yes we need to append it
if pos<Length(mon) then
Append(qmon,mon{[pos..Length(mon)]});
fi;
# reduce!
qmex:=[qmon,-ep[x+1]/lmc[y]];
qmon:=PolynomialByExtRepNC(fam,[qmon,ep[x+1]/lmc[y]]); #quotient monomial
quot[y]:=quot[y]+qmon;
qmex:=ZippedProduct(qmex,elist[y],
fam!.zeroCoefficient,fam!.zippedProduct);
qmex:=ZippedSum(ep,qmex,fam!.zeroCoefficient,fam!.zippedSum);
#poly:=PolynomialByExtRep(fam,ep);
#poly:=poly-qmon*plist[y]; # reduce
#ep:=ExtRepPolynomialRatFun(poly);
#if ep<>qmex then Error("RED!"); fi;
ep:=qmex;
fi;
od;
return [rem,quot];
end);
#############################################################################
##
#F PolynomialReducedRemainder(poly,plist,order)
##
InstallGlobalFunction(PolynomialReducedRemainder,function(poly,plist,order)
local fam, elist, lmp, lmo, lmc, ep, rem, noreduce, x, mon, y, mon2,
z, pos, qmon, di,qmex;
if IsMonomialOrdering(order) then
order:=MonomialExtrepComparisonFun(order);
fi;
fam:=FamilyObj(poly);
elist:=List(plist,ExtRepPolynomialRatFun);
lmp:=List(elist,i->LeadingMonomialPosExtRep(fam,i,order));
lmo:=List([1..Length(lmp)],i->elist[i][lmp[i]]);
lmc:=List([1..Length(lmp)],i->elist[i][lmp[i]+1]);
ep:=ExtRepPolynomialRatFun(poly);
rem:=Zero(poly);
while Length(ep)>0 do
# now try whether we can reduce at x
noreduce:=true;
x:=LeadingMonomialPosExtRep(fam,ep,order);
mon:=ep[x];
y:=1;
while y<=Length(plist) and noreduce do
mon2:=lmo[y];
#check whether the monomial at position x is a multiple of the
#y-th leading monomial
z:=1;
pos:=1;
qmon:=[]; # potential quotient
noreduce:=false;
while noreduce=false and z<=Length(mon2) and pos<=Length(mon) do
if mon[pos]>mon2[z] then
noreduce:=true; # indet in mon2 does not occur in mon -> does not
# divide
elif mon[pos]<mon2[z] then
Append(qmon,mon{[pos,pos+1]}); # indet only in mon
pos:=pos+2;
else
# the indets are the same
di:=mon[pos+1]-mon2[z+1];
if di>0 then
#divides and there is remainder
Append(qmon,[mon[pos],di]);
elif di<0 then
noreduce:=true; # exponent to small
fi;
pos:=pos+2;
z:=z+2;
fi;
od;
# if there is a tail in mon2 left, cannot divide
if z<=Length(mon2) then
noreduce:=true;
fi;
y:=y+1;
od;
if noreduce then
mon:=PolynomialByExtRepNC(fam,[ep[x],ep[x+1]]);
rem:=rem+mon;
ep:=ep{Difference([1..Length(ep)],[x,x+1])};
else
y:=y-1; # re-correct incremented numbers
# is there a tail in mon? if yes we need to append it
if pos<Length(mon) then
Append(qmon,mon{[pos..Length(mon)]});
fi;
# reduce!
qmex:=[qmon,-ep[x+1]/lmc[y]];
qmex:=ZippedProduct(qmex,elist[y],
fam!.zeroCoefficient,fam!.zippedProduct);
qmex:=ZippedSum(ep,qmex,fam!.zeroCoefficient,fam!.zippedSum);
#qmon:=PolynomialByExtRepNC(fam,[qmon,ep[x+1]/lmc[y]]); #quotient monomial
#poly:=PolynomialByExtRep(fam,ep);
#poly:=poly-qmon*plist[y]; # reduce
#ep:=ExtRepPolynomialRatFun(poly);
#if ep<>qmex then Error("RED0!"); fi;
ep:=qmex;
fi;
od;
return rem;
end);
#############################################################################
##
#F PolynomialDivisionAlgorithm(poly,plist,order)
##
InstallGlobalFunction(PolynomialDivisionAlgorithm,function(poly,plist,order)
local fam,quot,elist,lmp,lmo,lmc,x,y,z,mon,mon2,qmon,noreduce,
ep,pos,di,rem;
if IsMonomialOrdering(order) then
order:=MonomialExtrepComparisonFun(order);
fi;
fam:=FamilyObj(poly);
quot:=List(plist,i->Zero(poly));
elist:=List(plist,ExtRepPolynomialRatFun);
lmp:=List(elist,i->LeadingMonomialPosExtRep(fam,i,order));
lmo:=List([1..Length(lmp)],i->elist[i][lmp[i]]);
lmc:=List([1..Length(lmp)],i->elist[i][lmp[i]+1]);
rem:=Zero(poly);
ep:=ExtRepPolynomialRatFun(poly);
while Length(ep)>0 do
# now try whether we can reduce the leading monomial
x:=LeadingMonomialPosExtRep(fam,ep,order);
mon:=ep[x];
y:=1;
noreduce:=true;
while y<=Length(plist) and noreduce do
mon2:=lmo[y];
#check whether the monomial at position x (leading mon of poly)
#is a multiple of the y-th leading monomial
z:=1;
pos:=1;
qmon:=[]; # potential quotient
noreduce:=false;
while noreduce=false and z<=Length(mon2) and pos<=Length(mon) do
if mon[pos]>mon2[z] then
noreduce:=true; # indet in mon2 does not occur in mon -> does not
# divide
elif mon[pos]<mon2[z] then
Append(qmon,mon{[pos,pos+1]}); # indet only in mon
pos:=pos+2;
else
# the indets are the same
di:=mon[pos+1]-mon2[z+1];
if di>0 then
#divides and there is remainder
Append(qmon,[mon[pos],di]);
elif di<0 then
noreduce:=true; # exponent to small
fi;
pos:=pos+2;
z:=z+2;
fi;
od;
# if there is a tail in mon2 left, cannot divide
if z<=Length(mon2) then
noreduce:=true;
fi;
y:=y+1;
od;
if noreduce=false then
y:=y-1; # re-correct incremented numbers
# is there a tail in mon? if yes we need to append it
if pos<Length(mon) then
Append(qmon,mon{[pos..Length(mon)]});
fi;
# reduce!
qmon:=PolynomialByExtRepNC(fam,[qmon,ep[x+1]/lmc[y]]); #quotient monomial
quot[y]:=quot[y]+qmon;
poly:=poly-qmon*plist[y]; # reduce
else
# remove leading monomial
qmon:=PolynomialByExtRepNC(fam,[mon,ep[x+1]]);
rem:=rem+qmon;
poly:=poly-qmon;
fi;
ep:=ExtRepPolynomialRatFun(poly);
od;
return [rem,quot];
end);
BindGlobal("SyzygyCriterion",function(baslte,i,j,t,B)
local li, lj, lcm, a, b, k;
lcm:=false;
for k in [1..t] do
if k<>i and k<>j and
(not Set([i,k]) in B) and (not Set([j,k]) in B) then
if lcm=false then
li:=baslte[i];
lj:=baslte[j];
lcm:=[];
a:=1;
b:=1;
while a<Length(li) and b<Length(lj) do
if li[a]<lj[b] then
# li variable is smaller
Add(lcm,li[a]);
Add(lcm,li[a+1]);
a:=a+2;
elif li[a]>lj[b] then
# lj-variable is smaller
Add(lcm,lj[b]);
Add(lcm,lj[b+1]);
b:=b+2;
else
# same variable
Add(lcm,li[a]);
Add(lcm,Maximum(li[a+1],lj[b+1]));
a:=a+2;
b:=b+2;
fi;
od;
# add remaining tails. Only one of these loops is run
while a<Length(li) do
Add(lcm,li[a]);
Add(lcm,li[a+1]);
a:=a+2;
od;
while b<Length(lj) do
Add(lcm,lj[b]);
Add(lcm,lj[b+1]);
b:=b+2;
od;
fi;
# does baslte[k] divide?
li:=baslte[k];
a:=1;
b:=1;
while a<Length(li) and b<Length(lcm) do
# lcm has extra variables?
while b<Length(lcm) and li[a]>lcm[b] do
b:=b+2;
od;
# test whether variable OK and variable power not bigger
if b<Length(lcm) then
if li[a]=lcm[b] and li[a+1]<=lcm[b+1] then
a:=a+2;
b:=b+2;
else
# either a has an extra (smaller) variable or a bigger
# exponent. In both cases division fails which we indicate by
# only b running out
b:=Length(lcm)+1;
fi;
fi;
od;
# if b has variables left or a and b both ran out, we're fine
if b<Length(lcm) or (b>Length(lcm) and a>Length(li)) then
return true;
fi;
fi;
od;
return false;
end);
BindGlobal("GAPGBASIS",MakeImmutable(rec(
name:="naive GAP version of Buchberger's algorithm",
GroebnerBasis:=function(elms,order)
local orderext, bas, baslte, fam, t, B, i, j, s;
orderext:=MonomialExtrepComparisonFun(order);
bas:=Filtered(elms,x->not IsZero(x));
baslte:=List(bas,ExtRepPolynomialRatFun);
fam:=FamilyObj(bas[1]);
baslte:=List(baslte,i->i[LeadingMonomialPosExtRep(fam,i,orderext)]);
t:=Length(bas);
B:=Concatenation(List([1..t],i->List([1..i-1],j->[j,i])));
while Length(B)>0 do
j:=B[1][1];
i:=B[1][2];
# remove first entry of B
Remove(B, 1);
if Length(Intersection(baslte[i]{[1,3..Length(baslte[i])-1]},
baslte[j]{[1,3..Length(baslte[j])-1]}))=0 then
Info(InfoGroebner,2,"Pair (",i,",",j,") avoided by product criterion");
elif SyzygyCriterion(baslte,i,j,t,B) then
Info(InfoGroebner,2,"Pair (",i,",",j,") avoided by chain criterion");
else
s:=SPolynomial(bas[i],bas[j],order);
if InfoLevel(InfoGroebner)<3 then
Info(InfoGroebner,2,"Spol(",i,",",j,")");
else
Info(InfoGroebner,3,"Spol(",i,",",j,")=",s);
fi;
s:=PolynomialReducedRemainder(s,bas,orderext);
if not IsZero(s) then
Info(InfoGroebner,3,"reduces to ",s);
Add(bas,s);
s:=ExtRepPolynomialRatFun(s);
Add(baslte,s[LeadingMonomialPosExtRep(fam,s,orderext)]);
t:=t+1;
# add new pairs
for i in [1..t-1] do
Add(B,[i,t]);
od;
Info(InfoGroebner,1,"|bas|=",t,", ",Length(B)," pairs left");
fi;
fi;
od;
return bas;
end)
));
#############################################################################
##
#V GBASIS
##
## This variable is assigned to a record which determines which
## implementation is used to compute Groebner bases. By default it is
## assigned to `GAPGBASIS', a naive implementation of Buchberger's
## algorithm, which is only provided for educational purposes. For serious
## work a better implementation should be used.
GBASIS:=GAPGBASIS;
InstallGlobalFunction(GroebnerBasisNC,
function(elms,order)
return GBASIS.GroebnerBasis(elms,order);
end);
InstallMethod(GroebnerBasis,"polynomials",true,
[IsHomogeneousList and IsRationalFunctionCollection,IsMonomialOrdering],0,
function(elms,order)
local ov,i,d;
# do some tests
if not ForAll(elms,IsPolynomial) then
Error("<elms> must be a list of polynomials");
fi;
ov:=OccuringVariableIndices(order);
if ov<>true then
for i in elms do
d:=Difference(OccuringVariableIndices(i),ov);
if Length(d)>0 then
Error("Ordering is undefined for variables ",
List(d,j->Indeterminate(DefaultRing([FamilyObj(i)!.zeroCoefficient]),j)));
fi;
od;
fi;
return GroebnerBasisNC(elms,order);
end);
InstallOtherMethod(GroebnerBasis,"fix function",true,
[IsObject,IsFunction],0,
function(obj,f)
if f=MonomialLexOrdering or f=MonomialGrlexOrdering or
f=MonomialGrevlexOrdering then
return GroebnerBasis(obj,f());
fi;
TryNextMethod();
end);
InstallMethod(StoredGroebnerBasis,"ideal",true,
[IsPolynomialRingIdeal],0,function(I)
local ord;
ord:=MonomialGrlexOrdering();
return [ReducedGroebnerBasis(GeneratorsOfTwoSidedIdeal(I),ord),ord];
end);
InstallMethod(GroebnerBasis,"ideal",true,
[IsPolynomialRingIdeal,IsMonomialOrdering],0,
function(I,order)
local bas;
if HasStoredGroebnerBasis(I) then
bas:=StoredGroebnerBasis(I);
if IsIdenticalObj(bas[2],order) then return bas[1];fi;
fi;
# do some tests
bas:=ReducedGroebnerBasis(GeneratorsOfIdeal(I),order);
if not HasStoredGroebnerBasis(I) then
SetStoredGroebnerBasis(I,[bas,order]);
fi;
return bas;
end);
InstallMethod(GroebnerBasis,"ideal with stored GB",true,
[IsPolynomialRingIdeal and HasStoredGroebnerBasis,IsMonomialOrdering],0,
function(I,order)
# do some tests
if not IsIdenticalObj(StoredGroebnerBasis(I)[2],order) then
TryNextMethod();
fi;
return StoredGroebnerBasis(I)[1];
end);
InstallMethod(ReducedGroebnerBasis,"polynomials",true,
[IsHomogeneousList and IsRationalFunctionCollection,IsMonomialOrdering],0,
function(elms,order)
local orderext, bas, i, l, nomod, pol;
orderext:=MonomialExtrepComparisonFun(order);
# do some tests
bas:=ShallowCopy(GBASIS.GroebnerBasis(elms,order));
i:=1;
l:=Length(bas);
repeat
nomod:=true;
while i<=l do
pol:=PolynomialReducedRemainder(bas[i],
bas{Difference([1..l],[i])},orderext);
if pol<>bas[i] then nomod:=false;fi;
if IsZero(pol) then
bas[i]:=bas[l];
Unbind(bas[l]);
l:=l-1;
else
bas[i]:=pol;
i:=i+1;
fi;
od;
until nomod;
for i in [1..Length(bas)] do
bas[i]:=bas[i]/LeadingCoefficientOfPolynomial(bas[i],order);
od;
Sort(bas,MonomialComparisonFunction(order));
return bas;
end);
InstallMethod(ReducedGroebnerBasis,"ideal",true,
[IsPolynomialRingIdeal,IsMonomialOrdering],0,
function(I,order)
local bas;
# do some tests
if HasStoredGroebnerBasis(I)
and IsIdenticalObj(StoredGroebnerBasis(I)[2],order) then
bas:=StoredGroebnerBasis(I)[1];
else
bas:=GeneratorsOfIdeal(I);
fi;
# reduce
bas:=ReducedGroebnerBasis(bas,order);
if not HasStoredGroebnerBasis(I) then
SetStoredGroebnerBasis(I,[bas,order]);
fi;
return bas;
end);
InstallOtherMethod(ReducedGroebnerBasis,"fix function",true,
[IsObject,IsFunction],0,
function(obj,f)
if f=MonomialLexOrdering or f=MonomialGrlexOrdering or
f=MonomialGrevlexOrdering then
return ReducedGroebnerBasis(obj,f());
fi;
TryNextMethod();
end);
InstallMethod(\in,"polynomial ideal",true,
[IsPolynomial,IsPolynomialRingIdeal],0,function(p,I)
local g;
g:=StoredGroebnerBasis(I);
p:=PolynomialReducedRemainder(p,g[1],MonomialExtrepComparisonFun(g[2]));
return IsZero(p);
end);
# We only provide a single GcdOp method (taking three arguments) for the
# general case, and no specialized two argument version, which we usually
# provide to avoid costly construction of the polynomial ring object, simply
# because in this case the cost of Gcd computation dwarfs the cost for
# creating the polynomial ring.
InstallOtherMethod(GcdOp,"multivariate Gcd based on Groebner bases",
IsCollsElmsElms,[IsPolynomialRing,IsPolynomial,IsPolynomial],0,
# Input: f, g are multivariate polys in R=F[x1,...,xn]
# Output: a gcd of f,g in R
function(R,f,g)
local basis_elt,t,F,vars,vars2,GB,coeff_t;
F:=CoefficientsRing(R);
vars:=IndeterminatesOfPolynomialRing(R);
t:=Indeterminate(F,vars);
vars2:=Concatenation([t],vars);
GB:=ReducedGroebnerBasis([t*f,(1-t)*g],MonomialLexOrdering(vars2));
for basis_elt in GB do
coeff_t:=PolynomialCoefficientsOfPolynomial(basis_elt,t);
if Length(coeff_t)=1 then
return StandardAssociate(f*g/coeff_t[1]);
fi;
od;
return One(F);
end);
InstallOtherMethod(LcmOp,"multivariate Lcm based on Groebner bases",
IsCollsElmsElms,[IsPolynomialRing,IsPolynomial,IsPolynomial],0,
function(R,f,g)
return f*g/GcdOp(R,f,g);
end);
BindGlobal("NaiveBuchberger",function(elms,order)
local orderext, bas, baslte, fam, t, B, i, j, s;
orderext:=MonomialExtrepComparisonFun(order);
bas:=ShallowCopy(elms);
baslte:=List(bas,ExtRepPolynomialRatFun);
fam:=FamilyObj(bas[1]);
baslte:=List(baslte,i->i[LeadingMonomialPosExtRep(fam,i,orderext)]);
t:=Length(bas);
B:=Concatenation(List([1..t],i->List([1..i-1],j->[j,i])));
while Length(B)>0 do
i:=B[1]; # take one
j:=i[1];
i:=i[2];
s:=SPolynomial(bas[i],bas[j],order);
if InfoLevel(InfoGroebner)<3 then
Info(InfoGroebner,2,"Spol(",i,",",j,")");
else
Info(InfoGroebner,3,"Spol(",i,",",j,")=",s);
fi;
s:=PolynomialReducedRemainder(s,bas,orderext);
if not IsZero(s) then
Info(InfoGroebner,3,"reduces to ",s);
Add(bas,s);
s:=ExtRepPolynomialRatFun(s);
Add(baslte,s[LeadingMonomialPosExtRep(fam,s,orderext)]);
t:=t+1;
# add new pairs
for i in [1..t] do
Add(B,[i,t]);
od;
Info(InfoGroebner,1,"|bas|=",t,", ",Length(B)," pairs left");
fi;
# remove first entry of B
Remove(B,1);
od;
return bas;
end);
# setup for quotient rings of polynomial rings
#############################################################################
##
#M NaturalHomomorphismByIdeal(<R>,<I>)
##
InstallMethod( NaturalHomomorphismByIdeal,"polynomial rings",IsIdenticalObj,
[ IsPolynomialRing,IsRing],
function(R,I)
local ord,b,ind,num,c,corners,i,j,a,rem,bound,mon,n,monb,dim,sc,k,l,char,hom;
if not IsIdeal(R,I) then
Error("I is not an ideal!");
fi;
if not IsField(LeftActingDomain(R)) then
TryNextMethod();
fi;
char:=Characteristic(LeftActingDomain(R));
ord:=MonomialGrlexOrdering();
b:=GroebnerBasis(I,ord);
# special case: unit in ideal -- quotient is trivial
if ForAny(b,IsUnit) then
a:=Subring(R,[Zero(R)]);
hom:=MappingByFunction(R,a,x->Zero(a));
SetImagesSource(hom,a);
return hom;
fi;
# determine all monomials bounded
ind:=IndeterminatesOfPolynomialRing(R);
n:=Length(ind);
num:=List(ind,IndeterminateNumberOfUnivariateRationalFunction);
c:=List(b,x->ExtRepNumeratorRatFun(LeadingMonomialOfPolynomial(x,ord))[1]);
corners:=[];
bound:=List(ind,x->infinity);
for i in c do
a:=List(ind,x->0);
for j in [1,3..Length(i)-1] do
a[Position(num,i[j])]:=i[j+1];
od;
# single variable bound
if Number(a,x->x<>0)=1 then
bound[PositionNonZero(a)]:=Sum(a);
else
Add(corners,a); # extra corner
fi;
od;
if not ForAll(bound,IsInt) then
Info(InfoWarning,1,"quotient ring has infinite dimension");
TryNextMethod();
fi;
#run through bounded simplex
a:=List(ind,x->0);
mon:=[a];
i:=Length(a);
repeat
a[i]:=a[i]+1;
while i>0 and a[i]>=bound[i] do
a[i]:=0;
i:=i-1;
if i>0 then
a[i]:=a[i]+1;
fi;
od;
if i>0 then
if ForAny(corners,x->ForAll([1..n],j->x[j]<=a[j])) then
# set last coordinate to become 0 again
a[n]:=bound[n];
i:=n;
else
Add(mon,ShallowCopy(a));
i:=Length(a);
fi;
fi;
until i=0;
# basis of the quotient as vector space
mon:=List(mon,x->Product([1..n],y->ind[y]^x[y]));
monb:=Basis(VectorSpace(LeftActingDomain(R),mon));
# now form the quotient
dim:=Length(mon);
sc:=EmptySCTable(dim,0);
for i in [1..dim] do
for j in [1..dim] do
rem:=PolynomialReducedRemainder(mon[i]*mon[j],b,ord);
a:=rem;
if not IsZero(a) then
a:=Coefficients(monb,a);
if char>0 then
a:=List(a,Int);
fi;
l:=[];
for k in [1..dim] do
if not IsZero(a[k]) then
Add(l,a[k]);
Add(l,k);
fi;
od;
SetEntrySCTable(sc,i,j,l);
fi;
od;
od;
l:=List(mon,x->Concatenation("(",
Filtered(String(x),y->y<>'^' and y<>'*'),")"));
a:=PositionProperty(mon,IsOne);
if a<>fail then
l[a]:="(1)";
fi;
a:=AlgebraByStructureConstants(LeftActingDomain(R),sc,l);
if Length(l)<50 then
j:=Concatenation("<ring ",String(LeftActingDomain(R)));
for i in l do
Append(j,",");
Append(j,i);
od;
Append(j,">");
SetName(a,j);
fi;
#k:=Concatenation([One(R)],IndeterminatesOfPolynomialRing(R));
k:=IndeterminatesOfPolynomialRing(R);
l:=[];
for i in k do
j:=Position(mon,i);
if j<>fail then
Add(l,GeneratorsOfLeftOperatorRing(a)[j]);
else
rem:=PolynomialReducedRemainder(i,b,ord);
rem:=Coefficients(monb,rem);
Add(l,GeneratorsOfLeftOperatorRing(a)*rem);
fi;
od;
hom:=AlgebraWithOneHomomorphismByImagesNC(R,a,k,l);
SetIsSurjective(hom,true);
SetKernelOfAdditiveGeneralMapping(hom,I);
j:=AlgebraWithOneGeneralMappingByImages(a,R,l,k);
SetInverseGeneralMapping(hom,j);
return hom;
#x*y^3-x^2,x^3*y^2-y
end);
# BindGlobal("CANGB",function(elms,order)
#local orderext, bas, i, l, nomod, pol;
# orderext:=MonomialExtrepComparisonFun(order);
# # do some tests
# bas:=ShallowCopy(elms);
# i:=1;
# l:=Length(bas);
# repeat
# nomod:=true;
# while i<=l do
# pol:=PolynomialReducedRemainder(bas[i],
# bas{Difference([1..Length(bas)],[i])},orderext);
# if pol<>bas[i] then nomod:=false;fi;
# if IsZero(pol) then
# bas[i]:=bas[l];
# Unbind(bas[l]);
# l:=l-1;
# else
# bas[i]:=pol;
# i:=i+1;
# fi;
# od;
# until nomod;
# for i in [1..Length(bas)] do
# bas[i]:=bas[i]/LeadingCoefficientOfPolynomial(bas[i],order);
# od;
# Sort(bas,MonomialComparisonFunction(order));
# return bas;
#end);
#
#DoTest:=function(l,ord)
#local a,b,t1,t2;
# t1:=Runtime();
# a:=NaiveBuchberger(l,ord);
# t1:=Runtime()-t1;
# t2:=Runtime();
# b:=GroebnerBasis(l,ord);
# t2:=Runtime()-t2;
# a:=CANGB(a,ord);
# b:=CANGB(b,ord);
# return [t1,t2,a=b];
#end;
[ Verzeichnis aufwärts0.65unsichere Verbindung
Übersetzung europäischer Sprachen durch Browser
]
|
