PolyBoRi TutorialMichael Brickenstein |
The core of PolyBoRi is a C++ library. On top of it there exists a Python interface. Additionally to the Python interface a integration into SAGE was provided by Burcin Erocal. The main difference is, that PolyBoRi’s built-in Python interface makes use of the boost library, while the SAGE interface relies on Cython. However the wrappers for SAGE and the original Python interface are designed, that it is possible to run the same code under both bindings.
We provide an interactive shell for PolyBoRi using ipython for the SAGE interface (which is invoked the command sage) as well as for the built-in one, which can be accessed by typing ipbori at the command line prompt.
In ipbori a global ring is predefined and a set of variables called x(0), …, x(9999). The default ordering is lexicographical ordering (lp).
While in ipbori usually a standard ring is predefined,
it is possible to have more advanced ring declarations.
The declare_ring
-function takes two arguments.
The first argument is a list of variables or blocks or variables, and the second is a dictionary where the ring and the variable declarations are written to. The ring always get the name r in that dictionary (the best choice is to use the dictionary of global variables). In addition to that the ring is returned.
In [1]: declare_ring([Block("x",2),Block("y",3)],globals()) Out[1]: <polybori.dynamic.PyPolyBoRi.Ring object at 0x18436f0> In [2]: r Out[2]: <polybori.dynamic.PyPolyBoRi.Ring object at 0x18436f0> In [3]: [Variable(i) for i in xrange(r.nVars())] Out[3]: [x(0), x(1), y(0), y(1), y(2)] In [4]: declare_ring([AlternatingBlock(["x","y"],2)],globals()) Out[4]: <polybori.dynamic.PyPolyBoRi.Ring object at 0x2370b70> In [5]: [Variable(i) for i in xrange(r.nVars())] Out[5]: [x(0), y(0), x(1), y(1)] In [6]: [Variable(i) for i in xrange(r.nVars())] Out[6]: [x(0, 0), x(0, 1), x(0, 2), x(1, 0), x(1, 1), x(1, 2)] In [7]: declare_ring(["x","y","z"],globals()) Out[7]: <polybori.dynamic.PyPolyBoRi.Ring object at 0x2eb4630> In [8]: [Variable(i) for i in xrange(r.nVars())] Out[8]: [x, y, z]
The monomial ordering can be changed by calling change_ordering (Ordering.code), where code can be either lp (lexicographical ordering), dlex (degree lexicographical ordering), dp_asc (degree reverse lexicographical ordering with ascending variable ordering), block_dlex or block_dp_asc (for ordering composed out of blocks in the corresponding ordering). When using block ordering, after changing to that ordering, blocks have to be defined using the append_ring_block function.
In contrast to the lexicographical, degree lexicographical ordering, and the degree reverse lexicographical ordering in Singular, our degree reverse lexicographical ordering has a reverse variable order (the first ring variable is smaller than the second, the second smaller than the third). This is a result of the fact, that efficient implementation of monomial orderings using ZDD structures is quite difficult (and the performance depends on the ordering).
In [1]: r=declare_ring([Block("x",10),Block("y",10)],globals()) In [2]: x(0)>x(1) Out[2]: True In [3]: x(0)>x(1)*x(2) Out[3]: True In [4]: change_ordering(dlex) In [5]: x(0)>x(1) Out[5]: True In [6]: x(0)>x(1)*x(2) Out[6]: False In [7]: change_ordering(dp_asc) In [8]: x(0)>x(1) Out[8]: False In [9]: x(0)>x(1)*x(2) Out[9]: False In [10]: change_ordering(block_dlex) In [11]: append_ring_block(10) In [12]: x(0)>x(1) Out[12]: True In [13]: x(0)>x(1)*x(2) Out[13]: False In [14]: x(0)>y(0)*y(1)*y(2) Out[14]: True In [15]: change_ordering(block_dp_asc) In [16]: x(0)>x(1) Out[16]: False In [17]: x(0)>y(0) Out[17]: False In [18]: x(0)>x(1)*x(2) Out[18]: False In [19]: append_ring_block(10) In [20]: x(0)>y(0) Out[20]: True In [21]: x(0)>y(0)*y(1) Out[21]: True
In this example, we have an ordering composed of two blocks, the first with ten variables, the second contains the rest of variables y(0), … y(9) (per default indices start at 0).
Even, if there is a natural block structure, like in this example, we have to explicit call append_ring_block
in a block ordering to set the start indices of these blocks.
This can be simplified using the variable block_start_hints
created by our ring declaration.
In [1]: declare_ring([Block("x",2),Block("y",3),Block("z",2)],globals()) Out[1]: <polybori.dynamic.PyPolyBoRi.Ring object at 0x1848b10> In [2]: change_ordering(block_dp_asc) In [3]: for b in block_start_hints: ...: append_ring_block(b) ...: ...: In [4]: block_start_hints Out[4]: [2, 5]
Another important feature in PolyBoRi is the ability to iterate over a polynomial in the current monomial ordering.
In [1]: r=declare_ring([Block("x",10),Block("y",10)],globals()) In [2]: f=x(0)+x(1)*x(2)+x(2) In [3]: for t in f.terms(): print t x(0) x(1)*x(2) x(2) In [4]: change_ordering(dp_asc) In [5]: for t in f.terms(): print t x(1)*x(2) x(2) x(0)
This is a nontrivial functionality, as the natural order of paths in ZDDs is lexicographical.
Basic arithmetic is provided in the domain of Boolean polynomials. Boolean Polynomial polynomials are polynomials over ℤ2 where the maximal degree per variable is one. If exponents bigger than one per variable appear reduction by the field ideal (polynomials of the form x2+x) is done automatically.
In [1]: Polynomial(1)+Polynomial(1) Out[1]: 0 In [2]: x(1)*x(1) Out[2]: x(1) In [3]: (x(1)+x(2))*(x(1)+x(3)) Out[3]: x(1)*x(2) + x(1)*x(3) + x(1) + x(2)*x(3)
In addition to polynomials PolyBoRi implements a data type for sets of monomials, called BooleSet. This data type is also implemented on the top of ZDDs and allows to see polynomials from a different angle. Also, it makes high-level set operations possible, which are in most cases faster than operations handling individual terms, because the complexity of the algorithms depends only on the structure of the diagrams.
Polynomials can easily be converted to BooleSet s by using the
member function set()
.
In [1]: f=x(2)*x(3)+x(1)*x(3)+x(2)+x(4) In [2]: f Out[2]: x(1)*x(3) + x(2)*x(3) + x(2) + x(4) In [3]: f.set() Out[3]: {{x(1),x(3)}, {x(2),x(3)}, {x(2)}, {x(4)}}
One of the most common operations is to split the set into cofactors of a variable. This illustrates the following example.
In [4]: s0=f.set().subset0(x(2).index()) In [5]: s0 Out[5]: {{x(1),x(3)}, {x(4)}} In [6]: s1=f.set().subset1(x(2).index()) In [7]: s1 Out[7]: {{x(3)}, {}} In [8]: f==Polynomial(s1)*x(2)+Polynomial(s0) Out[8]: True
Even more low-level than operation with subset-methods is the use of navigators. Navigators provide an interface to diagram nodes, accessing their index as well as the corresponding then- and else-branches.
In [1]:f=x(1)*x(2)+x(2)*x(3)*x(4)+x(2)*x(4)+x(3)+x(4)+1 In [2]:s=f.set() In [3]:s Out[3]:{{x(1),x(2)}, {x(2),x(3),x(4)}, {x(2),x(4)}, {x(3)}, {x(4)}, {}} In [4]:x(1).index() Out[4]:1 In [5]:s.subset1(1) Out[5]:{{x(2)}} In [6]:s.subset0(1) Out[6]:{{x(2),x(3),x(4)}, {x(2),x(4)}, {x(3)}, {x(4)}, {}} In [7]:nav=s.navigation() In [8]:BooleSet(nav,s.ring()) Out[8]:{{x(1),x(2)}, {x(2),x(3),x(4)}, {x(2),x(4)}, {x(3)}, {x(4)}, {}} In [9]:nav.value() Out[9]:1 In [10]:nav_else=nav.elseBranch() In [11]:nav_else Out[11]:<polybori.dynamic.PyPolyBoRi.CCuddNavigator object at 0xb6e1e7d4> In [12]:BooleSet(nav_else,s.ring()) Out[12]:{{x(2),x(3),x(4)}, {x(2),x(4)}, {x(3)}, {x(4)}, {}} In [13]:nav_else.value() Out[13]:2
You should be very careful and always keep a reference to the original object, when dealing with navigators, as navigators contain only a raw pointer as data. For the same reason, it is necessary to supply the ring as argument, when constructing a set out of a navigator.
The opposite of navigation down a ZDD using navigators is to construct new ZDDs in the same way, namely giving their else- and then-branch as well as the index value of the new node.
In [1]:f0=x(2)*x(3)+x(3) In [2]:f1=x(4) In [3]:if_then_else(x(1),f0,f1) Out[3]:{{x(1),x(2),x(3)}, {x(1),x(3)}, {x(4)}} In [4]:if_then_else(x(1).index(),f0,f1) Out[4]:{{x(1),x(2),x(3)}, {x(1),x(3)}, {x(4)}} In [5]:if_then_else(x(5),f0,f1) --------------------------------------------------------------------------- <type 'exceptions.ValueError'> Traceback (most recent call last) /home/user/PolyBoRi/<ipython console> in <module>() <type 'exceptions.ValueError'>: Node index must be smaller than top indices of then- and else-branch.
It is strictly necessary that the index of the created node is smaller than the index of the branches. But also operations on higher levels are possible, like the calculation of the minimal terms (with respect to division) in a BooleSet:
In [1]: f=x(2)*x(3)+x(1)*x(3)+x(2)+x(4) In [2]: f.set() Out[2]: {{x(1),x(3)}, {x(2),x(3)}, {x(2)}, {x(4)}} In [3]: f.set().minimalElements() Out[3]: {{x(1),x(3)}, {x(2)}, {x(4)}}
Gröbner bases functionality is available using the function groebner_basis from polybori.gbcore. It has quite a lot of options and a exchangable heuristic. In principle, there exist standard settings, but – in case an option is not provided explicitely by the user – the active heuristic function may decide dynamically by taking the ideal, the ordering and the other options into account, which is the best configuration.
In [1]: groebner_basis([x(1)+x(2),(x(2)+x(1)+1)*x(3)]) Out[1]: [x(1) + x(2), x(3)]
There exists a set of default options for groebner_basis.
They can be seen, but not manipulated via accessing groebner_basis.options
.
A second layer of heuristics is incorporated into the groebner_basis-function, to choose dynamically the best options depending on the ordering and the given ideal.
Every option given explicitely by the user takes effect, but for the other options the default may be overwritten by the script.
This behaviour can be turned off by calling
groebner_basis(I,heuristic=False).
Important options are the following:
other_ordering_first
, possible values are False or an ordering code.
In practice, many Boolean examples have very few solutions and a very easy Gröbner basis. So, a complex walk algorithm (which cannot be implemented using the PolyBoRi data structures) seems unnecessary, as such Gröbner bases can be converted quite fast by the
normal Buchberger algorithm from one ordering into another ordering.
faugere
, turn off or on the linear algebra
linear_algebra_in_last_block
, this affects the last block of block orderings and degree orderings. If it is set to True linear algebra takes affect in this block.
selection_size
, maximum number of polynomials for parallel reductions
prot
, turn off or on the protocol
red_tail
, tail reductions in intermediate polynomials, this options affects mainly heuristics. The reducedness of the output polynomials can only be guaranteed by the option redsb
minsb
, return a minimal Gröbner basis
redsb
, return a reduced Gröbner basis (minimal and all elements are tail reduced)
Given Boolean generators G of an ideal I⊂ ℤ2[x1,…,xn,y1,…,ym]/⟨ x12+x1,…,xn2+xn,y1,…,ym ⟩ we would like to compute a generating system H of Boolean polynomials, where
⟨ H⟩={p∈ I| p can be represented by a (Boolean) polynomial in ℤ2[y1,…,ym]}. |
This can be done as in the classical case despite the field equations using an elimination ordering for x1, …, xn
In [1]: declare_ring([Block("x",3),Block("y",3)],globals()) Out[1]: <polybori.dynamic.PyPolyBoRi.Ring object at 0x1848b10> In [2]: change_ordering(block_dp_asc) In [3]: for b in block_start_hints: ...: append_ring_block(b) ...: ...: In [4]: G=[x(0)*x(2)*y(0)*y(1) + y(1)*y(2) + y(1), ...: x(1)*x(2)*y(0)*y(2) + x(0)*x(1)*y(2) + y(1), ...: x(0)*x(1)*x(2)*y(1) + x(1) + y(1)*y(2)] In [5]: H=groebner_basis(G) In [6]: H Out[6]: [x(0)*y(0)*y(1) + x(0)*y(1) + y(0)*y(1) + y(1), x(1) + y(1), x(2)*y(1) + x(0)*y(1) + y(1), y(1)*y(2) + y(1)] In [7]: [p for p in H if p.set().navigation().value()>=y(0).index()] Out[7]: [y(1)*y(2) + y(1)]
For special cases elimination (depending on the formulation of the equations) elimination of (auxiliary) variables can be done much faster as can be seen in 4.3.
The goal of this section is to explain how to get most performance out of PolyBoRi using the underlying ZDD structure. This awareness can be seen on several levels
PolyBoRi is implemented as layer over a decision diagram library (CUDD at the moment).
In CUDD every diagram node is unique: If two diagrams have the same structure, they are in fact identical (same position in memory). Another observation is, that CUDD tries to build a functional style API in the C programming language. This means, that no data is manipulated, only new nodes are created. Functional programming is a priori very suited for caching and multithreading (at the moment however threading is not possible in PolyBoRi). The ite-operator is the most central function in CUDD. It takes two nodes/diagrams t, e and an index i and creates a diagram with root variable i and then-branch t, else-branch e. It is necessary that the branches have root variables with bigger index (or are constant). It creates either exactly one node, or retrieves the correct node from the cache. Function calls which come essentially down to a single ite call are very cheap.
For example taking the product x1 ⋆2 (x2⋆2 (x3⋆2 (x4⋆2 x5))) is much cheaper than ((((x1 ⋆2 x2)⋆2 x3)⋆2 x4)⋆2 x5). In the first case, in each step a single not is prepended to the diagram, while in the second case, a completely new diagram is created. The same argument would apply for the addition of these variables. This example shows, that having a little bit background about the data structure, it is often possible to write code, that looks as well algebraic as provides good performance.
Often there is an alternative description in terms of set operations for algebraic operations, which is much faster.
An example for this behaviour is the calculation of power sets (sets of monomials/polynomials containing each term in the specified variables). The following code constructs such a power set very inefficiently for the first three variables:
sum([x(1)**i1*x(2)**i2*x(3)**i3 for i1 in (0,1) for i2 in (0,1) for i3 in (0,1)])
The algorithm has of course exponential complexity in the number of variables. The resulting ZDD however has only as many nodes as variables. In fact it can be constructed directly using the following function (from specialsets.py).
def power_set(variables): variables=sorted(set(variables),reverse=True,key=key) res=Polynomial(1).set() for v in variables: res=if_then_else(v,res,res) return res
Note, that we switched from polynomials to Boolean sets. We inverse the order of variable indices for iteration to make the computation compatible with the principles in 2.1 (simple ite operators instead of complex operations in each step).
def all_monomials_of_degree_d(d,variables): if d==0: return Polynomial(1).set() if len(variables)==0: return BooleSet() variables=sorted(set(variables),reverse=True,key=key) m=variables[-1] for v in variables[:-1]: m=v+m m=m.set() def do_all_monomials(d): if d==0: return Polynomial(1).set() else: prev=do_all_monomials(d-1) return prev.cartesianProduct(m).diff(prev) return do_all_monomials(d)
We use the set of all monomials of one degree lower using the cartesian product with the set of variables and remove every term, where the degree did not increase (boolean multiplication: x2=x).
In the following we will show five variants to implement a function, that computes the sum of all terms of degree d in a polynomial f.
def simple_graded(f, d): return sum((t for t in f.terms() if t.deg()==d))
This solution is obvious, but quite slow.
def friendly_graded(f, d): vec=BoolePolynomialVector() for t in f.terms: if t.deg()!=d: continue else: vec.append(t) return add_up_polynomials(vec)
We leave it to the heuristic of the add_up_polynomials function how to add up the monomials. For example a divide and conquer strategy is quite good here.
def highlevel_graded(f,d): return Polynomial(f.set().intersect(all_monomials_of_degree_d(d,f.varsAsMonomial())))
This solution build on the fast intersection algorithm and decomposes the task in just two set operations, which is very good.
However it can be quite inefficient, when f has many variables. This can increase the number of steps in the intersection algorithm (which takes with high probability the else branch of the second argument in each step).
The repeated unnecessary iteration over all variables in f (during the intersection call in the last section) can be avoided by taking just integers as second argument for the recursive algorithm (in the last section this was intersection).
def recursive_graded(f,d): def do_recursive_graded(f,d): if f.empty(): return f if d==0: if Monomial() in f: return Polynomial(1).set() else: return BooleSet() else: nav=f.navigation() if nav.constant(): return BooleSet() return if_then_else( nav.value(), do_recursive_graded(BooleSet(nav.thenBranch(),f.ring()),d-1), do_recursive_graded(BooleSet(nav.elseBranch(),f.ring()),d)) return Polynomial(do_recursive_graded(f.set(),d))
Recursive implementations are very compatible with our data structures, so are quite fast. However this implementation does not use any caching techniques. Cudd recursive caching requires functions to have one, two or three parameters, which are of ZDD structure (so no integers). Of course we can encode the degree d by the d-th Variable in the Polynomial ring.
The C++ implementation of the functionality in PolyBoRi is given in this section, which is recursive and uses caching techniques.
// determine the part of a polynomials of a given degree template <class CacheType, class NaviType, class DegType, class SetType> SetType dd_graded_part(const CacheType& cache, NaviType navi, DegType deg, SetType init) { if (deg == 0) { while(!navi.isConstant()) navi.incrementElse(); return SetType(navi); } if(navi.isConstant()) return SetType(); // Look whether result was cached before NaviType cached = cache.find(navi, deg); if (cached.isValid()) return SetType(cached); SetType result = SetType(*navi, dd_graded_part(cache, navi.thenBranch(), deg - 1, init), dd_graded_part(cache, navi.elseBranch(), deg, init) ); // store result for later reuse cache.insert(navi, deg, result.navigation()); return result; }
The encoding of integers for the degree as variable is done implicitely by our cache lookup functions.
Given a Boolean polynomial f, a variable x and a constant c, we want to plug in the constant c for the variable x.
The following code shows how to tackle the problem, by manipulating individual terms. While this is a very direct approach, it is quite slow. The method reducibleBy gives a test for divisibility.
def subst(f,x,c): if c==1: return sum([t/x for t in f.terms() if t.reducibleBy(x)])+\ sum([t for t in f.terms() if not t.reducibleBy(x)]) else: #c==0 return sum([t for t in f.terms() if not t.reducibleBy(x)])
In fact, the problem can be tackled quite efficiently using set operations.
def subst(f,x,c): i=x.index() c=Polynomial(c)#if c was int is now converted mod 2, so comparison to int(0) makes sense s=f.set() if c==0: #terms with x evaluate to zero return Polynomial(s.subset0(i)) else: #c==1 return Polynomial(s.subset1(i))+Polynomial(s.subset0(i))
On the other hand, this linear rewriting forms a rewriting problem and can be solved by calculating a normal form against a Gröbner basis. In this case the system is {x+c} ∪ {x12+x1,…,xn2+xn} (we assume that x=xi for some i). For this special case, that all Boolean polynomials have pairwise different linear leading terms w. r. t. lexicographical ordering, there exist special functions.
First, we encode the system {x+c} into one diagram
d=ll_encode([x+c])
This is a special format representing a set of such polynomials in one diagram, which is used by several procedures in PolyBoRi. Then we may reduce f by this rewriting system
ll_red_nf_noredsb(f,d)
This can be simplified in our special case in two ways.
ll_encode
does essentially a type conversion only (and much overhead).
This type conversion can be done implicitely (at least using the
boost::python
-based interface ipbori
).So you may call
ll_red_nf_noredsb(f,x+c)
In this case, there is no need for calling ll_encode
.
The second argument is converted implicitely to BooleSet.
ll_red_nf(f,x+c)As {x+c} is a reduced Boolean Gröbner basis (equivalently {x+c,x12+x1,…,xn2+xn}\ {x2+x} is a reduced Gröbner basis).
We want to a polynomial f(x1,…, xn) by xi↦ ci, where c1,…, cy are constants.
First, we show it in a naive way, similar to the first solution above.
def evaluate(f,m): res=0 for term in f.terms(): product=1 for variable in term.variables(): product=m[variable]*product res=res+product return Polynomial(res)
The following approach is faster, as it does not involve individual terms, but set operations
def evaluate(f,m): while not f.constant(): nav=f.navigation() i=nav.value() v=Variable(i) c=m[v] if c==0: #terms with x evaluate to zero f=Polynomial(nav.thenBranch()) else: #c==1 f=Polynomial(nav.thenBranch())+Polynomial(nav.elseBranch()) return f
For example, the call
evaluate(x(1)+x(2),{x(1).index():1,x(2).index():0})
results in 1
.
We deal here with navigators, which is dangerous, because
they do not increase the internal reference count on the represented polynomial
substructure. So, one has
to ensure, that f is still valid, as long as we use a navigator on f.
But it will show its value on optimized code (e. g. PyRex), where it causes
less overhead.
A second point, why it is desirable to use navigators is, that their
thenBranch
- and elseBranch
-methods immediately return (without
further calculations) the
results of the subset0
and subset1
-functions, when the latter are
called together with the top variable of the diagram f.
In this example, this is the crucial point in terms of performance.
But, since we already call the polynomial construction on the branches,
reference counting of the corresponding subpolynomials is done anyway.
This is quite fast, but suboptimal, because only the inner functions (additions) use caching. Furthermore, it contradicts the usual ZDD recursion and generates complex intermediate results.
The same problem can also be tackled by the linear-lead routines. In the case, when
all variables are substituted by constants, all intermediate results
(generated during ll_red_nf
/ll_red_nf_noredsb
) are constant.
In general, we consider the overhead of generating the encoding d as small,
since it consists of very few, tiny ZDD operations only (and some Python overhead in the quite general ll_encode
).
d=ll_encode([x+cx,y+cy]) ll_red_nf_noredsb(f,d)
Since the tails of the polynomials in the rewriting system consist of constants only, this forms also a reduced Gröbner basis. Therefore, you may just call
ll_red_nf(f,d)
This is assumed to be the fastest way.
We used ll_red_nf
/ll_red_nf_noredsb
functions on rewriting systems, where the tails of the polynomials was constant and the leading term linear.
They can be used in a more general setting (which allows to eliminate auxiliary variable).
We know that such a system forms together with the complete set of field equations a Gröbner basis w. r. t. lexicographical ordering.
In particular we can use ll_red_nf
to speedup substitution of a variable x by a value v also in the more general case, that the lexicographical leading term of x+v is equal to x.
This can be tested most efficiently by the expression
x.set().navigation().value()>v.set().navigation().value().
In many cases, we have a bigger equation system, where many variables have a linear leading term w. r. t. lexicographical ordering (at least one can optimize the formulation of the equations to fulfill these condition). These systems can be handled by the function eliminate in the module polybori.ll. I returns three results:
In [1]: from polybori.ll import eliminate In [2]: E=[x(1)+1,x(1)+x(2),x(2)+x(3)*x(4)] In [3]: (L,f,R)=eliminate(E) In [4]: L Out[4]: [x(1) + 1, x(2) + x(3)*x(4)] In [5]: R Out[5]: [x(3)*x(4) + 1] In [6]: f(x(1)+x(2)) Out[6]: x(3)*x(4) + 1
In PolyBoRi we have default file format and tools, which read the files and generate references for our test suite.
The format is a normal Python-file with a few exceptions:
declare_ring([Block("x",4, reverse=False)]) ideal=[ x(1)+x(3), x(0)+x(1)*x(2)]
The data file can be loaded using the following commands.
In [1]: from polybori.gbrefs import load_file In [2]: data=load_file("data-sample.py") In [3]: data.ideal Out[3]: [x(1) + x(3), x(0) + x(1)*x(2)]
Let S be a Boolean set. For example:
In [1]:S=BooleSet([x(0),x(1)*x(2)]) In [2]:S Out[2]:{{x(1),x(2)}, {x(0)}}
S is a set of sets of variables. Our usual interpretation is to identify it with a polynomial with corresponding terms:
In [3]:Polynomial(S) Out[3]:x(1)*x(2) + x(0)
Another interpretation is to map a set of variables m to a vector v in ℤ2n. The i-th entry of v is set to 1, if and only if the i-th variable occurs in m. So we can identify {x0} with (1,0,0) and {x1,x2} with (0,1,1) in ℤ23. Extending this identification from sets of variables to sets of set of variables we can identify S with {(1,0,0), (0,1,1)}. Note, that the choice of n as 3 was not determined by S. In fact every bigger n would have also been a candidate. For this reason, some procedures interpreting Boolean sets as subsets of ℤ2n taking the monomial ambient space as an additional parameter. The full vector space can be constructed by multiplying all needed variables and the set of divisors of the product.
In [4]:(x(1)*x(2)*x(3)).divisors() Out[4]:{{x(1),x(2),x(3)}, {x(1),x(2)}, {x(1),x(3)}, {x(1)}, {x(2),x(3)}, {x(2)}, {x(3)}, {}}
We distingish between procedures, which use subsets of the ambient space (like finding zeros of a polynomial), and such procedures, where only the dimension/involved unit vectors/variables matter. The first kind of procedures usually gets the ambient space itself, the second kind gets the monomial consisting of all involved variables.
In [1]:f=x(0)+x(1)+x(2) In [2]:ambient_space=(x(0)*x(1)*x(2)).divisors() In [3]:ambient_space Out[3]:{{x(0),x(1),x(2)}, {x(0),x(1)}, {x(0),x(2)}, {x(0)}, {x(1),x(2)}, {x(1)}, {x(2)}, {}} In [4]:f.zerosIn(ambient_space) Out[4]:{{x(0),x(1)}, {x(0),x(2)}, {x(1),x(2)}, {}} In [5]:S=BooleSet([x(0),x(1)*x(2)]) In [6]:f.zerosIn(S) Out[6]:{{x(1),x(2)}}
A example of the second kind, where only the full ambient space can be considered count lex_groebner_basis_points
In [1]:S=BooleSet([x(0),x(1)*x(2)]) In [2]:from polybori.interpolate import * In [3]:lex_groebner_basis_points(S,x(0)*x(1)*x(2)) Out[3]:[x(0) + x(2) + 1, x(1) + x(2)] In [4]:lex_groebner_basis_points(S,x(0)*x(1)*x(2)*x(3)) Out[4]:[x(0) + x(2) + 1, x(1) + x(2), x(3)]
This function calculates the reduced lexicographical Gröbner basis of the vanishing ideal of S. Here the ambient space matters, as an additional component would mean, that the corresponding entries are zero, so we would get an additional generator for the ideal x3.
Let V be a set of points in ℤ2n, f a Boolean polynomial. V can be encoded as a BooleSet as described in 6. Then we are interested in the normal form of f against the vanishing ideal of V: (V). It turns out, that the computation of the normal form can be done by the computation of a minimal interpolation polynomial, which takes the same values as f on V.
In [1]:from polybori.interpolate import * In [2]:V=(x(0)+x(1)+x(2)+x(3)+1).set()
We take V={e0,e1,e2,e3,0}, where ei describes the i-th unit vector. For our considerations it does not play any role, if we suppose V to be embedded in ℤ24 or a vector space of higher dimension.
In [3]:V Out[3]:{{x(0)}, {x(1)}, {x(2)}, {x(3)}, {}} In [4]:f=x(0)*x(1)+x(1)+x(2)+1 In [5]:nf_lex_points(f,V) Out[5]:x(1) + x(2) + 1
In this case, the normal form of f w. r. t. the vanishing ideal of V consists of all terms of f with degree smaller or equal to 1.
It can be easily seen, that this polynomial forms the same function on V as f.
In fact, our computation is equivalent to the direct call of the interpolation function interpolate_smallest_lex
, which has two arguments: the set of interpolation points mapped to zero and the set of interpolation points mapped to one.
In [6]:z=f.zerosIn(V) In [7]:z Out[7]:{{x(1)}, {x(2)}} In [8]:o=V.diff(z) In [9]:o Out[9]:{{x(0)}, {x(3)}, {}} In [11]:interpolate_smallest_lex(z,o) Out[11]:x(1) + x(2) + 1
A partial Boolean function f is given by two disjoint set of points O, Z. f is defined to have value 1 on O, 0 on Z and is undefined elsewhere. For encoding sets of Boolean vectors we use the encoding defined in 6.
If we identify 1 with the Boolean value True and 0 with False, operations from propositional logic get a meaning for Boolean functions.
We can apply operations like xor, and, or to partial Boolean functions, defined everywhere where the result is uniquely determined on extensions of these functions.
In [1]: from polybori.partial import PartialFunction In [2]: O=BooleSet([Monomial(),x(0)*x(1)]) In [3]: Z=BooleSet([x(2),x(0)*x(2)]) In [4]: f=PartialFunction(zeros=Z,ones=O) In [5]: f Out[5]: PartialFunction(zeros={{x(0),x(2)}, {x(2)}}, ones={{x(0),x(1)}, {}}) In [6]: O2=BooleSet([x(1),x(2)]) In [7]: Z2=Bool BooleConstant BoolePolynomialVector BooleRing BooleSet In [7]: Z2=BooleSet([x(0)*x(1),x(1)*x(2),x(0)*x(2)]) In [8]: g=PartialFunction(zeros=Z2,ones=O2) In [9]: f & g Out[9]: PartialFunction(zeros={{x(0),x(1)}, {x(0),x(2)}, {x(1),x(2)}, {x(2)}}, ones={}) In [10]: f|g Out[10]: PartialFunction(zeros={{x(0),x(2)}}, ones={{x(0),x(1)}, {x(1)}, {x(2)}, {}}) In [11]: f^g Out[11]: PartialFunction(zeros={{x(0),x(2)}}, ones={{x(0),x(1)}, {x(2)}})
Since addition of in ℤ2 is equivalent to xor (using this identification with Boolean logic), the operators &
and +
coincide.
In [12]: f+g Out[12]: PartialFunction(zeros={{x(0),x(2)}}, ones={{x(0),x(1)}, {x(2)}})
We have also build our interpolation functions as method for our PartialFunction class, which is a more convenient way to use it.
In [13]: f.interpolateSmallestLex() Out[13]: x(2) + 1 In [14]: g.interpolateSmallestLex() Out[14]: x(0) + x(1) + x(2)
The central class for writing your own Gröbner bases algorithm is called GroebnerStrategy It represents a system of generators (Boolean polynomials) and contains information about critical pairs as well as some extra information like the set of leading terms belonging to these generators.
The most important operations are:
After construction several options can be set, e. g. optRedTail for tail reductions (affects also the nf method). The GroebnerStrategy keeps track not only of the single generators, but also of properties of the whole system:
There are several methods for adding a generator to a GroebnerStrategy. It may not contain two generators with the same leading monomial. In this way generators can be accessed with both their index and their leading term.
In [1]: g=GroebnerStrategy() In [2]: g.addGenerator(x(1)) In [3]: g[x(1)] Out[3]: x(1) In [4]: g.addGenerator(x(1)+1) --------------------------------------------------------------------------- ValueError Traceback (most recent call last) /Users/michael/sing/PolyBoRi/<ipython console> in <module>() ValueError: strategy contains already a polynomial with same lead
An alternative is to push the generator to the (generalized) set of critical pairs instead of adding it directly
In [5]: g.addGeneratorDelayed(x(1)+1)
Due to the absence of other pairs, in this example the polynomial is on top of the pair queue
In [6]: g.nextSpoly() Out[6]: x(1) + 1
A alternative approach is to let PolyBoRi decide, if an generator is added to the system directly or not.
In [1]: g=GroebnerStrategy() In [2]: g.addAsYouWish(x(1))
The interred-function gives back a system generating the same ideal, where no two leading terms coincide. Also, using the parameter completely ensures that no leading term divides the terms in the tails of the other generators. Even more, it is easier than the Buchberger algorithm, because no critical pairs have to be handled (actually the GroebnerStrategy applies some criteria, when adding a generator, which produces some overhead). The algorithm works by passing the generators sorted to the strategy. If a generator is (lead) rewriteable, it is rewriteable by generators with smaller leading terms. So it will be rewritten by this procedure. The algorithm stops, when no generator is lead rewriteable any more (completely = False) or rewriteable (completely = True).
def interred(l,completely=False): """Computes a new generating system (g1, ...,gn), spanning the same ideal modulo field equations. The system is interreduced: For i!=j: gi.lead() does not divide any term of gj """ l=[Polynomial(p) for p in l if not p==0] l_old=None l=tuple(l) while l_old!=l: l_old=l l=sorted(l,key=Polynomial.lead) g=GroebnerStrategy() if completely: g.optRedTail=True for p in l: p=g.nf(p) if not p.isZero(): g.addGenerator(p) l=tuple(g) return l
In this section the buchberger function from the module simplebb is presented. Unlike PolyBoRi’s more sophisticated routines this procedure was developed for educational purposes only:
def buchberger(l): "calculates a (non minimal) Groebner basis" l=interred(l) #for making sure, that every polynomial has a different leading term #needed for addGenerator g=GroebnerStrategy() for p in l: g.addGenerator(p) while g.npairs()>0: g.cleanTopByChainCriterion() p=g.nextSpoly() p=g.nf(p) if not p.isZero(): g.addGenerator(p) return list(g)
The criteria are handled by the addGenerator-method for immediately applicable criteria and by the function cleanTopByChainCriterion for the chain criterion.
In this section, it is presented, how to use the building blocks for Buchberger algorithms for other tasks like estimating the number of solutions.
First, we observe the following:
This gives a break condition for the number Buchberger algorithm. It becomes clear at a certain point of the computations, that no more than n solutions exist. However, if there are more than n solutions, the full Gröbner basis is computed by this presented algorithm.
def less_than_n_solutions(ideal,n): l=interred(ideal) g=GroebnerStrategy() all_monomials=Monomial([Variable(i) for i in xrange(number_of_variables())]).divisors() monomials_not_in_leading_ideal=all_monomials for p in l: g.addGenerator(p) while g.npairs()>0: monomials_not_in_leading_ideal=monomials_not_in_leading_ideal%g.minimalLeadingTerms if len(monomials_not_in_leading_ideal)<n: return True g.cleanTopByChainCriterion() p=g.nextSpoly() p=g.nf(p) if not p.isZero(): g.addGenerator(p) monomials_not_in_leading_ideal=monomials_not_in_leading_ideal%g.minimalLeadingTerms if len(monomials_not_in_leading_ideal)<n: return True else: return False
This document was translated from LATEX by HEVEA.