LibPoly is a C library for computing with polynomials. It provides data structures to represent multivariate polynomials, and algorithms ranging from simple arithmetic and GCD computation, to root isolation and computation with algebraic numbers. The goal of the library is to be simple and extensible, and it is targeting tools that reason in nonlinear arithmetic. LibPoly is already successfully used in the Yices 2 nonlinear arithmetic solver. The library is freely available under a permissive open-source license. We present the basic functionality of LibPoly through the Python interface, and we describe a simple implementation of a classical cylindrical algebraic decomposition algorithm.
The interface to the library is designed to support both the traditional, purely symbolic algorithms, and the more recent model-driven methods. Therefore, LibPoly supports operations not usually considered in traditional symbolic computer algebra systems, such as a variable order that can change freely between various operations.
To present the functionality of LibPoly, we go through a simple implementation of a
cylindrical algebraic decomposition (CAD) algorithm [
We will go through a Python implementation of the CAD construction algorithm by relying on LibPoly's Python bindings. We rst demonstrate the basics of polynomial construction and manipulation. Then we show how to isolate and manipulate roots of univariate polynomials and construct the sign table of a univariate polynomial. Finally, we extend the sign-table construction to the multivariate case, by relying on CAD projections, resulting in a simple implementation of CAD construction.3 2.1
Variables. A variable in LibPoly is an object representing real-valued variables. Internally, it has an associated printable name and a unique numerical ID. The following example imports the LibPoly library in Python and creates three variables x, y, and z. 1 import polypy # Import the library 2 x = polypy . Variable ( 'x ') # Variable x 3 [y , z] = [ polypy . Variable (s) for s in [ 'y ', 'z ']] # Variables y and z
Variable Order. An important feature of LibPoly is a exible order over the variables. By default, variables are ordered by numerical ID's (i.e., in order of creation), but this order can be changed dynamically. The order is updated by manipulating the order list that changes the variable order according to the following semantics: variables in the order list are smaller than the variables not in the order list; variables in the order list are ordered according to the list order; and variables not in the order list are ordered according to their numerical ID's.
The following example illustrates the order manipulation. 1 order = polypy . variable_order # order = [] , x < y < z 2 order . push (z) # order = [z] , z < x < y 3 order . push (y) # order = [z , y], z < y < x 4 order . pop () # order = [z] , z < x < y 5 order . push (x) # order = [z , x], z < x < y
2The presentation focuses on the Python API for clarity, but the library is implemented in C and provides a C API with equivalent functionality.
3Complete examples are available at https://github.com/SRI-CSL/libpoly/tree/master/examples. Polynomials. Polynomials in LibPoly always have integer coe cients. They can be created and combined using the standard Python arithmetic operators over variables and integer constants. The following example creates two polynomials f = (x2 2x + 1)(x2 2) and g = z(x2 y2). 1 f = (x **2 - 2* x + 1) *( x **2 - 2) # Univariate polynomial 2 g = z *( x **2 - y **2) # Multivariate polynomial
The internal representation of polynomials is recursive: a polynomial is a list of coe cients in its topmost variable, where each coe cient is a polynomial in the remaining variables. This representation is updated dynamically to respect the variable order. For example, if the variable order is x < y < z, then g is a polynomial in Z[x; y; z], with z as top variable. Therefore, g is represented as g = ( y2 + x2) z. On the other hand, if the variable order is z < y < x, then x is the top variable, and g = z x2 + ( z) y2. 1 order . set ([x , y , z ]) # z is the top variable , i.e. g in Z[x ,y ,z] 2 print g # Out : ( -1* y **2 + (1* x **2) )*z 3 order . set ([z , y ]) # x is the top variable , i.e. g in Z[z ,y ,x] 4 print g # Out : (1* z)*x **2 + (( -1* z)*y **2)
Several operations give access to this recursive structure: The var() method returns the top variable of the polynomial in the current order, the degree() method returns the degree of the polynomial in its top variables, and the coefficients() method returns the list of polynomial coe cients with respect to its top variable (in decreasing degree order). The derivative () method computes the derivative of the polynomial with respect to its top variable. The factor_square_free() method returns a factorization of the polynomial where no factor is a square and the factors are pairwise co-prime.4 1 g. var () # Out : Variable ('x ') 2 g. degree () # Out : 2 3 g. coefficients () # Out : [( -1* z)*y **2 , 0, 1* z] 4 g. derivative () # Out : (2* z)*x 5 f. factor_square_free () # Out : [(1* x **2 - 2, 1) , (1* x - 1, 2) ] 6 g. factor_square_free () # Out : [(1*z , 1) , (1* x **2 + ( -1* y **2) , 1) ]
In the example above, f is factored (over Z) into irreducible polynomials (x2 2)(x polynomial g is factored into g = z:(x2 y2). This is not a full factorization as (x2 be further decomposed but the factorization is square free. 1)2. The y2) could 2.2
Root Isolation. When analyzing the behavior of a polynomial f 2 Z[~x; y], the basic operation
is nding the roots of f . A root of f , given an assignment of variables ~x to values ~a, is a
value b such that f (~a; b) = 0. The ordered list of roots of f can be obtained by calling the
roots_isolate() method on the polynomial f .
1 m = polypy . Assignment ()
2 r = f. roots_isolate (m)
3 print r [0] # Out : <1*x **2 + ( -2) , ( -3/2 , -5/4) >
4 print r [
In the example above, we rst create an empty Assignment object, and then obtain the roots of f = (x2 2)(x 1) 2 Z[x]. The result is the list p2; 1; p2 . LibPoly displays p2 and p2 as algebraic numbers as explained in the following paragraph.
4This is not a full factorization, but it is su cient for most purposes where factoring is useful. Complete
factorization, although potentially useful, is not necessary for CAD and is currently not available in LibPoly.
Working with Values. The real values (including roots) in LibPoly are represented as value
objects and are either integers, dyadic rationals5, or algebraic numbers such as p2. Algebraic
numbers are represented symbolically as pairs hf; Ii, where f is a univariate polynomial f , and
I is an real interval in which f has a unique root. The algebraic number hf; Ii is the root of f
that occurs in interval I. In the previous example, the algebraic number p2 is represented as
the pair hx2 2; ( 45 ; 32 )i.
1 r [0]. to_double () # Out : -1.41421356237
2 r [0] < r [
The snippet above illustrates the basic operations on real values. Real values can be converted to oating point with the to_double() method, and they can be compared to each other with the usual comparison operators in Python. Given two values v1 and v2, method get_value_between() returns the a value in the interval (v1; v2). The method picks the \simplest" possible value, namely a dyadic rational in the interval with the smallest value of m. In particular, the method will return an integer if one exists in (v1; v2). Method set_value() assigns a value to a variable, while unset_value() clears the value. Method sgn() computes the sign of a polynomial at a point given by a full assignment.
Constructing a Sign Table. We now have all the ingredients to compute sign tables of univariate polynomials. For a set of univariate polynomials F a sign table of F is a decomposition of R into intervals I1; : : : ; In, such that, in each interval Ik, the polynomials in F are sign-invariant on Ik (i.e., they don't change sign). To construct a sign table of F , we rst isolate the roots of all f 2 F , then we sort the computed roots by increasing value, and, nally, we evaluate the sign of the polynomials in the intervals delimited by the roots. In each such interval, Ik, the polynomials of F do not change sign, so we can evaluate the sign of each f 2 F on Ik by evaluating the sign of f at a single evaluation point of Ik.
A Python function that constructs a sign table using LibPoly is presented in Figure 1. Applying this function to polynomials f = x2 2 and g = x2 3 will result in a sign table that resembles the following.
f g ( 1; 1 1 p3) p3
In the univariate case, building a sign table is a complete method for determining if a set of polynomial constraints has a solution. For example, we can solve the constraint (f > 0)^(g < 0) by examining the sign table. It is clear that the solutions are in the set ( p3; p2) [ (p2; p3). Moreover, in order to perform such an analysis, we only need to keep track of one evaluation point per interval. Figure 2 plots the two polynomials f and g of our example, and shows the evaluation points (red dots). 2.3
A sign table that we computed in the previous section decomposes R into nitely many regions. Each region is either an open interval or a point, and in each region, the polynomials of F are
5Rationals of the form 2nm . 1 def sign_table (x , polys , m): 2 # Get the roots and print the header 3 roots = set () # Set of roots 4 output (" poly / int ") 5 for f in polys : 6 output (f) 7 f_roots = f.roots isolate(m) 8 roots . update ( f_roots ) 9 stdout . write ("\n") 10 # Sort the roots and add infinities 11 roots = [ polypy . INFINITY_NEG ] + sorted(roots) + [ polypy . INFINITY_POS ] 12 # Print intervals and signs in the intervals 13 root_i , root_j = itertools . tee ( roots ) 14 next ( root_j ) 15 for r1 , r2 in itertools . izip ( root_i , root_j ): 16 output (( r1 . to_double () , r2 . to_double () )) 17 # The interval (r1 , r2 ) 18 v = r1.get value between(r2); using the Python built-in sorted() function that compares values using value comparison functionality from LibPoly. x**2 - 2 x**2 - 3 2 1 0 -1 -2 -3 -3 -2 -1 1 2
3
0
x
sign-invariant. The sign table completely characterizes the behavior of F on R. A natural
question arises: can we extend the concept of sign table (with evaluation points) to the multivariate
case, so that we can solve polynomials constraints with several variables ? The answer to this
question is positive, and the generalization of sign tables to Rn is called cylindrical algebraic
decomposition (CAD) [
CAD Lifting. A nave attempt at constructing a CAD (sign table) would be to extend the sign-tabling as follows. Let F Z[x1; : : : ; xn] be a set of polynomials. We partition F according to the top variable of each polynomial:
F = F1 [ F2 [ [ Fn : Every polynomial f 2 Fk contains only variables x1; : : : ; xk, and xk is the top variable of f . With this setup, we can directly extend the sign_table() function of Figure 1 into a recursive procedure that builds the CAD of F one variable at a time. We do so by building a sign table T1 for polynomials in F1. Then, for each evaluation point in T1, we build a sign table T ;2 for polynomials in F2, and so on. In CAD terminology this procedure is called lifting.
The function lift() from Figure 3 takes as input a map poly_map that maps xk to Fk as above, and a list vars of variables [x1; : : : ; xn]. The function outputs the computed evaluation points. It shouldn't be too hard to see that this lifting procedure on its own does not produce 1.5 0.5 2 1 0 -0.5 -1.5 -1 -2 y 1
1 lift x1: Since F1 is empty, the set of roots for x1 is also empty, and we can choose any value in ( 1; +1) as evaluation point for x1. For example, let x1 7! 2 be the evaluation point for x1. lift x2: F2 = ff g, and we search for the roots of the polynomial f when x1 7! 2. This amounts for nding the roots of f (2; x2) = x22 + 3. Since this polynomial does not have roots in x2, we can, as before, pick any evaluation point. We choose x2 7! 0.
In the above example, the procedures computes a single evaluation point, namely, (x1; x2) = (2; 0). At this point, f has positive sign/ This single evaluation point clearly does not capture the behavior of f over R2, since f is not uniformly positive.
CAD Projection. The above example shows that building sign tables variable by variable is
not in itself su cient to construct a proper sign-invariant decomposition for a set of polynomials
F . The problem lies in the fact that the lower-level polynomial sets Fk do not capture enough
information about the behavior of the whole set F , and the sign-tabling can miss important
evaluation points. We can x this, while keeping the same overall procedure, by adding carefully
crafted polynomials to the sets Fk. These additional polynomials will \guide" the construction
of sign tables to cover all possible behaviors of the original set F . In the seminal paper on CAD
construction [
De nition 2.1 (Collins Projection). Given a set of polynomials F Z[~y; x], the Collins projector operator Pc(F ; x) is de ned as [
coe (f; x) [ f2F [ f2F g2R (f;x) psc(g; gx0; x) [
psc(gi; gj ; x) : [ i<j gi2R (fi;x) gj2R (fj;x)
The projection operator above, when applied to F , introduces enough additional
polynomials that capture the necessary missing behavior. The individual operations in Pc are advanced
polynomial operations (for more details see, e.g., [
Before de ning the projection function, we rst de ne two utility function. First, we de ne a function that, given a polynomial f , returns all non-constant reductums of f in variable x. 1 # Returns reductum closure of f , excluding constants 2 def get_reductums (f , x): 3 R = [] 4 while f.var() == x: R. append (f); f = f.reductum() 5 return R
Additionally, instead of adding polynomials directly to the sets Fk, we will add them trough the following proxy function that performs factorization rst. The square-free factorization simpli es the polynomials and ensures that some unnecessary complexity is removed. We are now ready to de ne our overall projection function project(). This function projects the sets Fn; : : : ; F1 in sequence by applying the projection operator Pc to each of them. 1 # Add polynomials to projection map 2 def add_polynomial ( poly_map , f): 3 # Factor the polynomial f 4 for ( f_factor , multiplicity ) in f.factor square free(): # Add non - constant polynomials to poly_map if (f factor.degree() > 0) : x = f factor.var() if (x not in poly_map ): poly_map [x] = set () poly_map [x ]. add ( f_factor ) # Add a collection of polynomials to projection map def add_polynomials ( poly_map , polys ):
for f in polys : add_polynomial ( poly_map , f)
1 # Project : go down the variable stack and project
2 def project ( poly_map , vars ):
3 for x in reversed ( vars ):
4 # Project variable x
5 x_polys = poly_map [x]
6 # [
add_polynomial ( poly_map , g1.psc(g2)) eAxpapmlypilneg, rtehseultpsroinjeFctb(e)infguenncltaiorgnedtotothFe =poflyxn21om1i;axl21s+etxF22 =1gf.xT12h+e xn22ewly1ga,ddfreodmpothlyenofamilieadl adds two critical roots to the sign table of x1, which allows the lifting to succeed in CAD construction.
CAD Construction. All the tools are now in place for a full CAD construction. The main CAD construction routine cad() takes a list of polynomials and outputs the evaluation points corresponding to the CAD of the polynomials. The CAD is constructed by rst projecting the input polynomials, and then constructing the CAD evaluation points over the projects set using the lifting procedure. 1 # Run the CAD construction 2 def cad ( polys , vars ): 3 polypy . variable_order . set ( vars ) 4 poly_map = {} 5 add_polynomials ( poly_map , polys ) 6 project ( poly_map , vars ) 7 lift ( poly_map , vars )
Applying the procedure above on our example set F = fx2 + y2 points depicted in Figure 4. 1g results in the evaluation 3
LibPoly is freely available on GitHub7 under the LGPL license, with more information
available at the LibPoly web page.8 The implementation is about 15KLOC of C code, with the
only external dependency being the widely used GMP library [
Although the library is e cient9, there are many parts of the library that could bene t
from algorithmic improvements or new features. For example, more e cient algorithms for
polynomial factorization and GCD computation, and more e cient computation with algebraic
numbers (and support of in nitesimal numbers [