PolyBoRi Tutorial

Michael Brickenstein

September 9, 2009
1

1 Introduction

1.1

1.1 Interfaces

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).

1.2

1.2 Ring Declarations

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.

Example

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]

1.3

1.3 Ordering

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).

Example

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.

1.4

1.4 Arithmetic

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)

1.5

1.5 Set operations

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)}}

1.6

1.6 Gröbner bases

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:

1.6.1

1.6.1 Elimination of variables

Given Boolean generators G of an ideal I 2[x1,,xn,y1,,ym]x12 + x1,,xn2 + xn,y1,,ymwe 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

Definition 1 (Elimination orderings) Let R = 2[x1,,xn,y1,ym]. An ordering “>” is called an elimination ordering of x1,,xn, if xi > t for every monomial t in K[y1,,ym] and every i = 1,,n.

Example

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.

2

2 How to program efficiently

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

2.1

2.1 Low level friendly programming

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 x12(x22(x32(x42x5))) is much cheaper than ((((x12x2)2x3)2x4)2x5). 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.

2.2

2.2 Replace algebra by set operations

Often there is an alternative description in terms of set operations for algebraic operations, which is much faster.

2.2.1

2.2.1 Construct power sets

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).

2.2.2

2.2.2 All monomials of degree d
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).

3

3 Case study: Graded part of a polynomial

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.

3.1

3.1 Simple, algebraic solution

def simple_graded(f, d):  
    return sum((t for t in f.terms() if t.deg()==d))

This solution is obvious, but quite slow.

3.2

3.2 Low level friendly, algebraic solution

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.

3.3

3.3 Highlevel with set operations

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).

3.4

3.4 Recursive

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.

3.5

3.5 Decision-diagram style recursive implementation in PolyBoRi

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.

4

4 Case study: Evaluation of a polynomial

4.1

4.1 Substitute a single variable x in a polynomial by a constant c

Given a Boolean polynomial f, a variable x and a constant c, we want to plug in the constant c for the variable x.

4.1.1

4.1.1 Naive approach

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)])  

4.1.2

4.1.2 Solution 1: Set operations

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))

4.1.3

4.1.3 Solution 2: Linear Lexicographical Lead rewriting systems

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}∪{x21 + x1,,x2n + 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.

  1. If our system consists of exactly one Boolean polynomial, 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.

  2. A second optimization is to call just
    ll_red_nf(f,x+c)

    As {x + c} is a reduced Boolean Gröbner basis (equivalently {x + c,x21 + x1,,x2n + xn}\{x2 + x} is a reduced Gröbner basis).

4.2

4.2 Evaluate a polynomial by plugging in a constant for each variable

We want to a polynomial f(x1,,xn) by xi↦→ci, where c1,,cy are constants.

4.2.1

4.2.1 Naive approach

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)

4.2.2

4.2.2 Solution 1: n set operations

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.

4.2.3

4.2.3 Solution 2: Linear Lexicographical Lead rewriting systems

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.

4.3

4.3 General Linear Lexicographical Lead Rewriting Systems

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).

Definition 2 Let L be a list of Boolean polynomials. If all elements p of L have pairwise different leading terms with respect to lexicographical ordering, then we call L a linear lexicographical lead rewriting system.

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:

  1. A maximal subset L of the equation system, which forms a linear lexicographical lexicographical rewriting system.
  2. A normal form algorithm f, s. th. f(p) forms a reduced normal form of p against the Gröbner basis consisting of L and the field equations.
  3. A list of polynomials R, which are in reduced normal form against L, s. th. L R spans modulo field equations the same ideal as the original equation system.
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

5

5 Storing polynomial data in a file

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)]

6

6 Reinterpretation of Boolean sets as subsets of the vector space 2n

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.

6.1

6.1 Examples

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.

7

7 Lexicographical normal form of a polynomial against a variety

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 : I(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

8

8 Partial Boolean functions

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)

9

9 Building your own Gröbner basis algorithm

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:

9.1

9.1 Adding a Generator

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))

9.2

9.2 Interreduction

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

9.3

9.3 A minimal Buchberger algorithm

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.

9.4

9.4 Estimating the number of solutions

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