1 Introduction

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

In order to follow the instructions of this tutorial a basic knowledge of IPython is necessary. For instance, when interactively entering programming blocks (loops, functional definitions), the latter are closed by typing two additional line breaks.

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)]) 
Out[1]: <polybori.dynamic.PyPolyBoRi.Ring object at 0x18436f0> 
 
In [2]: r 
Out[2]: <polybori.dynamic.PyPolyBoRi.Ring object at 0x18436f0> 
 
In [3]: [r.variable(i) for i in xrange(r.n_variables())] 
Out[3]: [x(0), x(1), y(0), y(1), y(2)] 
 
In [4]: declare_ring([AlternatingBlock(["x","y"],2)]) 
Out[4]: <polybori.dynamic.PyPolyBoRi.Ring object at 0x2370b70> 
 
In [5]: [r.variable(i) for i in xrange(r.n_variables())] 
Out[5]: [x(0), y(0), x(1), y(1)] 
 
In [6]: declare_ring([HigherOrderBlock("x", (2, 3))]) 
Out[6]: <polybori.dynamic.PyPolyBoRi.Ring object at 0x26a86d8> 
 
In [7]: [r.variable(i) for i in xrange(r.n_variables())] 
Out[7]: [x(0, 0), x(0, 1), x(0, 2), x(1, 0), x(1, 1), x(1, 2)] 
 
In [8]: declare_ring(["x","y","z"]) 
Out[8]: <polybori.dynamic.PyPolyBoRi.Ring object at 0x2eb4630> 
 
In [9]: [r.variable(i) for i in xrange(r.n_variables())] 
Out[9]: [x, y, z]

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)]) 
 
In [2]: x(0)>x(1) 
Out[2]: True 
 
In [3]: x(0)>x(1)*x(2) 
Out[3]: True 
 
In [4]: r = r.clone(ordering=dlex) 
 
In [5]: r(x(0)) > r(x(1)) 
Out[5]: True 
 
In [6]: r(x(0)) > r(x(1)*x(2)) 
Out[6]: False 
 
In [7]: r = r.clone(ordering=dp_asc) 
 
In [8]: r(x(0)) > r(x(1)) 
Out[8]: False 
 
In [9]: r(x(0)) > r(x(1)*x(2)) 
Out[9]: False 
 
In [10]: r = r.clone(ordering=block_dlex, blocks=[10]) 
 
In [11]: r(x(0)) > r(x(1)) 
Out[11]: True 
 
In [12]: r(x(0)) > r(x(1)*x(2)) 
Out[12]: False 
 
In [13]: r(x(0)) > r(y(0)*y(1)*y(2)) 
Out[13]: True 
 
In [14]: r = r.clone(ordering=block_dp_asc) 
 
In [15]: r(x(0)) > r(x(1)) 
Out[15]: False 
 
In [16]: r(x(0)) > r(y(0)) 
Out[16]: False 
 
In [17]: r(x(0)) > r(x(1)*x(2)) 
Out[17]: False 
 
In [18]: r = r.clone(ordering=block_dp_asc, blocks=[10]) 
 
In [20]: r(x(0)) > r(y(0)) 
Out[20]: True 
 
In [21]: r(x(0)) > r(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]: r.clone(ordering=block_dp_asc, blocks=block_start_hints) 
 
In [3]: block_start_hints 
Out[3]: [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)]) 
 
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]: r=r.clone(ordering=dp_asc) 
 
In [5]: for t in r(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 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, r)+Polynomial(1, r) 
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 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 BooleSets 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.else_branch() 
 
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 (...) 
 
/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().minimal_elements() 
Out[3]: {{x(1),x(3)}, {x(2)}, {x(4)}}

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 exchangeable heuristic. In principle, there exist standard settings, but – in case an option is not provided explicitly 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 explicitly 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 Parallelization

With a view towards parallel computing PolyBoRi also has a preliminary support for concurrent processing of different strategies. Note, this needs the pyprocessing module to be installed. For this purpose the function
groebner_basis_first_finished had been established. Very much like the command groebner_basis the first argument is assumed to be the system of equations, which is given as a list of Boolean polynomials. Further arguments are the options for the different strategies (in form of dictionaries with keyword arguments for groebner_basis).

It returns the result of the variant, which terminates first.

 
In [1]: from polybori.parallel import groebner_basis_first_finished 
 
In [2]: groebner_basis_first_finished([x(2)*x(1), x(1)+1], 
   {'heuristic':False}, {'heuristic': True}) 
Out[2]: [x(2), x(1) + 1]

The two option sets – with and without heuristic – present a reasonable choice, in the case that two CPUs are available, but there is not much known about the ideal. Usually the computation time for these settings differ very much: the heuristic is designed to select the apparently best settings in a very progressive way. Indeed, it is not always possible to decide without further experiments, which option setting would lead to smallest computation time. Hence, any a priori choice will lead to a suboptimal performance (like all attempts in this area). Because there is much experience behind those heuristic strategies, it will still help to achieve good results in many cases. Also, it yields the best out of the box experience (as far as it is possible here). Disabling the heuristics will usually result in an algorithm, which is much closer to the original Buchberger algorithm, but without the overhead of the heuristic functionality itself.

Another important set of alternative option sets consists of computing with and without (dense) linear algebra.

 
In [1]: from polybori.parallel import groebner_basis_first_finished 
 
In [2]: groebner_basis_first_finished([x(2)*x(1), x(1)+1], 
   ...:   {'faugere':False, 
   ...:   'linear_algebra_in_last_block': False}, 
   ...:   {'faugere':True, 
   ...:   'linear_algebra_in_last_block': True}) 
Out[2]: [x(2), x(1) + 1]

This parallelization, i. e. running different strategies in parallel, seems to be promising, because it potentially provides a nonlinear speedup (higher factor than number of CPUs). This can also be observed in other domains: for instance, current SAT-solvers make use of this by choosing different random seeds.

On the other hand, this approach does not contain any cooperation between the threads/processes. At best it is just as fast as the best variant. When studying a particular kind of systems, there is usually a good guess in advance for the best strategy. In this setting, it would be better to use all available CPUs in order to apply the guessed variant to suitable subproblems. The used data structure in PolyBoRi and the functional approach for handling them, seems to be suited to a cooperative shared memory parallelization. However, the underlying library CUDD does not provide any support for that. In this way, this kind of parallelization, remains to be a challenge for the future. While the ZDD operations cannot be parallelized easily, the linear algebra backend using M4RI could make use of several CPUs (when compiled with appropriate settings).

1.6.2 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)]) 
Out[1]: <polybori.dynamic.PyPolyBoRi.Ring object at 0x1848b10> 
 
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 [2]: r = r.clone(ordering=block_dp_asc, blocks=block_start_hints) 
 
In [4]: G=[r(poly) for poly in G] 
 
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 2.5.3.