1 Introduction↩︎

Many industrial applications in computer science and engineering reduce to deciding formulas over the real numbers that involve transcendental functions such as the sine and the exponential [1], [2]. Yet, deciding such formulas is notoriously difficult, and in many cases outright impossible:

• The problem is already undecidable for univariate formulas with addition, multiplication, and the sine function [3]; and only known to be decidable subject to Schanuel’s conjecture for the exponential [4].

• Even representing and checking solutions is sometimes not possible, in particular for formulas whose only solutions are transcendental numbers, i.e., numbers that are not roots of any polynomial with integer coefficients.

These facts motivate the design of principled but incomplete methods to tackle problems involving transcendental functions on a best-effort basis. In this paper, we present such a method.

Specifically, we address the Satisfiability Modulo Theories (SMT) problem for Transcendental Real Arithmetics (TRAs), that is, deciding the satisfiability of quantifier-free formulas from first-order theories extending Non-linear Real Arithmetic (NRA) with transcendental functions. Working in the Model Constructive Satisfiability (MCSAT) framework (which generalizes ideas from CDCL to the theory level), we design a procedure applicable to any such theory whose transcendental functions are total and interval-computable.

At the high level, our procedure abstracts an input TRA formula to NRA, and uses NRA lemmas to incrementally refine the abstraction. A Transcendental Real Arithmetic plugin (TRA plugin) acts as an intermediary between the core solver of MCSAT and a (already existing) plugin for NRA. The role of the TRA plugin is that of ensuring that the trail of MCSAT is “consistent” in the following sense: for every literal fully assigned in the trail, its truth value in the NRA abstraction matches the truth value in the TRA concretization. When the TRA plugin detects an inconsistency, it triggers a conflict, producing a valid NRA lemma that is conjoined to the NRA abstraction, effectively refining it.

In TRAs, evaluating the truth value of a fully-assigned literal is in general an undecidable problem. We adopt the practical strategy of relaxing the consistency check of the TRA plugin up to a tolerance \(\delta > 0\). This means that, in case it is not possible to disprove consistency using over-approximations of the transcendental functions that depend on \(\delta\), the TRA plugin considers the check passed, and marks the literal as \(\delta\)-consistent. The benefit of this approach is that \(\delta\)-consistency is decidable for interval-computable functions, and hence the TRA plugin does not get stuck during its consistency check.

The procedure can operate in two modes:

1. With tolerance \(\delta\) fixed throughout the procedure and provided by the user. In this mode, the procedure may return SAT, UNSAT or \(\delta\)-SAT. The first two cases carry their usual meaning. In the case of \(\delta\)-SAT, the procedure has found a solution to the NRA abstraction for which all literals were proven consistent or \(\delta\)-consistent by the TRA plugin, with at least one being \(\delta\)-consistent.

2. With tolerance not fixed, and the procedure may only return SAT or UNSAT. A tolerance \(\delta\) is still used internally, but rather than answering \(\delta\)-SAT, the procedure decreases \(\delta\) and backtracks to the empty trail, while preserving all NRA lemmas produced so far.

While our procedure works for any TRA in which the transcendental functions are interval-computable, we further devote special attention to the case of the sine and exponential functions (plus the constant \(\pi\)), providing more precise NRA lemmas based on Taylor and Padé approximations.

For \(\sin\), \(\exp\) and \(\pi\), we implemented a prototype¹ of our procedure within the Yices2 SMT solver [5]. We compared the prototype to MathSAT [6], cvc5 [7] and dReal [8] on the 2 512 instances from [6]. Our experimental evaluation shows that our tool outperforms other tools on both SAT and UNSAT instances (and, when \(\delta\) is fixed, also on \(\delta\)-SAT instances).

2 Preliminaries↩︎

In this section we introduce the necessary background material on arithmetic theories over the reals, and provide a high-level overview of the MCSAT framework. We assume familiarity with propositional satisfiability and the Conflict-Driven Clause Learning (CDCL) algorithm (a brief overview of CDCL is given in the appendix).

We recall that the Satisfiability Modulo Theory (SMT) problem asks for the satisfiability of a quantifier-free formula in conjunctive normal form (CNF) from a given first-order theory.

2.1 Non-linear Real Arithmetic (NRA)↩︎

NRA is formally defined as the first-order theory of the structure \({(\mathbb{R}, 0, 1, +, \cdot, <, =)}\). This means that, in NRA, variables range over real numbers, and wlog. atomic formulas are strict inequalities \(p(\boldsymbol{x}) < 0\) and equalities \(p(\boldsymbol{x}) = 0\), where the term \(p(\boldsymbol{x})\) is a polynomial with integer coefficients. We write \(p \leq 0\) and \(p \neq 0\) as shortcuts for \(\lnot (-p < 0)\) and \(\lnot (p = 0)\), respectively. Throughout the paper, we only consider quantifier-free formulas, given in CNF. The satisfiability problem for NRA is decidable, and in practice most decision procedures for this theory are based on Collin’s Cylindrical Algebraic Decomposition (CAD) [9].

A fundamental property of NRA is that if a formula has a solution, then it has an algebraic solution [10], i.e., a solution assigning algebraic numbers to all variables. A real number is algebraic if it is the root of a univariate polynomial with integer coefficients; else, it is transcendental. For instance, \(\sqrt{2}\) is algebraic since it is a solution of the equation \(x^2 - 2 = 0\), while \(\pi\) is transcendental. We write \(\overline{\mathbb{Q}}\) to denote the set of algebraic numbers. An algebraic number \(\alpha\) can be effectively represented with a triple \((p,\ell,u)\) where \(p\) is a univariate polynomial of which \(\alpha\) is a root, and \(\ell\) and \(u\) are rational numbers such that \(\alpha\) is the only root of \(p\) in the interval \([\ell,u]\). All standard decision procedures for NRA restrict the search of solution to algebraic solutions.

2.2 Real Arithmetics with Transcendental Functions↩︎

A real-valued function \(f \colon \mathbb{R}^n \to \mathbb{R}^m\) is called transcendental if its graph cannot be characterized by any formula in NRA; that is, there is no formula \(\varphi(\boldsymbol{x}, \boldsymbol{y})\) in NRA such that, for every \(\boldsymbol{r} \in \mathbb{R}^n\) and \(\boldsymbol{s} \in \mathbb{R}^m\), \(\varphi(\boldsymbol{r}, \boldsymbol{s})\) holds if and only if \(f(\boldsymbol{r}) = \boldsymbol{s}\). A Transcendental Real Arithmetic (TRA) is understood as any first-order theory of a structure that expands NRA with (computable) transcendental functions. Classical examples of TRAs are given by the first-order theories of the structures \({(\mathbb{R}, 0,1, +, \cdot, \sin, <,=)}\) and \({(\mathbb{R}, 0,1, +, \cdot, \exp, <,=)}\), i.e., the expansions of NRA by the sine function and the exponential function \(\exp(x) \mathrel{\vcenter{:}}= e^x\), respectively.

TRAs are very expressive. For example, even the univariate fragment of \({(\mathbb{R}, 0,1, +, \cdot, \sin, <,=)}\) is undecidable [3], while \({(\mathbb{R}, 0,1, +, \cdot, \exp, <,=)}\) is only known to be decidable subject to Schanuel’s conjecture [4]. Even worse, solutions are often unverifiable [11]: there is no general algorithm to decide if a given variable assignment satisfies a TRA formula. This stems from the nature of transcendental functions, whose outputs can typically only be computed to arbitrary precision, whereas the formula may be satisfied only at their exact value. For these reasons, except for procedures targeting very restricted theories (see, e.g., [12], [13]), all decision procedures for TRAs are incomplete. The procedure we describe in the next section is no exception. More precisely, our procedure is designed to solve extensions of NRA by interval-computable total functions, defined next.

Three observations are in order. First, the composition of two interval computable functions is interval computable. Second, addition, multiplication, and many transcendental functions of practical interest are interval computable, including the sine and exponential functions. For these natural functions, the algorithms \(I_f\) are effective and already available in several libraries; in particular, our tool uses the ARB library [14] for arbitrary-precision interval computation of the sine and exponential functions. Lastly, any computable constant \(c\) (e.g., \(\pi\) or \(e\)) can be treated as a unary interval computable function: given a box of some width \(\epsilon\), \(I_c\) computes an interval of width \(\epsilon\) around \(c\). Throughout the paper, we tacitly see transcendental constants as unary functions applied to the constant \(1\).

\(\delta\)-satisfiability↩︎

Because of undecidability and unverifiability, symbolic computation often offers very limited support in deciding TRAs, and must be complemented by numerical methods. A principled framework for combining these methods with formal soundness guarantees is given by the notion of \(\delta\)-satisfiability procedure from [15], which we now recall (in a slightly extended form that still allows SAT results).

Let \(\delta > 0\) be a rational (the tolerance). Let \(\varphi\) be a quantifier-free CNF formula from some expansion of NRA, with atomic formulas of the form \(\tau < 0\) or \(\tau = 0\). The \(\delta\)-weakening of \(\varphi\), denoted \(\varphi^\delta\), is the formula obtained from \(\varphi\) by replacing the literals in each clause following the rules: \(\tau < 0 \, \rightarrow\, \tau < \delta\), \(\lnot (\tau < 0) \, \rightarrow\, \lnot (\tau < -\delta)\), \(\tau = 0 \, \rightarrow\, -\delta \leq \tau \leq \delta\), and \(\lnot (\tau = 0) \, \rightarrow\, \top\). A formula \(\varphi\) is said to be \(\delta\)-satisfiable if and only if its \(\delta\)-weakening \(\varphi^\delta\) is satisfiable.

2.3 Background on MCSAT↩︎

We briefly introduce MCSAT, referring the reader to [16] and [17] for more detailed presentations.

The framework consists of a core solver that orchestrates a set of theory plugins. Each theory plugin implements a procedure for a specific logical theory. The control flow of the core solver, illustrated in 1, closely follows that of CDCL. As in CDCL, the central data structure is the trail, which is a list of four different types of assignments:

• Boolean decisions. Assignments \(\ell \mapsto b\) where \(\ell\) is a literal and \(b\) is a Boolean. These are added to the trail during the step.

• Semantic decisions. Assignments \(x \mapsto v\) where \(x\) is a first-order variable over some sort \(S\) (e.g., \(\mathbb{R}\)), and \(v\) is an element of \(S\). These decisions are also added during the step: calls the theory plugin associated to \(S\) (e.g., the NRA plugin), asking a value for \(x\).

• Boolean propagations. Assignments \(\ell \mapsto b\) where \(\ell\) is a literal and \(b\) is a Boolean, that are added during the step and derived by unit propagation.

• Semantic propagations. Assignments \(\ell \mapsto b\) added during the step by a theory plugin. Specifically, when the literal \(\ell\) contains variables already assigned in the trail, the theory plugin evaluates its truth value \(b\).

At each  step, exactly one Boolean or semantic decision is added to the trail; this decision can only feature variables or literals that are not yet in the trail. The other steps of MCSAT are described below.

Step ↩︎

This step saturates the trail with all available Boolean and Semantic propagations. Whenever the trail is updated, notifies all plugins of the update. The plugins use this information to refine a feasibility set for each variable of their competence. For example, if \(x\mapsto 2\) is in the trail, and notifies the NRA plugin that \((x^2-y \leq 0)\mapsto \top\) has been added, the plugin updates the feasibility set of \(y\) to exclude the interval \((-\infty,4]\). The feasibility sets are used to detect conflicts, as well as to perform semantic decisions.

Step ↩︎

A trail is propositionally consistent if it does not contain both \(\ell \mapsto \top\) and \(\ell \mapsto \bot\), for any literal \(\ell\). As semantics decisions do not involve variables in the trail, a propositionally consistent trail \(T\) can be viewed as a map (instead of a list), with \(T(t) = v\) if and only if \(t \mapsto v\) occurs in \(T\). We write \(\mathrm{dom}(T)\) for the domain of this map. A literal is said to be fully-assigned in \(T\) if \(\mathrm{vars}(\ell) \subseteq \mathrm{dom}(T)\), where \(\mathrm{vars}(\ell)\) is the set of first-order variables occurring in \(\ell\). For a fully-assigned literal \(\ell\), we write \(T \models \ell\) when the assignments of \(\mathrm{vars}(\ell)\) in \(T\) satisfy the literal. The notation \(T \models \psi\) extends to Boolean combinations \(\psi\) of fully-assigned literals.

After , checks the of the trail \(T\) by verifying that it is propositionally consistent, and by asking the plugins to check (1) that every unassigned variable has a non-empty feasibility set, and (2) that the value of fully-assigned literals in the trail is consistent with their satisfaction. Specifically, for each fully-assigned literal \(\ell\) with \(\ell \mapsto b\) in \(T\), it must hold that \(T \models \ell\) if and only if \(b = \top\).

Step ↩︎

If the check passes, verifies  of the trail, i.e., each clause in the input CNF formula contains a fully-assigned literal \(\ell\) such that \(\ell \mapsto \top\) or \((\lnot \ell) \mapsto \bot\) is present in the trail (where \(\lnot \ell\) is understood as \(\ell'\) if \(\ell = \lnot \ell'\)). If this is the case, returns  since the trail satisfies all clauses in the input formula; otherwise it moves to the  step.

Step ↩︎

When the  check fails, if the trail does not contain a decision, then the formula is . Otherwise, must backtrack and learn a new clause \(\psi\) that is implied by the input formula \(\varphi\) and explains the conflict. If the conflict stems from Boolean reasoning (e.g., it is due to the Boolean propagations performed after a Boolean decision), behaves as CDCL. Otherwise, it invokes the theory plugins. The plugins analyse the trail, and generate a suitable clause \(\psi\) explaining the conflict, along with the semantic or Boolean decision causing it. conjoins \(\psi\) to \(\varphi\), and then backtracks to before the problematic decision. When backtracking on a semantic decision \(x \mapsto v\), does not necessarily decide \(x\) first upon reaching the  step.

3 TRAs in MCSAT↩︎

We describe our MCSAT-based procedure for the first-order theory of a structure \((\mathbb{R}, 0, 1, +, \cdot, f_1,\dots,f_k, <, =)\) extending NRA with \(k\) interval computable transcendental functions \({f_j \colon \mathbb{R}^{n_j} \to \mathbb{R}}\). The only assumption we make is to have access, for each \(f_j\), to an algorithm \(I_{f_j}\) as described in 2. This general procedure is later tailored to the extension of NRA by the exponential and sine functions in 4.

3.1 High-level view of the procedure↩︎

As illustrated in 2, our procedure starts with a preprocessing phase that constructs from the input formula \(\varphi\) an abstraction \(\varphi^\alpha\) in NRA (defined in 3.2). The MCSAT framework is then executed on \(\varphi^\alpha\), with (the core solver of) MCSAT querying the NRA plugin. The main component of our procedure is given by a Transcendental Real Arithmetic plugin (TRA plugin) that acts as an intermediary between the core solver and its NRA plugin. The TRA plugin has access to both the formula \(\varphi\) and the abstraction \(\varphi^\alpha\) (as well as information to connect the two formulas, discussed later).

In a nutshell, the TRA plugin is limited to observing the communication between MCSAT and the NRA plugin, except for when the consistency of the trail must be checked or conflicts must be analysed. For instance, if the NRA plugin determines that the trail is consistent, the TRA plugin runs its own consistency check, but over constraints from \(\varphi\) rather than \(\varphi^{\alpha}\). If the check fails, the TRA plugin proposes clauses that are implied by \(\varphi\) but not by \(\varphi^\alpha\), effectively refining \(\varphi^\alpha\) enough to trigger a conflict in the NRA plugin.

Apart from assuming that MCSAT performs the loop described in the previous section, our procedure relies on the following assumptions on the NRA plugin. The existence of a plugin satisfying them is guaranteed by the decidability of NRA, and the NRA plugin of Yices2, which implements the CAD-based procedure from [18], is a concrete instance.

At its core, our procedure respects the specification of \(\delta\)-satisfiability procedure (3), but can also run in a standard SAT/UNSAT mode (as explained in 3.7).

In practical terms, literals in the trail can take three values: \(\top\), \(\bot\) and \(\top^\delta\). Except for the TRA plugin, all other theory plugins treat \(\top^\delta\) as \(\top\). The TRA plugin can weaken literals in the trail, replacing \(\top\) with \(\top^\delta\) while performing its consistency check (explained in 3.4). Roughly speaking, an assignment \({\ell \mapsto \top^\delta}\) in the trail should be understood as “\(\ell\) is a fully-assigned literal in the trail that is true with respect to the abstraction, but whose satisfaction in the original formula is only guaranteed for its \(\delta\)-weakening”.

In the core solver of MCSAT, the only difference concerns the check. When this check passes (seeing \(\top^\delta\) as \(\top\)), but a clause in the formula is only true because of a literal with value \(\top^\delta\), MCSAT returns \(\delta\)-SAT instead of SAT.

3.2 Preprocessing: the first NRA abstraction↩︎

We discuss how the first abstraction \(\varphi^\alpha\) is built. We tag each clause in the input formula \(\varphi\) as original, so that the TRA plugin can distinguish literals from these clauses among those in \(\varphi^\alpha\). These are referred to as original clauses and literals. Afterwards, we iteratively consider all terms \(f(t_1,\dots,t_r)\) in \(\varphi\) starting from the outermost ones, where \(f\) is a transcendental function, and perform the rewriting step \[\varphi\longrightarrow \varphi[f_{(t_1,\dots,t_r)}^{\text{out}}\,/\,f(t_1,\dots,t_r)] \land \textstyle\bigwedge_{i=1}^r f_{t_i}^{\text{in}} = t_i\,,\] where each \(f_{t_i}^{\text{in}}\) and \(f_{(t_1,\dots,t_r)}^{\text{out}}\) are fresh real-valued variables, and \(\varphi[t'\,/\,t'']\) stands for the formula obtained from \(\varphi\) by replacing every occurrence of \(t''\) with \(t'\). This rewriting step is applied until no transcendental function is left in the formula, and the resulting (NRA) formula is \(\varphi^\alpha\).

The TRA plugin is initialised by providing the formulas \(\varphi\) and \(\varphi^\alpha\), as well as maps to reconstruct terms \(f(t_1,\dots,t_r)\) from the variables \(f_{t_i}^{\text{in}}\) and \(f_{(t_1,\dots,t_r)}^{\text{out}}\). Given a literal \(\ell\) (resp. term \(\tau\)) occurring in \(\varphi^\alpha\), we write \(\ell^\gamma\) (resp. \(\tau^\gamma\)) for the literal (resp. term) obtained by replacing every \(f_{(t_1,\dots,t_r)}^{\text{out}}\) with \(f(t_1,\dots,t_r)\) and every \(f_{t_i}^{\text{in}}\) with \(t_i\). When the literal \(\ell\) is original, we write \(\mathrm{allvars}(\ell)\) for all the variables in \(\ell^\gamma\), plus all variables added when constructing \(\ell\) from \(\ell^\gamma\).

The following lemma, whose proof is immediate, states fundamental properties of \(\varphi^\alpha\). The TRA plugin preserves these properties as invariants throughout the procedure, while the NRA abstraction \(\varphi^\alpha\) is refined.

After constructing \(\varphi^\alpha\) and initializing the TRA plugin, MCSAT runs on \(\varphi^\alpha\) as described in 2.3.

We now detail the role of the TRA plugin at each step of the MCSAT framework.

3.3 Steps: and↩︎

In these two steps, the role of the TRA plugin is limited to forwarding data from MCSAT to the NRA plugin, and returning the results computed by the latter back to MCSAT. Note that, by 1, this implies that all semantic decisions assign algebraic values to variables.

3.4 Step:↩︎

Recall that during this step, the theory plugins verify (1) that the feasibility sets they maintain for the unassigned variables are non-empty, and (2) that the values of fully-assigned literals in the trail are consistent with their satisfaction.

The TRA plugin begins this step by forwarding the request to the NRA plugin, which perform its own consistency check. Should the trail be found inconsistent at this stage, the NRA plugin raises a conflict that the TRA plugin propagates back to MCSAT. Otherwise, the TRA plugin proceeds with an additional consistency check.

The additional check does not concern feasibility sets: the TRA plugin relies on the feasibility sets of the NRA plugin, which, by definition of \(\varphi^\alpha\), are overapproximations of the possible values that variables can take. While it is sometimes possible to refine these approximations beyond those of the NRA plugin (see, e.g., [12]), the undecidability status of even the univariate extension of NRA by the sine function [3] (which is interval computable) prevents the determination of exact feasibility sets in general.

The consistency check of the TRA plugin is therefore focused solely on the second of the points above: verifying that the values of fully-assigned literals in the trail are consistent with their satisfaction. Because this verification problem is undecidable (as explained in 2.2), we perform it only “up to \(\delta\)”. We call this the TRA \(\delta\)-consistency check:

3, which we discuss in a moment, implements the TRA \(\delta\)-consistency check. Following 4, we see that there are two cases in which the TRA plugin will not perform this check on a literal \(\ell\) occurring in the trail:

• When \(\ell\) is not original. The decision is pragmatic: clauses will only be added to improve the NRA abstraction and therefore soundness is still guaranteed.

• If some variable in the set \(\mathrm{allvars}(\ell)\) is unassigned in the trail. In this case, the TRA plugin places \(\ell\) in a watchlist and defers its consistency check until all variables in \(\mathrm{allvars}(\ell)\) are assigned (note: this must occur for the check to pass).

Consider thus an assignment \(\ell \mapsto b'\) in the trail \(T\), such that \(\ell\) is original, all variables in \(\mathrm{allvars}(\ell)\) are assigned in \(T\), and that passes the consistency check of the NRA plugin, i.e., we have \(T \models \ell\) if and only if \(b' = \top\). 3 works as follows. First, if \(\ell^\gamma\) contains no transcendental function, it returns \(b'\). Indeed, in this case \(\ell = \ell^\gamma\), and the NRA plugin has already verified \(\ell \mapsto b'\). Otherwise, let \(\tau\) be the term occurring in \(\ell\) (as in line [pseu:cons:tau]). Since interval computable functions are closed under composition, \(\tau^\gamma\) is interval computable. The algorithm computes an algebraic box \([i,j] \in \mathcal{B}_1(\overline{\mathbb{Q}})\) of width at most \(\delta\), such that \(T \models i \leq \tau^\gamma \leq j\). This interval is computed bottom-up, computing boxes \(I(\rho)\) for every subterm \(\rho\) of \(\tau^\gamma\), starting from the innermost subterms. Interval computability ensures that by choosing sufficiently small intervals around the values of variables in \(\tau^\gamma\) and the constants \(0\) and \(1\), the final interval \(I(\tau^\gamma) = [i,j]\) will have width less than \(\delta\). 3 terminates by performing interval analysis, following the table in line [pseu:cons:table]. Note that \(\top^\delta\) is only returned in cases where \(i\) and \(j\) are neither both positive, both negative, nor both zero.

Returning to the assignment \(\ell \mapsto b'\), let \(b\) denote the output of 3 on input \(\ell\) and \(T\). After performing the TRA \(\delta\)-consistency check, the TRA plugin completes its consistency step on \(\ell\) as follows:

3.5 Step:↩︎

Suppose the TRA plugin raises a conflict, causing MCSAT to enter the  step (when the trail to contain a decision) and query the TRA plugin in return. If the conflict originates from the NRA plugin, the TRA plugin lets the NRA plugin resolve it, again acting as an intermediary. Otherwise, the conflict arose following the TRA \(\delta\)-consistency check. In particular, there is an assignment \({\ell \mapsto b'}\) in the trail \(T\) with the following properties:

• \(\ell\) is original, and \(\mathrm{allvars}(\ell) \subseteq \mathrm{dom}(T)\).

• \(T \models \ell\), but \(T \not\models \ell^\gamma\). More precisely, 3 returned \(b \in \{\bot,\top\}\) on input \(\ell\) and \(T\), with \(b \neq b'\).

The TRA plugin uses \(\ell\) in order to refine \(\varphi^\alpha\) so that it continues to satisfy the properties in [lemma:abstractionsound] while also triggering a conflict in the NRA plugin.

The refinement consists in conjoining to \(\varphi^\alpha\) the formula \(\psi\) computed by 4. In a nutshell, this algorithm considers all subterms \(f(t_1,\dots,t_r)\) of \(\ell^\gamma\) where \(f\) is a transcendental function, and adds (line [refinement:psi]) entailments stating that if the value assigned to \(f_{t_j}^{\text{in}}\) lies in \(I(t_j)\) for every \(j \in [1..r]\), then the value assigned to \(f_{(t_1,\dots,t_r)}^{\text{out}}\) must lie in \(I(f(t_1,\dots,t_r))\), where \(I(\rho)\) denotes the interval computed for subterm \(\rho\) during the execution of 3 on \(\ell\) and \(T\).

After computing \(\psi\), the TRA plugin identifies a conjunct \((\bigwedge_{\rho \in \mathcal{T}} i_{\rho} \leq f_{\rho}^{\text{in}} \leq j_{\rho}) \Rightarrow i_\eta \leq f_{\boldsymbol{t}}^\text{out} \leq j_\eta\) of this formula that is not satisfied by \(T\). This reduces to evaluating inequalities between algebraic numbers, which can be done in polynomial time, as it corresponds to a query in the existential theory of the reals with a fixed number of variables [19].

Let \(x\) be the variable in the identified conjunct that was decided last in the trail \(T\). The TRA plugin returns to MCSAT the formula \(\psi\), to be conjoined to \(\varphi^\alpha\) (refining the abstraction), along with the variable \(x\) to backtrack. Before doing so, however, the TRA plugin performs some bookkeeping to ensure that the watchlist of literals waiting for the \(\delta\)-consistency check, and their values in the trail, will remain consistent with the trail after backtracking. More precisely, the plugin will re-add to the watchlist every literal \(\ell\) on the trail that, after backtracking \(x\), will have at least one variable in \(\mathrm{allvars}(\ell)\) become unassigned. Moreover, for any such literal, any trail assignment \(\ell \mapsto \top^\delta\) is set back to \(\ell \mapsto \top\).

3.6 Soundness of the procedure↩︎

We show that ours is a \(\delta\)-satisfiability procedure:

3.7 More on \(\delta\)-satisfiability↩︎

The notion of \(\delta\)-satisfiability procedure was introduced in [15] to provide theoretical grounding for the use of numerical methods in decision procedures. Different procedures may exploit the flexibility this notion offers in different ways. In our setting, we only employ \(\delta\)-weakening when checking the consistency of literals involving transcendental functions (3). This means, for example, that we never relax literals that are in NRA.

Note that even the definition of TRA\(\delta\)-consistency check (Def. 4) allows some flexibility in the output. For example, when the interval \([i,j]\) computed for a literal \(\tau < 0\) is included in \((0,\delta)\), a \(\delta\)-consistency check may return either \(\bot\) or \(\top^\delta\). However, the table that we use in 3, which returns \(\bot\) in this case, is the most precise possible. This guarantees that, within a single consistency check, an inconsistency output (that leans toward proving UNSAT) is always preferred over a \(\delta\)-consistency one (which leans toward proving \(\delta\)-SAT).

Our \(\delta\)-satisfiability procedure naturally extends to a standard SAT/UNSAT procedure. It suffices to iteratively run the \(\delta\)-satisfiability procedure starting with a large value of \(\delta\) (e.g., \(\delta = 1\)), halving \(\delta\) and restarting when finding \(\delta\)-SAT. Crucially, restarting can be done while preserving the current abstraction \(\varphi^\alpha\), since the learned lemmas are independent of \(\delta\). Decision heuristics (based, e.g., on caching [20]), can further improve efficiency during restarts. By never returning \(\delta\)-SAT, [theorem:soundness-delta] guarantees that this SAT/UNSAT mode is sound with respect to the standard notion of satisfiability.

4 Specialized treatment of the sine and exponential functions↩︎

The previous section outlined a general procedure for TRAs based solely on interval computability. Here, we show how additional assumptions can be used to refine this base procedure. We focus on the sine and exponential functions, though the approach extends in a modular way to other analytic functions. Formally, we now consider the structure \((\mathbb{R}, 0, 1, \pi, e, +, \cdot, \sin, \exp, f_1,\dots,f_k, <, =)\), where \(f_1,\dots,f_k\) are interval computable functions. (We recall that we see the constants \(\pi\) and \(e\) as interval-computable unary constant functions applied to the constant \(1\).)

4.1 Refining clause learning↩︎

The main improvement² to the procedure concerns refinements (i.e., the conjunctions added in line [refinement:psi] of 4). Let \(h\) be either \(\sin\) or \(\exp\), and consider a conjunct of the form \[\label{eq:gamma} \gamma(h^{\text{in}}_t,h^\text{out}_t) \mathrel{\vcenter{:}}= (i_{\text{in}} \leq h^{\text{in}}_t \leq j_{\text{in}} \Rightarrow i_{\text{out}} \leq h^{\text{out}}_t \leq j_{\text{out}}),\,\tag{1}\] conjoined to \(\psi\) at this line. Our goal is to strengthen \(\psi\) by replacing \(\gamma\) with a formula \(\chi\) that entails it, while preserving the specification of 4. To obtain a substantial strengthening of \(\gamma\), we furthermore require \(\chi\) to satisfy:

P phantomsection :

for every \([a,b] \subseteq [\round{i_{\text{in}}}-1,\round{j_{\text{in}}}+1]\) with \(b - a \leq j_{\text{in}} - i_{\text{in}}\), where \(\round{k}\) stands for the integer closest to \(k\), the possible values \(h_t^{\text{out}}\) can take in a solution to \(a \leq h_t^{\text{in}} \leq b \land \chi\) all lie in an interval of width at most \(j_{\text{out}} - i_{\text{out}}\).

Property [property:p] guarantees that \(\chi\) approximates \(h\) well throughout an interval of width (at least) \(2\). We now define the formula \(\chi\) for the sine and exponential functions.

4.2 Definition of \(\chi\) for the sine function↩︎

For the sine function, we define \(\chi\) through Taylor approximations. Let \({f \colon \mathbb{R}\to \mathbb{R}}\) be an infinitely differentiable function, and \({c \in \mathbb{R}}\). We recall that the Taylor expansion of \(f\) centred at \(c\) is the infinite sum \(\sum_{k=0}^\infty \frac{f^{(k)}(c)}{k!}(x-c)^k\), where \(f^{(k)}\) denotes the \(k\)th derivative of \(f\). The truncation to the first \(n\) terms of the sum is the \(n\)th Taylor polynomial of \(f\) centred at \(c\).

For \(r \in \mathbb{Q}\) and \(d \in \mathbb{N}\), define \(({\rm T}^d_r\sin)(x)\) as the term \(\sum_{k=0}^d \frac{\sin^{(k)}(r \cdot \pi)}{k!}(x- r \cdot \pi_1^{\text{out}})^k\). It has the following properties:

1. If \(r\) is of the form \(2 \cdot j + \frac{1}{2}\) or \(2 \cdot j + \frac{3}{2}\), for some \(j \in \mathbb{Z}\), then \(\sin^{(k)}(r \cdot \pi)\) is in \(\{-1,0,1\}\).

2. Assuming that \(\pi_1^{\text{out}}\) is assigned the value \(\pi\), then \[\begin{align} \textstyle ({\rm T}^{4 \cdot n+2}_{2 \cdot j + \frac{1}{2}}\sin)(x) &\leq \sin(x) \leq ({\rm T}^{4 \cdot n}_{2 \cdot j + \frac{1}{2}}\sin)(x), \text{ and}\\ \textstyle ({\rm T}^{4 \cdot n}_{2 \cdot j + \frac{3}{2}}\sin)(x) &\leq \sin(x) \leq ({\rm T}^{4 \cdot n+2}_{2 \cdot j + \frac{3}{2}}\sin)(x), \end{align}\] for all \(j \in \mathbb{Z}\), \(n \in \mathbb{N}\), and \(x \in \mathbb{R}\).

3. For every \(r \in \mathbb{Q}\), \(y = ({\rm T}^d_r\sin)(x)\) tends to the graph of \(y = \sin(x)\) as \(d \to \infty\) and \(\pi_1^{\text{out}} \to \pi\).

Consider \(\gamma(\sin_t^{\text{in}},\sin_t^{\text{out}})\) from 1 . Let \(p_\ell,p_u \in \mathbb{Z}\) be such that \([\round{i_{\text{in}}}-1,\round{j_{\text{in}}}+1] \subseteq [2 \cdot \pi \cdot p_\ell, 2 \cdot \pi \cdot p_u]\) and \(p_u-p_\ell\) is minimal. For \(\ell < u \in \mathbb{Q}\) with \(\pi \in [\ell,u]\) and \(n \in \mathbb{N}\), let \(\chi_{[\ell,u]}^n(\sin_t^{\text{in}},\sin_t^{\text{out}},\pi_1^{\text{out}})\) be the formula given by: \[\begin{align} & \ell < \pi_1^{\text{out}} < u\\ \land& \bigwedge\nolimits_{j = p_\ell}^{p_u} \textstyle ({\rm T}^{4 \cdot n+2}_{2 \cdot j + \frac{1}{2}}\sin)(\sin_{t}^{\text{in}}) \leq \sin_{t}^{\text{out}} \leq ({\rm T}^{4 \cdot n}_{2 \cdot j + \frac{1}{2}}\sin)(\sin_{t}^{\text{in}})\\ \land& \bigwedge\nolimits_{j = p_\ell}^{p_u} \textstyle ({\rm T}^{4 \cdot n}_{2 \cdot j + \frac{3}{2}}\sin)(\sin_{t}^{\text{in}}) \leq \sin_{t}^{\text{out}} \leq ({\rm T}^{4 \cdot n+2}_{2 \cdot j + \frac{3}{2}}\sin)(\sin_{t}^{\text{in}}). \end{align}\] By Items [prop:sin:1] and [prop:sin:2], \(\chi_{[\ell,u]}^n\) is a satisfiable formula of NRA. Moreover, Item [prop:sin:3] implies that, as \(n\) increases and the interval \([\ell,u]\) is taken as a tighter approximation of \(\pi\), the set of solutions to \(\chi_{[\ell,u]}^n\) converges to that of the TRA formula \({\pi_1^{\text{out}} = \pi \land \sin_t^{\text{out}} = \sin(\sin_t^{\text{in}})}\). We may thus define \(\chi \mathrel{\vcenter{:}}= \chi_{[\ell,u]}^n\) for \(n\) sufficiently large and \([\ell,u]\) sufficiently tight rational approximation of \(\pi\). This construction is effective:

4.3 Definition of \(\chi\) for the exponential function↩︎

To define \(\chi\) for the exponential function, we can still rely on Taylor approximations for lower bounds. For \(c \in \mathbb{Z}\) and \({d \in \mathbb{N}}\), define \((\rm{T}^d_c \exp)(x)\) as the term \({(e_1^{\text{out}})^c \cdot \sum_{k=0}^d \frac{1}{k!}(x-c)^k}\). When \(e_1^{\text{out}}\) is assigned the value \(e\), this term corresponds to the \(d\)-th Taylor polynomial of \(\exp\) centred at \(c\), and we have \({(\rm{T}^{2n+1}_c \exp)(x) \leq \exp(x)}\) for every \(n \in \mathbb{N}\), \(c \in \mathbb{Z}\) and \(x \in \mathbb{R}\).

Taylor approximations are inadequate for upper-bounding the exponential function, as no finite polynomial approximation can serve as a global upper bound for a super-polynomial function. For this purpose, we use Padé approximants instead, which are generalizations of Taylor polynomials by rational functions [21].

Given \(c \in \mathbb{Z}\) and \(d \in \mathbb{N}\), define the term \[\begin{align} ({\rm P}^d_c \exp)(x) \mathrel{\vcenter{:}}= (e_1^{\text{out}})^c \cdot \frac{Q_d(c,x)}{Q_d(x,c)}, \end{align}\] where \(Q_d(x,y) \mathrel{\vcenter{:}}= \textstyle\sum_{k=0}^{d}\frac{\left(2d-k\right)!}{\left(d-k\right)!d!k!}\left(y-x\right)^{k}\). When \(e_1^{\text{out}}\) is assigned the value \(e\), this term corresponds to the \(d\)-th Padé approximant of \(\exp\) centred at \(c\), and has the following well-known properties for every \(n \in \mathbb{N}\) and \(c \in \mathbb{Z}\):

1. The polynomial \(Q_{2n+1}(x,c)\) has a single root \(r > c\), and \(\exp(x) \leq ({\rm P}^{2n+1}_c \exp)(x)\) for every \(x \in [c,r)\).

2. \(\exp(x) \leq ({\rm P}^{2n}_c \exp)(x)\) for every \(x \in (-\infty,c]\).
   (In this case, \(Q_{2n}(x,c)\) has no roots.)

3. \(y = ({\rm P}^d_c \exp)(x)\) tends to the graph of \(y = \exp(x)\) as \(d \to \infty\) (and the root \(r\) tends to infinity).

Let \(\gamma(\exp_t^{\text{in}},\exp_t^{\text{out}})\) as in 1 , and \(c \mathrel{\vcenter{:}}= \round{j_{\text{in}}}+1\). For \(\ell < u \in \mathbb{Q}\) with \(e \in [\ell,u]\), and \(n \in \mathbb{N}\), let \(\chi_{[\ell,u]}^n\) be the formula in variables \(\exp_t^{\text{in}}\), \(\exp_t^{\text{out}}\) and \(e_1^{\text{out}}\) given by: \[\begin{align} & \ell < e_1^{\text{out}} < u \land (T_c^{2n+1}\exp)(\exp_t^{\text{in}}) \leq \exp_t^{\text{out}} \land {}\\ & (\exp_t^{\text{in}} \leq c \Rightarrow \exp_t^{\text{out}} \leq (P_c^{2n}\exp)(\exp_t^{\text{in}})) \land {}\\ &\exists r \big( Q_{2n+1}(r,c) = 0 \land \big(c \leq \exp_t^{\text{in}} < r \Rightarrow {}\\ && \exp_t^{\text{out}} \leq (P_c^{2n+1}\exp)(\exp_t^{\text{in}})\big)\big). \end{align}\] One can show that [lemma:sinechi] also holds for the formula \(\chi_{[\ell,u]}^n\) above, with respect to \(\gamma(\exp_t^{\text{in}},\exp_t^{\text{out}})\). It thus suffices to compute \(\ell\), \(u\) and \(n\) as given in [lemma:sinechi], and define \(\chi \mathrel{\vcenter{:}}= \chi_{[\ell,u]}^n\).

5 Implementation and experimental evaluation↩︎

We implemented our procedure, restricted to the structure \((\mathbb{R}, 0, 1, \pi, e, +, \cdot, \sin, \exp, <, =)\), as an open-source plugin for Yices2, relying on its MCSAT infrastructure and NRA plugin. We use the ARB library [14] for the arbitrary-precision interval arithmetic required by the \(\delta\)-consistency check. Below, we refer to the resulting tool as Yices-TRA.

For the benchmark, we use the same instances and time limit as [6]: 2 512 instances with a time limit of 1000 seconds, and no memory limit. We refer to [6] for further details on the instances. The results of our experiments are available in the accompanying artifact [22].

We performed two runs of Yices-TRA, one with \(\delta\) fixed at \(10^{-3}\) and one in the SAT/UNSAT mode with decreasing \(\delta\). The latter is denoted as Yices-TRA (dec\(\delta\)) in plots. We compared our tool with the SMT solvers MathSAT(v5.6.15) and cvc5(v1.3.2), and with the \(\delta\)-satisfiability tool dReal (v4.21.06.2), with tolerances \(\delta = 10^{-3}\) and \(\delta = 10^{-9}\). The two runs of dReal performed very similarly, suggesting that varying \(\delta\) would not affect the results significantly. 6 gives more information on the methods these tools implement.

The results are presented using survival plots in 5, and show the competitiveness of our approach. Runs taking less than \(0.1\) seconds are normalized to \(0.1\) seconds for clarity. With \(\delta\) fixed, our tool reported 101 SAT, 954 UNSAT and 230 \(\delta\)-SAT, outperforming all other tools in each category. The SAT/UNSAT mode classified a further 9 SAT and 8 UNSAT instances that were reported as \(\delta\)-SAT by the fixed \(\delta\) run. In total, our tool classified as SAT or UNSAT 218 instances that remained unsolved by the other three tools.

Among the other tools, cvc5 performs better on SAT instances and MathSAT on UNSAT ones. Interestingly, although Yices-TRA solves significantly more UNSAT instances than MathSAT within short time limits (e.g., less than 10s), the latter tool almost catches up at 1000s, falling short of just 44 benchmarks. This seems to suggest that the linear lemmas MathSAT uses (see incremental linearization in 6) converge slower than our NRA lemmas, but might scale better. Exploring how the degrees of the polynomial lemmas affects the performance and scalability of our procedure, and the development of hybrid approaches to automatically tune these degrees, are interesting directions for future work.

The last plot also counts \(\delta\)-SAT instances as solved. We observe that the ratio of \(\delta\)-SAT to solved instances differs markedly between Yices-TRA and dReal. Specifically, with \({\delta = 10^{-3}}\), dReal classifies 643 instances as UNSAT and 558 as \(\delta\)-SAT; that is, 46% of the instances solved by dReal are classified as \(\delta\)-SAT. By contrast, Yices-TRA classifies only 18% of the instances it solves as \(\delta\)-SAT. Looking more closely at the \(\delta\)-SAT instances, Yices-TRA returned UNSAT for 226 instances that dReal classified as \(\delta\)-SAT, while only 9 instances classified as \(\delta\)-SAT by Yices-TRA were classified as UNSAT by dReal. This suggests that, while \(\delta\)-SAT might intuitively give the impression that a formula is “probably satisfiable”, in practice it often is not.

6 Related Work↩︎

Incremental Linearization. The idea of incremental linearization [6] is to lazily approximates all non-linear functions with linear constraints, and solve the approximation with linear solvers. If the approximation is found unsatisfiable, so is the original formula. Otherwise, the model is either actual or spurious. Linear lemmas are added to exclude the spurious model. MathSAT and cvc5 implement this procedure. Ours is also an abstraction-refinement approach, but uses NRA lemmas with higher-degree polynomials. A benefit of this approach is that we can return irrational solutions, while linear solvers are limited to rational ones. Another distinction is that both MathSAT and cvc5 are based on CDLC(T). We instead work in the MCSAT framework, which is known to perform particularly well with NRA, as it allows CAD computations to be performed locally near conflicting assignments.

Interval Constraint Propagation (ICP). ICP is a branch-and-prune technique [23], [24]. Given a starting box, is uses interval arithmetic to determine constraint signs over the box. If all constraint signs are constant, ICP returns whether the formula is satisfied. Else, it splits the box into smaller boxes and recurs. The results of the recursive calls are computed lazily: as long as one answers true, so does ICP. It is implemented in iSAT3 [25] and (under the \(\delta\)-satisfiability framework) in dReal. In our tool, we also use interval arithmetic, but its use is limited to computing the over-approximation of a function over a single point.

\(\delta\)-satisfiability. We have already introduced \(\delta\)-satisfiability in 2.2. The tool dReal is tailored for this relaxed notion of satisfiability, and can only return UNSAT or \(\delta\)-SAT. This differs from the \(\delta\)-consistency mode of our tool, which can still answer SAT, and treats NRA constraints exactly, since they are not affected by the TRA consistency check. The kSMT calculus [26] combines \(\delta\)-sat with linearization (however, the tool does not support transcendental functions).

Deductive methods. The MetiTarski [27] theorem prover can prove unsatisfiability for some formulas in TRAs. It replaces transcendental functions with user-defined lower and upper bounds, and relies on decision procedures for NRA to prove unsatisfiability. Compared to our method, MetiTarski is not able to prove satisfiability, and it only relies on the initial abstraction provided by the user.

Complete methods. As stated in the introduction, even small fragments of TRAs are undecidable. A notable exception is the result by Chen and Xia [12], who prove decidability for the univariate fragment with the sine function, where atomic formulas are restricted to rational polynomial inequalities in \(x\) and \(\sin(x)\) (disallowing complex terms inside sine). In proving decidability, they devise a root isolation algorithm that could in principle be used in our procedure to tighten the feasibility sets computed by the NRA plugin. This, together with methods capable of proving satisfiability in the presence of transcendental solutions [28], [29], as well as local search based approaches  [30]–[32], are promising directions for further improving our procedure.

Acknowledgments. This work is part of a project co-funded by the European Union (GA 101154447) and by MCIN/AEI (A CEX2024-001471-M and GA PID2022-138072OB-I00). Views and opinions expressed are those of the authors only and do not necessarily reflect those of the European Union or European Commission. Neither the European Union nor the granting authority can be held responsible for them.

Additional material for Section 2.3↩︎

Formulas↩︎

Throughout the paper, we assume that all formulas are quantifier-free and in conjunctive normal form (CNF). We recall that a literal is defined as an atomic formula or its negation. A clause is a disjunction of literals, and a CNF formula is a conjunction of clauses.

The CDCL Algorithm↩︎

The main data structure of CDCL is given by the trail, a sequence \({\ell_1 \mapsto b_1}\), …, \(\ell_n \mapsto b_n\) assigning Booleans \(b_i\) to propositional literals \(\ell_i\). The trail is consistent whenever \(\bigwedge_{i=1}^n(\ell_i \leftrightarrow b_i)\) is satisfiable. A consistent trail is complete whenever substituting each \(\ell_i\) by \(b_i\) in the input formula yields, after the application of simple validities (e.g., \(\top \lor \psi \leftrightarrow \top\)), to the formula \(\top\). The algorithm reports SAT as soon as it constructs a complete trail.

Elements of the trail are annotated with information to implement a backtracking mechanism: each element \(\ell_i \mapsto b_i\) is annotated as either decided or propagated, along with its decision level (an integer). Decided literals are guessed by the algorithm and increment the decision level. Propagated literals are derived by unit propagation, and inherit the current decision level. A clause is unit if all but one of its literals are assigned false in the trail. Unit propagation consists of assigning the remaining literal to true.

When a trail becomes inconsistent, CDCL analyses the annotations on the trail to derive a learned clause \(\psi\) that explains the inconsistency (a.k.a. conflict). By construction, this clause is implied by the input formula \(\varphi\), and the algorithm conjoins \(\psi\) to \(\varphi\). CDCL then backtracks (non-chronologically) to a decision level at which the learned clause becomes unit. If no such level exists, i.e., the conflict occurs at decision level 0, the algorithm terminates and reports UNSAT.

Additional material for Section 3↩︎

The goal of this appendix is to show that our procedure is sound ([theorem:soundness-delta]). We start with some additional definitions and notation. Given a term \(t\) (from a TRA) and a map \(\nu\) assigning a value to each variable in \(t\), we write \(\llbracket t \rrbracket_\nu\) for the evaluation of \(t\) under \(\nu\). It is defined inductively as one expects:

When it is clear from the context, we omit the subscript \(\nu\) and simply write \(\llbracket t\rrbracket\).

The following lemma is mostly helpful to clarify the computation performed in line [pseu:cons:K] of 3.

6.1 An invariant for the abstraction \(\varphi^\alpha\)↩︎

Throughout its execution, our procedure maintains and updates an NRA abstraction \(\varphi^\alpha\) of the input formula \(\varphi\). Towards a proof of soundness of the procedure, we show that this abstraction satisfies the following property throughout the procedure:

Invariant (\(\ast\)):

(1) All variables in \(\varphi\) occur in \(\varphi^\alpha\) and (2) for every solution \(\nu\) to \(\varphi\), \(\nu^{\text{ext}}\) is a solution to \(\varphi^\alpha\).

Let us prove that Invariant (\(\ast\)) is satisfied by the first abstraction \(\varphi^\alpha\). This is essentially [lemma:abstractionsound]:

Throughout the execution of our procedure, the abstraction \(\varphi^\alpha\) is only updated by conjoining it to the formulas computed by 4. The following lemma shows that these formulas are satisfied by the canonical extension of every solution to \(\varphi\).

6.2 Correctness of 3 and 4↩︎

We now move to the proof of the correctness of 4. We will first show it assuming that the formula \(\varphi^\alpha\) satisfies the invariant mentioned at the beginning of this section.

We can now complete the proof of soundness of the procedure:

Additional material for 4↩︎

We prove [lemma:sinechi] in both the cases of \(\chi^n_{[\ell,u]}\) refining the lemmas for the sine and for the exponential functions.

We proceed first with the proof for the sine function.

Now we comment on some straightforward adaptations for the proof for the exponential function.

References↩︎

---

1. Available at https://github.com/arith-lab/yices-tra.↩︎

2. A further improvement, omitted for lack of space, consists of adding clauses (none marked as original) to the initial abstraction \(\varphi^\alpha\), limiting the values of “in” and “out” variables. For example, we add constraints such as \(-1 \leq \sin_t^\text{out} \leq 1\), \(1 + \exp_t^{\text{in}} \leq \exp_t^\text{out}\), and \(3.14 < \pi_1^{\text{out}} < 3.15\).↩︎