JBernstein: A Validity Checker for Generalized Polynomial Constraints Chih-Hong Cheng1 , Harald Ruess1 , and Natarajan Shankar2 1

2

1

fortiss GmbH, Guerickestr. 25, 80805 M¨unchen, Germany {cheng,ruess}@fortiss.org

SRI International, 333 Ravenswood, Menlo Park, CA 94025, United States [email protected]

Overview

Efficient and scalable verification of nonlinear real arithmetic constraints is essential in many automated verification and synthesis tasks for hybrid systems, control algorithms, digital signal processors, and mixed analog/ digital circuits. Despite substantial advances in verification technology, complexity issues with classical decision procedures for non-linear real arithmetic are still a major obstacle for formal verification of real-world applications. Recently, Mu˜noz and Narkawicz [3] proposed a procedure for deciding the validity of quantifierfree nonlinear real arithmetic based on the well-known transformation to Bernstein polynomials. Their Kodiak system outperforms tools based on cylindric algebraic decomposition, including QEPCAD [1] and REDLOG [2], in many cases. Starting from the algorithms by Mu˜noz and Narkawicz [3] we extended the approach and implemented the (little) verification engine JBernstein, which checks the validity of finite conjunctions of nonlinear arithmetic constraints of the form ∀x0 ∈ [l0 , u0 ], . . . , xmV∈ [lm , um ] : n ( j=1 Pj (x0 , . . . , xm ) ≺j cj ) → Q(x0 , . . . , xm ) ≺ d Pj and Q are real polynomials over the variables x0 through xm , each xi is interpreted over closed intervals [li , ui ] with real-valued lower and upper bounds, cj and d are real-valued constants, and the symbols ≺j and ≺ are arithmetic inequalities in {>, ≥, <, ≤}. These constraints support assume-guarantee style of reasoning about open systems, with Pi the assumptions on the environment and Q the corresponding guarantee of the system under consideration. The Java implementation JBernstein includes a number of algorithmic optimizations as described in Section 2, for example, for avoiding unnecessary case splits. In particular, JBernstein uses double-precision floating-point arithmetic of Java in a sound way via constraint strengthening (Section 3). The resulting runtimes of JBernstein are compared in Figure 1 with those reported by Mu˜noz and Narkawicz [3] for their Kodiak implementation, QEPCAD, and REDLOG. This comparison uses the PVS test suite as compiled by Mu˜noz and Narkawicz [3]. The experimental evaluation of these optimizations indicates that JBernstein is usually an order of magnitude faster than earlier results by Mu˜noz and Narkawicz or our un-optimized implementation, and it is two orders of magnitude faster than QEPCAD or REDLOG. For the experimental evaluation we use a similar hardware setting (Intel Core Duo 2.4 Ghz, MacOS, 8 GB RAM) as reported in [3]. In contrast to the other verification systems JBernstein is able to solve all benchmark problems, and it solves them one or two orders of magnitudes faster than the other systems. These initial experimental results are indeed promising, but, clearly, an extended and improved set of benchmarks is needed to obtain an indication about the asymptotic behavior of these solvers.

2

Chih-Hong Cheng, Harald Ruess and Natarajan Shankar

JBernstein is implemented in Java without further dependencies. It is freely available under the LGPL version 3 license at http://sourceforge.net/projects/jbernstein/

Problem

JBernstein Kodiak [3] REDLOGrlqe [3] REDLOGrlcad [3] QEPCAD [3] Metit [3]

Schwefel (∀) Schwefel (∃) Reaction Diffusion (∀) Reaction Diffusion (∃) Caprasse (∀) Caprasse (∃) Lotka-Volterra (∀) Lotka-Volterra (∃) Butcher (∀) Butcher (∃) Magnetism (∀) Magnetism (∃) Heart Dipole (∀) Heart Dipole (∃)

593§ 572‡ 7 3 23 8 5 4 19 65 48 50 271 231

§ ‡

940 280 < 10 < 10 290 310 100 < 10 200 200 73540 320 7360 3700

490 138900 340 340 1750 15060 360 350 420 360 670 420 > 300000 > 300000

> 300000 > 300000 370 350 > 300000 > 300000 450 400 > 300000 > 300000 360 360 > 300000 > 300000

840 910 10 10 6540 6540 10 10 (abort) (abort) 180 350 > 300000 > 300000

110 (n/a) 90 (n/a) 160 (n/a) 100 (n/a) (abort) (n/a) 540 (n/a) > 300000 (n/a)

For this example, we only ask > −0.000001, as precision can not be maintained for the original property. For this example, we only ask > 0.000001.

Table 1. Performance of JBernstein and other tools. Results from other tools are taken from [3]. The measured unit for execution time is in milliseconds.

2

Algorithmic Optimizations

We outline the algorithm of Mu˜noz and Narkawicz in [3] and the optimizations thereof in JBernstein using simplified constraints of the form ∀x0 ∈ [l0 , u0 ], . . . , xm ∈ [lm , um ] : φ(x0 , . . . , xm ) > c . The Bernstein approach first performs a range-preserving transformation from [a, b] to [0, 1] to obtain a constraint of the form ∀y0 ∈ [0, 1], . . . , ym ∈ [0, 1] : φ0 (y0 , . . . ym ) > c , such that xi = li + yi (ui − li ) for all i ∈ {1, . . . , m} and φ(x1 , . . . , xm ) = φ0 (y1 , . . . , ym ). The polynomial φ0 is then translated from polynomial basis into Bernstein basis. For example, φ0 (y) = 4y 2 − 4y + 1 has polynomial basis {y 2 , y, 1}. It is rewritten as 1(20 )(1 − y)2 -1(21 )y(1 − y) + 1(22 )(y)2 , with {(2k )y k (1 − y)2−k |k = 0, 1, 2} being the Bernstein basis. If all coefficients of the polynomial represented in Bernstein basis are greater than c, then the original constraint holds (true). Otherwise one checks if there exists a coefficient of an endindex (Bernstein basis vector where every term is either 0 or of highest degree) that is smaller or equal to c. For example, for ∀y ∈ [0, 1] : φ0 (y) > −3, as all Bernstein coefficients 1, −1, 1 are greater than −3, and the property holds. If the property is ∀y ∈ [0, 1] : φ0 (y) > 2, as for the first and the last coefficient (in {(2n )y k (1 − y)2−n |n = 0, 1, 2}, n = 0, 2 are endindices), we have 1 6> 2, and the property fails to hold at y = 0 and y = 1 (false). Lazy refinement. The checking process returns unknown if neither of these two conditions holds, i.e., there exists some non-endindex coefficients less or equal to c. In these cases, the algorithm does range splitting on some chosen variable from [0, 1] to [0, 21 ] and [ 12 , 1], it generates Bernstein

JBernstein: Checking Polynomial Constraints Variable index 0 1 2 Array index

0 0

1 1

0

2 1

0

3

offset[0] = 4 × (1 + 1) = 8

1

offset[1] = 1 × (3 + 1) = 4

0 1 2 3 0 1 2 3 0 1 2 3 0 1 2 3 0 1 2 3 0 1 2 3 offset[2] = 1 4 20 20 = 4 + (2 − 0)offset[0] Fig. 1. Index shifting by table look-up.

polynomials for each subspace, and it checks if the property does indeed hold for both polynomials. The Bernstein polynomial checker refines a subspace only when a decisive proof has not been found for values within this subspace. Inputs of the form ρ > d → ϕ > e are rewritten in JBernstein to its disjunctive form ρ ≤ d ∨ ϕ > e. The checker returns true if one of them is valid, and it returns false if there exists an endindex of both polynomials whose coefficients violate the constraint separately. Recursive variable exploration. The refinement process involves recursive calls of variable selection and domain partitioning to generate new Bernstein polynomials. For efficiency it is important to avoid variable selections which lead to the generating of an excessive number of polynomials for which the checker returns unknown. As the algorithm is recursive, a na¨ıve deterministic implementation can, given a fixed recursion depth, try to refine the first variable continuously until the end of the recursion depth. It turns out that this strategy often leads to a generation of numerous useless Bernstein polynomials without providing definite proofs. The reason is explained using the following extreme case. Continuously refine only one variable, say y0 , then eventually the generated Bernstein polynomial will have y0 with its domain converged from [0, 1] to one single point α ∈ [0, 1]. This means that in the original polynomial, the refinement tries to analyze the polynomial by setting x0 to l0 + α(u0 − l0 ), whose proving process is, intuitively, nearly as hard as working on the original polynomial. Given a certain budget for maximum allowed recursion, first trying deep recursion over one particular variable is therefore, unlikely to succeed in most cases. Therefore, the implementation of JBernstein uses a round-robin selection strategy on the variables during recursion. Avoiding superfluous computations. The third aspect concerning efficiency focuses on reducing the unit cost for generating every new Bernstein polynomial and checking the property. Generating a new Bernstein polynomial (domain partition) entails generating all of its coefficients, which themselves are derived from the coefficients of the original Bernstein polynomial (without domain partition). The new coefficients bL k for the left partitioned polynomial (for the right polynomial, the formula is similar), where yj is chosen for partition, is achieved using the following formula [3].   kj X 1 kj L bk = bk with [j:=r] (1) 2kj r r=0 In Eq. 1, for an m-variate polynomial, k is a vector of m tuples of positive integers or 0, where each term is smaller than the cardinality (i.e., every k is the unique signature of each Bernstein basis vector). kj is the j-th value of k, and bk with [j:=r] is the coefficient of the original Bernstein polynomial, where ”k with [j := r]” is an m-tuple that is equal to k on every index, except in index j where it has value r. The following optimizations are used by JBernstein: 1 2kj

from summation.
2. Replace the computation of krj by table look-up.

1. Factor out the constant

4

Chih-Hong Cheng, Harald Ruess and Natarajan Shankar

3. Compute bk with [j:=r] in an optimized way as follows: • Statically store k linearly in an array (prior to the range transformation). Figure 1 indicates how vectors are arranged. E.g., when k = (0, 1, 0), it is located in index 4. Store each bk linearly in an array following the above order, for every Bernstein polynomial. • When the solver iterates the array to create bL k , the index of k is known. • To find bk with [j:=r] , it amounts to finding the index of ”k with [j := r]”. Due to our formulation, it can now be translated to an index offset problem: given the index of k, what is the offset when replacing the j-th index with value r? The offset is (r − kj )offset[j], where offset[j] = (Degj+1 + 1) × . . . × (Degm−1 + 1), where Degs is the highest degree that appears for variable vs in the polynomial. The array offset can also be computed statically prior to the range-preserving transformation. Figure 1 illustrates the computation of the index of of (2, 1, 0) from k = (0, 1, 0). • Therefore, for bk with [j:=r] the solver uses three table look-ups (for offset[j], kj and the final value), one substraction, multiplication and addition. 4. Integrate the coefficient generation and the checking process. In each Bernstein polynomial, create an internal Boolean variable field isUnknown that is initially set to false. During the construction, check if a particular coefficient satisfies the property only when it belongs to the endindex or isUnknown equals false. Once if the polynomial is diagnosed as unknown, then there is no need to check for coefficients from non-endindices. Deciding whether a certain index is an endindex is also done statically once and is replaced by table look-up in later computations.

3

Sound Usage of double

We now justify the use of double (double-precision 64-bit IEEE 754 floating point) for data representation. In JBernstein, potential errors due to imprecision of double are handled by the following methodology (for the ease of explanation, we again set the property to be ∀x0 ∈ [l0 , u0 ], . . . , xm ∈ [lm , um ] : φ(x0 , . . . , xm ) > c): • Select a positive error-estimate  such as 10−5 . • To return true, in property checking all coefficients shall be greater than c + . • To return false, in property checking the solver needs to find a coefficient bk from an endindex such that bk ≤ c − . • If neither of the above two cases holds, the solver either proceeds with domain refinement (when recursion is still allowed) or returns unknown. The correctness relies on a crucial requirement that the accumulated error for each computed coefficient should never exceed . Instead of keeping track of the error during the computation, we apply static analysis on the algorithm to generate a safe error-estimate that holds for each computed coefficient, based on the polynomial constraint itself and the number of maximum refinement attempts. Due to space limits, we only review the key feature in refining a subspace to generate new Bernstein polynomials (a full text concerning the sound usage of double can be found in the appendix). Pkj 1 kj
For each refinement, recall in Eq. 1 where we have bL r=0 2kj r bk with [j:=r] . E.g., k =
1 kj 1 4 6 4 1 when kj = 4, 2kj r equals 16 , 16 , 16 , 16 , 16 (the sum of these values is 1) for r = 0, 1, 2, 3, 4. This means that each new coefficient is a weighted sum from old coefficients. If each original coefficient bk with [j:=r] has an error-estimate bounded by κ, in an idealized computation, the generated bL k also has an error-estimate bounded by κ. This gives an intuition that the growth of error should be very slow (nearly linear to the number of refinement attempts) within the refine-

JBernstein: Checking Polynomial Constraints

5

Fig. 2. Screenshot for running JBernstein.

ment process. The behavior of linear growth is the key ingredient that makes our methodology applicable3 , which is one of the nice properties for Bernstein polynomials.

4

Examples

Figure 2 demonstrates the usage of the graphical interface of JBernstein for the simple constraint ∀x0 ∈ [3, 5], x1 ∈ [2, 7] : 4x20 + 5x21 + 3x0 x1 + 2x0 + 5 > 85.

(2)

The statement VAR 2 in Figure 2 indicates that the constraint has two variables of name x0 and x1 , and statements such as "BOUND x0 [3,5]" are used to specify that x0 is interpreted over the closed interval [3, 5]. Finally, a statement "COEF (2,0) 4" specifies that the monomial x20 (i.e., x20 x01 ) has coefficient 4. The value ”Recursive expansion depth” specifies the maximum number box-refinement attempts for each variable to prove/disprove the constraint. Hitting the ”Analyze” button starts the checker, which provides the counter-example x0 = 3.0, x1 = 2.0 within 15 milliseconds. Options such as error-estimation (none, auto, user-selected) or printing precise derivations for counter-examples are also part of the user-interface. Assume-guarantee-style constraints such as ∀x0 ∈ [0, 1] : (6x0 − 1 > 0 ∧ 3x0 − 1 < 0) → 125x30 − 175x20 + 70x0 − 8 > 0 3

If every refinement computation brings κ twice as large, in static analysis, applying recursive expansion for small steps like 100 will make an initially small error-estimate prohibitively huge in the lastly generated Bernstein polynomial.

6

Chih-Hong Cheng, Harald Ruess and Natarajan Shankar

## FORALL x\in[0,1]: (6x-1>0 && 3x-1<0) -> (5x-1)(5x-2)(5x-4)>0 ## Specify to use one variables x0 VAR 1 ## Specify the NUMBER of conjunctions ## Usage: "CONJUNCTION NUMBER" CONJUNCTION 1 ## Specify the NUMBER of assumptions in the INDEX-th conjunction element ## Usage: "ASSUMP INDEX NUMBER" ASSUMP 0 2 ## Specify the coefficient of the polynomial (assumption) 5x+0>0 as ## E.g., A1_0 means the 1st assumption (indexed 0) in the 2nd conjunction (indexed 1) ## First assumption COEF A0_0 (1) 6 COEF A0_0 (0) -1 SIGN A0_0 GT VALUE A0_0 0 ## Second assumption COEF A0_1 (1) 3 COEF A0_1 (0) -1 SIGN A0_1 LT VALUE A0_1 0 ## Specify the coefficient of the polynomial (guarantee) ## (5x-1)(5x-2)(5x-4) = 125xˆ3-175xˆ2+70x-8 COEF G0 (3) 125 COEF G0 (2) -175 COEF G0 (1) 70 COEF G0 (0) -8 SIGN G0 GT VALUE G0 0 ## Specify the bound for each variable BOUND x0 [0, 1] ## Result: FALSE

Fig. 3. Assume-guarantee-style constraints in JBernstein.

are handled by the textual interface of JBernstein as described in Figure 4: • CONJUNCTION 1 specifies the use of one assume-guarantee rule. • ASSUMP 0 2 indicates that the first assume-guarantee rule (indexed 0) has two assumptions 6x0 − 1 > 0 and 3x0 − 1 < 0. • For the first assumption (indexed 0) of the first assume-guarantee rule (6x0 − 1 > 0), use ”COEF A0_0 (1) 6”, ”COEF A0_0 (0) -1”, ”SIGN A0_0 GT”, and ”VALUE A0_0 0” to specify the polynomial. Acknowledgements. We thank Dr. C´esar Mu˜noz (NASA Langley) for his support and helpful suggestions.

References 1. C. W. Brown. QEPCAD-B: a program for computing with semi-algebraic sets using CADs. SIGSAM Bull., 37(4):97–108, Dec. 2003. 2. A. Dolzmann and T. Sturm. REDLOG: computer algebra meets computer logic. SIGSAM Bull., 31(2):2– 9, June 1997. 3. C. Mu˜noz and A. Narkawicz. Formalization of a representation of Bernstein polynomials and applications to global optimization. Journal of Automated Reasoning, 2012. Accepted for publication.

JBernstein: Checking Polynomial Constraints

7

Appendix A

Further Remarks on the Use of double

In this section, we detail the discussion for error-estimation in the recursion process. Then, we outline how JBernstein handles error-estimates in (1) range-preserving transformation and (2) the conversion from polynomial basis to Bernstein basis. Notice that these methods may not be optimized, and other methodologies (such as precise program analysis) can derive more precise estimates. A.1

General methodology

In JBernstein, potential errors due to imprecision of double are handled by the following methodology (for the ease of explanation, we again set the property to be ∀x0 ∈ [l0 , u0 ], . . . , xm ∈ [lm , um ] : φ(x0 , . . . , xm ) > c): • Select a positive error-estimate  such as 10−5 . • To return true, in property checking all coefficients shall be greater than c + . • To return false, in property checking the solver needs to find a coefficient bk from an endindex such that bk ≤ c − . • If neither of the above two cases holds, the solver either proceeds with domain refinement (when recursion is still allowed) or returns unknown. A.2

Error due to generating new Bernstein polynomials

The correctness relies on a crucial requirement that the accumulated error for every computed coefficient should never exceed . Instead of keeping track of the error during the computation, we apply static analysis on the algorithm to generate a safe error-estimate that holds for every computed coefficient, based on the polynomial constraint itself and the number of maximum refinement attempts. We now provide some important observations and implementation strategies that help to establish the correctness of using double. When generating a new Bernstein polynomial from an old one via subspace partitioning, as can be seen from Eq. 1, the new coefficient is a sum of terms, where each term involves numerical computation from old coefficients. If all coefficients of the original polynomial have an errorestimate that is bounded by κ, we need to check the effect of κ (due to arithmetic computation) that is brought into the newly computed coefficient. E.g., if the computation brings κ twice as large, then applying recursive expansion for small steps like 100 will make an initially small error-estimate prohibitively huge in the lastly generated Bernstein polynomial. Fortunately, the overall computation in Bernstein polynomial only generates errors that are a weighted sum of errors from the parent polynomial. Recall in Eq. 1 where we have bL k =
Pkj 1 kj
1 kj 1 4 6 4 1 b . E.g., when k = 4, equals , , , , (the sum of j k k k with [j:=r] r=0 2 j r 16 16 16 16 16 2 j r these values is 1) for r = 0, 1, 2, 3, 4. If every bk with [j:=r] has an error-estimate bounded by κ, so does the generated bL k . This gives an intuition that the growth of error should be very slow (nearly linear to the number of refinement attempts) within refinement process. We illustrate how correctness is ensured in JBernstein. • Consider if in the first Bernstein polynomial (i.e., the one that is derived from converting polynomial basis into Bernstein basis), the largest coefficient (when estimating errors in range-preserving and basis transformation) is below 950 and the smallest is above −875.

8

Chih-Hong Cheng, Harald Ruess and Natarajan Shankar

One can fetch κ[−1024,1024] , an overapproximation on the largest possible truncation error when trying to generate a new Bernstein coefficient from original coefficients within [−1024, 1024]. As [−875, 950] ⊆ [−1024, 1024], in the first subspace partitioning, all arithmetic errors are bounded by κ[−1024,1024] . On polynomials with highest variable degree such as 4, κ[−1024,1024] is very small4 . In JBernstein, κ[−1024,1024] is computed statically for different degrees, together with κ[−128,128] , κ[−8192,8192] , etc. • In the second refinement, coefficients of the new Bernstein polynomial are weighted sum of original coefficients whose values are now bounded by −875 − κ[−1024,1024] and 950 + κ[−1024,1024] . If −875−κ[−1024,1024] > −1024 and 950+κ[−1024,1024] < 1024, κ[−1024,1024] can again act as a safe error-estimate for arithmetic truncation, when considering the second refinement. • In general, if the algorithm detects that the first Bernstein polynomial has coefficients whose range is [a, b] ⊆ [−1024, 1024], and if the solver only allows the refinement to proceed α times, the error that arises in the recursion will always be bounded by ακ[−1024,1024] . Let δ be the error-estimate on (1) range-preserving transformation and (2) conversion from polynomial basis to Bernstein basis (see later sections). Then implementing a run-time check on (a−ακ[−1024,1024] > −1024) ∧ (b+ακ[−1024,1024] < 1024) ∧ (δ +ακ[−1024,1024] < ) can guarantee correctness: If the check fails, the engine immediately returns unknown5 . The first two conjunctions ensure the validity of κ[−1024,1024] as an error-estimate for each refinement attempt, while the last conjunction is the guarantee stating that the accumulated error does not influence the correctness of property checking. The quantity of the errorestimate grows linearly with the number of refinement steps, not exponentially. The behavior of linear growth is the key ingredient that makes our methodology applicable (by returning true in most cases), which is one of the nice properties for Bernstein polynomials. For problems within our benchmark, coefficients of the first Bernstein polynomial are mainly within [−128, 128], and JBernstein only uses a small number (such as 100) for refinement attempts. The error accumulated to due arithmetic truncation in the refinement is very small, compared to the user-specified  which is in general between 10−5 to 10−8 . Lastly, to report a counter-example, as the refinement always splits intervals into two equally sized ones, a counterexample can be restored perfectly by recording the branching sequence for each variable. A.3

Range-preserving transformation

In range-preserving transformation, a new coefficient rk is generated via P the following Qm−1 rule
(which is based on linear transformation; details can be found in [3]): rk = k≤i≤n ci j=0 kijj (uj − i −k

lj )kj ljj j , where i and k are basis vectors, n is the largest basis vector in the polynomial (nj is the largest degree for variable indexed j), uj and lj are upperbound and lowerbound for variable xj , and ci is the coefficient from the original polynomial. Our strategy for error estimation tries 4

5

If the highest degree equals to 4, a weighted sum that generates a new coefficient can only contain at 4 6 4 1 1 , 16 , 16 , 16 , 16 . For each term, we most 5 terms from old coefficients. Each coefficient is paired with 16 only need to consider the error due to multiplication of 1, 4, or 6, as 16 is divided at the end due to our optimization, and in double it only decreases the exponent by 4. Conservatively, the largest error that can happen in addition-arithmetic appears when the first value is 6 × 1024 and the second is 2−52 times (double uses 52 bits for fraction) smaller than 1024. When this happens, the addition returns the first value, bring an error < 2−39 . With 4 additions, we derive 2−37 . Every element is created from coefficients that is multiplied by 1, 4 = 22 , or 6 = 22 + 2. These multiplication only cause at most 3 bit loss (1, 4 does not cause loss of precision). Therefore, we use 2−34 to be a safe error-estimate. This is conservative and can be further improved by checking the error for every generated Bernstein polynomial that returns true or false.

JBernstein: Checking Polynomial Constraints

9


to find the maximum possible value for each ci , kijj , (uj − lj ), kj , lj , and ij − kj to make the result of the multiplication largest. Then the error is conservatively estimated by considering the loss due to multiplication and addition. Notice that if for basis vector i where ci = 0, we directly omit addition (as thePvalue equals 0) and it is impossible to introduce error. Therefore, when considering the sum k≤i≤n , we only consider non-zero terms. Here we use examples to illustrate the method. Consider 4x20 − 5x21 + 3x0 x1 − 2x0 + 5, where x0 ∈ [−7, 2], x1 ∈ [3, 5]. The largest possible term that
can appear when constructing
a new coefficient, following the formulation in [3], is | − 5|( 21 (2 − (−7))2 (| − 7|)2 )( 21 (5 −
3)2 (|3|)2 ) = 2857680. The first term | − 5| is the largest quantity of a coefficient. 21 is the binomial coefficient that chooses half of the from 2 (which is the highest degree). (2 − (−7)) is the interval for the domain of x0 . The exponent 2 is the highest degree; for base that is smaller than 1, use 0 to ensure generating largest possible value. In (| − 7|)2 , −7 is the lowerbound of x0 . The relative error (the quantity of the error divides the actual value) for multiplication is bounded by 26(2−52 ), where 26 is the total number of elements that participate in the multiplication plus 1 (i.e., 12+1) times a factor of two6 . Representing 2857680 requires 22 bits. To create a new coefficient, one takes a sum which uses 5 non-zero coefficients (for basis vector whose coefficient is 0, 0 is added without posing any error). So the sum can be represented in 26 bits. Thus error can at most be 5(2−52+26 )dlog2 (26)e, where 2−52+26 is the single highest error that can be achieved in one addition operation, and it appears at most 5 times. dlog2 (26)e is the precision loss in multiplication. The above checking only involves properties of the polynomial such as the largest degree, the size of the domain for every variable, and the lowerbound. The error-estimate is implemented efficiently without decreasing the performance.

A.4

Transformation from polynomial basis into Bernstein basis

In basis transformation, the Bernstein coefficient indexed k is generated from the formula below: k P Qm−1 ( ijj ) p n i≤k i j=0 ( j ) , where pi is the coefficient for polynomial basis indexed i, i and k are basis ij vectors, and nj is the largest degree for variable indexed j. We consider the scenario where every pi is replaced by a common error-estimate κ. k k P P Qm−1 ( i j ) Qm−1 Qm−1 ( i j ) In i≤k κ j=0 njj = κ i≤k j=0 njj , the number of additions is j=0 (kj + 1). In ( ij ) ( ij ) k k Qm−1 ( i j ) Qm−1 k ! Qm−1 P Qm−1 ( ijj ) 1 each multiplication, j=0 njj = j=0 njj ! = j=0 (nj −k . Therefore, = nj i≤k j=0 j )! ( ij ) ( ij ) Qm−1 1 j=0 (kj +1) (nj −kj )! . For every kj 6= nj , each term is either smaller or equal to 1, so the multiplication is always less Q or equal to 1. Therefore, given a common error-estimate κ on coefficients m−1 of polynomial basis, κ j=0 (nj + 1) is an estimate on the resulting error due to the change of basis. 6

Consider two values a, b that perform multiplication. a, b have relative error  = 2−52 (double uses 52 bits for its fraction). Then the multiplication is a(1 + )b(1 + ) = ab(1 + 2 + 2 ). The value will also have a relative imprecision (1 + ) due to truncation. We have ab(1 + 2 + 2 )(1 + ) ≈ ab(1 + 3). The next multiplication with value c(1 + ) brings relative error ≈ abc(1 + 5). Therefore, the number of relative error after 12 multiplications is 25. We use 25+1 = 26, i.e., by introducing additional term rather than 25, for sound overestimation, i.e., the sum of all higher terms such as 2 , 3 , . . . are no larger than .

10

Chih-Hong Cheng, Harald Ruess and Natarajan Shankar

k Qm−1 ( i j ) Notice that κ is also a safe estimate for the lose of addition, as every term pi j=0 njj is at (i ) Qm−1 Qm−1 j most equal to pi . For at most j=0 (kj + 1) additions, the error is bounded by κ j=0 (nj + 1). Qm−1 Therefore, we use 2κ j=0 (nj + 1) as an error-estimate.

A.5

Error-estimate for the Heart Dipole problem

(Range-transformation) For the heart-dipole problem [3], the polynomial constraint is −x1 x36 + 3x1 x6 x27 − x3 x37 + 3x3 x7 x6 − x2 x35 + 3x2 x5 x28 − x4 x38 + 3x4 x8 x25 , There are 8 non-zero coefficients. Following the transformation the largest term is of the following.

|3|(

! 1 (0.4 − (−0.1))0 (| − 0.1|)0 1 ! 1 (1 − (−0.4))0 (|0.4|)0 1 !

1 (−0.4 − (−0.7))0 (| − 0.7|)0 1 ! 1 (0.4 − (−0.7))1 (| − 0.7|)0 1 ! 3 (0.2 − (0.1))0 (|0.1|)0 1 ! 3 (0.2 − (−0.1))0 (| − 0.1|)0 1 ! 3 (1.1 − (0.3))3 (|0.3|)0 1 ! 3 (−0.3 − (−1.1))0 (| − 1.1|)3 ) 1

(3)

This value is bounded by 27 , and within the generation, is summed with non-zero term for at most 8 times. Therefore, the coefficient can at most be 210 . The error thus can at most be 8(2−52+10 ) = 2−39 , where 2−52+10 is the single highest error that can be introduced by one addition operation, and it appears 8 times (recall that zero coefficients do not introduce any error). The number of multiplications constituting each element is 25. Therefore, the bit imprecision due to relative error is dlog2 ((25 + 1)2)e = 6, and the error can be bounded by 2−39+6 . (Basis-transformation) When considering the error due to transformation from polynomial baQm−1 sis into Bernstein basis, we have an error-estimate 2κ j=0 (nj + 1) = 2(2−39+6 )(1 + 1)4 (3 + 1)4 = 2−20 , where (1 + 1)4 are for x1 , x2 , x3 , x4 (they have highest degree 1), and (3 + 1)4 are for x5 , x6 , x7 , x8 (they have highest degree 3). This value (2−20 ) is smaller than 10−5 , which partially explains why the solver can give correct solution on the Heart-Dipole problem (whose property contains 4 digits below decimal) by setting  to values such as 10−5 .

JBernstein: Checking Polynomial Constraints

B

11

Demonstration Plans

As the tool requires nearly no dependencies except Java, in the demonstration, we will invite the audience to download the tool directly and play around with the tool themselves. When time permitted, we will also outline the underlying methodology with one or two slides, together with verification applications that uses JBernstein in our institute.

C

Images of JBernstein running PVS testsuites

We captured the screenshot (containing the analysis time) of JBernstein running on every example. These images can be downloaded from http://www6.in.tum.de/˜chengch/ efsmt/jbernstein.html.