A Formalization of Convex Polyhedra\\ based on the Simplex Method\thanks{The authors were partially supported by the programme ``Ingénierie Numérique \& Sécurité'' of ANR, project ``MALTHY'', number ANR-13-INSE-0003, by a public grant as part of the Investissement d’avenir project, reference ANR-11-LABX-0056-LMH, LabEx LMH and by the PGMO program of EDF and FMJH.
Xavier Allamigeon\inst{1
We present a formalization of convex polyhedra in the proof assistant~\Coq{. The cornerstone of our work is a complete implementation of the simplex method, together with the proof of its correctness and termination. This allows us to define the basic predicates over polyhedra in an effective way (\ie~as programs), and relate them with the corresponding usual logical counterparts. To this end, we make an extensive use of the Boolean reflection methodology. The benefit of this approach is that we can easily derive the proof of several essential results on polyhedra, such as Farkas Lemma, duality theorem of linear programming, and Minkowski Theorem.
1706.10269.tex
\documentclass[runningheads]{llncs}

\usepackage[utf8]{inputenc}
\usepackage{amsmath,amsfonts,amssymb}
\usepackage{mathtools}
\usepackage{color}
\usepackage{paralist}
\usepackage{tikz}
\usepackage[disable,colorinlistoftodos,bordercolor=orange,backgroundcolor=orange!20,linecolor=orange,textsize=scriptsize]{todonotes}
\usepackage{hyperref}

\usepackage{inconsolata}

\newcommand{\etc}{etc}
\newcommand{\ie}{i.e.}
\newcommand{\eg}{e.g.}
\newcommand{\cf}{cf.}
\newcommand{\eps}{\varepsilon}
\newcommand{\R}{\mathbb{R}}
\newcommand{\Pcal}{\mathcal{P}}
\newcommand{\Qcal}{\mathcal{Q}}
\newcommand{\Ccal}{\mathcal{C}}
\newcommand{\Hcal}{\mathcal{H}}
\newcommand{\Fcal}{\mathcal{F}}
\newcommand{\Coq}{\textsc{Coq}}
\newcommand{\baux}{\tilde{b}}
\newcommand{\xaux}{\tilde{x}}
\newcommand{\vecone}{{\bf 1}}
\newcommand{\geqlex}{\mathrel{\geq_\mathrm{lex}}}
\newcommand{\gtlex}{\mathrel{>_\mathrm{lex}}}
\newcommand{\leqlex}{\mathrel{\leq_\mathrm{lex}}}
\newcommand{\gap}{\mathrm{gap}}
\newcommand{\CoqPolyhedra}{\textrm{Coq-Polyhedra}}
\newcommand{\MathComp}{Mathematical Components}
\newcommand{\MathCompShort}{MathComp}
\newcommand{\scalar}[2]{\langle #1, #2 \rangle}
\newcommand{\LP}{\mathrm{LP}}
\newcommand{\DualLP}{\mathrm{DualLP}}
\newcommand{\trans}[1]{#1^T}
\newcommand{\transinv}[1]{#1^{-T}}
\newcommand{\Pos}{K}
\newcommand{\Neg}{L}
\newcommand{\Ext}{\text{ext}}
\newcommand{\Opt}{\text{opt}}
\newcommand{\Init}{\text{init}}

\definecolor{dkblue}{rgb}{0,0.1,0.5}
\definecolor{lightblue}{rgb}{0,0.5,0.5}
\definecolor{dkgreen}{rgb}{0,0.4,0}
\definecolor{dk2green}{rgb}{0.4,0,0}
\definecolor{dkviolet}{HTML}{932191}
\definecolor{dkpink}{rgb}{1.2,0,1.6}
\definecolor{iden}{HTML}{0332FF}
\definecolor{comment}{HTML}{B12122}

\usepackage{listings}
\lstdefinelanguage{SSR}{
mathescape=true,
texcl=false,
morekeywords=[1]{From, Section, Module, End, Require, Import, Export, Defensive, Function, Variable, Variables, Parameter, Parameters, Axiom, Hypothesis, Hypotheses, Notation, Local, Tactic, Reserved, Scope, Open, Close, Bind, Delimit, Definition, Let, Ltac, Fixpoint, CoFixpoint, Add, Morphism, Relation, Implicit, Arguments, Set, Unset, Contextual, Strict, Prenex, Implicits, Inductive, CoInductive, Record, Structure, Canonical, Coercion, Theorem, Lemma, Corollary, Proposition, Fact, Remark, Example, Proof, Goal, Save, Qed, Defined, Hint, Resolve, Rewrite, View, Search, Show, Print, Printing, All, Graph, Projections, inside, outside, Locate, Maximal},
morekeywords=[2]{forall, exists, exists2, fun, fix, cofix, struct, match, with, end, as, in, return, let, if, is, then, else, for, of, nosimpl},
morekeywords=[3]{Type, Prop},
morekeywords=[4]{pose, set, move, case, elim, apply, clear, hnf, intro, intros, generalize, rename, pattern, after, destruct, induction, using, refine, inversion, injection, rewrite, congr, unlock, compute, ring, field, replace, fold, unfold, change, cutrewrite, simpl, have, gen, generally, suff, wlog, suffices, without, loss, nat_norm, assert, cut, trivial, revert, bool_congr, nat_congr, abstract, symmetry, transitivity, auto, split, left, right, autorewrite},
morekeywords=[5]{by, done, exact, reflexivity, tauto, romega, omega, assumption, solve, contradiction, discriminate},
morekeywords=[6]{do, last, first, try, idtac, repeat},
literate={isn't }{{{\ttfamily\color{dkgreen} isn't }}}1,
showstringspaces=false,
morestring=[b]",
morestring=[d]',
tabsize=2,
extendedchars=true,
sensitive=true,
breaklines=true,
basicstyle=\ttfamily,
captionpos=b,
columns=[l]fullflexible,
identifierstyle={\ttfamily\color{black}},
keywordstyle=[1]{\ttfamily\color{dkviolet}},
keywordstyle=[2]{\ttfamily\color{dkgreen}},
keywordstyle=[3]{\ttfamily\color{lightblue}},
keywordstyle=[4]{\ttfamily\color{dkblue}},
keywordstyle=[5]{\ttfamily\color{red}},
keywordstyle=[6]{\ttfamily\color{dkpink}},
stringstyle=\ttfamily,
commentstyle=\rmfamily,
}
\lstset{language=SSR,basicstyle=\ttfamily\small,keepspaces=true,aboveskip=\smallskipamount,belowskip=\smallskipamount,escapechar=\?}
\lstset{moredelim=[is][\color{iden}\bfseries\ttfamily]{|*}{*|}}
\lstset{moredelim=[is][\color{comment}\ttfamily]{((*}{*))}}
\lstset{moredelim=[is][\bfseries\ttfamily]{!*}{*!}}
\newenvironment{coq}{\begin{lstlisting}}{\end{lstlisting}}
\let\C=\lstinline

\usepackage{etoolbox}
\expandafter\patchcmd\csname \string\lstinline\endcsname{%
  \leavevmode
  \bgroup
}{%
  \leavevmode
  \ifmmode\hbox\fi
  \bgroup
}{}{%
  \typeout{Patching of \string\lstinline\space failed!}%
}

\usetikzlibrary{calc}
\usetikzlibrary{arrows}
\usetikzlibrary{positioning}

\title{A Formalization of Convex Polyhedra\\ based on the Simplex Method\thanks{The  authors were partially supported by the programme ``Ingénierie Numérique \& Sécurité'' of ANR, project ``MALTHY'', number ANR-13-INSE-0003, by a public grant as part of the Investissement d’avenir project, reference ANR-11-LABX-0056-LMH, LabEx LMH and by the PGMO program of EDF and FMJH.}}

\author{Xavier Allamigeon\inst{1} \and Ricardo D.~Katz\inst{2}}

\institute{INRIA and CMAP, Ecole polytechnique, CNRS, Université de Paris--Saclay
\email{xavier.allamigeon@inria.fr}\\[1ex]
\and CIFASIS-CONICET\\[-.2ex]
\email{katz@cifasis-conicet.gov.ar}}

\begin{document}

\maketitle

\begin{abstract}
We present a formalization of convex polyhedra in the proof assistant~\Coq{}. The cornerstone of our work is a complete implementation of the simplex method, together with the proof of its correctness and termination. This allows us to define the basic predicates over polyhedra in an effective way (\ie~as programs), and relate them with the corresponding usual logical counterparts. To this end, we make an extensive use of the Boolean reflection methodology. The benefit of this approach is that we can easily derive the proof of several essential results on polyhedra, such as Farkas Lemma, duality theorem of linear programming, and Minkowski Theorem.
\end{abstract}

\section{Introduction}

Convex polyhedra play a major role in many different application areas of mathematics and computer science, including optimization and operations research, control theory, combinatorics, software verification, compilation and program optimization, constraint solving, \etc. Their success mainly comes from the fact that they provide a convenient tradeoff between expressivity (conjunction of linear inequalities) and tractability. As an illustration of the latter aspect, linear programming, \ie, the class of convex optimization problems over linear inequality constraints, can be solved in polynomial time~\cite{Khachiyan80}. 

Among the aforementioned applications of polyhedra, there are some which are critical. For instance, in software verification or control theory, polyhedra are used to provide guarantees on the safety of programs~\cite{CousotHalbwachs78} or the stability of dynamical systems~\cite{Guglielmi2017}. On the mathematical side, polyhedra are still a very active research subject. Let us mention Steve Smale's $9^\text{th}$ problem for the $21^\text{th}$ century (whether linear programming can be solved in strongly polynomial complexity)~\cite{Smale98}, or the open questions on the diameter of polytopes following the disproof of the Hirsch conjecture~\cite{Santos2012}. In particular, (informal) mathematical software play an increasing role in testing or disproving conjectures (see \eg{}~\cite{Bremner2013}). All this strongly motivates the need to formalize convex polyhedra in a proof assistant, in order to increase the level of trust in their applications. 

In this paper, we present the first steps of a formalization of the theory of convex polyhedra in the proof assistant~\Coq{}. A motivation for using~\Coq{} comes from the longer term objective of formally proving some mathematical results relying on large-scale computation (\eg{}, Santos' counterexample to the Hirsch conjecture~\cite{Santos2012}). The originality of our approach lies in the fact that our formalization is carried out in an effective way, in the sense that the basic predicates over polyhedra (emptiness, boundedness, membership, etc) are defined by means of \Coq{} programs. All these predicates are then proven to correspond to the usual logical statements. The latter take the form of the existence of certificates: for instance, the emptiness of a polyhedron is shown to be equivalent to the existence of a certificate \emph{a la} Farkas (see Corollary~\ref{cor:farkas} for the precise statement). This equivalence between Boolean predicates and formulas living in the kind \C$Prop$ is implemented by using the boolean reflection methodology, and the supporting tools provided by the \MathComp{} library and its tactic language~\cite{MathComp}. The benefit of the effective nature of our approach is demonstrated by the fact that we easily arrive at the proof of important results on polyhedra, such as several versions of Farkas Lemma, duality theorem of linear programming, separation from convex hulls, Minkowski Theorem, etc.

Our effective approach is made possible by implementing the simplex method inside~\Coq{}, and proving its correctness and termination. Recall that the simplex method is the first algorithm introduced to solve linear programming~\cite{Dantzig51}. Two difficulties need to be overcome to formalize it. On the one hand, we need to deal with its termination. More precisely, the simplex method iterates over the so-called bases. Its termination depends on the specification of a pivoting rule, whose aim is to determine, at each iteration, the next basis. In this work, we have focused on proving that the lexicographic rule~\cite{Dantzig1955} ensures termination. On the other hand, the simplex method is actually composed of two parts. The part that we previously described, called Phase II, requires an initial basis to start with. Finding such a basis is the purpose of Phase I. It consists in building an extended problem (having a trivial initial basis), and applying to it Phase~II. Both phases need to be formalized to obtain a fully functional algorithm. 

We point out that our goal here is \emph{not} to obtain a practically efficient implementation of the simplex method (\eg{}, via the code extraction facility of \Coq{}). Rather, we use the simplex method as a tool in our proofs and, in fact, it turns out to be the cornerstone of our approach, given the intuitionistic nature of the logic in \Coq{}. Thus, we adopt the opposite approach of most textbooks on linear programming where, firstly, theoretical results (like the ones mentioned above) are proven, and then the correctness of the simplex method is derived from them.

The formalization presented in this paper can be found in a library developed by the authors called~\CoqPolyhedra{}.\footnote{Available in a git repository at~{\small \url{https://github.com/nhojem/Coq-Polyhedra}}.} As mentioned above, our formalization is based on the \MathComp{} library (\MathCompShort{} for short). On top of providing a convenient way to use Boolean reflection, this library contains most of the mathematical tools needed to formalize the simplex method (linear algebra, advanced manipulations of matrices, \etc{}).

\vspace*{-1ex}
\subsubsection*{Related Work.} Our approach has been strongly influenced by the formalization of abstract linear algebra in the \MathComp{} library, which is done in an effective way by exploiting a variant of Gaussian elimination~\cite{Gonthier2011}. 

As far as we know, this is the first formalization of the simplex method in the Calculus of Constructions. In this paradigm, the only work concerning convex polyhedra we are aware of is the implementation of Fourier--Motzkin elimination on linear inequalities in \Coq{}, leading to a proof of Farkas Lemma~\cite{Sakaguchi2016}. Our work follows a different approach, relying on the theory of linear programming, which has the advantage of providing certificates for the basic predicates over polyhedra. Concerning other families of logics, HOL Light provides a very complete formalization of convex polyhedra, including several important results (Farkas Lemma, Minkowski Theorem, Euler--Poincaré formula, \etc)~\cite{Harrison2013}. The classical nature of the logic implemented in HOL Light makes it difficult to compare this work with ours. In Isabelle, an implementation of a simplex-based satisfiability procedure for linear arithmetics has been carried out~\cite{Spasic2012}. This is motivated by obtaining a practical and executable code for SMT solving purposes. Here, we are driven by using the simplex method for mathematical proving, which explains why we obtain a completely different kind of formalization.

Finally, the theory of convex polyhedra is widely used in the area of formal proving as an ``informal backend'' which helps to establish the validity of some linear inequalities. In more detail, such inequalities are proven by formally checking certificates which are built by untrusted oracles based on linear programming. As illustrations, this allows to automate the deduction of some linear inequalities in proof assistants (see \eg{}~\cite{Besson2007}), or to certify the computations made by static analysis tools~\cite{Fouilhe2014}.

\vspace*{-1ex}
\subsubsection*{Organization of the Paper.} In Sect.~\ref{sec:preliminaries}, we introduce basic concepts and results on polyhedra and linear programming. In Sect.~\ref{sec:bases}, we describe the main components of the simplex method, and start its formalization. The lexicographic rule is dealt with in Sect.~\ref{sec:lex}. The two phases of the simplex method are formalized in Sect.~\ref{sec:PhaseII} and~\ref{sec:Phase1}, along with some of the main mathematical results that can be derived from them. Finally, we discuss the outcome of our work in Sect.~\ref{sec:outcome}.

By convention, all \Coq{} definitions, functions, theorems, \etc{} introduced in our work are highlighted in blue. This is to distinguish them from the existing material, in particular, the ones brought from the \MathCompShort{} library. We inform the reader that the vast majority of the results described in this paper (especially the ones of Sect.~\ref{sec:bases} to~\ref{sec:Phase1}) are gathered in the file \C$simplex.v$ of \CoqPolyhedra.

\section{Polyhedra, Linear Programming and Duality}\label{sec:preliminaries}

A {\em (convex) polyhedron} is a set of the form $\Pcal(A,b) \coloneqq \{ x \in \R^n \mid A x \geq b \}$, where $A \in \R^{m \times n}$ and $b \in \R^m$. The notation $\geq$ stands for the partial ordering over vectors, meaning that $y \geq z$ when $y_i \geq z_i$ for all $i$. In geometric terms, a polyhedron corresponds to the intersection of finitely many halfspaces. A \emph{(affine) halfspace} refers to a set of the form $\{ x \in \R^n \mid \scalar{a}{x} \geq \beta \}$, where $a \in \R^n$, $\beta \in \R$, and $\scalar{\cdot}{\cdot}$ stands for the Euclidean scalar product, \ie{}, $\scalar{x}{y} \coloneqq \sum_i x_i y_i$.

More generally, convex polyhedra can be defined over any ordered field. This is why our formalization relies on a variable \C$R$ of the type $\C$realFieldType$$ of \MathCompShort{}, whose purpose is to represent an ordered field in which the inequality is decidable. Assume that \C$m$ and \C$n$ are variables of type \C$nat$. The types \C$'M[R]_(m,n)$ and \C$'cV[R]_m$ 
provided by \MathCompShort{} respectively represent matrices of size $\C$m$ \times \C$n$$ and column vectors of size \C$m$ with entries of type \C$R$. In this paper, we usually omit \C$R$ in the notation of these types, for the sake of readability. The polyhedron associated with the matrix \C$A:'M_(m,n)$ and the vector \C$b:'cV_m$ is then defined by means of a Boolean predicate, using the construction \C$pred$ of \MathCompShort{}:
\begin{lstlisting}
Definition |*polyhedron*| A b := [pred x:'cV_n | (A *m x) >=m b].
\end{lstlisting}
Here, \C$*m$ stands for the matrix product, and \C$>=m$ for the entrywise ordering of vectors: \C$y <=m z$ if and only if \C$y i 0 <= z i 0$ for all \C$i$, where \C$y i 0$ and \C$z i 0$ are respectively the $i$th entry of the vectors \C$y$ and \C$z$ (see \C$vector_order.v$).

\emph{Linear programming} consists in optimizing a linear map $x \in \R^n \mapsto \scalar{c}{x}$ over a polyhedron, such as:
\begin{equation}\tag{$\LP(A,b,c)$} \label{eq:lp}
\begin{array}{r@{\quad}l}
\text{minimize} & \scalar{c}{x} \\[\jot]
\text{subject to} & A x \geq b \, , \; x \in \R^n
\end{array}
\end{equation}
Let us introduce a bit of terminology. A problem of the form~\ref{eq:lp} is referred to as a \emph{linear program} (see Fig.~\ref{fig:lp} for an example). A vector $x \in \R^n$ satisfying the constraint $A x \geq b$ is a \emph{feasible point} of this linear program. The polyhedron $\Pcal(A,b)$, which consists of the feasible points, is called the \emph{feasible set}. The map $x \mapsto \scalar{c}{x}$ is the \emph{objective} function. The \emph{optimal value} is defined as the infimum of $\scalar{c}{x}$ for $x \in \Pcal(A,b)$. A point $x \in \Pcal(A,b)$ reaching this infimum is called \emph{optimal solution}. When $\Pcal(A,b)$ is not empty, the linear program~\ref{eq:lp} is said to be \emph{feasible}, and its optimal value is either finite, or $-\infty$ (when the quantity $\scalar{c}{x}$ is not bounded from below over $\Pcal(A,b)$). In the latter case, we say that the linear program is \emph{unbounded (from below)}. Finally, when $\Pcal(A,b)$ is empty, the linear program is \emph{infeasible}, and its value is defined to be $+\infty$. 

A fundamental result in linear programming relates the optimal value of \ref{eq:lp} with the one of another linear program which is dual to it. In more detail, the \emph{dual linear program} of~\ref{eq:lp} is the following linear program:
\begin{equation}\tag{$\DualLP(A,b,c)$} \label{eq:dual_lp}
\begin{array}{r@{\quad}l}
\text{maximize} & \scalar{b}{u} \\[\jot]
\text{subject to} & \trans{A} u = c \, , \; u \geq 0 \, , \; u \in \R^m
\end{array}
\end{equation}
where $\trans{A}$ stands for the transpose of $A$. Notice that \ref{eq:dual_lp} is a linear program as well. Indeed, its constraints can be rewritten into a block system $\biggl(\begin{smallmatrix} \trans{A} \\[-0.4ex] -\trans{A} \\ I_m \end{smallmatrix}\biggr) u \geq \begin{psmallmatrix} c \\ -c \\ 0 \end{psmallmatrix}$, where $I_m$ is the $m \times m$ identity matrix. Besides, the maximization problem can be turned into a minimization problem with objective function $x \mapsto \scalar{-b}{x}$. We denote by $\Qcal(A,c)$ the feasible set  of~\ref{eq:dual_lp}, and we refer to it as the \emph{dual polyhedron}. Assuming \C$c$ is a variable of type $\C$'cV_n$$ (\ie, representing a vector in $\R^n$), we adopt a specific formalization for this polyhedron, as follows:
\begin{lstlisting}
Definition |*dual_polyhedron*| A c := [pred u:'cV_m | A^T *m u == c & (u >=m 0)].
\end{lstlisting}
As opposed to the dual linear program, $\LP(A,b,c)$ is referred to as the \emph{primal linear program}. The interplay between the primal and dual linear programs is described by the following result:
\begin{theorem}[Strong duality]\label{th:strong_duality}
If one of the two linear programs \ref{eq:lp} or \ref{eq:dual_lp} is feasible, then they have the same optimal value.

In addition, when both are feasible, then the optimal value is attained by a primal feasible point $x^* \in \Pcal(A,b)$ and by a dual feasible point $u^* \in \Qcal(A,c)$.
\end{theorem}

\begin{figure}[t]
\begin{center}
\begin{tikzpicture}[convex/.style={draw=none,fill=lightgray,fill opacity=0.7},convexborder/.style={very thick},point/.style={orange},level_set/.style={blue!50, ultra thick, dotted},>=stealth']
\begin{scope}[scale=.5]
\clip (-0.5,-0.5) rectangle (9.5,8.5);

\draw[help lines,gray!40] (-1,-1) grid (10,10);
\draw[gray!70,->] (-1,0) -- (9.5,0) node[above left] {$x_1$};
\draw[gray!70,->] (0,-1) -- (0,8.5) node[below right] {$x_2$};

\coordinate (v1) at (3,1);
\coordinate (v2) at (6,1);
\coordinate (v3) at (8,5);
\coordinate (v4) at (2,7);
\coordinate (v5) at (1,3);

\fill[convex] (v1) -- (v2) -- (v3) -- (v4) -- (v5) -- cycle;
\draw[convexborder] (v1) -- (v2) -- (v3) -- (v4) -- (v5) -- cycle;

\filldraw[point] (v1) circle (4pt) node[below right=0.2ex and -.8ex] {$\{1,5\}$};
\filldraw[point] (v2) circle (4pt) node[below right=0ex and -.4ex] {$\{4,5\}$};
\filldraw[point] (v3) circle (4pt) node[above right=0ex and -.8ex] {$\{2,4\}$};
\filldraw[point] (v4) circle (4pt) node[above left=0ex and -.8ex] {$\{2,3\}$};
\filldraw[point,green!70!black] (v5) circle (4pt) node[above left=0ex and -.8ex] {$\{1,3\}$};

\coordinate (objective) at ($(3,1)$);
\path let \p1=(objective) in coordinate (level_dir) at ($10*(-\y1, \x1)$);
\draw[dashdotted, thick,<-, >=stealth'] ($(0,2) -0.16*(objective)$) -- ($(0,2)+5*(objective)$) ; 
\draw[level_set] (5,-1) -- + (level_dir);
\draw[level_set] (3,-1) -- + (level_dir);
\draw[level_set] (1,-1) -- + (level_dir);
\end{scope}
\begin{scope}[shift={(8,2)}]
\node[text width=6cm,font=\normalfont] {
\[
\begin{array}{r@{\quad}l}
\text{minimize} & \, 3 x_1 + x_2 \\[1ex]
\text{subject to} & 
\begin{aligned}[t]
x_1 + x_2 & \geq 4 \\ 
-x_1 - 3x_2 & \geq -23 \\ 
4 x_1 - x_2 & \geq 1 \\ 
- 2x_1 + x_2 & \geq -11 \\ 
x_2 & \geq 1
\end{aligned}
\end{array}
\]
};
\begin{scope}[shift={(2.5,0.73)}]
\node[draw,circle,font=\tiny,inner sep=1pt] at (0,0) {$1$};
\node[draw,circle,font=\tiny,inner sep=1pt] at (0,-0.5) {$2$};
\node[draw,circle,font=\tiny,inner sep=1pt] at (0,-1) {$3$};
\node[draw,circle,font=\tiny,inner sep=1pt] at (0,-1.5) {$4$};
\node[draw,circle,font=\tiny,inner sep=1pt] at (0,-2) {$5$};
\end{scope}
\end{scope}
\end{tikzpicture}
\end{center}
\caption{A linear program. The feasible set is depicted in gray. The direction in which the objective function decreases is represented by a dashdotted oriented line, and some level sets (\ie, sets of the form $\{ x\in \R^2 \mid 3 x_1 + x_2 = \alpha \}$, where $\alpha \in \R$) are drawn in blue dotted lines. The basic points are represented by dots, and are annotated with the corresponding bases. The optimal basic point is highlighted in green.}\label{fig:lp}
\end{figure}

In particular, when \ref{eq:dual_lp} is feasible and its optimal value is $+\infty$, the primal linear program \ref{eq:lp} is necessarily infeasible. This holds for any choice of the vector $c$, including $c = 0$. Observe that $\DualLP(A,b,0)$ obviously admits $u = 0$ as a feasible point. Hence, we readily obtain a characterization of the emptiness of the polyhedron $\Pcal(A,b)$:
\begin{corollary}[Farkas Lemma]\label{cor:farkas}
The polyhedron $\Pcal(A,b)$ is empty if, and only if, the optimal value of $\DualLP(A,b,0)$ is $+\infty$, or, equivalently, there exists $u \in \R^m$ such that $u \geq 0$, $\trans{A} u = 0$ and $\scalar{b}{u} > 0$.
\end{corollary}

The first part of Corollary~\ref{cor:farkas} shows a way to formalize the emptiness property of polyhedra in an effective fashion, \eg, as a program computing the value of $\DualLP(A,b,0)$ inside the proof assistant and comparing it to $+\infty$. This is precisely the approach that we have adopted in this work. As we shall see in Sect.~\ref{sec:outcome}, it also applies to other properties over polyhedra.

\section{The Three Ingredients of the Simplex Method}\label{sec:bases}

\subsubsection{Bases and Basic Points.} In order to solve the linear program~\ref{eq:lp}, the simplex method iterates over the feasible bases, up to reaching one corresponding to an optimal solution or concluding that the optimal value is $-\infty$. A \emph{basis} is a subset $I$ of $\{1,\ldots ,m\}$ with cardinality $n$ such that the square matrix $A_I$, formed by the rows $A_i$ of $A$ indexed by $i \in I$, is invertible. With each basis $I$, it is associated a \emph{basic point} defined as $x^I \coloneqq (A_I)^{-1} b_I$. The basis $I$ is said to be \emph{feasible} when the point $x^I$ is feasible. It is said to be \emph{optimal} when $x^I$ is an optimal solution of the linear program. We refer to Fig.~\ref{fig:lp} for an illustration.

In geometric terms, a basis corresponds to a set of $n$ hyperplanes $A_i x = b_i$ which intersect in a single point. The basis is feasible when this point belongs to the feasible set $\Pcal(A,b)$. It can be shown that feasible basic points precisely correspond to the vertices, \ie, the $0$-dimensional faces, of the polyhedron $\Pcal(A,b)$. 

Formalization of bases and feasible bases is performed by introducing three layers of types. We start with a type  corresponding to \emph{prebases}, \ie, subsets of $\{1,\ldots ,m\}$ with cardinality $n$. 
\begin{lstlisting}
Inductive |*prebasis*| := !*Prebasis*! (I: {set 'I_m}) of (#|I| == n).
\end{lstlisting}
Here, \C$'I_m$ stands for the finite subtype of \C$i:nat$ such that \C$i < m$ (\cf~Interface \C$finType$ of \MathCompShort{}). A term \C$I$ of type \C${set 'I_m}$ represents a finite set of elements of type \C$'I_m$, and \C$#|I|$ corresponds to its cardinality. 

Defining bases then requires us to deal with submatrices of the form $A_I$. This is the purpose of the library \C$row_submx.v$ of~\CoqPolyhedra{}, where we define:
\begin{lstlisting}
Definition |*row_submx*| (A:'M_(m,n)) (I:{set 'I_m}) := 
	(\matrix_(i < #|I|, j < n) A (enum_val i) j) : 'M_(#I, n).
\end{lstlisting}
In this definition, \C$\matrix_(i < p,j < q) Expr(i,j)$ is the matrix (of type \C$'M_(p,q)$) whose $(i,j)$ entry is \C$Expr(i,j)$. The function \C$enum_val$ retrieves the \C$i$th element of the set \C$I$. Even when \C$I$ has cardinality \C$n$, the submatrix \C$row_submx A I$ does not have type \C$'M_n$, \ie, that of square matrices of size $\C$n$ \times \C$n$$. Indeed, in \MathCompShort{}, matrices are defined using dependent types (depending on the size). Thus the two types \C$'M_n$ and \C$M_(#|I|,n)$ are distinct, and we use the function \C$castmx$ to explicitly do the glueing job. The square matrix $A_I$ is thus formalized as follows:
\begin{lstlisting}
Definition |*matrix_of_prebasis*| (A:'M_(m,n)) (I:prebasis) := 
	castmx (prebasis_card I, erefl n) (row_submx A I) : 'M_n.
\end{lstlisting}
where \C$|*prebasis_card*| I$ is a proof of the fact that \C$#|I| = n$ and \C$erefl n$ of the fact that \C$n = n$. Assuming the variables \C$A:'M_(m,n)$ and \C$b:'cV_m$ have been previously declared, the type representing bases is then defined by:
\begin{lstlisting}
Inductive |*basis*| := !*Basis*! (I:prebasis) of (matrix_of_prebasis A I) \in unitmx.
\end{lstlisting}
where the type \C$unitmx$ represents the set of invertible matrices. The basic point associated with a basis \C$I$ is determined by a function called \C$|*point_of_basis*|$: 
\begin{lstlisting}
Definition |*point_of_basis*| (I:basis) := 
  (invmx (matrix_of_prebasis A I)) *m (matrix_of_prebasis b I).
\end{lstlisting}
where \C$invmx Q$ returns the inverse of the matrix \C$Q$. From this, we can define the type of feasible bases:
\begin{lstlisting}
Inductive |*feasible_basis*| := 
  !*FeasibleBasis*! (I:basis) of point_of_basis I \in polyhedron A b.
\end{lstlisting}

\subsubsection{Reduced Costs.} The simplex method stops when the current feasible basic point is an optimal solution of $\LP(A,b,c)$. This is determined thanks to the so-called \emph{reduced cost vector}. The reduced cost vector associated with the basis $I$ is defined as $u \coloneqq\transinv{A_I} c$, where $\transinv{A_I}$ denotes the inverse of the transpose matrix of $A_I$. On the \Coq{} side, assuming \C$c$ is a variable of type \C$'cV_n$, this leads to:
\begin{lstlisting}
Definition |*reduced_cost_of_basis*| (I:basis) :=
	(invmx (matrix_of_prebasis A I)^T) *m c : 'cV_n.
\end{lstlisting}
where \C$Q^T$ stands for the transpose of the matrix \C$Q$. When $u \geq 0$ and $I$ is feasible, the associated basic point is optimal:
\begin{lstlisting}
Lemma |*optimal_cert_on_basis*| (I:feasible_basis) :
  (reduced_cost_of_basis I) >=m 0 ->
    forall y, y \in polyhedron A b -> '[c, point_of_basis I] <= '[c,y].
\end{lstlisting}
Here, the notation \C$'[.,.]$ corresponds to the scalar product $\scalar{\cdot}{\cdot}$ (see the file \C$inner_product.v$ in \CoqPolyhedra{}). 

Strong duality lies in the core of the simplex method. To see this, consider the \emph{extended reduced cost vector} $\bar u \in \R^m$, which is defined by $\bar u_i \coloneqq u_i$ if $i \in I$, and $\bar u_i \coloneqq 0$ otherwise. On the \Coq{} side, this extended vector is built by the function \C$|*ext_reduced_cost_of_basis*|: basis -> 'cV_m$. When $u \geq 0$, $\bar u$ is a feasible point of the dual linear program \ref{eq:dual_lp}, and it has the same objective value as the basic point $x^I$:
\begin{lstlisting}
Lemma |*ext_reduced_cost_dual_feasible*| (I:basis) :
  let: u := reduced_cost_of_basis I in
  u >=m 0 = (ext_reduced_cost_of_basis I \in dual_polyhedron A c).
Lemma |*eq_primal_dual_value*| (I:basis) :
  '[c, point_of_basis I] = '[b, ext_reduced_cost_of_basis I].
\end{lstlisting}
As a consequence, proving the termination of the simplex method is one of the possible ways to establish the duality theorem of linear programming. 

\subsubsection{Pivoting.} \emph{Pivoting} refers to the operation of moving from a feasible basis to a ``better'' one, chosen according to what is known as the \emph{pivoting rule}.
More precisely, when the reduced cost vector $u$ associated with the current feasible basis $I$ does not satisfy $u \geq 0$, the pivoting rule selects an index $i \in I$ such that $u_i < 0$, which is called the \emph{leaving variable}, and builds the direction vector $d \coloneqq (A_I)^{-1} e_i$ (where $e_i$ is the $i$th vector of the canonical base of $\R^n$):
\begin{lstlisting}
Definition |*direction*| (I:basis) (i:'I_n) :=
	let: ei := (delta_mx i 0):'cV_n in 
	(invmx (matrix_of_prebasis A I)) *m ei.
\end{lstlisting}
along which the objective function $x \mapsto \scalar{c}{x}$ decreases:
\begin{lstlisting}
Lemma |*direction_improvement*| c (I:basis) (i:'I_n) :
	let: u := reduced_cost_of_basis c I in
	u i 0 < 0 -> '[c, direction I i] < 0.
\end{lstlisting}
As a consequence, the simplex method moves along the halfline $\{x^I + \lambda d \mid \lambda \geq 0\}$ in order to decrease the value of the objective function. When $d$ is a \emph{feasible direction}, \ie, $A d \geq 0$, this halfline is entirely contained in the polyhedron $\Pcal(A,b)$. In this case, we can easily show that the linear program~\ref{eq:lp} is unbounded:
\begin{lstlisting}
Definition |*feasible_dir*| A := [pred d | (A *m d) >=m 0]. 
Lemma |*unbounded_cert_on_basis*| (I:feasible_basis) (i:'I_n) :
	let: u := reduced_cost_of_basis c I in
	let: d := direction I i in
	feasible_dir A d -> u i 0 < 0 ->
		forall M, exists x, (x \in polyhedron A b) /\ ('[c,x] < M).
\end{lstlisting}
In contrast, if $d$ is not a feasible direction, moving along the halfline $\{x^I + \lambda d \mid \lambda \geq 0\}$ makes the simplex method eventually hit the boundary of one of the halfspaces $\{ x \in \R^n \mid A_j x \geq b_j \}$ delimiting $\Pcal(A,b)$. This happens precisely when $\lambda$ reaches the threshold value $\bar \lambda$ defined by:
\begin{equation}
\bar \lambda = \min_j \biggl\{ \frac{b_j - A_j x^I}{A_j d} \mid A_j d < 0 \biggr\} \, . \label{eq:mingap}
\end{equation}
The indexes attaining the minimum in Eq.~\eqref{eq:mingap} correspond to the halfspaces which are hit. Then, the pivoting rule selects one of them, say $j$, which is called the \emph{entering variable}, and the next basis is defined as $J \coloneqq (I \setminus \{i\}) \cup \{j\}$. In this way, it can be shown that $J$ is a feasible basis, and that $\scalar{c}{x^J}\leq \scalar{c}{x^I}$.

The major difficulty arising in this scheme is the possibility that $\bar \lambda = 0$, or, equivalently, that several bases correspond to the same basic point. Such bases are said to be \emph{degenerate}, and constitute the only obstacle to the termination of the simplex method. In the presence of degenerate bases, the pivoting rule needs to choose carefully the entering and leaving variables in order to avoid cycling over them. Our formalization of the simplex method is based on a rule having this property, called the {\em lexicographic rule}~\cite{Dantzig1955}, which is described in the next section. 

\section{Lexicographic Pivoting Rule}\label{sec:lex}

In informal terms, the lexicographic rule acts as if the vector $b$ was replaced by a perturbed vector $\tilde b$ defined by $\tilde b_i \coloneqq b_i - \eps^i$, where $\eps$ is a small positive parameter (here $\eps^i$ is the usual exponentiation). The advantage of perturbing $b$ in such a way is that there is no degenerate basis anymore. However, as we shall see, the feasible bases of the polyhedron $\Pcal(A, \tilde b)$ only form a subset of the feasible bases of $\Pcal(A, b)$. The former are called \emph{lex-feasible bases}, and they constitute the set of bases over which the simplex method with the lexicographic rule iterates.

In the formalization, which is carried out in \C$Section |*Lexicographic_rule*|$ of \C$simplex.v$, we have chosen to use a symbolic perturbation scheme in order to avoid dealing with numerical values for $\eps$.\footnote{Finding how small $\eps$ must be chosen is tedious, and this would make proofs unnecessarily complicated.} In this symbolic perturbation scheme, a row vector $v = (v_0, \dots, v_m) \in \R^{1 \times (1+m)}$ encodes the perturbed quantity $v_0 + \sum_{i = 1}^m v_i \eps^i$. The vector $\tilde b$ is then implemented as a row block matrix built from the column vector \C$b$ and the opposite of the identity matrix \C$-(1%:M)$ of order $m$:
\begin{lstlisting}
Definition |*b_pert*| := (row_mx b -(1%:M)):'M_(m,1+m).
\end{lstlisting}
In this way, the matrix \C$b_pert$ can be thought of as a column vector whose $i$th entry is the row vector $(b_i, 0, \dots, 0, -1, 0, \dots, 0)$, representing the quantity $b_i - \eps^i$, as desired.
Given a basis, the associated ``perturbed'' basic point is then: 
\begin{lstlisting}
Definition |*point_of_basis_pert*| (I:basis) :=
  (invmx (matrix_of_prebasis A I)) *m (matrix_of_prebasis b_pert I).
\end{lstlisting}
Now we can define the type of lex-feasible bases: 
\begin{lstlisting}
Definition |*is_lex_feasible*| (I:basis) :=   
  [forall i, ((row i A) *m (point_of_basis_pert I)) >=lex (row i b_pert)].
Inductive |*lex_feasible_basis*| := 
  !*LexFeasibleBasis*! (I:basis) of is_lex_feasible I.
\end{lstlisting}
where \C$>=lex$ is the lexicographic ordering over row vectors (see \C$vector_order.v$ in \CoqPolyhedra). 
We first observe that any lex-feasible basis is feasible:
\begin{lstlisting}
Lemma |*lex_feasible_basis_is_feasible*| (I:lex_feasible_basis): is_feasible I.
\end{lstlisting}

Following the description of the pivoting step in Sect.~\ref{sec:bases}, we now assume that the variables \C$I:lex_feasible_basis$ and \C$i:'I_n$ have been declared, and we make the following assumptions: 
\begin{lstlisting}
Hypothesis |*leaving*|: (reduced_cost_vector_of_basis I) i 0 < 0. 
Hypothesis |*infeas_dir*|: ~~(feasible_dir A (direction I i)).
\end{lstlisting}
where \C$~~b$ stands for the negation of the Boolean $b$. Our aim is to determine an entering variable $j$. In the symbolic perturbation scheme, every ratio appearing in Eq.~\eqref{eq:mingap} turns out to be a row vector encoding a perturbed quantity:
\begin{lstlisting}
Definition |*lex_gap*| (d:'cV_n) (j:'I_m) :=
  let: x_pert := point_of_basis_pert I in
	 ((A *m d) j 0)^-1 *: ((row j b_pert) - ((row j A) *m x_pert)) : 'rV_(1+m).
\end{lstlisting}
In order to obtain in the perturbed setting the analog of the threshold value $\bar \lambda$ defined in Eq.~\eqref{eq:mingap}, we determine the minimum of these ratios in the lexicographic sense, using the function \C$|*lex_min_seq*| S$ introduced in \C$vector_order.v$. The entering variable is then computed as follows:
\begin{lstlisting}
Definition |*lex_ent_var_nat*| :=
  let: d := direction I i in
  let: J := [seq j <- (enum 'I_m) | (A *m d) j 0 < 0] in
  let: min_gap := lex_min_seq [seq lex_gap d j | j <- J] in
    find (fun j => (j \in J) && (min_gap == lex_gap d j)) (enum 'I_m).
\end{lstlisting}
where the \MathCompShort{} function \C$find p S$ returns the index of the first item in the sequence \C$S$ for which the predicate \C$p$ holds, if any. Next, we prove that the result (of type \C$nat$) returned by \C$lex_ent_var_nat$ is strictly less than \C$m$, which allows us to convert it into an element of type \C$'I_m$ called \C$|*lex_ent_var*|$. We are finally ready to build the next basis:
\begin{lstlisting}
Definition |*lex_rule_set*| := lex_ent_var |: (I :\ (enum_val [...] i))).
\end{lstlisting}
where \C$k |: S$ and \C$S :\ k$ respectively adds and removes the element \C$k$ from the set \C$S$. With this definition, we show that the lexicographic rule provides a lex-feasible basis called \C$|*lex_rule_lex_bas*|$, by proving the following successive results:
\begin{lstlisting}
Lemma |*lex_rule_card*| : #|lex_rule_set| == n.
Lemma |*lex_rule_is_basis*| : is_basis (!*Prebasis*! lex_rule_card).
Lemma |*lex_rule_lex_feasibility*| : is_lex_feasible (!*Basis*! lex_rule_is_basis).
Definition |*lex_rule_lex_bas*| := !*LexFeasibleBasis*! lex_rule_lex_feasibility.
\end{lstlisting}
We finally prove that the analog of the objective function in the perturbed setting is \emph{strictly} decreasing in the lexicographic sense: 
\begin{lstlisting}
Lemma |*lex_rule_dec*| : let: J := lex_rule_lex_bas in
  (c^T *m point_of_basis_pert I) >lex (c^T *m point_of_basis_pert J).
\end{lstlisting}
As mentioned above, this comes from the fact that the analog of the threshold $\bar \lambda$ in this setting is nonzero, thanks to the absence of degenerate bases:
\begin{lstlisting}
Lemma |*eq_pert_point_imp_eq_bas*| (I I':basis) :
  point_of_basis_pert I = point_of_basis_pert I' -> I == I'.
\end{lstlisting}
Let us sketch the proof of this key result. Recall that \C$point_of_basis_pert I$ is a $n \times (1+m)$-matrix. Given \C$j:'I_m$, we can show that the \C$(1+j)$th column of this matrix is nonzero if, and only if, \C$j$ belongs to \C$I$ (we refer to \C$Lemma |*col_point_of_basis_pert*|$ in \C$simplex.v$). Indeed, since the matrix $A_I$ is invertible, the \C$(1+j)$th column of \C$point_of_basis_pert I$ is nonzero if, and only if, the \C$(1+j)$th column of \C$matrix_of_prebasis b_pert I$ is. By construction of \C$b_pert$, the latter column vector has only zero entries, except in the case where \C$j \in I$ (in this case, the entry corresponding to the index of \C$j$ in \C$I$ is $-1$).

\section{Phase II of the Simplex Method, and Farkas Lemma}\label{sec:PhaseII}

\subsubsection*{Phase II.} In this section, we present our formalization of \emph{Phase II} of the simplex method. We do it before the one of Phase I because as we will explain in Sect.~\ref{sec:Phase1}, Phase II is used in Phase I. Phase II of the simplex method determines the optimal value of the linear program~\ref{eq:lp}, supposing that an initial feasible basis \C$bas0:feasible_basis$ is known. \emph{De facto}, this makes the underlying assumption that the linear program is feasible. 

Our implementation of Phase II, which is developed in \C$Section |*Phase2*|$ of \C$simplex.v$, consists in iterating the function \C$lex_rule_lex_bas$ until finding an optimal basis (\ie~identifying that the associated reduced cost vector is nonnegative), or determining that the linear program is unbounded (\ie~identifying that the direction vector is feasible). Termination is expected to be guaranteed by $\C$Lemma lex_rule_dec$$ and the fact that the number of bases is finite. In addition, it looks reasonable to start the iteration of \C$lex_rule_lex_bas$ from the basis \C$bas0$. However, albeit feasible, the basis \C$bas0$ has no reason to be lex-feasible. Fortunately, it can be shown that, up to reordering the inequalities defining $\Pcal(A,b)$, we can make \C$bas0$ be lex-feasible. Instead of applying permutations on the rows of $A$ and $b$, we choose to apply the inverse permutation on the symbolic perturbation components of \C$b_pert$, and leave the initial problem~\ref{eq:lp} unchanged. As a consequence, we modify the previous definition of \C$b_pert$ as follows:
\begin{lstlisting}
Definition |*b_pert*| := (row_mx b (-(perm_mx s))).
\end{lstlisting}  
where \C$s:'S_m$ represents a permutation of the set $\{1, \ldots ,m\}$, and \C$perm_mx$ builds the corresponding permutation matrix (see the libraries \C$perm$ and \C$matrix$ of \MathCompShort{}). All the previous results remain valid under this change. The only difference is that they are now additionally parametrized by the permutation \C$s$, appearing as a global variable in \C$Section |*Lexicographic_rule*|$. For reason of space, we omit the description of the construction of the permutation \C$s0$ associated with \C$bas0$. We only mention that it satisfies the expected result:
\begin{lstlisting}
Lemma |*feasible_to_lex_feasible*| : is_lex_feasible s0 bas0.
\end{lstlisting}

The function performing one iteration of the Phase II algorithm with the lexicographic rule is built as follows:
\begin{lstlisting}
Definition |*basic_step*| (bas: lex_feasible_basis) :=
  let u := reduced_cost_of_basis c bas in
  if [pick i | u i 0 < 0] is Some i ((*(* picks i such that u i 0 < 0, if any *)*))
    then let d := direction bas i in
      if (@idPn (feasible_dir A d)) is ReflectT infeas_dir
    	then !*Lex_next_basis*! (lex_rule_lex_bas infeas_dir)
    	else !*Lex_final*! (!*Lex_res_unbounded*! (bas, i))
    else !*Lex_final*! (!*Lex_res_optimal_basis*! bas).
\end{lstlisting}
where \C$@idPn (feasible_dir A d)$ returns a proof \C$infeas_dir$ of the fact that the direction vector \C$d$ is not feasible, when this is the case. As a consequence, the function \C$basic_step$ returns either a next basis (constructor \C$!*Lex_next_basis*!$), or indicates that the method should stop (constructor \C$!*Lex_final*!$). 

The recursive function which iterates the function \C$basic_step$ is the following:
\begin{lstlisting}
Function |*lex_phase2*| bas {measure basis_height bas} :=
  match basic_step bas with
  | !*Lex_final*! final_res => final_res
  | !*Lex_next_basis*! bas' => lex_phase2 bas'
  end.
\end{lstlisting}
It is defined in the framework provided by the library \C$RecDef$ of \Coq{}, see~\cite{Barthe2006}. More precisely, its termination (and subsequently, the fact that \Coq{} accepts the definition) is established by identifying an integer quantity which is strictly decreased every time the function \C$basic_step$ returns a next basis:
\begin{lstlisting}
Definition |*basis_height*| bas :=  #| [ set bas':(lex_feasible_basis s0) | 
	(c^T *m (point_of_basis_pert s0 bas')) <lex 
			(c^T *m (point_of_basis_pert s0 bas)) ] |.
\end{lstlisting}
This quantity represents the number of lex-feasible bases for which the value of the ``perturbed'' objective function is (lexicographically) strictly less than the value of this function at the current lex-feasible basis. The fact that \C$basis_height$ decreases at every iteration is a consequence of \C$Lemma lex_rule_dec$. 

Gathering all these components, we finally arrive at the definition of the function implementing Phase~II:
\begin{lstlisting}
Definition |*phase2*| :=
  let: lex_bas0 := !*LexFeasibleBasis*! feasible_to_lex_feasible in
  lex_to_phase2_final_result ((@lex_phase2 s0) c lex_bas0).
\end{lstlisting} 
We present the correctness specification of this function by means of an adhoc inductive predicate. Such a presentation is idiomatic in the \MathComp{} library. The advantage is that it provides a convenient way to perform case analysis on the result of \C$phase2$. 
\begin{lstlisting}
Inductive |*phase2_spec*| : phase2_final_result -> Type :=
| !*Phase2_unbounded*! (p: feasible_basis * 'I_n) of 
  (reduced_cost_of_basis c p.1) p.2 0 < 0 /\ feasible_dir A 
    (direction p.1 p.2) : phase2_spec (!*Phase2_res_unbounded*! p)
| !*Phase2_optimal_basis*! (bas: feasible_basis) of 
  (reduced_cost_of_basis c bas) >=m 0 : 
    phase2_spec (!*Phase2_res_optimal_basis*! bas).?\smallskip?
Lemma |*phase2P*| : phase2_spec phase2.
\end{lstlisting} 
More precisely, \C$Lemma phase2P$ states that when the function \C$phase2$ returns a result of the form \C$!*Phase2_res_unbounded*! (bas, i)$, the pair \C$(bas, i)$ satisfies \C$(reduced_cost_of_basis c bas) i 0 < 0$ and \C$feasible_dir A (direction bas i)$. It precisely corresponds to the hypotheses of \C$Lemma unbounded_cert_on_basis$, and indicates that~\ref{eq:lp} is unbounded. Similarly, if the result of \C$phase2$ is of the form \C$!*Phase2_res_optimal_basis*! bas$, we have \C$(reduced_cost_of_basis c bas) >=m 0$, \ie, the basis \C$bas$ is an optimal basis (see \C$Lemma optimal_cert_on_basis$). 

\subsubsection*{Effective Definition of Feasibility, and Farkas Lemma.}

We can now formalize the notion of feasibility, \ie, the property that the polyhedron $\Pcal(A,b)$ is empty or not, as a Boolean predicate.\footnote{We make a slight abuse of language, since feasibility usually applies to linear programs. By extension, we apply it to polyhedra: $\Pcal(A,b)$ is \emph{feasible} if it is nonempty.}
We still assume that the variables \C$A$ and \C$b$ are declared. Following the discussion at the end of Sect.~\ref{sec:preliminaries}, the predicate is defined by means of the function \C$phase2$ executed on the dual problem $\DualLP(A, b, 0)$. To this end, we first build a feasible basis \C$|*dual_feasible_bas0*|$ for this problem, whose associated basic point is the vector $0 \in \R^m$. Feasibility of the polyhedron $\Pcal(A, b)$ is then defined as follows:
\begin{lstlisting}
Definition |*feasible*| :=
  if phase2 dual_feasible_bas0 (-b) is !*Phase2_res_optimal_basis*! _ then true else false.
\end{lstlisting}
Note that \C$-b$ corresponds to the objective function of the dual linear program (when written as a minimization problem). The correctness of our definition is established by showing that the predicate \C$feasible$ is equivalent to the existence of a point $x \in \Pcal(A, b)$. This is presented by means of Boolean reflection, using the \C$reflect$ relation of \MathCompShort{}:
\begin{lstlisting}
Lemma |*feasibleP*| : reflect (exists x, x \in polyhedron A b) feasible.
\end{lstlisting}
We point out that the feasibility certificate \C$x$ is constructed from the extended reduced cost vector of the optimal basis of $\DualLP(A, b, 0)$ returned by \C$phase2$. 

In a similar way, we prove the following characterization of the emptiness of~$\Pcal(A,b)$, which precisely corresponds to Farkas Lemma: 
\begin{lstlisting}
Definition |*dual_feasible_dir*| := [pred d | (A^T *m d == 0) && (d >=m 0)].
Lemma |*infeasibleP*| :
  reflect (exists d, dual_feasible_dir A d /\ '[b,d] > 0) (~~feasible). 
\end{lstlisting}
Indeed, \C$~~feasible$ amounts to the fact that \C$phase2$ returns an unboundedness certificate \C$!*Phase2_res_unbounded*! (bas,i)$ for $\DualLP(A, b, 0)$. The emptiness certificate \C$d$ of $\Pcal(A,b)$ is then obtained from the dual feasible direction \C$direction bas i$. 

\section{Complete Implementation of the Simplex Method}\label{sec:Phase1}

\subsubsection{The Pointed Case.} In order to obtain a full formalization of the simplex method, it remains to implement a \emph{Phase~I} algorithm. Its purpose is twofold: 
\begin{inparaenum}[(i)] 
\item determine whether the linear program~\ref{eq:lp} is feasible or not,
\item in the former case, return an initial feasible basis for Phase~II. 
\end{inparaenum}
There is one obstacle to the definition of such a Phase~I algorithm: even if a linear program is feasible, it may not have any feasible basis. For instance, consider the linear program over the variables $x_1, x_2$ which aims at minimizing $x_2$ subject to $-1 \leq x_2 \leq 1$. The feasible set is a cylinder around the $x_1$-axis, and it does not have any vertex, or, equivalently, basic point. A necessary and sufficient condition for the existence of a feasible basis is that the rank of $A$ is $n$. When this condition is fulfilled, the feasible set $\Pcal(A, b)$ is said to \emph{pointed}. We now describe the Phase~I algorithm under this assumption. This is developed in \C$Section |*Pointed_simplex*|$ of \C$simplex.v$.

From the hypothesis on the rank of $A$, we can extract an invertible square submatrix of $A$, which provides an initial basis \C$bas0$ of~\ref{eq:lp}. Beware that this basis is not necessarily a feasible one. As a consequence, we split the inequalities in the system $A x \geq b$ into two complementary groups, $A_\Pos x \leq b_\Pos$ and $A_\Neg x \leq b_\Neg$, where the $\Pos$ is the set of indexes $i \in \{1,\ldots ,m\}$ for which the basic point \C$point_of_basis bas0$ does not satisfy the inequality $A_i x \leq b_i$, and $\Neg \coloneqq \{1,\ldots ,m\} \setminus \Pos$. We denote by $p$ the cardinality of the set $K$. 
Phase~I is based on applying Phase~II algorithm to the following ``extended'' problem over the vector $z = (x,y) \in \R^{n+p}$:
\begin{equation}\tag{$\LP_\text{Phase~I}$} \label{eq:phaseI}
\begin{array}{r@{\quad}l}
\text{minimize} & \scalar{e}{y - A_\Pos x} \\[\jot]
\text{subject to} & A_\Pos x \leq b_\Pos + y \, , \; A_\Neg x \geq b_\Neg \, , \; y \geq 0 \, , \; (x,y) \in \R^{n+p} 
\end{array}
\end{equation}
where $e \in \R^p$ stands for the all-$1$-vector. The constraints defining~\ref{eq:phaseI} are gathered into a single system $A_\Ext z \geq b_\Ext$. Similarly, the objective function of~\ref{eq:phaseI} can be rewritten as a sole linear function $z = (x,y) \mapsto \scalar{c_\Ext}{z}$. 

The linear program~\ref{eq:phaseI} has two important properties. On the one hand, its optimal value can be bounded (from below) by the quantity $M_\Ext \coloneqq \scalar{e}{-b_\Pos}$: 
\begin{lstlisting}
Definition |*Mext*| := '[const_mx 1, - (row_submx b K)].
Lemma |*cext_min_value*| z : (z \in polyhedron Aext bext) -> '[cext, z] >= Mext.
\end{lstlisting}
On the other hand, the optimal value of~\ref{eq:phaseI} is equal to $M_\Ext$ if, and only if, the original problem~\ref{eq:lp} is feasible. The ``only if'' implication follows from the following lemma, which also provides a feasibility witness of~\ref{eq:lp}:
\begin{lstlisting}
Lemma |*feasible_cext_eq_min_active*| z :
  z \in polyhedron Aext bext -> '[cext,z] = Mext ->
    (usubmx z \in polyhedron A b).
\end{lstlisting}
where the \MathCompShort{} function \C$usubmx$ returns the upper subvector $x$ of a block vector of the form $z = \begin{psmallmatrix} x \\ y
\end{psmallmatrix}$. 
Regarding the ``if'' implication, an infeasibility certificate of~\ref{eq:lp} can be constructed by means of a feasible point $\bar u \in \R^{m+p}$ of the dual of~\ref{eq:phaseI} whose objective value $\scalar{b_\Ext}{\bar u}$ is strictly greater than $M_\Ext$. This certificate is built by the following function:
\begin{lstlisting}
Definition |*dual_from_ext*| (u:'cV[R]_(m+p)) :=
  \col_i (if i \in K then 1 - (usubmx u) i 0 else (usubmx u) i 0).
\end{lstlisting}
where \C$\col_i Expr(i)$ is the column vector whose $i$th entry is \C$Expr(i)$. As expected, this certificate satisfies:
\begin{lstlisting}
Lemma |*dual_polyhedron_from_ext*| u :
  (u \in dual_polyhedron Aext cext) -> dual_feasible_dir A (dual_from_ext u).
Lemma |*dual_from_ext_obj*| u : '[bext, u] > Mext -> '[b, dual_from_ext u] > 0.
\end{lstlisting}
In this way, we readily obtain a proof that~\ref{eq:lp} is infeasible, by using \C$Lemma infeasibleP$.

Finally, we can build an initial feasible basis \C$|*feasible_bas0_ext*|$ for~\ref{eq:phaseI} by considering the union of the basis \C$bas0$ with the set $\{m+1, \dots, m+p\}$ of the indexes of the last $p$ constraints $y \geq 0$ of~\ref{eq:phaseI}.\footnote{We let the reader check that the associated basic point is $\begin{psmallmatrix} x \\ 0 \end{psmallmatrix} \in \R^{n+p}$, where $x$ is the basic point associated with the basis \C$bas0$, and that this point is feasible.}
As a consequence, we can apply \C$phase2$ to solve~\ref{eq:phaseI}, starting from the basis \C$feasible_bas0_ext$. 
In this way, we obtain an optimal basis \C$bas$ of~\ref{eq:phaseI}. If the associated basic point \C$z$ satisfies \C$'[cext,z] > Mext$, we build an infeasibility certificate of~\ref{eq:lp} using the function \C$dual_from_ext$, as described above. Otherwise, we construct a feasible basis \C$bas'$ of~\ref{eq:lp}. This is performed by the function \C$|*extract_feasible_basis*|$ which we do not describe here for the sake of concision. Then, we use \C$bas'$ to execute \C$phase2$ on~\ref{eq:lp} and finally obtain its optimal value. 

The previous discussion precisely describes the way we have implemented the function \C$|*pointed_simplex*|$, which completely solves the linear program~\ref{eq:lp} under the pointedness assumption.

\subsubsection{The General Case.} In general, we can always reduce to the pointed case by showing that \ref{eq:lp} is equivalent to the following linear program in which the original variable $x \in \R^n$ is substituted by $v-w$ with $v,w \geq 0$:
\begin{equation}\tag{$\LP_\text{Pointed}$}\label{eq:lpEquiv}
\begin{array}{r@{\quad}l}
\text{minimize} & \scalar{c}{(v-w)}  \\[\jot]
\text{subject to} & A (v - w) \geq b \, , \; v \geq 0 \, , \; w \geq 0 \, , \;(v,w) \in \R^{n+n}
\end{array}
\end{equation}
The feasible set of~\ref{eq:lpEquiv} is pointed because of the constraints $v,w \geq 0$. Thus, we can apply to it the function \C$pointed_simplex$ of the previous section. In this way, we define the function \C$|*simplex*|$, which is able to solve any linear program~\ref{eq:lp}. It is implemented in \C$Section |*General_simplex*|$ of \C$simplex.v$. Its correctness proof is formalized by means of the following inductive type:
\begin{lstlisting}
Inductive |*simplex_spec*| : simplex_final_result -> Type :=
| !*Infeasible*! d of (dual_feasible_dir A d /\ '[b, d] > 0): 
		simplex_spec (!*Simplex_infeasible*! d)
| !*Unbounded*! p of [/\ (p.1 \in polyhedron A b), (feasible_dir A p.2) & 
		('[c,p.2] < 0)] : simplex_spec (!*Simplex_unbounded*! p)
| !*Optimal_point*! p of [/\ (p.1 \in polyhedron A b), 
		(p.2 \in dual_polyhedron A c) & '[c,p.1] = '[b, p.2]] : 
		simplex_spec (!*Simplex_optimal_point*! p).?\smallskip?
Lemma |*simplexP*|: simplex_spec simplex.
\end{lstlisting}
In other words, when \C$simplex$ returns a result of the form \C$!*Simplex_infeasible*! d$, then \C$d$ is a certificate of infeasibility of $\LP(A,b,c)$, see \C$Lemma infeasibleP$. 

Similarly, the unboundedness of the linear program $\LP(A,b,c)$ is characterized by the fact that \C$simplex$ returns a result of the form \C$!*Simplex_unbounded*! (x,d)$. Equivalently, we can define a predicate corresponding to this situation, and prove that it is correct, as follows:
\begin{lstlisting}
Definition |*unbounded*| := 
	if simplex is (Simplex_unbounded _) then true else false.
Lemma |*unboundedP*| : 
  reflect (forall M, exists y, y \in polyhedron A b /\ '[c,y] < M) unbounded.
\end{lstlisting}
Given any \C$M$, the certificate \C$y$ is built by taking a point of the form $\C$x$ + \lambda \C$d$$, where $\lambda \geq 0$ is sufficiently large.

Finally, when \C$simplex$ returns \C$!*Simplex_optimal_point*! (x,u)$, this means that \C$x$ is an optimal solution of $\LP(A,b,c)$, and \C$u$ is a dual feasible element which certificates its optimality (\ie, $\scalar{c}{x} = \scalar{b}{u}$). Thanks to this, we can define in an effective way the fact that $\LP(A,b,c)$ admits an optimal solution (we say that the linear program is \emph{bounded}), and, in this case, deal with the optimal value:
\begin{lstlisting}
Definition |*bounded*| :=
  if simplex is (Simplex_optimal_point _) then true else false.
Definition |*opt_value*| :=
  if simplex is (Simplex_optimal_point (x,_)) then '[c,x] else 0.
Lemma |*boundedP*| :
  reflect ((exists x, x \in polyhedron A b /\ '[c,x] = opt_value) /\ (forall y, y \in polyhedron A b -> opt_value <= '[c,y])) bounded.
\end{lstlisting}

\section{Outcome of the Effective Approach}\label{sec:outcome}

Duality results immediately follow from the correctness statements  of the simplex method and the resulting predicates \C$feasible$, \C$unbounded$ and \C$bounded$. For instance, when~\ref{eq:lp} and~\ref{eq:dual_lp} are both feasible, we have:
\begin{lstlisting}
Theorem |*strong_duality*| : feasible A b -> dual_feasible A c -> 
  exists p, [/\ p.1 \in polyhedron A b, p.2 \in dual_polyhedron A c 
              & '[c,p.1] = '[b,p.2]].
\end{lstlisting}
which corresponds to the second part of Th.~\ref{th:strong_duality}. The remaining cases of Th.~\ref{th:strong_duality} (when one of the two linear programs is infeasible) are dealt with in the file \C$duality.v$. All these statements are obtained in a few lines of proof. We also obtain another well-known form of Farkas Lemma, characterizing the logical implication between linear inequalities (\C$Lemma |*farkas_lemma_on_inequalities*|$).

The membership to the convex hull of a finite set of points is another  property which can be defined in an effective way in our framework. Recall that a point $x \in \R^n$ belongs to the convex hull of a (finite) set $V = \{v^i\}_{1 \leq i \leq p} \subset \R^n$ if there exists $\lambda \in \R^p$ such that $x = \sum_{i = 1}^p \lambda_i v^i$, $\lambda \geq 0$ and $\sum_i \lambda_i = 1$. The latter constraints define a polyhedron over $\lambda \in \R^p$, and the membership of $x$ amounts to fact that this polyhedron is feasible. This is how we arrive at the definition of a Boolean predicate \C$|*is_in_convex_hull*|$, see the file \C$minkowski.v$. The  \emph{separation result} states that if $x$ does not belong to the convex hull of $V$, then there is a hyperplane \emph{separating} $x$ from $V$. This means that $x$ is located on one side of the hyperplane, while the points of $V$ are on the other side. Formalizing $V$ as the matrix of size $n \times p$ with columns $v^i$, we establish this result as follows:
\begin{lstlisting}
Theorem |*separation*| (x: 'cV_n) : 
  ~~ (is_in_convex_hull x) -> exists c, [forall i, '[c, col i V] > '[c, x]].
\end{lstlisting}
The certificate \C$c$ can be built directly from the infeasibility certificate of the underlying polyhedron over $\lambda \in \R^p$. Our proof of the separation result reduces to the technical manipulations of block matrices performing this conversion.

Finally, Minkowski Theorem states that every bounded polyhedron equals the convex hull of its vertices. We recover this result as the extensional equality of the predicates \C$polyhedron A b$ and \C$is_in_convex_hull matrix_of_vertices$, where \C$|*matrix_of_vertices*|$ is the matrix whose columns are the basic points of $\Pcal(A,b)$:
\begin{lstlisting}
Theorem |*minkowski*| : bounded_polyhedron A b ->
  polyhedron A b =i is_in_convex_hull matrix_of_vertices.
\end{lstlisting}
The most difficult part of the statement is proven in a few lines: if $x \in \Pcal(A,b)$ does not belong to the convex hull of the basic points, \C$Lemma separation$ exhibits a separating hyperplane $c$ such that $\scalar{c}{x} < \scalar{c}{x^I}$ for all feasible bases $I$ of $\Pcal(A,b)$. However, the program \C$pointed_simplex$ is able to provide an optimal feasible basis $I^*$, \ie, which satisfies $\scalar{c}{x^{I^*}} \leq \scalar{c}{x}$. This yields a contradiction.

\section{Conclusion}

We have presented a formalization of convex polyhedra in~\Coq{}. Its main feature is that it is based on an implementation of the simplex method, leading to an effective formalization of the basic predicates over polyhedra. We have illustrated the outcome of this approach with several results of the theory of convex polyhedra. As a future work, we plan to deal with faces, which are a central notion in the combinatorial theory of polyhedra (the early steps of an effective definition of faces are already available in the file \C$face.v$ of \CoqPolyhedra{}). The simplex method should also greatly help us to prove adjacency properties on faces, in particular, properties related with the connectivity of the (vertex-edge) graph of polyhedra. Another direction of work is to exploit our library to certify computational results on polyhedra, possibly on large-scale instances. A basic problem is to formally check that a certain polyhedron (defined by inequalities) is precisely the convex hull of a certain set of points. This is again a problem in which the simplex method plays an important role~\cite{Avis1992}. To cope with the computational aspects, we plan to investigate how to  translate our formally proven statements to lower-level data structures, like in~\cite{Cohen2013}. 

\vspace*{-1ex}
\subsubsection*{Acknowledgments.}
The authors are very grateful to A.~Mahboubi for her help to improve the presentation of this paper, and to G.~Gonthier, F.~Hivert and P.-Y.~Strub  for fruitful discussions. The second author is also grateful to M.~Cristi\'a for introducing him to the topic of automated theorem proving. The authors finally thank the anonymous reviewers for their suggestions and remarks.

\bibliographystyle{splncs}
\begin{thebibliography}{10}
\providecommand{\url}[1]{\texttt{#1}}
\providecommand{\urlprefix}{URL }

\bibitem{Avis1992}
Avis, D., Fukuda, K.: A pivoting algorithm for convex hulls and vertex
  enumeration of arrangements and polyhedra. Discrete {\&} Computational
  Geometry  8(3) (1992)

\bibitem{Barthe2006}
Barthe, G., Forest, J., Pichardie, D., Rusu, V.: Defining and reasoning about
  recursive functions: A practical tool for the {C}oq proof assistant. In:
  Hagiya, M., Wadler, P. (eds.) Proceedings of FLOPS 2006. Springer Berlin
  Heidelberg (2006)

\bibitem{Besson2007}
Besson, F.: Fast reflexive arithmetic tactics the linear case and beyond. In:
  Altenkirch, T., McBride, C. (eds.) Proceedings of TYPES 2006. Springer (2007)

\bibitem{Bremner2013}
Bremner, D., Deza, A., Hua, W., Schewe, L.: More bounds on the diameters of
  convex polytopes. Optimization Methods and Software  28(3) (2013)

\bibitem{Cohen2013}
Cohen, C., D{\'e}n{\`e}s, M., M{\"o}rtberg, A.: Refinements for free! In:
  Gonthier, G., Norrish, M. (eds.) Proceedongs of CPP 2013. Springer (2013)

\bibitem{CousotHalbwachs78}
Cousot, P., Halbwachs, N.: Automatic discovery of linear restraints among
  variables of a program. In: Proceedings of POPL 1978. ACM Press, Tucson,
  Arizona (1978)

\bibitem{Dantzig51}
Dantzig, G.B.: {Maximization of a Linear Function of Variables Subject to
  Linear Inequalities, in Activity Analysis of Production and Allocation}.
  Wiley (1951)

\bibitem{Dantzig1955}
Dantzig, G.B., Orden, A., Wolfe, P.: The generalized simplex method for
  minimizing a linear form under linear inequality restraints. Pacific J. Math.
   5(2) (1955)

\bibitem{Fouilhe2014}
Fouilhe, A., Boulm{\'e}, S.: A certifying frontend for (sub)polyhedral abstract
  domains. In: Giannakopoulou, D., Kroening, D. (eds.) Proceedings of VSTTE
  2014. Springer (2014)

\bibitem{Gonthier2011}
Gonthier, G.: Point-free, set-free concrete linear algebra. In: van Eekelen,
  M., Geuvers, H., Schmaltz, J., Wiedijk, F. (eds.) Proceedings of ITP 2011.
  Springer (2011)

\bibitem{MathComp}
Gonthier, G., Mahboubi, A., Tassi, E.: {A Small Scale Reflection Extension for
  the Coq system}. Research Report RR-6455, {Inria Saclay Ile de France} (2016)

\bibitem{Guglielmi2017}
Guglielmi, N., Laglia, L., Protasov, V.: Polytope {L}yapunov functions for
  stable and for stabilizable {LSS}. Foundations of Computational Mathematics
  17(2) (2017)

\bibitem{Harrison2013}
Harrison, J.: The {HOL} {L}ight theory of {E}uclidean space. Journal of
  Automated Reasoning  50 (2013)

\bibitem{Khachiyan80}
Khachiyan, L.: Polynomial algorithms in linear programming. USSR Computational
  Mathematics and Mathematical Physics  20(1) (1980)

\bibitem{Sakaguchi2016}
Sakaguchi, K.: Vass. \url{https://github.com/pi8027/vass} (2016)

\bibitem{Santos2012}
{Santos}, F.: {A counterexample to the Hirsch conjecture.} {Ann. Math.}  176(1)
  (2012)

\bibitem{Smale98}
Smale, S.: Mathematical problems for the next century. Mathematical
  Intelligencer  20 (1998)

\bibitem{Spasic2012}
Spasi{\'{c}}, M., Mari{\'{c}}, F.: Formalization of incremental simplex
  algorithm by stepwise refinement. In: Giannakopoulou, D., M{\'e}ry, D. (eds.)
  Proceedings of FM 2012. Springer (2012)

\end{thebibliography}


\newpage

\end{document}


1706.10269.extracted_bib
{Avis1992 Avis, D., Fukuda, K.: A pivoting algorithm for convex hulls and vertex enumeration of arrangements and polyhedra. Discrete {\& Computational Geometry 8(3) (1992)
1706.10269.extracted_bib
{Barthe2006 Barthe, G., Forest, J., Pichardie, D., Rusu, V.: Defining and reasoning about recursive functions: A practical tool for the {Coq proof assistant. In: Hagiya, M., Wadler, P. (eds.) Proceedings of FLOPS 2006. Springer Berlin Heidelberg (2006)
1706.10269.extracted_bib
{Besson2007 Besson, F.: Fast reflexive arithmetic tactics the linear case and beyond. In: Altenkirch, T., McBride, C. (eds.) Proceedings of TYPES 2006. Springer (2007)
1706.10269.extracted_bib
{Bremner2013 Bremner, D., Deza, A., Hua, W., Schewe, L.: More bounds on the diameters of convex polytopes. Optimization Methods and Software 28(3) (2013)
1706.10269.extracted_bib
{Cohen2013 Cohen, C., D{\'en{\`es, M., M{\"ortberg, A.: Refinements for free! In: Gonthier, G., Norrish, M. (eds.) Proceedongs of CPP 2013. Springer (2013)
1706.10269.extracted_bib
{CousotHalbwachs78 Cousot, P., Halbwachs, N.: Automatic discovery of linear restraints among variables of a program. In: Proceedings of POPL 1978. ACM Press, Tucson, Arizona (1978)
1706.10269.extracted_bib
{Dantzig51 Dantzig, G.B.: {Maximization of a Linear Function of Variables Subject to Linear Inequalities, in Activity Analysis of Production and Allocation. Wiley (1951)
1706.10269.extracted_bib
{Dantzig1955 Dantzig, G.B., Orden, A., Wolfe, P.: The generalized simplex method for minimizing a linear form under linear inequality restraints. Pacific J. Math. 5(2) (1955)
1706.10269.extracted_bib
{Fouilhe2014 Fouilhe, A., Boulm{\'e, S.: A certifying frontend for (sub)polyhedral abstract domains. In: Giannakopoulou, D., Kroening, D. (eds.) Proceedings of VSTTE 2014. Springer (2014)
1706.10269.extracted_bib
{Gonthier2011 Gonthier, G.: Point-free, set-free concrete linear algebra. In: van Eekelen, M., Geuvers, H., Schmaltz, J., Wiedijk, F. (eds.) Proceedings of ITP 2011. Springer (2011)
1706.10269.extracted_bib
{MathComp Gonthier, G., Mahboubi, A., Tassi, E.: {A Small Scale Reflection Extension for the Coq system. Research Report RR-6455, {Inria Saclay Ile de France (2016)
1706.10269.extracted_bib
{Guglielmi2017 Guglielmi, N., Laglia, L., Protasov, V.: Polytope {Lyapunov functions for stable and for stabilizable {LSS. Foundations of Computational Mathematics 17(2) (2017)
1706.10269.extracted_bib
{Harrison2013 Harrison, J.: The {HOL {Light theory of {Euclidean space. Journal of Automated Reasoning 50 (2013)
1706.10269.extracted_bib
{Khachiyan80 Khachiyan, L.: Polynomial algorithms in linear programming. USSR Computational Mathematics and Mathematical Physics 20(1) (1980)
1706.10269.extracted_bib
{Sakaguchi2016 Sakaguchi, K.: Vass. \url{https://github.com/pi8027/vass (2016)
1706.10269.extracted_bib
{Santos2012 {Santos, F.: {A counterexample to the Hirsch conjecture. {Ann. Math. 176(1) (2012)
1706.10269.extracted_bib
{Smale98 Smale, S.: Mathematical problems for the next century. Mathematical Intelligencer 20 (1998)
1706.10269.extracted_bib
{Spasic2012 Spasi{\'{c, M., Mari{\'{c, F.: Formalization of incremental simplex algorithm by stepwise refinement. In: Giannakopoulou, D., M{\'ery, D. (eds.) Proceedings of FM 2012. Springer (2012)