"""
Faugere's F_5 algorithm
This implementation is along the lines of John Perry's pseudocode
and Singular implementation. It was inspired by Martin Albrecht's
implementation and discussions with him and Ludovic Perret.
See http://www.math.usm.edu/perry/Research/ for details.
AUTHOR:
-- 20081022 Simon King
EXAMPLE:
sage: R. = PolynomialRing(GF(29))
sage: I = R* [3*x^4*y + 18*x*y^4 + 4*x^3*y*z + 20*x*y^3*z + 3*x^2*z^3, \
3*x^3*y^2 + 7*x^2*y^3 + 24*y^2*z^3,
12*x*y^4 + 17*x^4*z + 27*y^4*z + 11*x^3*z^2]
sage: J = I.homogenize()
sage: f5 = F5() # original F5
sage: gb = f5(J)
Generator 2/3
1 critical pairs in degree 7
1 critical pairs in degree 9
1 critical pairs in degree 11
Generator 3/3
1 critical pairs in degree 6
2 critical pairs in degree 7
5 critical pairs in degree 8
7 critical pairs in degree 9
11 critical pairs in degree 10
10 critical pairs in degree 11
WARNING: top-reduction to zero
WARNING: reduction to zero
WARNING: reduction to zero
3 critical pairs in degree 12
5 critical pairs in degree 13
2 critical pairs in degree 14
1 critical pairs in degree 16
Note that we encountered 3 S-polynomials reducing to zero
sage: f5.zeroes(), len(gb)
(3, 18)
"""
import sage
import sage.all
from sage.all import copy
import sys
from sage.structure.sequence import Sequence
from sage.rings.ideal import Ideal
#from sage.rings.polynomial.multi_polynomial_libsingular cimport MPolynomial_libsingular
class CounterClass:
def __init__(self):
self.Gen = 0
self.Del = 0
def __call__(self,i):
if i>0:
self.Gen += 1
else:
self.Del += 1
Counter = CounterClass()
CounterB = CounterClass()
def dehomogenize(I,R):
"dehomogenize(I,R): dehomogenize an ideal I that was obtained by homogenizing an ideal of ring R"
h = I[0].parent()('h')
return ([R(X(h=1)) for X in I]*R).reduced_basis()
cdef class DecoratedPolynomial:
"""
A decorated polynomial is a polynomial p in an ideal I = of a polynomial
ring R, together with a 'signature' (u,i), where u is a monomial in R and 1 <= i <= m
is an integer. The signature has the following meaning:
* One can express p as an ideal combination $\sum_{k=1}^i t_k\cdot f_k$, and u is the leading
term of $t_i$ of f_1,...,f_i.
The signature allows for avoiding some useless critical pairs in the Buchberger algorithm.
"""
cdef int nu
cdef object mu # the multiplier
cdef tuple Emu # exponents of the multiplier
#cdef MPolynomial_libsingular p
cdef object p
def __init__(self,mu, int nu, p):
self.mu = mu
self.Emu = mu.degrees()
self.nu = nu
self.p = p
#self.rule = -1
#self.hist = []
#def __init__(self,mu, int nu, object p):
# pass
def sig(DecoratedPolynomial self):
return (self.mu,self.nu)
def poly(DecoratedPolynomial self):
return self.p
#def info(self):
# return self.hist
def __richcmp__(DecoratedPolynomial self, DecoratedPolynomial other, op):
# < 0, <= 1, == 2, != 3, > 4, >= 5
# Idea for applications: When considering a pair of decorated monomials and
# reducing them, we always reduce the "bigger" one,
if op == 0:
if self.nuother.nu:
return True
if (self.nu==other.nu) and (self.mu>other.mu):
return True
return False
if op == 5:
if self.nu>other.nu:
return True
if (self.nu==other.nu) and (self.mu>=other.mu):
return True
return False
cdef class F5:
cdef dict Rules # rewriting rules (indexed by leading monomial, yielding a monomial ideal)
cdef object Rgb # reduced GB of the first few input generators
cdef object RgbMon # their leading ideal
cdef list Inputgen # the ideal generators we got, interreduced and increasingly sorted
cdef list Basis # list of decorated polynomials that yields self.Rgb
cdef object Ring # the polynomial ring we started with
cdef object HRing # homogenized ring
cdef object One # "One" in the ring we work with (may be homogenisation of self.Ring
cdef int nvars # number of variables of HRing
cdef int homog # 1, if the input was homogeneous; 0 otherwise (then we homogenize/dehomogenize)
cdef int Zeroes # Number of S-polynomials reducing to zero, to test the efficiency of
# our implementation of the F5 criterion.
# We are not keeping track of the corresponding critical pairs.
def __init__(F5 self):
"F = F5() set up the machinery to compute Groebner bases with the F5 algorithm"
self.Rules={}
self.Rgb = None
self.RgbMon = None
self.Inputgen = []
self.Basis = []
self.Ring = None
self.HRing = None
self.One = None
self.nvars = 0
self.homog = 0
self.Zeroes = 0
#def __init__(self):
# pass
#def __dealloc__(self):
# self.Rules={}
# self.Rgb = None
# self.RgbMon = None
# self.Inputgen = []
# self.Basis = []
# self.Ring = None
# self.HRing = None
# self.One = None
# self.nvars = 0
# self.homog = 0
# self.Zeroes = 0
# print "Deleted"
def __call__(F5 self, I):
"""
F(I): Apply the F5 algorithm to an ideal I and return its reduced Groebner basis
"""
if not isinstance(I, sage.rings.polynomial.multi_polynomial_ideal.MPolynomialIdeal):
I = Ideal(I)
self.Ring = I.ring()
self.Zeroes = 0 # S-polynomials reducing to zero.
# We do not keep track of the corresponding critical pair.
self.Rgb = Sequence([I.ring()(0)],I.ring()) # (reduced) GB of the first few input generators
self.RgbMon = self.Rgb
# TODO: is homogenization really needed?
# Certainly it is needed in the F5 matrix version. But here??
if not I.is_homogeneous():
print "homogenizing the input"
J = I.homogenize()
self.homog = 0
else:
J = I
self.homog = 1
if J.gens()==(0,):
self.Inputgen = [J.ring()(0)] # reduced_basis would yield a segfault
else:
self.Inputgen = list(J.reduced_basis()) # the ideal generators we got, sorted increasingly
self.HRing = J.ring()
self.Inputgen.sort()
self.Rules = {} # rewriting rules
self.Basis = [] # List of decorated polynomials that will eventually form a
# a Groebner basis of J
self.One = J.ring()(1) # self.One.parent() =self.HRing != self.Ring, if there was a homogenisation
self.nvars = len(J.ring().gens())
return self.basis()
## Some methods that reveal the C-attributes of self.
def rules(self):
return self.Rules
def zeroes(self):
return self.Zeroes
#############
## Main procedure
## Return the previously computed basis, or do the computation when first called
def basis(slf, maxgen=None):
"""
F.basis(): Compute a reduced Groebner basis using the F5 algorithm.
"""
cdef F5 self = slf
if self.Basis:
if self.homog:
return self.Rgb # or [X.poly for X in self.Basis]*self.HRing, currently without reduction
else:
return dehomogenize(self.Rgb,self.Ring)
if not self.homog:
h = self.HRing('h')
if self.Inputgen[0] == 0:
return Sequence([self.HRing(0)], self.HRing)
cdef DecoratedPolynomial newDP = DecoratedPolynomial(self.One,0,self.Inputgen[0]/self.Inputgen[0].lc())
self.Rules[0]=[]
self.Basis.append(newDP)
cdef int i
cdef int UpTo
if maxgen is None:
UpTo = len(self.Inputgen)
else:
UpTo = maxgen
cdef int laststop = 0
cdef DecoratedPolynomial X,Y
for i in range(1,UpTo):
print "Generator %d/%d"%(i+1,UpTo)
# by induction, we have a reduced GB for the generators 1,...,i-1.
self.Rgb = ([Y.p for Y in self.Basis]*self.HRing).reduced_basis()
self.RgbMon = ([Y.p.lm() for Y in self.Basis]*self.HRing).reduced_basis()
self.insert(i)
for X in self.Basis[laststop:]:
# it it contains a constant polynomial, the ideal is trivial
if self.homog:
if X.p.degree()==0:
self.Rgb = Sequence([self.One], self.HRing)
return Sequence([self.Ring(1)], self.Ring)
else:
if X.p(h=1).degree()==0:
self.Rgb = Sequence([self.One], self.HRing)
return Sequence([self.Ring(1)], self.Ring)
laststop = len(self.Basis)
# Now, there only remains to toss out the reduced GB
if self.Zeroes:
print "Note that we encountered %d S-polynomials reducing to zero"%(self.Zeroes)
self.Rgb = ([Y.p for Y in self.Basis]*self.HRing).reduced_basis()
if self.homog:
return self.Rgb
else:
return dehomogenize(self.Rgb,self.Ring)
cdef insert(F5 self,int i):
"""F.insert(i): Compute a Groebner basis for the ideal spanned by generators #0,...,i,
when a Groebner basis for generators #0,...,i-1 is already known.
"""
# We may reduce the i-th generator f_i by everything that is in the "old" ideal (f_0,...,f_{i-1})
f = self.Inputgen[i].reduce(self.Rgb)
self.Rules[i]=[]
if f.lc()==0:
return
cdef DecoratedPolynomial newDP = DecoratedPolynomial(self.One, i, f/f.lc())
cdef int newindices = len(self.Basis) # self.Basis[:newindices] was our starting point
self.Basis.append(newDP)
cdef int j,d
cdef list P
cdef list Pd,tmpL,X
cdef DecoratedPolynomial Y,Z
# generate new S-polynomials
# The argument i is used for the F5 criterion (it appears in the signature of
# polynomials we want to reduce with)
P=self.crit_pairs(newDP,newindices,i,newindices) # only crit pairs involving newDP and old GB elements
if not self.homog:
h = self.HRing('h')
while P:
#print i
#print "len P =",len(P)
d = min([X[0].degree() for X in P])
tmpL,Pd = [],[]
for X in P:
if X[0].degree() == d:
Pd.append(X)
else:
tmpL.append(X)
print "%d critical pairs in degree %d"%(len(Pd),d)
P = tmpL
S = self.SPolys(Pd)
R = self.reduction(S, i, newindices)
for Y in R:
self.Basis.append(Y)
if self.homog:
if Y.p.degree()==0:
return
else:
if Y.p(h=1).degree()==0:
return
P.extend(self.crit_pairs(Y,len(self.Basis),i,newindices)) # crit pairs for all of self.Basis (including the new ones)
cdef list crit_pair(self, DecoratedPolynomial X, DecoratedPolynomial Y, int i, int newindices):
"""
F.crit_pair(X,Y,i,newindices) returns (t,u_X,X,u_Y,Y) corresponding to a critical pair
X,Y when necessary for the computation of a Groebner basis of (f_1,...,f_i). Here,
the F5 criterion for label i is used. By switching X and Y, if necessary, we will
have XY.nu) or ((X.nu==Y.nu) and (umuX>umuY)):
return [[t,Y,u_Y,X,u_X]]
return [[t,X,u_X,Y,u_Y]]
cdef list crit_pairs(self, DecoratedPolynomial X, int ub, int i, int newindices):
"compare with crit_pair, but here Y runs on self.Basis[:ub]"
cdef list CP = []
cdef DecoratedPolynomial Y
for Y in self.Basis[:ub]:
CP.extend(self.crit_pair(X,Y,i,newindices))
return CP
cdef list SPolys(self,list P):
"""
F.SPolys(P), P a list of critical pairs, given by a 5-tuple of
1) a monomial t (the least common multiple of the leading monomials of X.poly() and Y.poly()),
2) the cofactor u_X = t/lm(X.poly()) for X,
3) a decorated polynomial X,
4) the cofactor u_Y = t/lm(X.poly()) for Y,
5) a decorated polynomial Y.
Moreover, (u_X*X)<(u_Y*Y).
The S-polynomials are computed and tested for rewritability. The non-rewritable S-polynomials
are added to the rewrite rules and are returned in a list.
"""
# P is a list of tuples (lcm, u_X,X, u_Y,Y), with u_X*lm(X) = u_Y*lm(Y) = lcm
P.sort() # which means increasing by lcm
cdef list S = []
cdef DecoratedPolynomial X,Y,Z
for (t,X,u_X,Y,u_Y) in P:
if (not self.is_rewritable(u_X,X)) and (not self.is_rewritable(u_Y,Y)):
s = (u_X*X.p*Y.p.lc() - u_Y*Y.p*X.p.lc()) # it will be reduced soon, in a different method...
if s!=0:
# We know that u_X*X is on top of u_Y*Y.
# Hence, thinking abouz the F5 matrix version, we would
# create rows labeled (X.mu*u_X,X.nu) and (Y.mu*u_Y,Y.nu),
# and a reduction from top to bottom would yield
# the S-polynomial at (Y.mu*u_Y,Y.nu).
# TODO: Is this thinking correct?
Z = DecoratedPolynomial(Y.mu*u_Y,Y.nu, s/s.lc())
#Z.hist.append(('S-poly',X,X.sig(),Y,Y.sig()))
self.add_rule(Z)
S.append(Z)
S.sort() # by increasing signature, hence, from top to bottom in the F5 matrix version
return S
cdef list reduction(self,list S, int i, int newindices):
"""
F.reduction(S,i,newindices): S a list of decorated polynomials, and F.Basis[:newindices] (a private
attribute) provides a Groebner basis for the ideal generatd by the generators #0,...,i-1.
It is assumed that the list is sorted by increasing signature.
We will reduce each element of S by the GB and the elements of S of lower signature.
Returns the list of completely reduced decorated polynomials.
"""
cdef DecoratedPolynomial X,Y
cdef list truncGB = list(self.Rgb)
# We use self.Rgb, since it is a reduced GB for the ideal of the first i-1 indices.
# Later on, we add to it new reduced S-polynomials.
cdef list completed = []
cdef list newly_completed = []
cdef list redo = []
#sys.stdout.write('r')
#sys.stdout.flush()
while S:
#sys.stdout.write('.')
#sys.stdout.flush()
X = S.pop(0) # this is the one with the smallest signature
# By construction, the signature of all our input polynomials has index i:
#sys.stdout.write('r')
#sys.stdout.flush()
X.p = X.p.reduce(truncGB)
#sys.stdout.write('.')
#sys.stdout.flush()
#X.hist.append('reduction')
if X.p.lc()!=0:
X.p = X.p/X.p.lc()
newly_completed, redo = self.top_reduction(X,newindices,completed)
completed.extend(newly_completed)
truncGB.extend([Y.p for Y in newly_completed]) # hence basis is "old basis" plus "completed"
S.extend(redo)
S.sort() # by increasing signature, since reductions must be done from top to bottom.
#print "S",len(S)
#print "rule",len(self.Rules[i])
#print "Basis", len(self.Basis)
#sys.stdout.write('\n')
#sys.stdout.flush()
return completed
cdef tuple top_reduction(self, DecoratedPolynomial X, int newindices,list completed):
if X.p.lc()==0:
print "WARNING: reduction to zero"
self.Zeroes = self.Zeroes+1
return ([],[])
#X.hist.append(('test',X.p))
#sys.stdout.write('t')
#sys.stdout.flush()
J = self.find_reductor(X,newindices,completed)
#sys.stdout.write('.')
#sys.stdout.flush()
if not J:
return ([X],[])
cdef DecoratedPolynomial Y = J[1]
u = J[0]
p = X.p-u*Y.p
#sys.stdout.write('\n')
#sys.stdout.flush()
if p.lc()!=0:
p = p/p.lc()
else:
print "WARNING: top-reduction to zero"
self.Zeroes = self.Zeroes+1
return ([],[])
newmono = Y.mu*u
# TODO: is it newmono