\section{Introduction}
\textit{Quantifier elimination (QE)} is the process of transforming quantified formulas into logically equivalent quantifier-free formulas.
In this paper, we consider QE for the first order logic of real numbers ($\text{FOL}_{\mathbb{R}}$), so quantifiers range over the real numbers.
The Tarski-Seidenberg theorem proves that QE is admissible for the theory of real-closed fields \cite{Tarski, Seidenberg}.
Real quantified statements arise in a number of application domains, including geometry, chemistry, life sciences, and the verification of cyber-physical systems (CPS) \cite{DBLP:journals/mics/Sturm17}.
Many of the applications which require QE are safety-critical \cite{platzer2018logical}; accordingly, it is crucial to have both efficient and formally verified support for QE to trust the resulting decisions.

Unfortunately, QE algorithms are quite intricate, which makes it difficult to formally verify their correctness.
In practice, this necessitates the use of unverified tools.
For example, the theorem prover KeYmaera~X \cite{fulton2015keymaera}, which is designed to formally verify models of CPS (such as planes and surgical robots) uses Mathematica and/or Z3 as unverified blackbox solvers for QE.
Although these are admirable tools, they are unverified, and their use introduces a weak link \cite{Misfortunes} into what would otherwise be a (fully verified \cite{DBLP:conf/cpp/BohrerRVVP17}) trustworthy proof.

To help fill this gap, we formally verify linear and quadratic \textit{virtual substitution (VS)} due to Weispfenning \cite{weispfenning1988complexity, weispfenning1997quantifier}.
VS applies to QE queries that involve a quantified variable $x$ that occurs in a polynomial $f(x)$ with low degree (in theory, this means degree at most 4; in practice, it means degree 2, or at most 3).
VS is widely considered to be practical; it serves to improve SMT solvers and is central to REDLOG with good practical results \cite{sturm2018thirty}.
To our knowledge, ours is the first formally verified algorithm that handles VS with quadratic inequalities.

As we focus on correct and practical VS, we export our verified Isabelle/HOL code to SML for experimentation.
We test our exported formalization of the equality VS algorithm (\rref{sec:eqVsubst}) and of the general VS algorithm (\rref{sec:genVsubst}).
We compare to REDLOG, SMT-RAT, Z3, and Mathematica.
For SMT-RAT, we test both the general case of the tool and its restriction to VS.
With 291 examples, we solve fewer examples than existing tools (which solve at least 338).
This is to be expected given that other tools have been optimized and fine-tuned (some for decades) and use efficient general-purpose fallback QE algorithms when VS does not succeed.
But we provide the fastest VS implementation that is correct (assuming the orthogonal challenge of correct code generation from Isabelle \cite{DBLP:conf/esop/HupelN18}).

Our formalization is approximately 21,000 lines in Isabelle/HOL, available at: \url{https://pi314mm.github.io/assets/data/VirtualSubstitution.zip}.

\section{Related Work}
The fastest known QE algorithm is Cylindrical Algebraic Decomposition (CAD)\cite{collins1975quantifier}, which has not yet been fully formally verified.
There are few general-purpose formally verified QE algorithms, and there appears to exist a tradeoff between the practicality of an algorithm and the ease of formalization.
Mahboubi and Cohen verified Tarski's original QE algorithm \cite{AssiaQE} and McLaughlin and Harrison have a proof-producing QE procedure based on Cohen-H\"ormander \cite{mclaughlin2005proof}; unfortunately, Tarski's algorithm and Cohen-H\"ormander both have non-elementary complexity which limits the computational feasibility of these formalizations.

There has already been some work on formally verified VS: Nipkow \cite{nipkow2010linear} formally verified a VS procedure for \emph{linear} equations and inequalities.
The building blocks of $\text{FOL}_{\mathbb{R}}$~formulas, or ``atoms," in Nipkow's work only allow for polynomials of the form $\sum_i a_i x_i\sim c$ to be represented, where $\sim\ \in \{=,<\}$, the $x_i$'s are quantified variables and $c$ and the $a_i$'s are real valued coefficients.
These restrictions ensure that linear QE can always be performed, and they also simplify the substitution procedure and associated proofs.

Nipkow additionally provides a generic framework that can be applied to several different kinds of atoms.
Each new atom requires implementing several new code theorems in order to create an exportable algorithm, as some functions were defined generically.
While this is an excellent theoretical framework---we utilize several similar constructs in our formulation---we create an independent formalization that is specific to general $\text{FOL}_{\mathbb{R}}$~formulas, as our main focus is to provide an efficient algorithm in this domain.
Specializing to one type of atom allows us to implement several optimizations, such as our modified DNF algorithm, that would be unwieldy to develop in a generic setting.

Chaieb \cite{chaieb2008automated} discusses an extension of Nipkow's work to handle VS with quadratic equalities.
Chaieb never made his formalizations for the quadratic case publicly available, but in correspondence he generously provided us with his code in case we wanted to extend it directly.
While this was very helpful for reference, we chose to not build directly on his code.
The main reason is that Chaieb represented multivariate polynomials with a custom polynomial implementation, and we wanted to build our implementation on one of Isabelle/HOL's newer polynomial libraries, for purposes of consistency.
Another difference is that Chaieb explicitly designs VS to be used in conjunction with an auxiliary QE procedure, and we envision/develop VS as an exportable standalone procedure.

Some more specialized works are of interest.
For example, some work has been done in constraint solving with falsification: RSolver \cite{DBLP:conf/adhs/RatschanS06} was designed for hybrid systems verification and can find concrete counterexamples for fully quantified existential QE problems on \emph{compact} domains.
dReal \cite{DBLP:conf/cade/GaoKC13} is based on similar ideas and slightly relaxes the notion of satisfiability to $\delta$-satisfiability.
Constraint solving has also been considered in the SMT-solving with Z3's nlsat \cite{DBLP:conf/cade/JovanovicM12}, which uses CDCL to decide systems of nonlinear inequalities and equations.

\section{The Virtual Substitution Algorithm}
Informally (and broadly) speaking, VS discretizes the QE problem by solving for the roots of one or more low-degree polynomials $f_1(x), \dots, f_n(x)$.
When it applies, VS focuses on these roots and the intervals around them to identify and substitute appropriate representative ``sample points'' for $x$ in the rest of the formula.
However, these sample points may contain fractions, square roots, and/or other extensions of our language, and so they must be substituted ``virtually'': That is, VS creates a formula \textit{in $\text{FOL}_{\mathbb{R}}$\ proper} that behaves exactly as if the points had been directly substituted.
VS applies in two cases: an equality case and a general case.
We formalize both, and discuss each in turn.

\begin{remark}\label{rem:trick} Since the overall correctness proof for VS of multivariate polynomials proceeds per valuation, it turns out to suffice to state most (but not all) of our correctness lemmas for \emph{univariate} polynomials with \emph{real coefficients}.
We utilize this trick in the general VS case to simplify difficult proofs.
The virtual substitution algorithms, however, crucially need to work for \emph{multivariate} polynomials.
\end{remark}

\subsection{Equality Virtual Substitution Algorithm}\label{sec:eqVsubst}
Let $a, b$ and $c$ be arbitrary polynomials with real coefficients that do not mention the variable $x$.
Consider the formula $\exists x. (ax^2+bx+c=0 \land F).$
There are three possible cases: Either $a$ is nonzero, or $a = 0$ and $b$ is nonzero, or all of $a, b, c$ are zero (so $ax^2+bx+c=0$ uninformative).
Letting $F_x^r$ denote the substitution of $x=r$ for $x$ in $F$, and solving for the roots of $ax^2 + bx + c$, we have the following:
\\ \\
\begin{align*}
&\exists x. (ax^2+bx+c=0 \land F) \longleftrightarrow \\
&\Big( (a=0\land b=0 \land c=0 \land \exists x. F)~\lor \\
&(a=0 \land b\neq0 \land F_x^{-c/b})~\lor \\
&(a\neq 0\land b^2-4ac\geq 0\land (F_x^{(-b+\sqrt{b^2-4ac})/(2a)}\lor F_x^{(-b-\sqrt{b^2-4ac})/(2a)}))  \Big).
\end{align*}
If we were to substitute $F_x^{-c/b}$, $F_x^{(-b+\sqrt{b^2-4ac})/(2a)}$, and $F_x^{(-b-\sqrt{b^2-4ac})/(2a)}$ directly, the resulting formula would no longer be in $\text{FOL}_{\mathbb{R}}$.
Instead, VS avoids directly dividing polynomials or taking square roots by rewriting them equivalently in $\text{FOL}_{\mathbb{R}}$.
This involves two procedures: for fractions and for square roots.

To virtually substitute a fraction $p/q$ of polynomials where $q \neq 0$ into the atom $\sum_{i=0}^n a_ix^i\sim0$, where $\sim\ \in \{=, <, \leq, \neq\}$, it suffices to normalize the denominator of the LHS, with the caveat that for $<$ and $\leq$ atoms, we must ensure that we are not flipping the direction of the inequality by normalizing by a negative value.
When $n$ is even, $q^n\geq 0$ under any possible valuation, so normalizing by $q^n$ does not flip the inequality.
Alternatively, if $n$ is odd, $q^{n+1}\geq 0$.
We formalize this in our \isa{linear\isacharunderscore substitution} function (see \rref{app:linsubst}).

What remains is the problem of substituting $y = \sqrt{c}$ into an atom, where $c\geq 0$.
Its direct substitution can be separated out into even and odd exponents:
\begin{equation*}
\sum_{i=0}^n p_i(\sqrt{c})^i = \sum_{i=0}^{n/2}p_{2i}c^{i}+\sum_{i=0}^{n/2}p_{2i+1}c^{i}\sqrt{c}
\end{equation*}
Now our polynomial has the form $A+B\sqrt{c}$.
Then, we have the following cases:
\begin{align*}
A+B\sqrt{c}=0&\longleftrightarrow AB\leq 0 \land A^2-B^2c=0\\
A+B\sqrt{c}<0&\longleftrightarrow (A<0\land B^2c-A^2<0) \lor (B\leq 0 \land (A<0\lor A^2-B^2c<0))\\
A+B\sqrt{c}\leq 0&\longleftrightarrow (A\leq 0\land B^2c-A^2\leq 0) \lor (B\leq 0 \land A^2-B^2c\leq 0)\\
A+B\sqrt{c}\neq 0&\longleftrightarrow -AB<0 \lor A^2-B^2c\neq 0
\end{align*}
The equivalences for $=$ and $\neq$ atoms are derived from the observation that if $B \neq 0$, $A+B\sqrt{c}=0$ can be solved to find $\sqrt{c} = -A/B$, which holds iff $A^2 = B^2c$ and $-A/B \geq 0$.
The inequality cases involve casework to determine when polynomial $A$ is negative and dominates $B\sqrt{C}$ as $A^2>B^2C$, and when $B$ is negative and $B\sqrt{C}$ dominates $A$ as $B^2C>A$.
We formalize the virtual substitution procedure for quadratic roots in \isa{quadratic\isacharunderscore sub} (see \rref{app:quadsubst}).

\subsection{General Virtual Substitution Algorithm}\label{sec:genVsubst}
As we have seen, QE very naturally leads to finitely many cases (discretizes) for formulas that involve quadratic equality atoms (we call this the \textit{equality case}).
The VS algorithm for the \textit{general case}, which also handles inequality atoms, is more involved, because, unlike equalities, inequalities may have uncountably many solutions.
As a simple example, consider the formula $\exists x. (p<0 \land q<0)$, where $p$ and $q$ are the univariate quadratic polynomials (in variable $x$) depicted in Fig. \ref{graph}.
Noting that the roots of $p$ and $q$ cannot possibly satisfy the strict inequalities, we partition the number line into ranges in between these zeros.

\begin{wrapfigure}[9]{R}{0.44\textwidth}
\vspace{-1\baselineskip}
\includegraphics[width=0.4\textwidth]{graph.png}
\vspace{-0.5\baselineskip}
\caption{\label{graph}Two quadratics, their roots (black dots) and off-roots (red \textsf{x})}
\end{wrapfigure}

We recognize a key property: In each of the ranges between the roots of $p,q$, the signs of both $p$ and $q$ remain consistent.
Since the ranges cover all roots ($p,q$), the truth value of the formula at a single point in a range is representative of the truth value of the formula on the entire range.
To discretize the QE problem, we need only pick one sample point for each range.

However, we want to pick appropriate sample points \textit{for any} possible $p$ and $q$.
The points we pick as representatives are called the off-roots, which occur $\epsilon$ units away from the roots, where $\epsilon > 0$ is arbitrarily small.
We additionally need a representative for the leftmost range, which we represent with the point $-\infty$, where $-\infty$ is arbitrarily negative.
Of course, we cannot directly substitute $\epsilon$ and $-\infty$: they are not real numbers!
However, we can \textit{virtually} substitute them.


\subsubsection{Negative Infinity.}
Given any formula $F$, the VS of $-\infty$ should satisfy the equivalence $F_x^{-\infty} \longleftrightarrow \exists y.~\forall x{<}y.~F(x)$ (where $y$ does not occur in $F$).
Intuitively, this says that $-\infty$ acts as if it is arbitrarily negative (so less than the $x$ component of all roots of the polynomials in $F$) and captures information for the leftmost range on the real number line in any valuation of the non-$x$ variables. 

If formula $\exists y.~\forall x{<}y.~ax^2 + bx + c = 0$ is true, where $a, b, c$ are polynomials that do not involve $x$, then $ax^2 + bx + c = 0$ holds at infinitely many $x$; since nonzero polynomials have finitely many roots, this can only happen if $ax^2+bx+c$ is the zero polynomial in $x$, i.e., it holds that:
\begin{equation}\label{eq:eqatoms}
(ax^2+bx+c=0)_x^{-\infty}  \longleftrightarrow a=0\land b=0 \land c=0
\end{equation}

The negation of \rref{eq:eqatoms} captures the behavior of $\neq$ atoms. For $<$ atoms, note that the sign value at $-\infty$ is dominated by the leading coefficient, so:
\begin{equation*}
(ax^2+bx+c<0)_x^{-\infty} \longleftrightarrow a<0 \lor (a=0\land (b>0\lor (b=0\land c<0)))
\end{equation*}
Finally, \({(ax^2+bx+c\leq 0)_x^{-\infty}} {\longleftrightarrow} {(ax^2+bx+c= 0)_x^{-\infty}} \lor {(ax^2+bx+c< 0)_x^{-\infty}}\).

In Isabelle/HOL, we formalize that our virtual substitution of $-\infty$ satisfies the desired equivalence (on $\mathbb{R}$ using \rref{rem:trick}) in the following lemma:
\begin{isabelle}
\infinityevalLem
\end{isabelle}

To explain this lemma, we need to take a slight detour and discuss a few structural details of our framework (which is discussed in greater detail in \rref{sec:framework}).
The datatype \isa{atom'} contains a triple of real numbers (which represent the coefficients of a univariate quadratic polynomial) and a sign condition:
\begin{isabelle}
\datatypeatomprime
\end{isabelle}
The \isa{aEval'} function has type \isa{atom' $\Rightarrow$ real $\Rightarrow$ bool}; that is, it takes a sign condition, a triple of real numbers $(a,b,c)$ and a real number $x$ and evaluates whether $ax^2 + bx + c$ satisfies the sign condition.
The \isa{eval'} function has type \isa{atom' fm' $\Rightarrow$ real $\Rightarrow$ bool}, where an \isa{atom' fm'} is a formula that involves conjunctions and disjunctions of elements of type \isa{atom'} (and ``True'' and ``False'').
That is, the  \isa{eval'} function takes such a formula and a real number and evaluates whether the formula is true at the real number.
Thus, \isa{infinity\isacharunderscore eval'} states that, given \isa{At} of type \isa{atom'}, with tuple $(a,b,c)$ and sign condition $\sim\ \in \{<,=,\leq, \neq\}$, $\isa{At}_x^{-\infty}$ holds iff $\exists y. \forall x{<}y. a*x^2 + b*x + c \sim 0$.
This captures the desired equivalence.

Note that these definitions are set up for univariate polynomials (coefficients are assumed to be real numbers).
This is deliberate.
In \rref{sec:eval}, we will discuss how our framework reduces correctness results for multivariate polynomials to univariate lemmas like \isa{infinity\isacharunderscore eval'}.


\subsubsection{Infinitesimals.}\label{sec:Infinitesimals}Given arbitrary $r$ (not containing $x$), VS for $r + \epsilon$ for variable $x$ should capture the equivalence
$F_x^{r + \epsilon} \longleftrightarrow \exists y{>}r. \forall x{\in}(r, y].\ F(x)$, where $F$ does not contain $y$.
Intuitively, this says that (in any valuation of the non-$x$ variables) $r + \epsilon$ captures information for the interval between $r$ and the next greatest $x$-root.

For $=$ and $\neq$ atoms, we proceed in the same manner as we did with $-\infty$, as $(r,y]$ contains infinitely many points and only the zero polynomial has infinitely many solutions.
As before, $\leq$ atoms turn into disjunctions of the inequality and equality representations at $r + \epsilon$.
We are left only to consider $<$ atoms.

Consider $ (p{<}0)_x^{r + \epsilon}$ where $p = ax^2 + bx + c$ with polynomials $a,b,c$ not containing $x$, and an arbitrary $r$ not containing $x$.
Notice that if $(p{<}0)_x^r$, then because polynomials are continuous, we can choose a small enough $y$ so that $\forall x{\in}(r, y].\ p{<}0$.
If instead $(p=0)_x^r$, then consider the partial derivative of $p$ evaluated at $r$.
If $\frac{\partial p}{\partial x}(r)$ is negative, then $\exists y{>}r. \forall x{\in}(r, y].\ p{<}0$  holds, because $p$ is decreasing in $x$ after $x{=}r$.
If $\frac{\partial p}{\partial x}(r)$ is positive, then $ \exists y{>}r. \forall x{\in}(r, y].\ p{<}0$  \textit{does not} hold, because $p$ is increasing in $x$ after $x{=}r$.
If $\frac{\partial p}{\partial x}(r)$ is zero, then to ascertain whether $ \exists y{>}r. \forall x{\in}(r, y].\ p{<}0$, we will need to check higher  derivatives.

This pattern forms the following recurrence, with the base case $(p<0)_x^{r + \epsilon} = (p<0)_x^r$ for polynomials $p$ of degree zero:
\begin{equation*}
(p<0)_x^{r + \epsilon} ~\stackrel{\text{def}}{=}~ (p<0)_x^r \lor \big((p=0)_x^r\land (({\partial p}/{\partial x})<0)_x^{r + \epsilon}\big)
\end{equation*}
We use the VS algorithm from Section~\ref{sec:eqVsubst} to characterize $(p<0)_x^r$ and $(p=0)_x^r$.

In Isabelle/HOL, we show that given a quadratic root $r$, the virtual substitution of $r + \epsilon$ satisfies the desired equivalence in the following theorem (on $\mathbb{R}$ using Remark \ref{rem:trick}; we have an analogous lemma for linear roots $r$): \begin{isabelle}
\infinitesimalquad
\end{isabelle}

Note that \isa{{\isacharbraceleft}{\kern0pt}{r}{\kern0pt}{\isacharless}{\kern0pt}{\isachardot}{\kern0pt}{\isachardot}{\kern0pt}y{\kern0pt}{\isacharbraceright}} in Isabelle stands for the range $(r,y]$. This says that, given \isa{At} of type \isa{atom'}, with tuple $(a,b,c)$ and sign condition $\sim\ \in \{<,=,\leq, \neq\}$, $\isa{At}_x^{r+\epsilon}$ holds iff $\exists y > r.\forall x \in (r, y]. ax^2 + bx + c \sim 0$, which is the desired equivalence.

\subsubsection{The General VS Theorem.}\label{sec:genVS}
Now that we have explained virtually substituting $-\infty$ and infinitesimals, we are ready to state the general VS theorem.

Let $F$ be a formula of the following shape, where each $a_i,b_i,c_i$, and $d_i$ is a polynomial that is at most quadratic in variable $x$:
\begin{equation*}
    F = \left(\bigwedge a_i =0\right) \land \left(\bigwedge b_i < 0\right) \land \left(\bigwedge c_i \leq 0\right) \land \left(\bigwedge d_i \neq 0\right).
\end{equation*}
Let $R(p)$ denote the set of roots of an arbitrary polynomial $p$ in $x$, returning the empty set for the zero polynomial.
We now define
\begin{equation*}
A = \bigcup R(a_i)\quad B = \bigcup R(b_i)\quad C = \bigcup R(c_i)\quad D = \bigcup R(d_i)
\end{equation*}
Then we obtain the following QE equivalence:
\begin{equation}\label{eq:QEequiv1}
(\exists x. F) \longleftrightarrow F_x^{-\infty} \lor \bigvee_{r\in A\cup C}F_x^{r}\lor \bigvee_{r\in B\cup C\cup D}F_x^{r+\epsilon}
\end{equation}

Intuitively, this formula states that if there is a particular $x$ that satisfies $F$, then it must be the case that $x$ is one of the equality roots from $A \cup C$, or that $x$ falls in one of the particular ranges (including $-\infty$ as a range) obtained by partitioning the number line by the roots in $B\cup C \cup D$.

Equation \rref{eq:QEequiv1} can be further optimized to:
\begin{equation}\label{eq:QEequiv2}
(\exists x. F) \longleftrightarrow F^{-\infty} \lor \bigvee_{r\in A\cup  C}F_x^{r}\lor \bigvee_{r\in B\cup D}F_x^{r+\epsilon}.
\end{equation}
Intuitively, this optimization holds because polynomials are continuous.
More precisely, if $F$ has the shape $F = (p{\leq}0 \land G)$, and if $r$ is an $x$-root of $p$, then $r$ already satisfies $p{\leq}0$ in any valuation of the non-$x$ variables, so including $r+\epsilon$ as a sample point on account of $p{\leq}0$ is redundant.
It is possible that $G$ contains some atom $q < 0$ or $q \neq 0$ where $r$ is an $x$-root of $q$.
In this case, $r + \epsilon$ will already be a sample point on account of $q$, and we do not need to add it in on account of $p$.
Alternatively, if $G$ does not contain such a $q$, then, in any valuation of the non-$x$ variables, it is impossible for $G$ to be satisfied by $r + \epsilon$ and not $r$, meaning that it is redundant to include $r+\epsilon$ as a sample point on account of $G$.

The general QE theorem is proved in Isabelle/HOL as the following, using Remark \ref{rem:trick} to restrict to the univariate case:
\begin{isabelle}
\generalqe
\end{isabelle}

Here, \texttt{`} is the Isabelle/HOL syntax for mapping a function over a set.
This theorem says that if a finite-length formula \isa{F} is of the requisite shape, then there exists an $x$ satisfying \isa{F} iff \isa{F} is satisfied at $-\infty$ (captured by $\isa{F}_{\isa{inf}}$), or there is a root \isa{r} of one of the $=$ or $\leq$ atoms where \isa{F} \isa{r} holds, or if there is a root \isa{r} of one of the $<$ or $\neq$ atoms where \isa{F}$_{\epsilon}$ \isa{r} holds.
The proof is quite lengthy and involves a significant amount of casework; however, because we are working with univariate polynomials thanks to \rref{rem:trick}, this casework mostly reduces to arithmetic computations and basic real analysis for univariate polynomials, and some of what we need, such as properties of discriminants and continuity properties of polynomials, is already formalized in Isabelle/HOL's standard library.

\subsection{Top Level Algorithms}\label{sec:toplevel}
We develop several top-level algorithms that perform these VS procedures on multivariate polynomials; these are described in more detail in \rref{app:exported}.
Crucially, each features its own proof of correctness.
For example, for the \isa{VSEquality} algorithm, which performs equality VS repeatedly, we have:
\begin{isabelle}
  \correctnessEquality
\end{isabelle}

Here, the \isa{eval} function expresses the truth value of the (multivariate) input formula given a valuation \isa{xs}, represented as a list of real numbers.
Since we quantify over all possible valuations and express that they are the same before and after running the algorithm, we prove the soundness of \isa{VSEquality}.
The correctness of this theorem only relies on Isabelle/HOL's trusted core.

As our algorithms are general enough to handle formulas with high degree polynomials where VS does not apply, we cannot assert that the result is quantifier free (it might not be).
To demonstrate the practical usefulness of these algorithms, we export our code to SML and experimentally show that these algorithms solve many benchmarks.
The code exports rely on the correctness of Isabelle/HOL's code export, which ongoing work is attempting to establish \cite{DBLP:conf/esop/HupelN18}.

\section{Framework}\label{sec:framework}
We turn to a discussion of our framework, which is designed with two key goals in mind: First, perform VS as many times as possible on any given formula.
Second, reduce unwieldy multivariate proofs to more manageable univariate ones.

\subsection{Representation of Formulas}\label{sec:representation}
We define our type for formulas in the canonical datatype \isa{fm}:
\begin{isabelle}
\datatypefm
\end{isabelle}

As in Nipkow's previous work\cite{nipkow2010linear}, we use De Brujin indices to express the variables: That is, the 0th variable represents the innermost quantifier, and variables greater than the number of quantifiers represent the free variables.

We utilize the multivariate polynomial library\cite{sternagel2010executable} to define our atoms:
\begin{isabelle}
\datatypeatom
\end{isabelle}
Each atom is normalized without loss of generality so that the atom \isa{Less p} means $p<0$, \isa{Eq p} means $p=0$, and so on.

For example, the $\text{FOL}_{\mathbb{R}}$~formula $\forall x.((\exists y. x*a=y^2*b) \land \neg(\forall y. 5x^2\leq y))$ is represented in our framework as follows:
\begin{verbatim}
AllQ (And (ExQ (Atom (Eq (Var 1 * Var 2 - (Var 0)^2 * Var 3))))
         (Neg (AllQ (Atom (Leq (Const 5 * (Var 1)^2 - Var 0)))))),
\end{verbatim}
\isa{Const n} represents the constant $n \in \mathbb{R}$, and \isa{Var i} represents the $i$th variable.

Note that we could restrict ourselves to the $\top, \neg, \lor, \exists$ connectives and normalize $\leq$ and $\neq$ atoms to combinations of $<$ and $=$ atoms, and we could still express all of $\text{FOL}_{\mathbb{R}}$.
We avoid this for two reasons because it would linearly increase the size of the formula, and because we want to handle $\leq$ atoms in the optimized way discussed in \rref{sec:genVS} (see \rref{eq:QEequiv2}).
We do, however, allow for the normalization of $p=q$ into $p-q=0$.
This does not affect the size of the formula, and can afford simplifications: For example, $x^3+x^2+x+1=x^3$ becomes $x^2+x+1=0$.

\subsection{Modified Disjunctive Normal Form}\label{sec:dnf}
Nipkow's prior work \cite{nipkow2010linear} avoided incurring cases where linear VS does not apply by constraining atoms to be linear.
In order to develop a general-purpose VS method which can be used, e.g., as a preprocessing method for CAD, we must reason about cases where VS fails to perform QE for a specific quantifier, and still continue the execution of the algorithm to the remaining quantifiers to simplify the formula as much as possible.
To help with this, we implement a modified disjunctive normal form (DNF) that allows expressions to involve quantifiers.

\subsubsection{Contextual Awareness.}
Let us analyze how to increase the informational content in a formula with respect to a quantified variable of interest.

Say we wish to perform VS to eliminate variable $x$ in the formula $\exists x. F$, where $F$ is not necessarily quantifier free.
In linear time, we remove all negations from the formula by converting it into negation normal form.
We can then normalize $\exists x. F$ into the following form, where the $A_i$'s are atoms (and thus quantifier free):
\begin{equation*}
\exists x. \bigvee \left(\left(\bigwedge A_i\right) \land \left(\bigwedge \forall y. F_j\right) \land \left(\bigwedge\exists z. F_k\right) \right).
\end{equation*}

This normalization procedure is similar to standard DNF, handling quantified formulas as if they were atomic formulas.
We are allowed to distribute the existential quantifier across the disjuncts, which results in the equivalent formula:
\begin{equation}\label{dnf}
    \bigvee \exists x.\left(\left(\bigwedge A_i\right) \land \left(\bigwedge \forall y. F_j\right) \land \left(\bigwedge\exists z. F_k\right) \right).
\end{equation}
Now we run the VS algorithm, i.e. the input to VS is a conjunction of atomic formulas and quantified formulas in the shape of (\ref{dnf}).
Notice that if equality VS applies for atom $A_i$, then the relevant roots can be substituted into the quantified formulas $F_i$ and $F_k$, but roots from $F_i$ or $F_k$ cannot be substituted into $A_i$ if they feature quantified variables which are undefined in the broader context.
As such, our informational content is maximized when the number of $A_i$ atoms is maximized and the sizes of the $F_j$ and $F_k$ are minimized.

\subsubsection{Innermost Quantifier Elimination.}\label{sec:innermosttooutermost}
The innermost quantifier has an associated formula which is entirely quantifier free (and thus has no $F_j$ and $F_k$).
As such, we opt to perform VS recursively, starting with the innermost quantifier and moving outwards, hoping that VS is successful and the quantifier-free property is maintained.
This is not always optimal.
Consider the following formula:
\begin{equation*}
   \exists x. (x=0\land \exists y.\ xy^3+y=0).
\end{equation*}
If we attempt to perform quadratic VS on the innermost $y$ quantifier, it will fail.
However, performing VS on the $x$ quantifier first fixes $x=0$, which converts the $xy^3+y=0$ equality into the linear $y=0$. 
So, an (unoptimized) run of inside-out VS would produce $\exists y. y = 0$, and we could completely resolve the QE query by running VS again.
To account for these potential simplifications, it may sometimes be beneficial to run the general VS strategy several times, working from inwards to outwards each time (though as reflected in \rref{sec:experiments}, this does not currently make a difference on the benchmarks).

\subsubsection{Reaching Under Quantifiers.}\label{sec:reach}
We would like to recover usable information from the $F_k$ formulas to increase the informational content going into our QE algorithm.
It would be ideal if we could ``reach underneath" the existential binders and ``pull out'' the atoms from the formulas.
We can achieve this through a series of transformations.
If we pull out each existential quantifier one by one, we get the following formula, which is equivalent to formula (\ref{dnf}):
\begin{equation*}
    \bigvee \exists z_0.\cdots \exists z_n.\exists x. \left(\left(\bigwedge A_i\right) \land \left(\bigwedge \forall y. F_j\right) \land \left(\bigwedge F_k\right) \right)
\end{equation*}
This works because the rest of the conjuncts do not mention the quantified variable $z_i$ and adjacent existential quantifiers can be swapped freely (without changing the logical meaning of the formula).

We can then recursively unravel the formulas $F_k$, moving as many existential quantifiers as possible to the front.
Our implementation does this via a bottom-up procedure, starting underneath the innermost existential quantifier and building upwards, normalizing the formula into the form:
\begin{equation*}
    \bigvee \exists z_0.\cdots \exists z_n.\exists x. \left(\left(\bigwedge A_i\right) \land \left(\bigwedge \forall y. F_j\right) \right)
\end{equation*}

On paper, these transformations are simple as they involve named quantified variables; however, because our implementation uses a locally nameless form for quantifiers with DeBrujin indices, shifting an existential quantifier $A \land \exists z. F \longleftrightarrow \exists z. (A \land\ F)$ requires a ``lifting'' procedure $A\uparrow$ which increments all the variable indices by one.
This allows for the following conversion: $A\land \exists. F \longleftrightarrow \exists. (A\uparrow \land F)$.

We are unable to eliminate the $\forall$ quantifiers in a similar fashion, as one cannot exchange the order of $\forall$ and $\exists$ quantifiers without ruining logical equivalence.
However, we can still shrink the size of the quantified formula and expose some of the content it holds in some specific cases  described in \rref{app:pushforall}.

\subsection{Evaluation}\label{sec:eval}

Our proofs show that the input formula and the output formula (after VS) are \textit{logically equivalent}, i.e., have the same truth value under any valuation.
This needs a method of ``plugging in'' the real-valued valuation into the variables of the polynomials.
Towards this, we define the \isa{eval} function, which accumulates new values into the valuation as we go underneath quantifiers, and the \isa{aEval} function, which perform the substitution of the valuations into the polynomials.

When proving correctness, we focus our attention on one quantifier at a time.
By \rref{rem:trick} correctness of general VS case follows when considering a formula $F$ with a single quantifier, where $F$ contains only polynomials of at most degree two (otherwise general VS does not apply).
With these restrictions, we can substitute a valuation into the non-quantified variables, which transforms the multivariate polynomials into univariate polynomials.
For example, let $a, b,$ and $c$ be arbitrary multivariate polynomials that do not mention variable $x$.
Let $\hat{p} = \gamma(p)$ denote the insertion of the valuation $\gamma$ into polynomial $p$ ($\hat{p}$ is a real number).
We obtain the following conversion between multivariate and univariate polynomials:
\begin{equation*}
\textbf{eval}\ (\exists x. ax^2+bx+c = 0)\ \gamma \longleftrightarrow \textbf{eval'}\ (\exists x. \hat{a}x^2+\hat{b}x+\hat{c} = 0)\ \hat{x}
\end{equation*}

As such, we develop an alternative VS algorithm for univariate polynomials, where atoms are represented as triples of real-valued coefficients (as seen in \rref{sec:genVsubst}), and show that under this specific valuation, the multivariate output is equivalent per valuation with the output of the univariate case. Thus we finish the proof of the multivariate case by lifting the proof for the univariate case.

\subsection{Polynomial Contributions}
We build on the polynomials library \cite{sternagel2010executable}, which was designed to support executable multivariate polynomial operations.
This choice naturally comes with tradeoffs, and a number of functions and lemmas that we needed were missing from the library.
For example, we added a partial derivative function (needed to formalize VS for infinitesimals, as discussed in \rref{sec:Infinitesimals}); and prove its correctness for polynomials of degree at most two.

Additionally, we needed an efficient way to isolate the coefficient of a variable within a polynomial, which we define in the \isa{isolate\isacharunderscore variable\isacharunderscore sparse} function.
Similarly, we needed to extend code theorems for the generically defined \isa{insertion} function, which inserts a valuation into a polynomial, in order to be able to compute the result of valuations and export our program.
This required a number of lemmas that rely on the \isa{monomials} function, which separates a polynomial into a sum of monomials, and the \isa{degree} function, which computes the degree of a polynomial with respect to a single variable.

To make these new functions usable within proofs, we formulated a large collection of lemmas for polynomials with real-valued coefficients.
This involved several simplification lemmas for the interaction of the \isa{isolate\isacharunderscore variable\isacharunderscore sparse}, \isa{insertion}, and \isa{degree} functions across summations of polynomials and products of polynomials; for example, we proved that if the degrees of polynomials $p$ and $q$ in $x$ are both less than $d$, then the degree of $p+q$ in $x$ is also less than $d$.

One particularly critical lemma that we added allows us to take a multivariate polynomial in $\mathbb{R}[a_1,\dots, a_n, x]$ and rewrite it as a nested polynomial, i.e., a univariate polynomial in x with coefficients that are polynomials in $\mathbb{R}[a_1, \dots, a_n]$:
\begin{isabelle}
\sumOverDegree
\end{isabelle}

This is needed rather frequently within VS, as we often seek to re-express polynomials with respect to a single quantified variable of interest.

To utilize the variables within polynomials as DeBrujin indices, we implemented various lifting and substitution operations.
Our \isa{liftPoly} function takes in a lower limit $d$, an amount $a$, and a polynomial $p$ and returns a polynomial which reindexes variables within $p$ such that every variable greater or equal to $d$ is increased by $a$.
This is commonly denoted by $p\uparrow^a_d$.
This \isa{liftPoly} is needed in cases where we want to increase the number of quantifiers surrounding a formula.
For example, $\forall. ((\exists A) \land B)$ is equivalent to $\forall. \exists. (A \land (B\uparrow^1_1))$.

We also need an inverse function to \isa{liftPoly}, which we call \isa{lowerPoly}.
Whenever we have eliminated a quantified variable with our QE procedure, we can drop that quantifier and use \isa{lowerPoly} to reindex all the other variables accordingly.
We implement this in the \isa{clearQuantifiers} procedure.

\section{Experiments}\label{sec:experiments}

The benchmark suite consists of 378 QE problems in category CADE09\xspace collected from 94 examples \cite{DBLP:conf/cade/PlatzerQR09}, and category Economics\xspace with 45 examples \cite{DBLP:journals/corr/abs-1804-10037}.
CADE09\xspace and Economics\xspace examples were converted into decision problems, powers were flattened to multiplications, and CADE09\xspace were additionally rewritten to avoid polynomial division.
For sanity checking, we also negated the CADE09\xspace examples \cite{DBLP:conf/cade/PlatzerQR09}.
We run on commodity hardware.\footnote{\label{hardware} MacBook Pro 2019 with 2.6GHz Intel Core i7 (model 9750H) and 32GB memory (2667MHz DDR4 SDRAM).}


\paragraph{Tools.}
We compare the performance of our VSEquality\xspace (E\xspace), VSGeneral\xspace (G\xspace), VSGeneral\_3\_times\xspace (G3\xspace), VSLucky\xspace (L\xspace), and VSLEG\xspace (LEG\xspace) algorithms (\rref{app:exported}) to REDLOG~\cite{DBLP:journals/cca/Dolzmann097} (R\xspace), SMT-RAT~\cite{DBLP:conf/sat/CorziliusKJSA15} with CAD (S-CAD\xspace) and restricted to VS (S-VS\xspace), the SMT solver Z3 (Z3\xspace), and Mathematica (M\xspace).
All tools, save Mathematica, were run in Docker containers on Ubuntu 18.04 with 8GB of memory and 6 CPU cores.
\begin{description}
  \item[VS] Algorithms \isa{VSEquality\xspace}, \isa{VSGeneral\xspace}, \isa{VSGeneral\_3\_times\xspace}, \isa{VSLucky\xspace}, and \isa{VSLEG\xspace} (\rref{app:exported}). 
  Examples translated from SMT-LIB to our SML data structures, compiled with MLton\footnote{\tiny\url{http://mlton.org/}}; as a result, measurements do not include parsing
  \item[REDLOG] Snapshot 2021-04-13\footnote{\tiny\url{https://sourceforge.net/projects/reduce-algebra/files/snapshot_2021-04-13/}}, using the Reduce SMT REPL command ``read''.
  \item[SMT-RAT] Release 19.06\footnote{\tiny\url{https://github.com/ths-rwth/smtrat/releases/tag/19.06}} with default strategy ``SMTComp2016'' and virtual substitution strategy ``JustVS''
  \item[Z3] Release 4.8.10\footnote{\tiny\url{https://github.com/Z3Prover/z3/releases/tag/z3-4.8.10}}
  \item[Mathematica] Version 12.1.0, executed natively on MacOS, calls to \texttt{Resolve} restricted to 8GB of memory using \texttt{MemoryLimit} but no CPU core limit
\end{description}


\paragraph{Benchmark Results.}
Each example has a timeout of 30s.
Figure~\ref{fig:evaluationsummary} summarizes the performance on the CADE09\xspace and Economics\xspace examples in terms of the cumulative time needed to solve (return ``true'', ``false'', or ``unsat'') the fastest $n$ problems with a logarithmic time axis: more problems solved and a flatter curve is better.

\begin{figure}[t!b!p]  
\centering
\includegraphics[scale=1]{Figure2.pdf}
\caption{Cumulative time to solve fastest $n$ problems
(flatter and more is better)
}
\label{fig:evaluationsummary}

\end{figure}

Mathematica solves all problems in the CADE09\xspace category, closely trailed by Z3, Redlog.
The near constant computation time offset of Redlog in comparison to Z3, SMT-RAT, and Mathematica might be attributable to the additional step of entering an SMT REPL.
Our verified VSEquality\xspace (E\xspace), VSGeneral\xspace (G\xspace), VSLucky\xspace (L\xspace), and VSLEG\xspace (LEG\xspace) algorithms rank in performance between the basic VS implementation in SMT-RAT (S-VS\xspace),  full SMT approaches (S-CAD\xspace, Z3\xspace), and combined virtual substitution plus CAD implementations (R\xspace, M\xspace).
The reduced startup time of our algorithms is attributable to the omitted parsing step.
Overall, VSEquality\xspace and VSLucky\xspace solve examples fast, but the wider applicability of VSGeneral\xspace and VSLEG\xspace allows them to solve considerably more examples.
Though we have already implemented a number of optimizations for VS, which we discuss in \rref{app:optimiz}, we do not expect to outperform the tools at this stage, as many of them have been optimized over a period of many years.

A comparison of duration per problem is in \rref{fig:detailedevaluation}.
Though there is considerable overlap between VSEquality\xspace, VSGeneral\xspace, and VSLucky\xspace, mutually exclusive sets of solved examples (and considerable performance differences on a number of examples) foreshadow the performance achievable with the combined VSLEG\xspace algorithm.
The simplifications of earlier methods in the combined algorithm even allows VSLEG\xspace to solve additional examples in the given maximum duration.

\begin{figure}[thbp]
\centering
\includegraphics[scale=1.05]{Figure3.pdf}
\caption{CADE09\xspace duration per problem (color indicates duration, lighter is better)}
\label{fig:detailedevaluation}
\end{figure}

In \rref{fig:falseevaluation}, we compare the CADE09\xspace results to the results on negated CADE09\xspace examples to highlight \emph{contradictions} between answers (e.g., both $A$ and $\lnot A$ are claimed to be true).
Mathematica and Z3 answer consistently on both formula sets, and solve (almost) all examples.
REDLOG, the main VS implementation, in general does not perform well on the negated formulas and reports 96 contradictory answers, while SMT-RAT performs slightly better on the negated formulas, but still contradicts itself on 38 examples for both original and negated formula.
Our VSLEG\xspace algorithm has similar performance for proving and disproving in terms of absolute number of solved examples. As for S-VS\xspace, combining proving and disproving solves more examples than just one question individually.

\begin{figure}[b!tp]
\centering
\includegraphics[scale=1.01]{Figure4.pdf}
\caption{CADE09\xspace consistency comparison between original and negated formula: color indicates discrepancies within tools (green \captioncolorsquare{3}{Paired}: answer on original and negated formula agree, dark-blue \captioncolorsquare{1}{Paired}: only original solved, light-blue \captioncolorsquare{0}{Paired}: only negated solved, red+long \captioncolorsquare{5}{Paired}: \textbf{contradictory} answers (both formulas unsat/proved or both sat/disproved), empty: both timeout/unknown)}
\label{fig:falseevaluation}
\end{figure}

\section{Conclusion and Future Work}\label{sec:future}
We verify linear and quadratic virtual substitution; our algorithms are \textit{provably correct} up to Isabelle/HOL's trusted core and up to Isabelle's code export.
Developing the groundwork for practical verified real VS in Isabelle/HOL required significantly many low-level improvements and extensions to Isabelle's multivariate polynomials library.
Our experiments confirm that, even if verification were not a virtue in and of itself, QE is so subtle that our approach is the most competitive to get a practical correct VS implementation.
They also reveal the benefits of our current optimizations and indicate room for future improvements.
Further optimizations to the polynomial libraries, such as efficient coefficient lookup for polynomials using red black trees, would be welcome.
It would also be beneficial to treat homogeneous blocks of quantifiers as single entities, to support variable ordering heuristics.
Finally, we plan to eventually support cubic VS.

\paragraph{Acknowledgment.}
We wish to thank Fabian Immler for his substantial contributions at CMU to the polynomial Isabelle/HOL theories and regret that his current industry position precludes our ability to include him as a coauthor.

This material is based upon work supported by the National Science Foundation under Grant No. CNS-1739629, a National Science Foundation Graduate Research Fellowship under Grants Nos. DGE1252522 and DGE1745016, and by the AFOSR under grant number FA9550-16-1-0288. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the National Science Foundation or of AFOSR.

\let\oldbibliography\thebibliography
\renewcommand{\thebibliography}[1]{%
  \oldbibliography{#1}%
  \addtolength{\itemsep}{-6pt}%
}


