OOn Sound Compilation of Reals
Eva Darulova
EPFL eva.darulova@epfl.ch
Viktor Kuncak
EPFL viktor.kuncak@epfl.ch
Abstract
Writing accurate numerical software is hard because of manysources of unavoidable uncertainties, including finite numericalprecision of implementations. We present a programming modelwhere the user writes a program in a real-valued implementationand specification language that explicitly includes different typesof uncertainties. We then present a compilation algorithm that gen-erates a conventional implementation that is guaranteed to meetthe desired precision with respect to real numbers. Our verificationstep generates verification conditions that treat different uncertain-ties in a unified way and encode reasoning about floating-pointroundoff errors into reasoning about real numbers. Such verifica-tion conditions can be used as a standardized format for verifyingthe precision and the correctness of numerical programs. Due totheir often non-linear nature, precise reasoning about such verifica-tion conditions remains difficult. We show that current state-of-theart SMT solvers do not scale well to solving such verification con-ditions. We propose a new procedure that combines exact SMTsolving over reals with approximate and sound affine and intervalarithmetic. We show that this approach overcomes scalability lim-itations of SMT solvers while providing improved precision overaffine and interval arithmetic. Using our initial implementation weshow the usefullness and effectiveness of our approach on severalexamples, including those containing non-linear computation.
1. Introduction
Writing numerical programs is difficult, in part because the pro-grammer needs to deal not only with the correctness of the al-gorithm but also with unavoidable uncertainties. Program inputsmay not be exactly known because they come from physical exper-iments or were measured by an embedded sensor. The computa-tion itself suffers from roundoff errors at each step, because of theuse of finite-precision arithmetic. In addition, resources like energymay be scarce so that only a certain number of bits are availablefor the numerical datatype. At the same time, the computed resultscan have far-reaching consequences if used to control, for example,a vehicle or a nuclear power plant. It is therefore of great impor-tance to improve confidence in numerical code [34]. One of the firstchallenges in doing this is that most of our automated reasoningtools work with real arithmetic, whereas the code is implemented infinite-precision (typically, floating point) arithmetic. Many currentapproaches to verifying numerical programs start with the floating-point implementation and then try to verify the absense of (runtime)errors. However, the absence of errors by itself does not guaranteethat program behavior matches the desired specification expressedusing real numbers. Fundamentally, the source code semantics isexpressed in terms of data types such as floating points. This isfurther problematic for compiler optimizations, because, e.g., theassociativity law is unsound with respect to such source code se-mantics. In this paper we advocate a natural but ambitious alternative:source code programs should be expressed in terms of ideal, math-ematical real numbers. In our system, the programmer writes a pro-gram using a
Real data type and states the desired postconditions,then specifies uncertainties explicitly in pre- and postconditions aswell as the desired target precision on the return values. It is thenup to our trustworthy compiler to check, taking into account all un-certainties and their propagation, that the desired precision can besoundly realized in a finite precision implemenation. If so, the com-piler chooses and emits one such finite-precision implementation.Following this model, we derive verification conditions that en-code reasoning about floating-point roundoff errors into reason-ing about real numbers. Our verification conditions explicitly sep-arate the ideal program without external uncertainties and round-offs from the actual program, which is actually executed in finiteprecision with possibly noisy inputs. Additionally, our constraintsare fully parametric in the floating-point precision and can thus beused with different floating-point hardware configurations, as wellas used for determining the minimum precision needed.To summarize, the view of source code as functions over realnumbers has several advantages: • Programmers can, for the most part, reason in real arithmeticand not floating-point arithmetic. We thus achieve separation ofthe design of algorithms (which may still be approximate) fromtheir realization using finite precision computations. • We can verify the ideal meaning of programs using techniquesdeveloped to reason over real numbers, which are more scalablethan techniques that directly deal with floating point arithmetic. • The approach allows us to quantify the deviation of implemen-tation outputs from ideal ones, instead of merely proving e.g.range bounds of floating-point variables which is used in sim-pler static analyses. • The compiler for reals is free to do optimizations as long as theypreserve the precision requirements. This allows the compilerto apply, for example, associativity of arithmetic, or even selectdifferent approximation schemes for transcendental functions. • In addition to roundoff errors, the approach also allows thedeveloper to quantify program behavior in the face of externaluncertainties such as input measurement errors.Using our verification condition generation approach, the cor-rectness and the precision of compilation for small programs canbe directly verified using an SMT solver such as Z3 [16] (see Sec-tion 5.3); this capability will likely continue to improve as thesolvers advance. However, the complexity of the generated veri-fication conditions for larger programs is still out of reach of suchsolvers and we believe that specialized techniques are and will con-tinue to be necessary for this task. This paper presents two special-ized techniques that improve the feasibility of the verification task.The first technique performs local approximation and is effective a r X i v : . [ c s . P L ] S e p ven in benchmarks containing nonlinear arithmetic; the secondtechnique specifically handles conditional expressions. Nonlinear arithmetic poses a significant challenge for verificationbecause it cannot directly be handled using Simplex-like algorithmsembedded inside SMT solvers. Although interesting relevant frag-ments are decidable and are supported by solvers, their complex-ity is much higher. Unfortunately, non-linear arithmetic is ubiq-uitous in numerical software. Furthermore, our verification condi-tions add roundoff errors to arithmetic expressions, so the resultingconstraints grow further in complexity, often becoming out of reachof solvers. An alternative to encoding into SMT solver input is touse a sound and overapproximating arithmetic model such as in-terval or affine arithmetic [15]. However, when used by itself onnon-linear code, these approaches yield too pessimistic results tobe useful.We show that we can combine range arithmetic computationwith SMT solving to overcome the limitations of each of the indi-vidual techniques when applied to non-linear arithmetic. We obtaina sound, precise, and somewhat more scalable procedure. Duringrange computation, our technique also checks for common prob-lems such as overflow, division by zero or square root of a negativenumber, emitting the corresponding warnings. Additionally, be-cause the procedure is a forward computation, it is suitable for auto-matically generating function summaries containing output rangesand errors of a function. From the point of view of the logicalencoding of the problem, range arithmetic becomes a specializedmethod to perform approximate quantifier elimination of boundedvariables that describe the uncertainty.
In the presence of uncertainties, conditional braches become an-other verification challenge. The challenge is that the ideal execu-tion may follow one branch, but, because of input or roundoff er-rors, the actual execution follows another. This behaviour may beacceptable, however, if we can show that the error on the outputremains within required bounds. Therefore, our approach woulddirectly benefit from automated analysis of continuity, which wasadvocated previously [48]. For such continuous functions, our anal-ysis can be done separately for each program path. In the absenceof knowledge of continuity, we present a new method to check thatdifferent paths taken by real-valued and floating point versions ofthe program still preserves the desired precision specification. Ourcheck does not require continuity (which becomes difficult to provefor non-linear code). Instead, it directly checks that the differencebetween the two values on different branches meets the requiredprecision. This technique extends our method for handling non-linear arithmetic, so it benefits from the combination of range arith-metic and SMT solving.
We have implemented our compilation and verification procedure,including the verification condition generation, analysis of possi-bly non-linear expressions, and the handling of conditionals. Oursystem is implemented as an extension of a verifier for functionalScala programs. The implementation relies on a range arithmeticimplementation for Scala as well as on the Z3 SMT solver.We have evaluated the system on a number of diverse bench-marks, obtaining promising results. We are releasing our bench-marks (and making them available as supplementary material forreviewers). We hope that the benchmarks can be used to com-pare future tools for error quantification, help the development ofnonlinear solvers, and also present challenges for more aggressivecompilation schemes, with different number representations and different selection of numerical algorithms. To support program-ming of larger code fragments our system also supports a modu-lar verification technique, which handles functions through inliningfunction bodies or using their postconditions. We thus expect thatour technique is applicable to larger code bases as well, possiblythrough refactoring code into multiple smaller and annotated func-tions. Even on the benchmarks that we release, we are aware of noother available system that would provide the same guarantees withour level of automation.
Our overall contribution is an approach for sound compilation ofreal numbers. Specifically: • We present a real-valued implementation and specification lan-guage for numerical programs with uncertainties; we define itssemantics in terms of verification constraints that they induce.We believe that such verification conditions can be used as astandardized format for verifying the precision and the correct-ness of numerical programs. • We develop an approximation procedure for computing preciserange and error bounds for nonlinear expressions which com-bines SMT solving with interval arithmetic. We show that suchan approach significantly improves computed range and errorbounds compared to standard interval arithmetic, and scalesbetter than SMT solving alone. Our procedure can also be usedindependently as a more precise alternative to interval arith-metic, and thus can perform forward computation without hav-ing the desired postconditions. • We describe an approach for soundly computing error boundsin the presence of branches and uncertainties, which ensuressoundness of compilation in case the function defined by aprogram is not known to be continuous. • We have implemented our framework and report our experienceon a set of diverse benchmarks, including benchmarks fromphysics, biology, chemistry, and control systems. The resultsshow that our technique is effective and that it achieves a syn-ergy of the techniques on which it relies.
2. Example
We demonstrate some aspects of our system on the example in Fig-ure 1. The methods triangle and triangleSorted compute the area ofa triangle with side lengths a, b and c . The notation a.in(1.0, 9.0) is short for < a && a < . We consider a particular applica-tion where the user may have two side lengths given, and mayvary the third. She has two functions available to do the compu-tation and wants to determine whether either or both satisfy theprecision requirement of 1e-11 on line 9. Our tool determines thatsuch requirement needs at least double floating point precision; thechallenge then is to establish that this precision is sufficient to en-sure these bounds, given that errors in floating code accumulate andgrow without an a priori bound.Our tool verifies fully automatically that the method triangleSorted indeed satisfies the postcondition and generates the source codewith the Double datatype which also includes a more precise andcomplete postcondition: ≤ res ∧ res ≤ ∧ res +/ − − To achieve this result, our tool first checks that the preconditionof the function call is satisfied using the Z3 solver. Then, it inlinesthe body of the function triangleSorted and computes a sound boundon the result’s uncertainty with our approximation procedure anduses it to show that the postcondition is satisfied. The error compu- ef mainFunction(a: Real, b: Real, c: Real): Real = { require (4.500005 < = a && a < = 6.5) val b = 4.0 val c = 8.5 val area = triangleSorted(a, b, c) //val area = triangleUnstable(a, b, c) area } ensuring (res = > res +/ − − def triangle(a: Real, b: Real, c: Real): Real = { require (a.in(1.0, 9.0) && b.in(1.0, 9.0) && c.in(1.0, 9.0) &&a + b > c + 1e − > b + 1e − > a + 1e − val s = (a + b + c)/2.0 sqrt(s ∗ (s − a) ∗ (s − b) ∗ (s − c)) } def triangleSorted(a: Real, b: Real, c: Real): Real = { require (a.in(1.0, 9.0) && b.in(1.0, 9.0) && c.in(1.0, 9.0) &&a + b > c + 1e − > b + 1e − > a + 1e − a < c && b < c) if (a < b) { sqrt((c+(b+a)) ∗ (a − (c − b)) ∗ (a+(c − b)) ∗ (c+(b − a)))/4.0 } else { sqrt((c+(a+b)) ∗ (b − (c − a)) ∗ (b+(c − a)) ∗ (c+(a − b))) / 4.0 }} Figure 1.
Computing the area of a triangle.tation takes into account in a sound way the input uncertainty (herean initial roundoff error on the inputs), its propagation and roundofferrors commited at each arithmetic operation. Additionally, due tothe initial roundoff error the comparison on line 24 is not exact, sothat some floating-point computations will take a different branchthan their corresponding real-valued computation. More precisely,the total error when computing the condition is . e − , ascomputed by our tool. That is, floating-point values that satisfy a < b + 7 . e − may take the else branch, even though thecorresponding real values would follow the then branch, and simi-larly in the opposite direction. Our tool verifies that the differencein the computed result in two branches, due to this divergence be-tween real arithmetic and floating point arithmetic code, remainswithin the precision requirement.Finally, our tool uses our novel range computation procedureto also compute a more precise output range than we could haveobtained in interval arithmetic. It then includes this more com-plete postcondition in the generated floating-point code. In fact, in-terval arithmetic alone computes the ranges [0 . , . and [0 . , . for using the methods triangle and triangleSorted respectively, seemingly suggesting that two methods perform en-tirely different computations. With our technique, the tool com-putes the same range [0 . , . for both methods, but shows adifference in the absolute error of the computation. For the method triangle , the verification fails, as desired, because the computed er-ror ( . e − ) exceeds the required precision bound. This resultis expected—the textbook formula for triangles is known to sufferfrom imprecision for flat triangles [31], which is somewhat rectifiedin the method triangleSorted .
3. Programs with Reals
Each program to be compiled consists of one top-level object withmethods written in a functional subset of the Scala programminglanguage [41]. All methods are functions over the
Real datatype and the user annotates them with pre- and postconditions that ex-plicitly talk about uncertainties.
Real represents ideal real num-bers without any uncertainty. We allow arithmetic expressions over
Real s with the standard arithmetic operators { + , − , ∗ , /, √} , andtogether with conditionals and function calls they form the bodyof methods. Our tool also supports immutable variable declarationsas val x = ... . This language allows the user to define a computationover real numbers. Note that this specification language is not exe-cutable.The precondition allows the user to provide a specification ofthe environment. A complete environment specification consists oflower and upper bounds for all method parameters and an upperbound on the uncertainty or noise. Range bounds are expressedwith regular comparison operators. Uncertainty is expressed withthe predicate such as x +/ − − , which denotes that the variable x is only known up to e − . Alternatively, the programmer canspecify the relative error as x +/ − − ∗ x . If no noise except forroundoff is present, roundoff errors are automatically added toinput variables.The postcondition can specify the range and the maximum un-certainty accepted on the output. In addition to the language al-lowed in the precondition, the postcondition may reference the er-rors on inputs directly in the following way: res +/ − ∗ !x , whichsays that the maximum acceptable error on the output variable res is bounded from above by . times the initial error on x . Whereasthe precondition may only talk about the ideal values, the postcon-dition can also reference the actual value directly via ∼ x . This al-lows us to assert that runtime values will not exceed a certain range,for instance. Floating-point arithmetic
Our tool and technique support in thegenerated target code any floating-point precision and in partic-ular, single and double floating-point precision as defined by theIEEE 754 floating-point standard [46]. It is also straightforward tocombine it with techniques to generate fixed-point arithmetic [13],which similarly relies on knowing ranges of variables. We assumerounding-to-nearest rounding mode and that basic arithmetic oper-ations { + , − , ∗ , /, √} are rounded correctly, which means that theresult from any such operation must be the closest representablefloating-point number. Hence, provided there is no overflow, theresult of a binary operation ◦ F satisfies x ◦ F y = ( x ◦ R y )(1 + δ ) , | δ | ≤ (cid:15) M , ◦ ∈ { + , − , ∗ , / } (1)where ◦ R is the ideal operation in real numbers and (cid:15) M is themachine epsilon that determines the upper bound on the relativeerror. This model provides a basis for our roundoff error estimates.
4. Compiling Reals to Finite Precision
Given a specification or program over reals and a possible targetfloating-point precision, our tool generates code over floating-pointnumbers that satisfy the given pre- and postconditions. Figure 2presents a high-level view of our compilation algorithm. Our toolfirst analyses the entire specification and generates one verificationcondition for each postcondition to be proven. To obtain a modularalgorithm, the tool also generates verification conditions that checkthat at each function call the precondition of the called functionis satisfied. The methods are then sorted by occuring functioncalls. This allows us to re-use already computed postconditionsof function calls in a modular analysis. If the user specified atarget precision, the remaining part of the compilation process isperfomed with respect to this precision, if not or in the case theuser specified several possible precisions, our tool will iterativelyselect the next more precise precision and attempt to compile theentire program. If compilation fails due to unsatisfied assertions,the next precision is selected. Thus, one task of our algorithm is to nput : spec: specification over Reals, prec: candidate precisions for fnc ← spec.fncsfnc.vcs = generateVCs(fnc)spec.fncs.sortBy((f1, f2) = > f1 ⊆ f2.fncCalls) while prec (cid:54) = ∅ and notProven(spec.fncs)precision = prec.nextPrecise for fnc ← spec.fncs for vc ← fnc.vcs while vc.hasNextApproximation ∧ notProven(vc)approx = getNextApproximation(vc, precision)vc.status = checkWithZ3(approx)generateSpec(fnc)generateCode(spec) Output : floating − point code with generated postconditions Figure 2.
Compilation algorithm.automatically determine the necessary least floating-point precisionto ensure the specification is met. Currently, each precision iterationis computed separately, which is not a big issue due to a smallnumber alternative floating-point targets. We did identify certainshared computations between iterations; we can exploit them in thefuture for more efficient compilation. In order for the compilationprocess to succeed, the specification has to be met with respectto a given floating-point precision, thus the principal part of ouralgorithm is spent in verification, which we describe in Section 5.We envision that in the future the compilation task will alsoinclude automatic precision-preserving code optimizations, but inthis paper we concentrate on the challenging groundwork of veri-fying the precision of code.Our tool currently generates Scala code over single, dou-ble, double-double and quad-double floating-point arithmetic. Fordouble-double and quad-double precision, which were imple-mented in software by [3], we provide a Scala interface with thegenerated code. In case the verification part of compilation fails,we nonetheless generate code (together with a failure report) withthe best postconditions our tool was able to compute. The user canthen use the generated specifications to gain insight why and whereher program does not satify requirements.Our techniques are not restricted to these floating-point preci-sions, however, and will work for any precision that follows theabstraction given in Equation (1). Furthermore, while we have im-plemented our tool to accept specifications in a domain specific lan-guage embedded in Scala and generate code in Scala, all our tech-niques apply equally to all programming languages and hardwarethat follow the floating-point abstraction we assume (Equation 1).
5. Verifying Real Programs
We will now describe the verification part of our compilation al-gorithm. In the following we will call the ideal computation thecomputation in the absence of any uncertainties and implementedin a real arithmetic, and the actual computation the computationthat will finally be executed in finite-precision and with potentiallyuncertain inputs.
Each method with its precondition P and postcondition Q impliesthe following verification condition: ∀ (cid:126)x, (cid:126)res, (cid:126)y. P ( (cid:126)x ) ∧ body ( (cid:126)x, (cid:126)y, (cid:126)res ) → Q ( (cid:126)x, (cid:126)res ) (*) a < = x && x < = b x ∈ [ a, b ] x +/ − k x ◦ = x + err x ∧ err x ∈ [ − k, k ] x +/ − m ∗ x x ◦ = x + err x ∧ err x ∈ [ −| mx | , | mx | ] ∼ x x ◦ !x err x x (cid:5) y ( x (cid:5) y ) ∧ ( x ◦ (cid:5) y ◦ )(1 + δ ) sqrt(x) sqrt ( x ) ∧ sqrt ( x ◦ )(1 + δ ) val z = x z = x ∧ z ◦ = x ◦ if (c(x)) e1(x) (( c ( x ) ∧ e x )) ∨ ( ¬ c ( x ) ∧ e x ))) ∧ else e2(x) (( c ◦ ( x ◦ ) ∧ e ◦ ( x ◦ )) ∨ ( ¬ c ◦ ( x ◦ ) ∧ e ◦ ( x ◦ ))) g(x) g ( x ) ∧ g ◦ ( x ◦ ) (cid:5) ∈ { + , − , ∗ , / }− (cid:15) m ≤ δ i ∧ δ i ≤ (cid:15) m , all δ are fresh cond ◦ and e ◦ denote functions with roundoff errors at each step Table 1.
Semantics of our specification language.where (cid:126)x, (cid:126)res, (cid:126)y denote the input, output and local variables respec-tively.Table 1 summarizes how verification constraints are generatedfrom our specification language. Each variable x in the specificationcorresponds to two real-valued variables x, x ◦ , the ideal one in theabsence of uncertainties and roundoff errors and the actual one,computed by the compiled program. Note that the ideal and actualvariables are related only through the error bounds in the pre- andpostconditions, which allows for the ideal and actual executions totake different paths.In the method body we have to take into account roundofferrors from arithmetic operations and the propagation of existingerrors. Our system currently supports operations { + , − , ∗ , /, √} ,but these can be in principle extended to elementary functions, forinstance by encoding them via Taylor expansions [37].Note that the resulting verification conditions are parametric inthe machine epsilon. In order to give feedback to developers and to facilitate automaticmodular analysis, our tool also provides automatic specificationgeneration. By this we mean that the programmer still needs toprovide the environment specification in form of preconditions, butour tool automatically computes a precise postcondition.Formally, we can rewrite the constraint (*) as ∀ (cid:126)x, res. ( ∃ (cid:126)y. P ( (cid:126)x ) ∧ body ( (cid:126)x, (cid:126)y, res )) → Q ( (cid:126)x, res ) where Q is now unknown. We obtain the most precise post-condition Q by applying quantifier elimination (QE) to P ( (cid:126)x ) ∧ body ( (cid:126)x, (cid:126)y, res ) and eliminate (cid:126)y . The theory of arithmetic over re-als admits QE so it is theoretically possible to use this approach.We do not currently use a full QE procedure for specification gen-eration, as it is expensive and it is not clear whether the returnedexpressions would be of a suitable format. Instead, we use our ap-proximation approach which computes ranges and maximum errorsin a forward fashion and allows us to compute an (over) approxi-mation of a postcondition of the form res ∈ [ a, b ] ∧ res ± u . For small functions we can already prove interesting properties byusing the exact encoding of the problem just described and dis-charging the verification constraints with Z3. Consider the follow- ef getNextApproximation(vc, precision): paths path-wise merging functions inlining postcondition inlining full function arithmetic approx. error approx. error & range first approximation getRange ! evalWithError ! getPathError ! Figure 3.
Approximation pipeline.ing code a programmer may write to implement the third B-splinebasic function which is commonly used in signal processing [29]. def bspline3(u: Real): Real = { require (0 ≤ u && u ≤ ± − − u ∗ u ∗ u / 6.0 } ensuring (res ≥ − ≤ res && res ≤ ± − Functions and the correspoding verification conditions of this com-plexity are already within the possibilities of the nonlinear solverwithin Z3. For more complex functions however, Z3 does not(yet) provide an answer in a reasonable time, or returns unknown.Whether alternative techniques in SMT solvers can help in suchcases remains to be seen [7, 30]. We here provide an approachbased on step-wise approximation that addresses the difficulty ofgeneral-purpose constraint solving.
In order to nontheless verify interesting programs, we have de-veloped an approximation procedure that computes a sound over-approximation of the range of an expression and of the uncertaintyon the output. This procedure is a forward computation and wealso use it to generate specifications automatically. We describe theapproximation procedure in detail in Section 6, for now we will as-sume that it exists and, given a precondition P and an expression expr , computes a sound bound on the output range and its associ-ated uncertainty: ([ a, b ] , err ) = evalWithError ( P, expr ) ⇔∀ (cid:126)x, (cid:126)x ◦ , res, res ◦ .P ( (cid:126)x, (cid:126)x ◦ ) ∧ res = expr ( (cid:126)x ) ∧ res ◦ = expr ◦ ( (cid:126)x ◦ ) → ¬ ( res < a ∨ b < res ) ∧ | res − res ◦ | < err We have identified three possibiities for approximation: nonlin-ear arithmetic, function calls, and paths and each can be approxi-mated at different levels. We have observed in our experiments, that“one size does not fit all” and a combination of different approxi-mations is most successful in proving the verification conditions weencountered. For each verification condition we thus construct ap-proximations until Z3 is able to prove one, or until we run out of ap-proximations where we report the verification as failed. We can thusview verification as a stream of approximations to be proven. We il-lustrate the pipeline that computes the different approximations inFigure 3. The routines getPathError, getRange and evalWithError aredescribed in the following sections in more detail.The first approximation (indicated by the long arrow in Fig-ure 3) is to use Z3 alone on the entire constraint constructed by therules in Table 1. This is indeed an approximation, as all functioncalls are treated as uninterpreted functions in this case. As notedbefore, this approach only works in very simple cases or when nouncertainties and no functions are present.
Function calls
If the verification constraint contains functioncalls and the first approximation failed, our tool will attempt to in-line postconditions and pass on the resulting constraint down the approximation pipeline. We support inlining of both user-providedpostconditions and postconditions computed by our own specifi-cation generation procedure. If this still is not precise enough, weinline the entire function body.Postcondition inlining is implemented by replacing the functioncall with a fresh variable and constraining it with the postcondi-tion. Thus, if verification suceeds with inlining the postcondition,we avoid having to consider each path of the inlined function sep-arately and can perform modular verification avoiding a potentialpath explosion problem. Such modular verification is not feasiblewhen postconditions are too imprecise and we plan to explore thegeneration of more precise postconditions in the future. One stepin this direction is to allow postconditions that are parametric inthe initial errors, for example with the operator ! x introduced inSection 3. While our tool currently supports postcondition inliningwith such postconditions, we do not yet generate these automati-cally. Arithmetic
The arithmetic part of the verification constraints gen-erated by Table 1 can be essentially divided into the ideal part andthe actual part, which includes roundoff errors at each computationstep. The ideal part determines whether the ideal range constraintsin the postcondition are satisfied and the actual part determineswhether the uncertainty part of the postcondition is satisfied. Wecan use our procedure presented in Section 6 to compute a soundapproximation of both the result’s range as well as its uncertainty.Based on this, our tool constructs two approximations. For the first,the ideal part of the constraint is left untouched and the actual partis replaced by the computed uncertainty bound. This effectively re-moves a large number of variables and many times suffiently sim-plifies the constraint for Z3 to succeed. If this fails, our tool addi-tionally replaces the ideal part by the computed range constraint.Note that the second approximation may not have enough informa-tion to prove a more complex postconditions, as correlation infor-mation is lost. We note that the computation of ranges and errorsis the same for both approximations and thus trying both does notaffect efficiency significantly.
Paths
In the case of several paths through the program, we havethe option to consider each path separately or to merge results ateach join in the control flow graph. This introduces a tradeoff be-tween efficiency and precision, since on one hand, considering eachpath separately leads to an exponential number of paths to consider.On the other hand, merging at each join looses correlation informa-tion between variables which may be necessary to prove certainproperties. Our approximation pipeline chooses merging first, be-fore resorting to a path-by-path verification in case of failure. Webelieve that other techniques for exploring the path space could alsobe integrated into our tool [9, 32]. Another possible improvementare heuristics that select a different order of approximations de-pending on particular characteristics of the verification condition.
Example
We illustrate the verification algorithm on the exam-ple in Figure 4. The functions sineTaylor and sineOrder3 are veri-fied first since they do not contain function calls. Verification withthe full verification constraint fails. Next, our tool computes the er-rors on the output and Z3 suceeds to prove the resulting constraintwith the ideal part untouched. From this approximation our tooldirectly computes a new, more precise postcondition, in particularit can narrow the resulting error to . Next, our tool con-siders the comparison function. Inlining only the postcondition isnot enough in this case, but computing the error approximation onthe inlined functions suceeds in verifying the postcondition. Z3 isagain able to verify the real-valued portion independently. Finally,the tool verifies that the preconditions of the function calls are satis-fied by using Z3 alone. Verification of the function comparisonInvalid fails with all approximations so that our tool reports unknown . It also ef comparisonValid(x: Real): Real = { require ( − < x && x < val z1 = sineTaylor(x) val z2 = sineOrder3(x)z1 − z2 } ensuring (res = > ∼ res < = 0.1) def comparisonInvalid(x: Real): Real = { require ( − < x && x < val z1 = sineTaylor(x) val z2 = sineOrder3(x)z1 − z2 } ensuring (res = > ∼ res < = 0.01) // counterexample: 1.0 def sineTaylor(x: Real): Real = { require ( − < x && x < − (x ∗ x ∗ x)/6.0 + (x ∗ x ∗ x ∗ x ∗ x)/120.0 − (x ∗ x ∗ x ∗ x ∗ x ∗ x ∗ x)/5040.0 } ensuring (res = > − < res && res < − − def sineOrder3(x: Real): Real = { require ( − < x && x < ∗ x − ∗ (x ∗ x ∗ x) } ensuring (res = > − < res && res < − − Figure 4.
Different polynomial approximations of sine.provides the counterexample it obtains from Z3 ( x = 1 . ), whichin this case is a valid counterexample. Our procedure is sound because our constraints overapproximatethe actual errors. Furthemore, even in the full constraint as gener-ated from Table 1, roundoff errors are overapproximated since weassume the worst-case error bound at each step. While this ensuressoundness, it also introduces incompleteness, as we may fail to val-idate a specification because our overapproximation is too large.This implies that counterexamples reported by Z3 are in generalonly valid, if they disprove the ideal real-valued part of the verifica-tion constraint. In general, this is easy to distinguish from the casewhere verification fails due to too large error bounds, as Z3 preferssimpler counterexamples over complex ones and thus chooses allerror variables to be zero if the ideal part is being violated.
6. Solving Nonlinear Constraints
Having given an overview of the approximation pipeline, we nowdescribe the computation of the approximation for nonlinear arith-metic, which corresponds to the last box in Figure 3. For complete-ness of presentation, we first review interval and affine arithmeticwhich are common choices for performing sound arithmetic com-putations and which we also use as part of our technique. We thenpresent our novel procedure for computing the output range of anonlinear expression given ranges for its inputs that can be a moreprecise substitute for interval or affine arithmetic. Finally, we con-tinue with a procedure that computes a sound overapproximationof the uncertainty on the result of a nonlinear expression.One possibility to perform guaranteed computations is to usestandard interval arithmetic [39]. Interval arithmetic computes abounding interval for each basic operation as x ◦ y = [ min ( x ◦ y ) , max ( x ◦ y )] ◦ ∈ { + , − , ∗ , / } and analogously for square root. Affine arithmetic was originallyintroduced in [15] and addresses the difficulty of interval arithmeticin handling correlations between variables. Affine arithmetic repre- sents possible values of variables as affine forms ˆ x = x + n (cid:88) i =1 x i (cid:15) i where x denotes the central value (of the represented interval) andeach noise symbol (cid:15) i is a formal variable denoting a deviation fromthe central value, intended to range over [ − , . The maximummagnitude of each noise term is given by the corresponding x i .Note that the sign of x i does not matter in isolation, it does, how-ever, reflect the relative dependence between values. For example,take x = x + x (cid:15) , then x − x = x + x (cid:15) − ( x + x (cid:15) )= x − x + x (cid:15) − x (cid:15) = 0 If we subtracted x = x − x (cid:15) instead, the resulting interval wouldhave width ∗ x and not zero.The range represented by an affine form is computed as [ˆ x ] = [ x − rad (ˆ x ) , x + rad (ˆ x )] , rad (ˆ x ) = n (cid:88) i =1 | x i | A general affine operation α ˆ x + β ˆ y + ζ consists of addition, sub-traction, addition of a constant ( ζ ) or multiplication by a constant( α, β ). Expanding the affine forms ˆ x and ˆ y we get α ˆ x + β ˆ y + ζ = ( αx + βy + ζ ) + n (cid:88) i =1 ( αx i + βy i ) (cid:15) i An additional motivation for using affine arithmetic is that differ-ent contributions to the range it represents remain, at least partly,separated. This information can be used for instance to help iden-tify the major contributor of a result’s uncertainty or to separatecontributions from external uncertainties from roundoff errors.
The goal of this procedure is to perform a forward-computation todetermine the real-valued range of a nonlinear arithmetic expres-sion given ranges for its inputs. Two common possibilities are in-terval and affine arithmetic, but they tend to overapproximate theresulting range, especially if the input intervals are not sufficientlysmall (order 1). Affine arithmetic improves over interval arithmeticsomewhat by tracking linear correlations, but in the case of non-linear expressions the results can become actually worse than forinterval arithmetic.
Observation:
A nonlinear theorem prover such as the one thatcomes with Z3 can decide with relatively good precision whethera given bound is sound or not. That is we can check with a proverwhether for an expression e the range [ a, b ] is a sound interval en-closure. This observation is the basis of our range computation.The input to our algorithm is a nonlinear expression expr and aprecondition P on its inputs, which specifies, among possibly otherconstraints, ranges on all input variables (cid:126)x . The output is an interval [ a, b ] which satisfies the following: [ a, b ] = getRange ( P, expr ) ⇔∀ (cid:126)x, res.P ( (cid:126)x ) ∧ res = expr ( (cid:126)x ) → ¬ ( res < a ∨ b < res ) The algorithm for computing the lower bound of a range is givenin Figure 5. The computation for the upper bound is symmetric. Foreach range to be computed, our tool first computes an initial soundestimate of the range with interval arithmetic. It then performs aninitial quick check to test whether the computed first approximationbounds are already tight. If not, it uses the first approximation as thestarting point and then narrows down the lower and upper bounds ef getRange(expr, precondition, precision, maxIterations):z3.assertConstraint(precondition)[aInit, bInit] = evalInterval(expr, precondition.ranges); //lower bound if z3.checkSat(expr a + precision) == UNSATa = aInitb = bInitnumIterations = 0 while (b − a) < precision ∧ numIterations < maxIterationsmid = a + (b − a) / 2numIterations++z3.checkSat(expr mid) matchcase SAT ⇒ b = mid case UNSAT ⇒ a = mid case Unknown ⇒ breakaNew = a else aNew = aInitbNew = ... //upper bound symmetrically return : [aNew, bNew] Figure 5.
Algorithm for computing the range of an expression.using a binary search. At each step of the binary search our tooluses Z3 to confirm or reject the newly proposed bound.The search stops when either Z3 fails, i.e. returns unknown fora query or cannot answer within a given timeout, the differencebetween subsequent bounds is smaller than a precision threshold,or the maximum number of iterations is reached. This stoppingcriterion can be set dynamically.
Additional constraints
In addition to the input ranges, the pre-condition may also contain further constraints on the variables. Forexample consider again the method triangle in Figure 1. The pre-condition bounds the inputs as a, b, c ∈ [1 , , but the formula isuseful only for valid triangles, i.e. when every two sides togetherare longer than the third. If not, we will get an error at the lat-est when we try to take the square root of a negative number. Ininterval-based approaches we can only consider input intervals thatsatisfy this constraint for all values, and thus have to check sev-eral (and possibly many) cases. In our approach, since we are usingZ3 to check the soundness of bounds, we can assert the additionalconstraints up-front and then all subsequent checks are performedwith respect to all additional and initial contraints. This allows us toavoid interval subdivisions due to imprecisions or problem specificconstraints such as those in the triangle example. This becomes es-pecially valuable in the presence of multiple variables, where wemay need an exponential number of subdivisions. We now describe our approximation procedure which, for a givenexpression expr and a precondition P on the inputs, computesthe range and error on the output. More formally, our proceduresatisfies the following: ([ a, b ] , err ) = evalWithError ( P, expr ) ⇔∀ (cid:126)x, (cid:126)x ◦ , res, res ◦ .P ( (cid:126)x, (cid:126)x ◦ ) ∧ res = expr ( (cid:126)x ) ∧ res ◦ = expr ◦ ( (cid:126)x ◦ ) → ¬ ( res < a ∨ b < res ) ∧ | res − res ◦ | < err where expr ◦ represents the expression evaluated in floating-pointarithmetic and (cid:126)x, (cid:126)x ◦ are the ideal and actual variables. The precon-dition specifies the ranges and uncertainties of initial variables and other additional constraints on the ideal variables. The uncertaintyspecification is necessary, as it related the ideal and actual variables.The idea of our procedure is to “execute” a computation whilekeeping track of the output range of the current expression andits associated errors. At each arithmetic operation, we propagateexisting errors, compute an upper bound on the roundoff errorand add it to the overall errors. Since the roundoff error dependsproportionally on the range of values, we also need to keep track ofthe ranges as precisely as possible.Our procedure is build on the abstraction that a computation isan ideal computation plus or minus some uncertainty. The abstrac-tion of floating-point roundoff errors that we chose also follows thisseparation: fl ( x (cid:5) y ) = ( x (cid:5) y )(1 + δ ) = ( x (cid:5) y ) + ( x (cid:5) y ) δ for δ ∈ [ − (cid:15) m , (cid:15) m ] and (cid:5) ∈ { + , − , ∗ , / } . This allows us to treat alluncertainties in a unified manner.Our procedure builds on the idea of the SmartFloat datatype [12],which uses affine arithmetic to track both the range and the errors.For nonlinear operations, however, the so computed ranges becomevery pessimistic quickly and the error computation may also sufferfrom this imprecision. We observed that since the errors tend to berelatively small, this imprecision does not affect the error propaga-tion itself to such an extent. If the initial errors are small (less thanone), multiplied nonlinear terms tend to be even smaller, whereasif the affine terms are larger than one, the nonlinear terms grow. Wethus concentrate on improving the ideal range of values and use ournovel range computation procedure for this part and leave the errorpropagation with affine arithmetic as in [12].In our adaptation, we represent every variable and intermediatecomputation result as a datatype with the following components: x : ( range : Interval, ˆ err : AffineF orm ) where range is the range of this variable, computed as described inSection 6.1 and ˆ err is the affine form representing the errors. The(overapproximation) of the actual range including all uncertaintiesis then given by totalRange = range + [ ˆ err ] , where ˆ err denotesthe interval represented by the affine form. Roundoff error computation
Roundoff errors are computed ateach computation step as ρ = δ ∗ maxAbs ( totalRange ) where δ is the machine epsilon, and added to ˆ err as a fresh noiseterm. Note that this roundoff error computation makes our errorcomputation parametric in the floating-point precision. Error propagation
For affine operations addition, subtraction,and multiplication by a constant factor the propagated errors arecomputed term-wise and thus as for standard affine arithmetic. Werefer the reader to [12, 15] for further details and describe hereonly the propagation for nonlinear arithmetic. For multiplication,division and square root, the magnitude of errors also dependson the ranges of variables. Since our ranges are not affine termsthemselves, propagation has to be adjusted. In the following, wedenote the range of a variable x by [ x ] and its associated error bythe affine form ˆ err x . When we write [ x ] ∗ ˆ err y we mean that theinterval [ x ] is converted into an affine form and the multiplicationis performed in affine arithmetic.Multiplication is computed as x ∗ y = ([ x ] + ˆ err x )([ y ] + ˆ err y )= [ x ] ∗ [ y ] + [ x ] ∗ ˆ err y + [ y ] ∗ ˆ err x + ˆ err x ∗ ˆ err y + ρ where ρ is the new roundoff error. Thus the first term contributesto the ideal range and the remaining three to the error affine form.The larger the factors [ x ] and [ y ] are, the larger the finally computed rrors will be. In order to keep the overapproximation as small aspossible, we evaluate [ x ] and [ y ] with our new range computation.Division is computed as xy = x ∗ y = ([ x ] + ˆ err x )([1 /y ] + ˆ err /y )= [ x ] ∗ [ 1 y ] + [ x ] ∗ ˆ err y + [ 1 y ] ∗ ˆ err x + ˆ err x ∗ ˆ err y + ρ For square root, we first compute an affine approximation of squareroot as in [12] √ x = α ∗ x + ζ + θ and then perform the affine multiplication term wise. Overflows and NaN
Our procedure allows us to detect potentialoverflows, division by zero and square root of a negative value, asour tool computes ranges of all intermediate values. We currentlyreport these issues as warnings to the user.
The limitation of this approach is clearly the ability of Z3 to checkour constraints. We found its capabilities satisfactory, although weexpect the performance to still significantly improve. To emphasizethe difference to the constraints that are defined by Table 1, theconstraints we use here do not add errors at each step and thusthe number of variables is reduced significantly. We also foundseveral transformations helpful, such as rewriting powers (e.g. x ∗ x ∗ x to x ), multiplying out products and avoiding non-strictcomparisons in the precondition, although the benefits were notentirely consistent. Note that at each step of our error computation,our tool computes the current range. Thus, even if Z3 fails to tightenthe bound for some expressions, we still compute more precisebounds than interval arithmetic overall in most cases, as the rangesof the remaining subexpressions have already been computed moreprecisely.
7. Conditional Statements
In this Section we consider the difference between the ideal andactual computation due to uncertainties on computing branch con-ditions and the resulting different paths taken. We note that thefull constraint constructed according to Section 3 automatically in-cludes this error. Recall that the ideal and actual computations areindependent except for the initial conditions, so that it is possiblethat they follow different paths through the program.In the case of approximation, however, we compute the erroron individual paths and have to consider the error due to divergingpaths separately. If a method encodes a continous function in theusual mathematical sense then we note that we only need to quan-tify errors for each path separately. Thus, if have a method [9] todetermine whether a function with conditional statements is con-tinuous, then our approach described so far is sufficient to providesound error bounds. For the case where such a procedure does notexist or fails to provide an answer, for example due to nonlinear-ity, or the function simply is not continuous, we propose the fol-lowing algorithm to explicitly compute the difference between theideal and the actual computation across paths. Note that we do notassume continuity, i.e. the algorithm allows us to compute errorbounds even in the case on non-continous functions.For simplicity, we present here the algorithm for the case of oneconditional statement: if (c(x) <
0) f1(x) else f2(x)
It generalizes readily to more complex expressions. W.l.o.g. weassume that the condition is of the form c ( x ) < . Indeed, any def getPathError: Input : pre ( x ∈ [ a, b ] ∧ x ± n )program ( if (cond(x) <
0) f1(x) else f2(x)) val pathError1 = computePathError(pre, cond, f1, f2) val pathError2 = computePathError(pre, ¬ cond, f2, f1) return max (pathError1, pathError2) def computePathError(pre, c, f1, f2):([c], err c ) = evalWithError(pre, c) ([f2] float , err float ) =evalWithError(pre ∧ c(x) ∈ [0, err c ], f2) [f1] real =getRange(pre ∧ c(x) ∈ [ − err c , 0], f1) return : max | [f1] real − ([f2] float + err float ) | Figure 6.
Computing error due to diverging paths.conditional of the form c ( x ) == 0 would yield different resultsfor the ideal and actual computation for nearly any input, so we donot allow it in our specification language.The actual computation commits a certain error when comput-ing the condition of the branch and it is this error that causes someexecutions to follow a different branch than the corresponding idealone would. Consider the case where the ideal computation evalu-ates f1, but the actual one evaluates f2. Algorithm 6 gives the com-putation of the path error in this case. The idea is to compute theranges of f1 and f2, but only for the inputs that could be diverg-ing. The final error is then the maximum difference of these value.The algorithm extends naturally to several variables. In the case ofseveral paths through the program, this error has to be, in princi-ple, computed for each combination of paths. We use Z3 to rule outinfeasible paths up front so that the path error computation is onlyperformed for those paths that are actually feasible.We have currently implemented this approach in our tool forthe case when we use merging to handle paths in order to avoidhaving to consider an exponential number of path combinations.We also use a higher default precision and number of iterationsthreshold during the binary search in the range computation as thiscomputation requires in general very tight intervals for each path.We identify two challenges for performing this computation:1. As soon as the program has multiple variables, the inputs for thedifferent branches are not two-dimensional intervals anymore,which makes an accurate evaluation of the individual pathsdifficult in standard interval arithmetic.2. The inputs for the two branches are inter-dependent. Thus,simply evaluating the two branches with inputs that are in thecorrect ranges, but are not correlated, yields pessimistic resultswhen computing the final difference (line 16).We overcome the first challenge with our range computationwhich takes into account additional constraints. For the secondchallenge, we use our range computation as well, however unfortu-nately Z3 fails to tighten the final range to a satisfactory precisiondue to timeouts. We still obtain much better error estimates thanwith interval arithmetic alone, as the ranges of values for the indi-vidual paths are already computed much more precisely. We reportin Section 8 on the type of programs whose verification is alreadyin our reach today.
8. Experiments
The examples in Figure 1 and 4 and Section 5.3 provide an ideaof the type of programs our tool is currently able to verify fullyautomatically. The B-spline example from Section 5.3 is the largest enchmark Our error (IA only) Simulated errordoppler1* 2.36e-6 5.97e-7doppler1 4.92e-13 (4.95e-13) 7.11e-14doppler2* 6.21e-5 1.85-5doppler2 1.29e-12 1.14e-13doppler3* 1.23e-4 5.96e-5doppler3 2.03e-13 (2.05e-13) 4.27e-14rigidBody1* 9.21e-7 8.24e-7rigidBody1 5.08e-13 2.28e-13rigidBody2* 1.51e-4 1.25e-4rigidBody2 6.48e-11 2.19e-11jetEngine* 0.15 (-) 3.58e-5jetEngine 1.62e-8 (-) 5.46e-12turbine1* 4.86e-6 3.71e-7turbine1 1.25e-13 (1.38e-13) 1.07e-14turbine2* 8.05e-6 7.66e-7turbine2 1.76e-13 (1.96e-13) 1.43e-14turbine3* 3.35e-6 1.04e-6turbine3 8.50e-14 (9.47e-14) 5.33e-15verhulst* 2.82e-4 2.40e-4verhulst 6.82e-16 2.23e-16predatorPrey* 9.22e-5 8.61e-5predatorPrey 2.94e-16 (2.96e-16) 1.12e-16carbonGas* 2114297.84 168874.70carbonGas 4.64e-8 (5.04e-8) 3.73e-9Sine (single) 1.03e-6 (1.57e-6) 1.79e-7Sine 9.57e-16 (1.46e-15) 4.45e-16Sqrt (single) 9.03e-7 (9.52e-7) 2.45e-7Sqrt 8.41e-16 (8.87e-16) 4.45e-16Sine, order 3 (single) 1.19e-6 (1.55e-6) 2.12e-7Sine, order 3 1.11e-15 (1.44e-15) 3.34e-16 Table 3.
Comparison of errors computed with our procedureagainst simulated errors. Simulations were performed with ran-dom inputs. (*) indicates that inputs have external uncertainties as-sociated.meaningful example we were able to find that Z3 alone could verifyin the presence of uncertainties. For all other cases, it was necessaryto use our approximation methods. To evaluate our range and error computation technique we havechosen several nonlinear expressions commonly used in physics,biology and chemistry [40, 43, 49] as benchmark functions, as wellas benchmarks used in control systems [1] and suitable benchmarksfrom [18].Experiments were performed on a desktop computerrunning Ubuntu 12.04.1 with a 3.5GHz i7 processor and 16GB ofRAM. Running times highly depend on the timeout used for Z3.Our default setting is 1 second; we did not find much improvementin the success rate above this threshold.
Range computation
Stepwise estimation of errors crucially de-pends on the estimate of the ranges of variables. The strength ofusing a constraint solver such as Z3 is that it can perform such esti-mation while taking into account the precise dependencies between variables in preconditions and path conditions. Table 2 comparesresults of our range computation procedure described in Section 6against ranges obtained with standard interval arithmetic. Intervalarithmetic is one of the methods used for step-wise range estima-tion; an alternative being affine arithmetic. We have also experi-mented with an affine arithmetic implementation [12]. However, wefound that affine arithmetic gives more pessimistic results for com-puting ranges for non-linear benchmarks. We believe that this isdue to imprecision in computing nonlinear operations. Note, how-ever, that we still use affine arithmetic to estimate errors given thecomputed ranges.We set the default precision threshold to − and maximumnumber of iterations for the binary search to 50. To obtain anidea about the ranges of our functions, we have also computeda lower bound on the range using simulations with randominputs and with exact rational arithmetic evaluation of expressions.We observe that our range computation can significantly improveover standard interval bounds. The jetEngine benchmark is a notableexample, where interval arithmetic yields the bound [ −∞ , ∞ ] ,but our procedure can still provide bounds that are quite close tothe true range. Running times are below 7 seconds for the mostcomplex benchmarks, except for jetEngine which runs in about 1minute due to timeouts from Z3 for some intermediate ranges. Error computation
Table 3 compares uncertainties computed byour tool against maximum uncertainties obtained through exten-sive simulation with random inputs. We ran the simulation inparallel with rational and their corresponding floating-point valueand obtained the error by taking the difference in the result. Bench-marks marked with (*) have added initial uncertainties. Unless oth-erwise indicated, we used double floating-point precision. To ourknowledge this is the first quantitative comparison of an error com-putation precision with (an approximation) of the true errors onsuch benchmarks. Except for the benchmarks jetEngine ∗ our com-puted uncertainties are within an order and many times even closerto the underapproximation of the true errors provided by simula-tion. In the case of the jetEngine ∗ benchmark, we believe that theimprecision is mainly due to its complexity and subsequent failuresof Z3. The values in parentheses in the second column indicate er-rors computed if ranges at each arithmetic operation are computedusing interval arithmetic alone. While we have not attempted to im-prove the affine arithmetic-based error computation from [12], wecan see that in some cases a more precise range computation cangain us improvements. The full effect of the imprecision of stan-dard range computation appears when, due to this imprecision, weobtain possible errors such as division-by-zero or square root of anegative number errors. The first case happens in the case of thenon-linear jetEngine benchmark, so with interval arithmetic alonewe would therefore not obtain any meaningful result. Similarly, forthe triangle example from Section 2, without being able to con-strain the inputs to form valid triangles, we cannot compute anyerror bound, because the radicand becomes possibly negative.Table 4 presents another relevant experiment, evaluating theability to use additional constraints during our range computation.We use the triangle example from Section 2 with additional con-straints allowing increasingly flat triangles by setting the thresholdon line 13 ( a + b > c + 1e − ) to the different values given in the firstcolumn. As the triangles become flatter, we observe an expected in-crease in uncertainty on the input since the formula becomes moreprone to roundoff errors. At threshold − our range computationfails to provide the necessary precision and the radicand becomespossibly negative. (We used double precision in this example aswell.) enchmark Our range interval arithmetic Simulated rangedoppler1 [-137.639, -0.033951] [-158.720, -0.029442] [-136.346, -0.035273]doppler2 [-230.991, -0.022729] [-276.077, -0.019017] [-227.841,-0.023235]doppler3 [-83.066, -0.50744] [-96.295, -0.43773] [-82.624, -0.51570]rigidBody1 [-705.0, 705.0] [-705.0, 705.0] [-697.132, 694.508]rigidBody2 [-56010.1, 58740.0] [-58740.0, 58740.0] [-54997.635, 57938.052]jetEngine [-1987.022, 5099.243] [ −∞ , ∞ ] [-1779.551, 4813.564]turbine1 [-18.526, -1.9916] [-58.330, -1.5505] [-18.284, -1.9946]turbine2 [-28.555, 3.8223] [-29.437, 80.993] [-28.528, 3.8107]turbine3 [0.57172, 11.428] [0.46610, 40.376] [0.61170, 11.380]verhulst [0.31489, 1.1009] [0.31489, 1.1009] [0.36685,0.94492]predatorPrey [0.039677, 0.33550] [0.037277, 0.35711] [0.039669,0.33558]carbonGas [4.3032 e6, 1.6740 e7] [2.0974 e6, 3.4344 e7] [4.1508 e6, 1.69074 e7]Sine [-1.0093, 1.0093] [-2.3012, 2.3012] [-1.0093, 1.0093]Sqrt [1.0, 1.3985] [0.83593, 1.5625] [1.0, 1.3985]Sine (order 3 approx.) [-1.0001, 1.0001] [-2.9420, 2.9420] [-1.0, 1.0] Table 2.
Comparison of ranges computed with out procedure against interval arithmetic and simulation. Simulations were performed with random inputs. Ranges are rounded outwards.Benchmark Range Max. abs. errortriangle1 (0.1) [0.29432, 35.0741] 2.72e-11triangle2 (1e-2) [0.099375, 35.0741] 8.04e-11triangle3 (1e-3) [3.16031e-2, 35.0741] 2.53e-10triangle4 (1e-4) [9.9993e-3, 35.0741] 7.99e-10triangle5 (1e-5) [3.1622e-3, 35.0741] 2.53e-9triangle6 (1e-6) [9.9988e-4, 35.0741] 7.99e-9triangle7 (1e-7) [3.1567e-4, 35.0741] 2.54e-8triangle8 (1e-8) [9.8888e-5, 35.0741] 8.08e-8triangle9 (1e-9) [3.0517e-5, 35.0741] 2.62e-7triangle10 (1e-10) - - Table 4.
Ranges and errors for increasingly flat triangles. All val-ues are rounded outwards. Interval arithmetic alone fails to provideany result.
Figure 7 presents several examples to evaluate our error compu-tation procedure across different paths from Section 7. The firstmethod cav10 [22] has been used before as a benchmark functionfor computing the output range. Our tool can verify the given post-condition immediately. Note that the error on the result is actuallyas large as the result itself, since the method is non-continuous,an aspect that has been ignored in previous work, but that ourtool detects automatically. The method squareRoot3 is also an non-continous function that computes the square root of x using anapproximation for small values and the regular library method oth-erwise. Note the additional uncertainty on the input, which couldoccur for instance if this method is used in an embedded controller.Our tool can verify the given spefication. If we change the conditionon line 10 to x < − however, verification fails. In this fashion,we can use our tool to determine the appropriate branch conditionto meet the precision requirement. The above examples verify allin under 5 seconds. Finally, the smartRoot method computes oneroot of a quadratic equation using the well-known more precise def cav10(x: Real): Real = { require (x.in(0, 10)) if (x ∗ x − x > = 0) x/10 else x ∗ x + 2 } ensuring (res = > < = res && res < = 3.0 && res +/ − def squareRoot3(x: Real): Real = { require ( x.in(0,10) && x +/ − −
10 ) if (x < − ∗ x else sqrt(1 + x) } ensuring ( res = > res +/ − − def smartRoot(a: Real, b: Real, c: Real): Real = { require (3 < = a && a < = 3 && 3.5 < = b && b < = 3.5 &&c.in( −
2, 2) && b ∗ b − a ∗ c ∗ > val discr = b ∗ b − a ∗ c ∗ if (b ∗ b − a ∗ c > { if (b > ∗ − b − sqrt(discr)) else if (b < − b + sqrt(discr))/(a ∗ else ( − b + sqrt(discr))/(a ∗ } else { ( − b + sqrt(discr))/(a ∗ } } ensuring (res = > res +/ − − Figure 7.
Path error computation examples.method from [23]. We are currently not aware of an automatic toolto prove programs continuous in the presence of nonlinear arith-metic, so that we need to compute the error across different pathsas well. Our tool succeeds in verifying the postcondition in about25s. The rather long running time is due to the complexity of theconditions when computing the error across paths, and thus Z3’slonger response time, and a number of Z3 timeouts (Z3 timeouthere means merely that some ranges have not been tightend to thebest precision). In the future, we envision that optimizations that elect alternative approximations can be performed automaticallyby a trustworthy compiler.
9. Related work
Current approaches for verifying floating-point code include ab-stract interpretation, interactive theorem proving and decision pro-cedures, which we survey in this section. We are not aware of workthat would automatically integrate reasoning about uncertainties.
Abstract interpretation (AI)
Abstract domains that are soundwith respect to floating-point computations can prove bounds onthe ranges of variables [5, 11, 20, 28, 38]. The only work in this areathat can also quantify roundoff errors is the tool Fluctuat[17, 24].These techniques use interval or affine arithmetic and together withthe required join and meet operations may yield too pessimisticresults. [42] improves the precision of Fluctuat by refining the in-put domains with a constraint solver. Our approach can be viewedas approaching the problem from a different end, starting with anexact constraint and then using approximation until the solver suc-ceeds. Unlike AI tools in general, our system currently handlesonly functional code, in particular it does not handle loops. If theuser can provide inductive postconditions, then we can still provethe code correct, but we do not in general discover these ourselves.Our focus lies on proving precise bounds on the ranges in the pres-ence of nonlinear computations and the quantification of roundofferrors and other uncertainties.
Theorem proving
The Gappa tool [14, 35] generates a proofcheckable by the interactive theorem prover Coq from source codewith specifications. It can reason about properties that can be re-duced to reasoning about ranges and errors, but targets, very preciseproperties of specialized functions, such as software implementa-tions of elementary functions. The specification itself requires ex-pertize and the proofs human intervention. A similar approach istaken by [2] which generate verification conditions that are dis-charged by various theorem provers. Harisson has also done signif-icant work on proving floating-point programs in the HOL Lighttheorem prover [26].Our approach makes a different compromise on the precisionvs. automation tradeoff, by being less precise, but automatic. TheGappa approach can be used complementary to ours, in that if wedetect that more precision is needed, Gappa is employed by anexpert user on selected methods, and the results are then used byour tool instead of automatically computed specifications.
Range computation
The Gappa tool and most constraint solversinternally use interval arithmetic for sound range computations,whose limitations are well-known. [19] describes an arithmeticbased on function enclosures and [37] use an arithmetic basedon taylor series as an alternative. This approach is useful whenchecking a constraint, but is not suitable for a forward computation of ranges and errors.
Decision procedures
An alternative approach to verification viarange computation are floating-point decision procedures. Bit-precise constraints, however, become very large quickly. [8] ad-dresses this problem by using a combination of over- and underap-proximations. [25] present an alternative approach in combininginterval constraint solving with a CDCL algorithm and [21] is adecision procedure for nonlinear real arithmetic combining inter-val constraint solving with an SMT solver for linear arithmetic.[44]formalizes the floating-points for the SMT-LIB format. While theseapproach can check ranges on numeric variables, they do not handleroundoff errors or other uncertainties and cannot compute specifi-cations automatically.Our techniques rely on the performance of Z3. We hope thatan integration of the recent new improved solver for nonlinear arithmetic [30] will make many more verification problems feasiblewith our techniques. An alternative to this approach is using linearapproximations to solve polynomial constraints [7]. We believe thatsuch advances are largely orthogonal to our use of range arithmeticand complement each other.
Testing
Symbolic execution is a well-known technique for gener-ating test inputs. [6] use a combination of meta-heuristic search andinterval constraint solving to solve the floating-point constraintsthat arise, whereas [33] combine random search and evolutionarytechniques. [47] test numerical code for precision by perturbinglow-order bits of values and rewriting expressions. The idea is toexagerate initial errors and thus make imprecisions more visible.Probabilistic arithmetic [45] is a similar approach but it does theperturbation by using different rounding modes. [4] also proposea testing produce to detect accuracy problems by instrumentringcode to perform a higher-precision computation side by side withthe regular computations. While these approaches are sound withrespect to floating-point arithmetic, they only generate or can checkindividual inputs and are thus not able to verify or compute outputranges or their roundoff errors.
Robustness analysis [27] combines abstract interpretation withmodel checking to check programs for stability by tracking the evo-lution of the width of the interval representing a single input. [36]use concolic execution to find inputs which, given maximum devi-ations on inputs, maximize the deviation on the outputs. These twoworks however, use a testing approach and cannot provide soundguarantees. [9] presents a framework for continuity analysis of pro-grams along the mathematical (cid:15) − δ definition of continuity and [10]builds on this work and presents a sound robustness analysis. Thisframework provides a syntactic proof of robustness for programsover reals and thus does not consider floating-points. Our approachdescribes a quantitative measure of robustness for nonlinear pro-grams with floating-point numbers and other uncertainties, and webelieve that it can complement the cited framework.
10. Conclusion
We have presented a programming model for numerical programsthat decouples the mathematical problem description from its re-alization in finite precision. The model uses a
Real data type thatcorresponds to mathematical real numbers. The developer specifiesthe program using reals and indicates the target precision; the com-piler chooses a floating point representation while checking that thedesired precision targets are met. We have described the soundnesscriteria by translating programs with precision requirements intoverification conditions over mathematical reals. The resulting veri-fication conditions, while a natural description of the problem beingsolved, are difficult to solve using a state-of-the art SMT solver Z3.We therefore developed an algorithm that combines SMT solvingwith range computation. Our notion of soundness incorporates fullinput/output behavior of functions, taking into account that, due toconditionals, small differences in values can lead to different pathsbeing taken in the program. For such cases our approach estimatesa sound upper bound on the total error of the computation.We have evaluated our techniques on a number of benchmarksfrom the literature, including benchmarks from physics, biology,chemistry, and control systems. We have found that invocation ofSMT solver alone is not sufficient to handle these benchmarks dueto scalability issues, whereas the use of range arithmetic by itselfis not precise enough. By combining these two techniques we wereable to show that a floating point version of the code conforms tothe real-valued version with reasonable precision requirements.We believe that our results indicate that it is reasonable to in-troduce
Real s as a data type, following a list of previously intro-duced mathematical abstractions in programming languages, such s unbounded integers, rationals, and algebraic data types. The fea-sibility of verified compilation of our benchmarks suggests that itis realistic to decouple the verification of executable mathematicalmodels over reals from their sound compilation. We therefore ex-pect that this methodology will help advance rigorous formal veri-fication of numerical software and enable us to focus more on high-level correctness properties as opposed to run-time errors alone. References [1] A. Anta and P. Tabuada. To Sample or not to Sample: Self-TriggeredControl for Nonlinear Systems.
IEEE Transactions on AutomaticControl , 55(9), 2010.[2] A. Ayad and C. March´e. Multi-prover verification of floating-pointprograms. In
IJCAR , 2010.[3] D. H. Bailey, Y. Hida, X. S. Li, and B. Thompson. C++/Fortran-90 double-double and quad-double package. http://crd-legacy.lbl.gov/˜dhbailey/mpdist/, 2013.[4] F. Benz, A. Hildebrandt, and S. Hack. A dynamic program analysis tofind floating-point accuracy problems. In
PLDI , 2012.[5] B. Blanchet, P. Cousot, R. Cousot, J. Feret, L. Mauborgne, A. Min´e,D. Monniaux, and X. Rival. A static analyzer for large safety-criticalsoftware. In
PLDI , pages 196–207, 2003.[6] M. Borges, M. d’Amorim, S. Anand, D. Bushnell, and C. S. Pasareanu.Symbolic Execution with Interval Solving and Meta-heuristic Search.In
ICST , 2012.[7] C. Borralleras, S. Lucas, A. Oliveras, E. Rodr´ıguez-Carbonell, andA. Rubio. Sat modulo linear arithmetic for solving polynomial con-straints.
J. Automated Reasoning , 48(1), 2012.[8] A. Brillout, D. Kroening, and T. Wahl. Mixed abstractions for floating-point arithmetic. In
FMCAD , pages 69–76, 2009.[9] S. Chaudhuri, S. Gulwani, and R. Lublinerman. Continuity analysisof programs. In
POPL , 2010.[10] S. Chaudhuri, S. Gulwani, R. Lublinerman, and S. Navidpour. ProvingPrograms Robust. In
ESEC/FSE , 2011.[11] L. Chen, A. Min´e, J. Wang, and P. Cousot. Interval Polyhedra: AnAbstract Domain to Infer Interval Linear Relationships. In
SAS , 2009.[12] E. Darulova and V. Kuncak. Trustworthy numerical computation inScala. In
OOPSLA , 2011.[13] E. Darulova, V. Kuncak, R. Majumdar, and I. Saha. Synthesis of fixed-point arithmetic code. In
EMSOFT , 2013.[14] F. de Dinechin, C. Lauter, and G. Melquiond. Certifying the Floating-Point Implementation of an Elementary Function Using Gappa.
IEEETrans. Comput. , 2011.[15] L. H. de Figueiredo and J. Stolfi.
Self-Validated Numerical Methodsand Applications . IMPA/CNPq, Brazil, 1997.[16] L. De Moura and N. Bjørner. Z3: an efficient SMT solver. In
TACAS ,2008.[17] D. Delmas, E. Goubault, S. Putot, J. Souyris, K. Tekkal, andF. Vedrine. Towards an industrial use of FLUCTUAT on safety-criticalavionics software. In
FMICS , 2009.[18] V. D’Silva, L. Haller, D. Kroening, and M. Tautschnig. NumericBounds Analysis with Conflict-driven Learning. In
TACAS , 2012.[19] J. A. Duracz and M. Konecny. Polynomial function enclosures andfloating point software verification.
Constraints in Formal Verifica-tion/IJCAR , 2008.[20] J. Feret. Static Analysis of Digital Filters. In the 13th EuropeanSymposium on Programming - ESOP , 2004.[21] S. Gao, M. Ganai, F. Ivancic, A. Gupta, S. Sankaranarayanan, andE. Clarke. Integrating ICP and LRA solvers for deciding nonlinearreal arithmetic problems. In
FMCAD , 2010.[22] K. Ghorbal, E. Goubault, and S. Putot. A Logical Product Approachto Zonotope Intersection. In
CAV , 2010.[23] D. Goldberg. What every computer scientist should know aboutfloating-point arithmetic.
ACM Comput. Surv. , 23(1), 1991. [24] E. Goubault, S. Putot, and F. V´edrine. Modular static analysis withzonotopes. In
SAS , 2012.[25] L. Haller, A. Griggio, M. Brain, and D. Kroening. Deciding floating-point logic with systematic abstraction. In
FMCAD , 2012.[26] J. Harrison. Floating-Point Verification using Theorem Proving. In
Formal Methods for Hardware Verification, 6th International Schoolon Formal Methods for the Design of Computer, Communication, andSoftware Systems, SFM , 2006.[27] F. Ivancic, M. Ganai, S. Sankaranarayanan, and A. Gupta. Numericalstability analysis of floating-point computations using software modelchecking. In
MEMOCODE , 2010.[28] B. Jeannet and A. Min´e. Apron: A Library of Numerical AbstractDomains for Static Analysis. In
CAV , 2009.[29] J. Jiang, W. Luk, and D. Rueckert. FPGA-Based Computation of Free-Form Deformations. In
Field-Programmable Logic and Applications ,2003.[30] D. Jovanovi´c and L. de Moura. Solving Non-linear Arithmetic. In
IJCAR 2012 , 2012.[31] W. Kahan. Miscalculating Area and Angles of a Needle-like Triangle.Technical report, University of California Berkeley, 2000.[32] V. Kuznetsov, J. Kinder, S. Bucur, and G. Candea. Efficient StateMerging in Symbolic Execution. In
PLDI , 2012.[33] K. Lakhotia, N. Tillmann, M. Harman, and J. de Halleux. FloPSy -Search-Based Floating Point Constraint Solving for Symbolic Execu-tion. In
Testing Software and Systems . Springer Berlin / Heidelberg,2010.[34] X. Leroy. Verified squared: does critical software deserve verifiedtools? In
POPL , 2011.[35] M. D. Linderman, M. Ho, D. L. Dill, T. H. Meng, and G. P. Nolan. To-wards program optimization through automated analysis of numericalprecision. In
CGO , 2010.[36] R. Majumdar, I. Saha, and Z. Wang. Systematic Testing for ControlApplications. In
MEMOCODE , 2010.[37] K. Makino and M. Berz. Taylor Models and Other Validated Func-tional Inclusion Methods.
International Journal of Pure and AppliedMathematics , 4, 2003.[38] A. Min´e. Relational Abstract Domains for the Detection of Floating-Point Run-Time Errors. In
ESOP , 2004.[39] R. Moore.
Interval Analysis . Prentice-Hall, 1966.[40] J. D. Murray.
Mathematical Biology, I. An Introduction . Springer,2002.[41] M. Odersky, L. Spoon, and B. Venners.
Programming in Scala: AComprehensive Step-by-step Guide . Artima Incorporation, 2008.[42] O. Ponsini, C. Michel, and M. Rueher. Refining Abstract Interpreta-tion Based Value Analysis with Constraint Programming Techniques.In CP , 2012.[43] A. Quarteroni, F. Saleri, and P. Gervasio. Scientific Computing withMATLAB and Octave . Springer, 3rd edition, 2010.[44] P. R¨ummer and T. Wahl. An SMT-LIB Theory of Binary Floating-Point Arithmetic. In
Informal proceedings of 8th International Work-shop on Satisfiability Modulo Theories (SMT) at FLoC , 2010.[45] N. Scott, F. J´ez´equel, C. Denis, and J.-M. Chesneaux. Numerical‘health check’ for scientific codes: the CADNA approach.
ComputerPhysics Communications , 2007.[46] I. C. Society. IEEE Standard for Floating-Point Arithmetic.
IEEE Std754-2008 , 2008.[47] E. Tang, E. Barr, X. Li, and Z. Su. Perturbing numerical calculationsfor statistical analysis of floating-point program (in)stability. In
ISSTA ,2010.[48] E. M. Westbrook and S. Chaudhuri. A Semantics for ApproximateProgram Transformations.
CoRR , 2013.[49] C. Woodford and C. Phillips.
Numerical Methods with Worked Exam-ples , volume 2nd. Springer, 2012.12