A New Look at the Automatic Synthesis of Linear Ranking Functions^{}

Roberto Bagnara^{a,b}, Fred Mesnard^{c}, Andrea Pescetti^{d}, Enea Zaffanella^{a,b}

aDipartimento di Matematica, Università di Parma, Italy

                        bBUGSENG srl, http://bugseng.com

cIREMIA, LIM, Université de la Réunion, France dUniversità di Parma, Italy

Abstract

The classical technique for proving termination of a generic sequential computer program involves the synthesis of a ranking function for each loop of the program. Linear ranking functions are particularly interesting because many terminating loops admit one and algorithms exist to automatically synthesize it. In this paper we present two such algorithms: one based on work dated 1991 by Sohn and Van Gelder; the other, due to Podelski and Rybalchenko, dated 2004. Remarkably, while the two algorithms will synthesize a linear ranking function under exactly the same set of conditions, the former is mostly unknown to the community of termination analysis and its general applicability has never been put forward before the present paper. In this paper we thoroughly justify both algorithms, we prove their correctness, we compare their worst-case complexity and experimentally evaluate their efficiency, and we present an open-source implementation of them that will make it very easy to include termination-analysis capabilities in automatic program verifiers.

Keywords: Static analysis, computer-aided verification, termination analysis.

1. Introduction

Termination analysis of computer programs (a term that here we interpret in its broadest sense) consists in attempting to determine whether execution of a given program will definitely terminate for a class of its possible inputs. The ability to anticipate the termination behavior of programs (or fragments thereof) is essential to turn assertions of partial correctness (if the program reaches a certain control point, then its state satisfies some requirements) into assertions of total correctness (the program will reach that point and its state will satisfy those requirements). It is worth observing that the property of termination of a program fragment is not less important than, say, properties concerning the absence of run-time errors. For instance, critical reactive systems (such as fly-by-wire avionics systems) must maintain a continuous interaction with the environment: failure to terminate of some program components can stop the interaction the same way as if an unexpected, unrecoverable run-time error occurred.

This work has been partly supported by PRIN project “AIDA—Abstract Interpretation:

Design and Applications” and by the Faculty of Sciences of the University of Parma.

Email addresses: bagnara@cs.unipr.it (Roberto Bagnara), Frederic.Mesnard@univ-reunion.fr (Fred Mesnard), andrea.pescetti@unipr.it (Andrea Pescetti), zaffanella@cs.unipr.it (Enea Zaffanella)

Preprint submitted to Information and Computation

Developing termination proofs by hand is, as any other program verification task, tedious, error-prone and, to keep it short, virtually impossible to conduct reliably on programs longer than a few dozens of lines. For this reason, automated termination analysis has been a hot research topic for more than two decades. Of course, due to well-known limitative results of computation theory, any automatic termination analysis can only be expected to give the correct answer (“the program does —or does not— terminate on these inputs”) for some of the analyzed programs and inputs: for the other programs and inputs the analysis will be inconclusive (“don’t know”). It is worth noticing that there is no need to resort to the halting problem to see how hard proving termination can be. A classical example is the 3x + 1 problem,^{1} whose termination for any n has been a conjecture for more than 70 years:

while n > 1 do

if (n mod 2) ≠ 0 then n := 3n + 1 else n := n div 2

The classical technique for proving termination of a generic sequential computer program consists in selecting, for each loop w of the program:

1. a set S_{w} that is well-founded with respect to a relation R_{w} ⊆ S_{w} × S_{w}; namely, for each U ⊆ S_{w} such that U ≠ ∅, there exists v ∈ U such that (u, v) ∉ R_{w} for each u ∈ U;

2. a function f_{w} from the set of program states that are relevant for w (e.g., those concerning the head of the loop and that are reachable from a designated set of initial states) to the set S_{w}, such that the values of f_{w} computed at any two subsequent iterations of w are in relation R_{w}.

The function f_{w} is called ranking function, since it ranks program states according to their “proximity” to the final states. Let us focus on deterministic programs, and consider a loop w and a set of initial states Σ_{w}^{I} for w. Assume further that the body of w always terminates when w is initiated in a state σ ∈ Σ_{w}^{I} and that Σ_{w}^{F} is a set of final states for w, that is, w immediately terminates when it is initiated in a state σ ∈ Σ_{w}^{F}. If we fix any enumeration of Σ_{w}^{I} = {σ_{0}^{0}, σ_{1}^{0}, …}, then the computations of w we are interested in can be represented by the (possibly infinite) sequence of (possibly infinite) sequences

1Also known as the Collatz problem, the Syracuse problem, Kakutani’s problem, Hasse’s algorithm, and Ulam’s problem: see, e.g., [1].

\begin{matrix}σ_{0}^{0} & σ_{0}^{1} & … \\ ⋮ & ⋮ & ... \\ σ_{i}^{0} & σ_{i}^{1} & … \\ ⋮ & ⋮ & ... \end{matrix}  (1)

Let Σ_{w} be the set of all states that occur in (1). Suppose that we succeed in finding a ranking function f_{w}: Σ_{w} → S_{w}, where S_{w} is well-founded with respect to R_{w} and, for each m, n ∈ ℕ, if σ_{m}^{n} and σ_{m}^{n+1} occur in (1), then (f_{w}(σ_{m}^{n+1}), f_{w}(σ_{m}^{n})) ∈ R_{w}. In this case we know that all the sequences in (1), and hence all the computations they represent, are finite.

Example 1.1. Consider the following loop, where x takes values in ℤ:

\begin{matrix}while x ≠ 0 do \\ x := x − 1 \end{matrix}

Here the state at the loop head can be simply characterized by an integer number: the value of x. If we take Σ^{I} := ℕ then the computation sequences of interest are

0

1 0
⋮  ⋮  ...
n n − 1 … 0
⋮  ⋮  ...

We thus have Σ = ℕ and Σ^{F} = {0}. If we define S := ℕ, f as the identity function over ℕ, and R := { (h, k) | h, k ∈ ℕ, h < k }, then S is well founded with respect to R and f is a ranking function (with respect to Σ, S and R).

We have seen that, if there exists a ranking function, then all computations summarized by (1) terminate. What is interesting is that the argument works also the other way around: if all the computations summarized by (1) do terminate, then there exists a ranking function (actually, there exists an infinite number of them). In fact, suppose all the sequences in (1) are finite. Since the program is deterministic, any state occurs only once in every sequence. Moreover, if a state σ occurs in more than one sequence, then the suffixes of these sequences that immediately follow σ are all identical (since the future of any computation is completely determined by its current state). The function mapping each σ in Σ_{w} to the natural number representing the length of such suffixes is thus well defined and is a ranking function with respect to Σ_{w} and ℕ with the well-founded ordering given by the ‘<’ relation.

It is worth observing that the above argument implies that if any ranking function exists, then there exists a ranking function over (ℕ, <). This observation can be generalized to programs having bounded nondeterminism [2]: therefore, ranking functions on the naturals are sufficient, for instance, when modeling the input of values for commonly available built-in data types. However, the use of more general well-founded orderings can simplify the search for a ranking function (see, e.g., [3, Example 1.2]). Moreover, such a generalization is mandatory when dealing with unbounded nondeterminism [2] (see also [4, Section 10]).

The termination of a set of computations and the existence of a ranking function for such a set are thus completely equivalent. On the one hand, this means that trying to prove that a ranking function exists is, at least in principle, not less powerful than any other method we may use to prove termination. On the other hand, undecidability of the termination problem implies that the existence of a ranking function is also undecidable. An obvious way to prove the existence of a ranking function is to synthesize one from the program text and a description of the initial states: because of undecidability, there exists no algorithm that can do that in general.

The use of ranking functions as a tool to reason about termination can be traced back to the seminal work of R. W. Floyd in [5], where they are introduced under the name of W -functions. Since then, several variations of the method have been proposed so as to extend its applicability from the realm of classical sequential programs to more general constructs (e.g., concurrency). In particular, in [4], seven different ‘à la Floyd’ induction principles for nondeterministic transition systems are formally shown to be sound, semantically complete and equivalent. For instance, it is shown that it is sufficient to consider a single, global ranking function, instead of a different ranking function for each program control point, as originally proposed in [5]; and that the decrease of such a global ranking function need not be verified at all program control points, but it is enough to consider a minimal set of loop cut-points; moreover, when trying to prove properties that only depend on the current state of the system (e.g., termination of a deterministic program), it is always possible to find a ranking function depending on the current state only, i.e., independent of the initial state of the system. Note that these results have been implicitly exploited in the examples above so as to simplify the presentation of the method.

In this paper we present, in very general terms so as to encompass any programming paradigm, the approach to termination analysis based on the explicit search of ranking functions. We then restrict attention to linear ranking functions obtained from linear approximations of the program’s semantics. For this restriction, we present and fully justify two methods to prove the existence of linear ranking functions: one, based on work dated 1991 by Sohn and Van Gelder, that is almost unknown outside the field of logic programming even though, as we demonstrate in the present paper, it is completely general; the other, due to Podelski and Rybalchenko, dated 2004, was proved correct by the authors but the reasons why it works were never presented. We then provide a proof of equivalence of the two methods, thus providing an independent assessment of their correctness and relative completeness. We also compare their theoretical complexity and practical efficiency on three related problems:

1. proving that one linear ranking function exists;

2. exhibiting one such function;

3. computing the space of all linear ranking functions.

The experimental evaluation is based on the implementation of the two methods provided by the Parma Polyhedra Library [6], a free software library of numerical abstraction targeted at software/hardware analysis and verification. These implementations are, to the best of our knowledge, the first ones that are being made available, in source form, to the community. For this reason, the implementations should be regarded as complementary to the present paper in the common aim of making the automatic synthesis of linear ranking functions known outside programming language barriers, understandable and accessible.

The plan of the paper is as follows: Section 2 recalls preliminary notions and introduces the notation used throughout the paper; Section 3 introduces the problem of automatic termination analysis of individual loops and its solution technique based on the synthesis of ranking functions; Section 4 presents a simple generalization of the approach of [7] that is generally applicable to termination analysis of any language; Section 5 shows and fully justifies the approach of [8]; Section 6 proves the two methods are equivalent and compares them from the point of view of computational complexity; Section 7 presents the implementation of the two approaches offered by the Parma Polyhedra Library and the corresponding experimental evaluation, providing a comparison of their practical efficiency; Section 8 concludes.

Readers who are interested in a more detailed exposition of the themes treated in this paper are referred to its technical report version [3].

2. Preliminaries

Set Theory. The set of all finite sequences of elements of S is denoted by S^{∗}. The empty sequence is denoted by ε and the length of a sequence w is denoted by |w|. The set of non-negative integers, rationals and reals are denoted by ℕ, ℚ_{+} and ℝ_{+}, respectively.

Linear Algebra. For each i ∈ {1, …, n}, v_{i} denotes the i-th component of the real (column) vector v = ⟨v_{1}, …, v_{n}⟩ ∈ ℝ^{n}. A vector v ∈ ℝ^{n} can also be interpreted as a matrix in ℝ^{n×1} and manipulated accordingly with the usual definitions for addition, multiplication (both by a scalar and by another matrix), and transposition, which is denoted by v^{T}, so that ⟨v_{1}, …, v_{n}⟩ = (v_{1}, …, v_{n})^{T}. If v ∈ ℝ^{n} and w ∈ ℝ^{m}, we will write ⟨v, w⟩ to denote the column vector in ℝ^{n+m} obtained by “concatenating” v and w, so that ⟨v, w⟩ = ⟨v_{1}, …, v_{n}, w_{1}, …, w_{m}⟩. The scalar product of v, w ∈ ℝ^{n} is the real number v^{T}w = ∑_{i=1}^{n} v_{i}w_{i}. The identity matrix in ℝ^{n×n} is denoted by I_{n}. We write 0 to denote a matrix in ℝ^{n×m} having all of its components equal to zero; the dimensions n and m will be clear from context. We sometimes treat scalars as vectors in ℝ^{1} or matrices in ℝ^{1×1}. For any relational operator ▷◁ ∈ {<, ≤, =, ≥, >}, we write v ▷◁ w to denote the conjunctive proposition ⋀_{i=1}^{n}(v_{i} ▷◁ w_{i}). Moreover, v ≠ w will denote the proposition ¬(v = w). We will sometimes use the convenient notation a ▷◁_{1} b ▷◁_{2} c to denote the conjunction a ▷◁_{1} b ∧ b ▷◁_{2} c and we will not distinguish conjunctions of propositions from sets of propositions. The same notation applies to vectors defined over other numeric fields and, for the supported operations, to vectors defined over numeric sets such as ℕ and ℚ_{+}.

First-Order Logic. Let L be a first-order language with variables in X. We will routinely confuse a tuple of variables with the set of its components. So, if ϕ is an L-formula, we will write ϕ[x̄] to denote ϕ itself, yet emphasizing that the set of free variables in ϕ is included in x̄. Let x̄, ȳ ∈ X^{∗} be of the same length and let ϕ be a L-formula: then ϕ[ȳ/x̄] denotes the formula obtained by simultaneous renaming of each free occurrence in ϕ of a variable in x̄ with the corresponding variable in ȳ, possibly renaming bound variable occurrences as needed to avoid variable capture. A formula with no free variable occurrences is termed closed or called a sentence. The universal closure of a formula ϕ is denoted by ∀(ϕ). If ϕ is a closed L-formula and A is an L-structure, we write A ⊧ ϕ if ϕ is satisfied in A. A set T of closed L-formulas is called an L-theory. We write A ⊧ T if A ⊧ ϕ for each ϕ ∈ T . If ϕ is a closed L-formula and T is an L-theory, we write T ⊧ ϕ if, for each L-structure A, A ⊧ T implies A ⊧ ϕ. In this case we say that ϕ is a logical consequence of T .

3. Termination Analysis of Individual Loops

We will start by restricting our attention to individual loops of the form

{ I } while B do C  (2)

where

• I is a loop invariant that a previous analysis phase has determined to hold just before any evaluation of B;

• B is a Boolean guard expressing the condition on the state upon which iteration continues;

• C is a command that, in the context set by (2), is known to always terminate.

Notice that, for maximum generality, we do not impose any syntactic restriction on I, B and C and will only observe their interaction with the program state: I and B express conditions on the state, and C is seen as a state transformer, that is, a condition constraining the program states that correspond to its initial and final states. We assume that such conditions are expressed in a fragment of some first-order language L that is closed under finite conjunction and implication (indeed a limited form of implication is often enough). We assume further that the meaning of the sentences in L is given by some theory T for which we are given a sound inference procedure denoted by ‘⊢’, that is, for each sentence ϕ ∈ L, if T ⊢ ϕ then T ⊧ ϕ. Finally, we fix an L-structure D such that D ⊧ T , which captures the domain over which computation and program reasoning take place. Let x̄ be the tuple of variables containing (among possible others) all the free variables of (2). The effect of C within the loop can be captured by stipulating that x̄ characterizes the state before execution of C, introducing a tuple of new variables x̄′ that characterizes the state after C’s execution, and by imposing restrictions on the combined tuple x̄x̄′. Our last assumption is that we are given formulas of L that correctly express the semantics of I, B, and C: let us call these formulas ϕ_{I}, ϕ_{B} and ϕ_{C}, respectively. With these definitions and assumptions, the semantics of loop (2) is correctly approximated as follows:

1. whenever the loop guard B is evaluated, ϕ_{I}[x̄] holds;

2. if ϕ_{I}[x̄] ∧ ϕ_{B}[x̄] is inconsistent, iteration of the loop terminates;

3. just before execution of C, ϕ_{I}[x̄] ∧ ϕ_{B}[x̄] holds;

4. just after execution of C, ϕ_{I}[x̄] ∧ ϕ_{B}[x̄] ∧ ϕ_{C}[x̄x̄′] holds.

It is worth observing that the presence of the externally-generated invariant I is not restrictive: on the one hand, ϕ_{I}[x̄] can simply be the “true” formula, when nothing better is available; on the other hand, non trivial invariants are usually a decisive factor for the precision of termination analysis. As observed in [9], the requirement that I must hold before any evaluation of B can be relaxed by allowing I not to hold finitely many times.^{2} The same kind of approximation can be applied to ϕ_{I}, ϕ_{B} and ϕ_{C} by only requesting that they eventually hold.

We would like to stress that, at this stage, we have not lost generality. While the formalization of basic iteration units in terms of while loops has an unmistakable imperative flavor, it is general enough to capture iteration in other programming paradigms [3].

The approach to termination analysis based on ranking functions requires that:

1. a set O and a binary relation ≺ ⊆ O × O are selected so that O is wellfounded with respect to ‘≺’;

2. a term δ[ȳ] of L is found such that

T ⊢ ∀((ϕ_{I}[x̄] ∧ ϕ_{B}[x̄] ∧ ϕ_{C}[x̄x̄′]) → ω(δ[x̄′/ȳ], δ[x̄/ȳ])),  (3)

where the interpretation of ω over D corresponds to ‘≺’; the function associated to δ in D is called ranking function for the loop (2).

Termination of (2) follows by the correctness of ϕ_{I}, ϕ_{B}, ϕ_{C} and ‘⊢’, and by well-foundedness of O with respect to ‘≺’. To see this, suppose, towards a contradiction, that loop (2) does not terminate. The mentioned soundness conditions would imply the existence of an infinite sequence of elements of O

o_{0} ≻ o_{1} ≻ o_{2} ≻ ⋯  (4)

2Such an invariant is called tail invariant in [9].

Let U ⊆ O be the (nonempty) set of elements in the sequence. Since O is well founded with respect to ‘≺’, there exists j ∈ ℕ such that, for each i ∈ ℕ, o_{i} ⊀ o_{j}. But this is impossible, as, for each j ∈ ℕ, o_{j+1} ≺ o_{j}. This means that the infinite chain (4) cannot exist and loop (2) terminates.

This general view of the ranking functions approach to termination analysis allows us to compare the methods in the literature on a common ground and focusing on what, besides mere presentation artifacts, really distinguishes them from one another. Real differences have to do with:

• the choice of the well-founded ordering (O, ≺);

• the class of functions in which the method “searches” for the ranking functions;

• the choice of the first-order languae L, the domain D and theory T ; this has to accommodate the programming formalism at hand, the semantic characterization upon which termination reasoning has to be based, the axiomatization of (O, ≺), and the representation of ranking functions;

• the class of algorithms that the method uses to conduct such a search.

We now briefly review these aspects.

The most natural well-founded ordering is, of course, (ℕ, <). This is especially indicated when the termination arguments are based on quantities that can be expressed by natural numbers. This is the case, for instance, of the work by Sohn and Van Gelder for termination analysis of logic programs [7, 10]. Orderings based on ℚ_{+} or ℝ_{+} can be obtained by imposing over them relations like those defined, for each ϵ > 0, by <_{ϵ} := { (h, k) ∈ 𝕊_{+}^{2} | h + ϵ ≤ k }, where ϵ ∈ 𝕊_{+} and 𝕊_{+} = ℚ_{+} or 𝕊_{+} = ℝ_{+}, respectively. Of course, this is simply a matter of convenience: a ranking function f with codomain (ℝ_{+}, <_{ϵ}) can always be converted into a ranking function g with codomain (ℕ, <) by taking g(ȳ) = ⌊f(ȳ)ϵ^{−1}⌋. Similarly, any ranking function over (ℝ_{+}, <_{ϵ}) can be converted into a ranking function over (ℝ_{+}, <_{1}). On tuples, the lexicographic ordering is the most common choice for a well-founded relation: given a finite number of well-founded relations ≺_{i} for i = 1, …, n over a set 𝕊, the lexicographic ordering over 𝕊^{n} is induced by saying that s ≺ t if and only if s_{i} ≺_{i} t_{i} for an index i and s_{j} = t_{j} for all indices j < i. The termination analyzer of the Mercury programming language [11, 12] first attempts an analysis using the (ℕ, <) ordering; if that fails then it resorts to lexicographic orderings. Lexicographic orderings on Cartesian products of (ℝ_{+}, <_{ϵ}) are also used in [13].

The synthesis of ranking functions is easily seen to be a search problem. All techniques impose limits upon the universe of functions that is the domain of the search. For instance, in the logic programming community, the works in [11, 12, 14, 15] use ranking functions of the form f(x_{1}, …, x_{n}) = ∑_{i=1}^{n} μ_{i}x_{i}, where, for i = 1, …, n, μ_{i} ∈ {0, 1} and the variable x_{i} takes values in ℕ. The method of Sohn and Van Gelder [7, 10] is restricted to linear functions of the form f(x_{1}, …, x_{n}) = ∑_{i=1}^{n} μ_{i}x_{i}, where, for i = 1, …, n, μ_{i} ∈ ℕ and the variable x_{i} takes values in ℕ. Its generalization to ℚ_{+} was proposed in [16] and further generalized by Mesnard and Serebrenik [17, 18] to obtain affine functions of the form f(x_{1}, …, x_{n}) = μ_{0} + ∑_{i=1}^{n} μ_{i}x_{i}, where μ_{i} ∈ ℤ and x_{i} take values in ℚ or ℝ, for i = 0, …, n. Use of the method of Podelski and Rybalchenko [8] was presented in [19] and is a component of Terminator, a termination prover of C systems code [20]. Nguyen and De Schreye [21] proposed, in the context of logic programming and following a thread of work in termination of term rewrite systems that can be traced back to [22], to use polynomial ranking functions. These are of the basic form f(x_{1}, …, x_{n}) = μ_{0}+∑_{j=1}^{m} μ_{j} ∏_{i=1}^{n} x_{i}^{k_{ij}} where μ_{0} ∈ ℤ and, for i = 1, …, n and j = 1, …, m, μ_{j} ∈ ℤ, k_{ij} ∈ ℕ and the variable x_{i} takes values in ℤ [23]. Several further restriction are usually imposed: first a domain A ⊆ ℕ is selected; then it is demanded that, for each x_{1}, …, x_{n} ∈ A, f(x_{1}, …, x_{n}) ∈ A and that f is strictly monotone over A on all its arguments. The set of all such polynomials is itself well-founded with respect to ‘<_{A}’: f <_{A} g if and only if, for each x_{1}, …, x_{n} ∈ A, f(x_{1}, …, x_{n}) < g(x_{1}, …, x_{n}). The condition of strict monotonicity, namely, for each x_{1}, …, x_{n} ∈ A, each i = 1, …, n, and each y, z ∈ A with y < z, f(x_{1}, …, x_{i−1}, y, x_{i+1}, …, x_{n}) < f(x_{1}, …, x_{i−1}, z, x_{i+1}, …, x_{n}), is ensured if, for each j = 1, …, m, we have μ_{j} ∈ ℕ and, for each i = 0, …, n, there exists j such that μ_{j} ≠ 0 and k_{ij} ≠ 0. Choosing A ≠ ℕ brings some advantages. For example, if A ⊆ { n ∈ ℕ | n ≥ 2 } then multiplication of polynomials is strictly monotone on both its arguments (i.e., f <_{A} f ·g and g <_{A} f ·g). Additional restrictions are often imposed in order to make the search of ranking functions tractable: both the maximum degree of polynomials and their coefficients —the μ_{j}’s above— can be severely limited (an upper bound of 2 both on degrees and on coefficients is typical). Quadratic ranking functions of the form f(x_{1}, …, x_{n}) = ⟨x_{1}, …, x_{n}, 1⟩^{T}M⟨x_{1}, …, x_{n}, 1⟩ are considered in [24], where the variables x_{i} and the unknown coefficients μ_{ij} of the (n + 1) × (n + 1) symmetric matrix M take values in ℝ. [13] considers a search space of tuples of (up to a fixed number of) linear functions.

The logic used in most papers about the synthesis of linear (or affine) ranking functions (such as [7, 8, 25]) is restricted to finite conjunctions of linear equalities or inequalities and simple implications (e.g., of a single inequality by a conjunction). In [13] this logic is extended to include disjunction, so as to capture precisely the effect of the loop body.

Concerning algorithms, the restriction to conjunctions of linear equalities or inequalities allows the use of the simplex algorithm (or other algorithms for linear programming) to prove the existence of linear ranking functions in [7, 8] or to synthesize one of them. When a space of ranking functions is sought, these can be obtained by projecting the systems of constraints onto a designated set of variables using, for instance, Fourier-Motzkin elimination. In these approaches, standard algorithms from linear programming work directly on an abstraction of the loop to be analyzed and are able to decide the existence of linear ranking functions for that abstraction. The algorithms used in other approaches belong to the category of “generate and test” algorithms: the “generate” phase consists in the selection, possibly guided by suitable heuristics, of candidate functions, while the “test” phase amounts to prove that a candidate is indeed a ranking function. This is the case, for instance, of [13], where generation consists in the instantiation of template functions and testing employs an algorithm based on a variant of Farkas’ Lemma. Non-linear constraints generated by the method described in [24] are handled by first resorting to semidefinite programming solvers and then validating the obtained results by using some other tools, since these solvers are typically based on interior point algorithms and hence may incur into unsafe rounding errors. Note that, in principle, the very same observation would apply to the case of linear constraints, if the corresponding linear programming problem is solved using an interior point method or even a floating-point based implementation of the simplex algorithm; however, there exist implementations of the simplex algorithm based on exact arithmetic, so that linear programming problems can be numerically solved incurring no rounding errors at all and with a computational overhead that is often acceptable.^{3}

4. The Approach of Sohn and Van Gelder, Generalized

As far as we know, the first approach to the automatic synthesis of ranking functions is due to Kirack Sohn and Allen Van Gelder [7, 10]. Even though Sohn and Van Gelder’s work concerned termination of logic programs, we will show that the key ideas of their approach can be applied, with only rather simple modifications, to the synthesis of ranking functions for any programming paradigm, thus going beyond what subsequent authors acknowledged.

4.1. Ranking Functions for Binary, Directly Recursive CLP(ℕ) Programs

We refer the interested reader to [3] for a more complete reconstruction of the approach by Sohn and Van Gelder. For the other readers, suffices it to say that termination analysis of logic programs can be mapped onto termination analysis of CLP(ℕ) programs, the termination of which implies the termination of the original program.^{4} In turn, the termination analysis of general CLP(ℕ) programs can be approximated by termination analysis of directly recursive, binary CLP(ℕ) clauses. In order to show how ranking functions can be computed from these, we deal first with a single clause

p(x̄) :− c[x̄, x̄′], p(x̄′),  (5)

where p is a predicate symbol, x̄ and x̄′ are disjoint n-tuples of variables, and c[x̄, x̄′] is a linear constraint involving variables in x̄ ∪ x̄′.^{5} The meaning of such a clause is that, if p is called on some tuple of integers x̄, then there are two cases:

3In contrast, an exact solver for non-linear constraints would probably require a truly symbolic computation, incurring a much more significant computational overhead.

4See, e.g., [26] for an introduction to CLP —constraint logic programming— languages, and [16] for the mapping of logic programs to CLP(ℕ) programs.

5As usual, we abuse notation by confusing a tuple with the set of its elements.

• c[x̄, x̄′] is unsatisfiable (i.e., there does not exist a tuple of integers x̄′ that, together with x̄, satisfies it), in which case the computation will fail, and thus terminate;

• there exists x̄′ such that c[x̄, x̄′] holds, in which case the computation proceeds with the (recursive) calls p(x̄′), for each x̄′ such that c[x̄, x̄′].

The question is now to see whether that recursive procedure is terminating, that is whether, for each x̄ ∈ ℕ^{n}, the call p(x̄) will only give rise to chains of recursive calls of finite length. The approach of Sohn and Van Gelder allows to synthesize a function f_{p}: ℕ^{n} → ℕ such that

∀x̄, x̄′ ∈ ℕ^{n} : c[x̄, x̄′] ⟹ f_{p}(x̄) > f_{p}(x̄′).  (6)

This means that the measure induced by f_{p} strictly decreases when passing from a call of p to its recursive call. Since the naturals are well founded with respect to ‘<’, this entails that p, as defined in (5), is terminating.

A very important contribution of Sohn and Van Gelder consists in the algorithm they give to construct a class of functions that satisfy (6). The class is constituted by linear functions of the form

f_{p}(y_{1}, …, y_{n}) = ∑_{i=1}^{n} μ_{i}y_{i},  (7)

where μ_{i} ∈ ℕ, for i = 1, …, n. For this class of functions and by letting μ̄ = (μ_{1}, …, μ_{n}), condition (6) can be rewritten as

∃μ̄ ∈ ℕ^{n} . ∀x̄, x̄′ ∈ ℕ^{n} : c[x̄, x̄′] ⟹ ∑_{i=1}^{n} μ_{i}x_{i} − ∑_{i=1}^{n} μ_{i}x′_{i} ≥ 1.  (8)

Given that c[x̄, x̄′] is a linear constraint, for any choice of μ̄ ∈ ℕ^{n} we can easily express (8) as an optimization problem over the naturals. In order to move from tuple notation to the more convenient vector notation, assume without loss of generality that, for some m ∈ ℕ, A_{c} ∈ ℤ^{m×2n} and b_{c} ∈ ℤ^{m} are such that A_{c}⟨x, x′⟩ ≥ b_{c} is logically equivalent to c[x̄, x̄′] under the obvious, respective interpretations. Then, for any candidate choice of μ ∈ ℕ^{n}, condition (8) is equivalent to imposing that the optimization problem

minimize θ = ⟨μ, −μ⟩^{T}⟨x, x′⟩
subject to A_{c}⟨x, x′⟩ ≥ b_{c}
x, x′ ∈ ℕ^{n}

(9)

is either unsolvable or has an optimal solution whose cost θ̂ is such that θ̂ ≥ 1. If this is the case, then μ induces, according to (7), a function f_{p} satisfying (6). Notice that, for any fixed choice of μ ∈ ℕ^{n}, θ is a linear expression and hence (9) is an integer linear programming (ILP) problem. This gives us an expensive way (since ILP is an NP-complete problem [27]) to test whether a certain μ ∈ ℕ^{n} is a witness for termination of (5), but gives us no indication about where to look for such a tuple of naturals.

A first step forward consists in considering the relaxation of (9) obtained by replacing the integrality constraints x, x′ ∈ ℕ^{n} with x, x′ ∈ ℚ_{+}^{n}. This amounts to trading precision for efficiency. In fact, since any feasible solution of (9) is also feasible for the relaxed problem, if the optimum solution of the latter has a cost greater than or equal to 1, then either (9) is unfeasible or θ̂ ≥ 1. However, we may have θ̂ ≥ 1 even if the optimum of the relaxation is less than 1.^{6} On the other hand, the relaxed problem is a linear problem: so by giving up completeness we have passed from an NP-complete problem to a problem in P for which we have, in addition, quite efficient algorithms.^{7} Furthermore, we observe that although the parameters μ are naturals in (8), this condition can be relaxed as well: if μ ∈ ℚ_{+}^{n} gives a relaxed problem with optimum greater than 1, then we can multiply this vector by a positive natural so as to obtain a tuple of naturals satisfying (8). The relaxation can now be written using the standard linear programming (LP) notation:

minimize ⟨μ, −μ⟩^{T}⟨x, x′⟩
subject to A_{c}⟨x, x′⟩ ≥ b_{c}
⟨x, x′⟩ ≥ 0.

(10)

We still do not know how to determine the vector of parameters μ so that the optimum of (10) is at least 1, but here comes one of the brilliant ideas of Sohn and Van Gelder: passing to the dual. It is a classical result of LP theory that every LP problem can be converted into an equivalent dual problem. The dual of (10) is

maximize bc^{T}y
subject to Ac^{T}y ≤ ⟨μ, −μ⟩
y ≥ 0,

(11)

where y is an m-column vector of (dual) unknowns. Duality theory ensures that if both (10) and (11) have bounded feasible solutions, then both of them have optimal solutions and these solutions have the same cost. More formally, for every choice of the parameters μ ∈ ℚ_{+}^{n}, if ⟨x̂, x̂′⟩ ∈ ℚ^{2n} is an optimal solution for (10) and ŷ ∈ ℚ^{m} is an optimal solution for (11), then ⟨μ, −μ⟩^{T}⟨x̂, x̂′⟩ = bc^{T}ŷ. Moreover, if one of (10) and (11) is unfeasible, then the other is either unbounded or unfeasible. In contrast, if one of (10) and (11) is unbounded, then the other is definitely unfeasible.

Thus, thanks to duality theory, the LP problems (10) and (11) are equivalent for our purposes and we can consider any one of them. Suppose we analyze the dual problem (11):

6Let us consider the clause: p(x) :− 2x ≥ 2x′ + 1, p(x′) with μ = 1. The optimization over the integers leads to θ̂ = 1, whereas the optimization for the relaxation has θ̂ = \frac{1}{2} .

7We denote by P the class of problems solvable in weakly polynomial time. For a formal definition of P and the notion of NP-completeness we refer the reader to, e.g., [28].

• If (11) is unfeasible then either (10) is unfeasible, which implies trivial termination of (5), or (10) is unbounded, in which case —since we are working on relaxations— nothing can be concluded about whether μ defines a ranking function for (5).

• If (11) is feasible and unbounded then (10) is unfeasible and (5) trivially terminates.

• If (11) is feasible and bounded, then we have proved termination (μ induces a ranking function) if the cost of the optimal solution is at least 1 (actually, any positive rational could be used instead of 1). The analysis is inconclusive otherwise.

The crucial point is that, in (11), the parameters μ occur linearly, whereas in (10) they are multiplied by ⟨x, x′⟩. So we can treat μ as a vector of variables and transform (11) into the new LP problem in m + n variables

maximize ⟨b_{c}, 0⟩^{T}⟨y, μ⟩
subject to (Ac^{T} \begin{matrix}−I_{n} \\ I & n \end{matrix} ) ⟨y, μ⟩ ≤ 0
⟨y, μ⟩ ≥ 0.

(12)

The requirement that, in order to guarantee termination of (5), the optimal solutions of (10) and (11) should not be less than 1 can now be captured by incorporating bc^{T}y ≥ 1 into the constraints of (12), yielding

maximize ⟨b_{c}, 0⟩^{T}⟨y, μ⟩
subject to (\begin{matrix}Ac^{T} & \begin{matrix}−I_{n} \\ I & n \end{matrix} \\ −bc^{T} & 0 \end{matrix} ) ⟨y, μ⟩ ≤ (\begin{matrix}0 \\ −1 \end{matrix})
⟨y, μ⟩ ≥ 0.

(13)

There are several possibilities:

1. If (13) is unfeasible, then:

(a) if (12) is unfeasible, then, for each μ ∈ ℚ_{+}^{n}, (11) is unfeasible and: i. if (10) is unfeasible, then (5) trivially terminates; ii. otherwise (10) is unbounded and we can conclude nothing about

the termination of (5).

(b) If (12) is feasible, then it is bounded by a rational number q < 1. Thus, for each μ̌ ∈ ℚ_{+}^{n} extracted from a feasible solution ⟨y̌, μ̌⟩ ∈ ℚ_{+}^{m+n} of (12), the corresponding LP problem (11) is also feasible, bounded, and its optimum q′ ∈ ℚ is such that q′ ≤ q < 1. Moreover, we must have q′ ≤ 0. In fact, if q′ > 0, problem (10) instantiated over μ̌′ := μ̌/q′ would have an optimal solution of cost 1; the same would hold for the corresponding dual (11), but this would contradict the hypothesis that (12) is bounded by q < 1. Hence q′ ≤ 0. Since by duality the optimum of problem (10) is q′, the analysis is inconclusive.

2. If (13) is feasible, let ⟨y̌, μ̌⟩ ∈ ℚ^{m+n} be any of its feasible solutions. Choosing μ̌ for the values of the parameters, (11) is feasible. There are two further possibilities:

(a) either (11) is unbounded, so (5) trivially terminates;

(b) or it is bounded by a rational q ≥ 1 and the same holds for its dual (10).

In both cases, μ̌, possibly multiplied by a positive natural in order to get a tuple of naturals, defines, via (7), a ranking function for (5).

The above case analysis boils down to the following algorithm:

1. Use the simplex algorithm to determine the feasibility of (13), ignoring the objective function. If it is feasible, then any feasible solution induces a linear ranking function for (5); exit with success.

2. If (13) is unfeasible, then try to determine the feasibility of (9) (e.g., by using the simplex algorithm again to test whether the relaxation (10) is feasible). If (9) is unfeasible then (5) trivially terminates; exit with success.

3. Exit with failure (the analysis is inconclusive).

An example should serve to better clarify the methodology we have employed.

Example 4.1. In the CLP(ℕ) program

p(x_{1}, x_{2}) :− x_{1} ≤ 1 ∧ x_{2} = 0,
p(x_{1}, x_{2}) :− x_{1} ≥ 2 ∧ 2x′_{1} + 1 ≥ x_{1} ∧ 2x′_{1} ≤ x_{1} ∧ x′_{2} + 1 = x_{2}, p(x′_{1}, x′_{2}),

p(x_{1}, x_{2}) is equivalent to

x_{2} = \begin{cases}⌊log_{2}(x_{1})⌋, & if x_{1} ≠ 0; \\ 0, & otherwise. \end{cases}

The relaxed optimization problem in LP notation (10) is^{8}

minimize ⟨μ_{1}, μ_{2}, −μ_{1}, −μ_{2}⟩^{T}⟨x_{1}, x_{2}, x′_{1}, x′_{2}⟩

subject to (\begin{matrix}1 & 0 & 0 & 0 \\ −1 & 0 & 2 & 0 \\ 1 & 0 & −2 & 0 \\ 0 & 1 & 0 & −1 \\ 0 & −1 & 0 & 1 \end{matrix}) (\begin{matrix}x_{1} \\ x_{2} \\ x′_{1} \\ x′_{2} \end{matrix}) ≥ (\begin{matrix}2 \\ −1 \\ 0 \\ 1 \\ −1 \end{matrix})

⟨x_{1}, x_{2}, x′_{1}, x′_{2}⟩ ≥ 0,

8We will tacitly replace an equality in the form α = β by the equivalent pair of inequalities α ≥ β and −α ≥ −β whenever the substitution is necessary to fit our framework.

and the dual optimization problem (11) is

maximize ⟨2, −1, 0, 1, −1⟩^{T}⟨y_{1}, y_{2}, y_{3}, y_{4}, y_{5}⟩

subject to (\begin{matrix}1 & −1 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 & −1 \\ 0 & 2 & −2 & 0 & 0 \\ 0 & 0 & 0 & −1 & 1 \end{matrix}) (\begin{matrix}y_{1} \\ y_{2} \\ y_{3} \\ y_{4} \\ y_{5} \end{matrix}) ≤ (\begin{matrix}μ_{1} \\ μ_{2} \\ −μ_{1} \\ −μ_{2} \end{matrix})

⟨y_{1}, y_{2}, y_{3}, y_{4}, y_{5}⟩ ≥ 0.

Incorporation of the unknown coefficients of μ among the problem variables finally yields as the transformed problem (13):

maximize ⟨2, −1, 0, 1, −1, 0, 0⟩^{T}⟨y_{1}, y_{2}, y_{3}, y_{4}, y_{5}, μ_{1}, μ_{2}⟩

subject to (\begin{matrix}1 & −1 & 1 & 0 & 0 & −1 & 0 \\ 0 & 0 & 0 & 1 & −1 & 0 & −1 \\ 0 & 2 & −2 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & −1 & 1 & 0 & 1 \\ −2 & 1 & 0 & −1 & 1 & 0 & 0 \end{matrix}) (\begin{matrix}y_{1} \\ y_{2} \\ y_{3} \\ y_{4} \\ y_{5} \\ μ_{1} \\ μ_{2} \end{matrix}) ≤ (\begin{matrix}0 \\ 0 \\ 0 \\ 0 \\ −1 \end{matrix})

⟨y_{1}, y_{2}, y_{3}, y_{4}, y_{5}, μ_{1}, μ_{2}⟩ ≥ 0

(14)

This problem is feasible so this CLP(ℕ) program terminates. Projecting the constraints of (14) onto μ we obtain, in addition, the knowledge that every μ with μ_{1} + μ_{2} ≥ 1 gives a ranking function. In other words, μ_{1}x_{1} + μ_{2}x_{2} is a ranking function if the non-negative numbers μ_{1} and μ_{2} satisfy μ_{1} + μ_{2} ≥ 1.

The following result illustrates the strength of the method:

Theorem 4.2. Let C be the binary CLP(ℚ_{+}) clause p(x̄) :− c[x̄, x̄′], p(x̄′), where p is an n-ary predicate and c[x̄, x̄′] is a linear satisfiable constraint. Let plrf(C) be the set of positive linear ranking functions for C and svg(C) be the set of solutions of (13) projected onto μ, that is,

plrf(C) := { μ ∈ ℚ_{+}^{n} | ∀x̄, x̄′ ∈ ℚ_{+}^{n} : c[x̄, x̄′] ⟹ ∑_{i=1}^{n} μ_{i}x_{i} − ∑_{i=1}^{n} μ_{i}x′_{i} ≥ 1 },
svg(C) := { μ̌ ∈ ℚ_{+}^{n} | ⟨y̌, μ̌⟩ is a solution of (13) }.

Then plrf(C) = svg(C).

Proof. As c[x̄, x̄′] is satisfiable, problem (10) is feasible. We prove each inclusion separately.

svg(C) ⊆ plrf(C). Assume that (13) is feasible and let ⟨y̌, μ̌⟩ be a solution of (13). For this choice of μ̌, the corresponding LP problems (10) and (11) are bounded by q ≥ 1 (case 2b of the discussion above). So μ̌ ∈ plrf(C). plrf(C) ⊆ svg(C). Let us pick μ ∈ plrf(C). For this choice, the corresponding LP problem (10) is bounded by r ≥ 1, so is its dual (11). Let ŷ be an optimal solution for (11). Thus ⟨ŷ, μ⟩ is a feasible solution of (12) and (13). Hence μ ∈ svg(C).

As an immediate consequence, the question “does a given binary recursive clause with linear constraint admit a positive linear mapping?” can be solved in weakly polynomial time.

Corollary 4.3. Let C be the binary CLP(ℚ_{+}) clause p(x̄) :− c[x̄, x̄′], p(x̄′), where c[x̄, x̄′] is a linear satisfiable constraint. The decision problem plrf(C) = ∅ is in P.

Proof. By Theorem 4.2 the problems plrf(C) = ∅ and svg(C) = ∅ are equivalent. So, if (13) is feasible then the answer is no: as c[x̄, x̄′] is satisfiable, we are in case (2)(b). Otherwise, again because of the satisfiability of c[x̄, x̄′], either (10) is unbounded (case (1)(a)ii. or it is bounded by q′ < 0 (case (1)(b)). In both cases, the answer is yes. Finally, testing the satisfiability of a linear system, as well as computing one of its solutions —and thus computing one concrete linear ranking function—, is in P (see, e.g., [28]).

For the case where we have more than one directly recursive binary CLP(ℕ) clauses, C_{1}, …, C_{n}, the set of global positive linear ranking functions, i.e., that ensure termination whichever clause is selected at each computation step, is given by ⋂_{i=1}^{n} svg(C_{i}). This can be computed by taking the conjunction of the constraints obtained, for each clause, from the projection of the constraints of the corresponding linear problem (13) onto μ.

To summarize, the main contribution of Sohn and Van Gelder lies in their encoding of the ranking function search problem into linear programming and their use of the duality theorem. As we will see, this idea is amenable to a generalization that makes it widely applicable to any programming paradigm, not just (constraint) logic programming.

4.2. The Generalization by Mesnard and Serebrenik

Fred Mesnard and Alexandre Serebrenik have generalized the method of Sohn and Van Gelder from the analysis of logic programs to the analysis of CLP(ℚ) and CLP(ℝ) programs in [17, 18]. In the following, for presentation purposes and without loss of generality, we consider the case of rational-valued variables. They use a class of affine ranking functions of the form

f_{p}(y_{1}, …, y_{n}) = μ_{0} + ∑_{i=1}^{n} μ_{i}y_{i},  (15)

where μ_{i} ∈ ℚ, for i = 0, …, n. Allowing for rational-valued coefficients μ_{i} and variables y_{i} (both the μ_{i}’s and the y_{i}’s were naturals in [7]) implies that (15) does not necessarily define a nonnegative function and that Zeno sequences^{9} are not automatically excluded. Consequently, to avoid these two problems, condition (6) is strengthened to^{10}

∀x̄, x̄′ ∈ ℚ^{n} : c[x̄, x̄′] ⟹ (f_{p}(x̄) ≥ 1 + f_{p}(x̄′) ∧ f_{p}(x̄) ≥ 0).  (16)

Note that the choice of the numbers 1 and 0 in the right hand side of the above implication preserves generality: the general form of the former condition, i.e., f_{p}(x̄) ≥ ϵ + f_{p}(x̄′) for a fixed and strictly positive ϵ ∈ ℚ_{+}, can be transformed as shown in Section 3, and the general form of the latter, i.e., f_{p}(x̄) ≥ b for a fixed b ∈ ℚ, can be transformed into f_{p}(x̄) ≥ 0 by a suitable choice of μ_{0}. Condition (16) can be rewritten as

∀x̄, x̄′ ∈ ℚ^{n} : c[x̄, x̄′] ⟹ (∑_{i=1}^{n} μ_{i}x_{i} − ∑_{i=1}^{n} μ_{i}x′_{i} ≥ 1 ∧ μ_{0} + ∑_{i=1}^{n} μ_{i}x_{i} ≥ 0). (17)

Using the same notation chosen for (9), the existence of a ranking function can now be equivalently expressed as the existence of a solution of cost at least 1 to the former and a solution of cost at least 0 to the latter of the following optimization problems:

minimize ⟨μ, −μ⟩^{T}⟨x, x′⟩ subject to A_{c}⟨x, x′⟩ ≥ b_{c}  minimize ⟨μ̃, 0⟩^{T}⟨x̃, x′⟩ subject to Ã_{c}⟨x̃, x′⟩ ≥ b̃_{c} (18)

where the extended vectors μ̃ := ⟨μ_{0}, μ⟩ and x̃ := ⟨x_{0}, x⟩ include the parameter μ_{0} and the new variable x_{0}, respectively, and the extended matrix and vector

Ã_{c} := (\begin{matrix}1 & 0 \\ −1 & 0 \\ 0 & A_{c} \end{matrix}) and b̃_{c} := ⟨1, −1, b_{c}⟩

encode the additional constraint x_{0} = 1.

Reasoning as in Section 4.1, the problems (18) can then be transformed, applying the suitable form of the duality theorem, into the following dual problems over new vectors of variables y and z, ranging over ℚ^{m} and ℚ^{m+2}, respectively:

maximize bc^{T}y subject to Ac^{T}y = ⟨μ, −μ⟩ y ≥ 0  maximize b̃c^{T}z subject to Ãc^{T}z = ⟨μ̃, 0⟩ z ≥ 0  (19)

Now the condition that the optimal solution is at least 1 (resp., 0) can be added to the constraints, thus reducing the optimization problems (18) to testing the satisfiability of the system:

9Such as 1, \frac{1}{2} , \frac{1}{4} , \frac{1}{8} , ….

10Our presentation is strictly more general than the formulation in [17, 18], which imposes that f_{p}(x̄) ≥ 1 + f_{p}(x̄′) ∧ f_{p}(x̄′) ≥ 0.

\{ \begin{matrix}bc^{T}y ≥ 1 \\ Ac^{T}y = ⟨μ, −μ⟩ \\ y ≥ 0 \\ b̃c^{T}z ≥ 0 \\ Ãc^{T}z = ⟨μ̃, 0⟩ \\ z ≥ 0 \end{matrix}

or equivalently, after incorporating the parameters μ (resp., μ̃) into the variables, to the generalization to ℚ of problem (13):

( \begin{matrix}^{Ac^{T}} & \begin{matrix}−I_{n} \\ I & n \end{matrix} \\ −Ac^{T} & \begin{matrix}I & n \\ −I_{n} \end{matrix} \\ −I_{m} & 0 \\ −bc^{T} & 0 \end{matrix} ) ⟨y, μ⟩ ≤ ( \begin{matrix}0 \\ 0 \\ 0 \\ −1 \end{matrix} ) ⋀ ( \begin{matrix}^{Ãc^{T}} & \begin{matrix}−I_{n+1} \\ 0 \end{matrix} \\ −Ãc^{T} & \begin{aligned}I & n+1 \\ & 0 \end{aligned} \\ − I_{m+2} & 0 \\ −b̃c^{T} & 0 \end{matrix} ) ⟨z, μ̃⟩ ≤ 0. (20)

The following completeness result generalizes Theorem 4.2:

Theorem 4.4. Let C be the binary CLP(ℚ) clause p(x̄) :− c[x̄, x̄′], p(x̄′), where p is an n-ary predicate and c[x̄, x̄′] is a linear satisfiable constraint. Let lrf(C) be the set of linear ranking functions for C and ms(C) be the set of solutions of (20) projected onto μ̃, that is,

lrf(C) := \{ μ̃ ∈ ℚ^{n+1} | \begin{matrix}∀x̄, x̄′ ∈ ℚ^{n} : c[x̄, x̄′] ⟹ \\ ∑_{i=1}^{n} μ_{i}x_{i} − ∑_{i=1}^{n} μ_{i}x′_{i} ≥ 1 ∧ μ_{0} + ∑_{i=1}^{n} μ_{i}x_{i} ≥ 0 \end{matrix} \},
ms(C) := { μ̃ ∈ ℚ^{n+1} | ⟨y, μ⟩ and ⟨z, μ̃⟩ are solutions of the problems (20) }.

Then lrf(C) = ms(C).

Proof. We use l and r as subscripts of our references to the LP problems (18), (19), and (20) to denote the LP problems on the left and the LP problems on the right.

ms(C) ⊆ lrf(C). Assume that (20) is feasible and let ⟨y̌, μ̌⟩ be a solution of (20)_{l} and ⟨ž, μ̌⟩ be a solution of (20)_{r}. For this choice of μ̌, the corresponding LP problems (19)_{l} and (18)_{l} are bounded by 1 while the corresponding LP problems (19)_{r} and (18)_{r} are bounded by 0. Hence we have

∀x̄, x̄′ ∈ ℚ^{n} : c[x̄, x̄′] ⟹ ∑_{i=1}^{n} μ_{i}x_{i} − ∑_{i=1}^{n} μ_{i}x′_{i} ≥ 1

and

∀x̄, x̄′ ∈ ℚ^{n} : c[x̄, x̄′] ⟹ μ_{0} + ∑_{i=1}^{n} μ_{i}x_{i} ≥ 0.

Thus

∀x̄, x̄′ ∈ ℚ^{n} : c[x̄, x̄′] ⟹ ∑_{i=1}^{n} μ_{i}x_{i} − ∑_{i=1}^{n} μ_{i}x′_{i} ≥ 1 ∧ μ_{0} + ∑_{i=1}^{n} μ_{i}x_{i} ≥ 0,

so that μ̌ ∈ lrf(C).

lrf(C) ⊆ ms(C). Let us pick μ̃ ∈ lrf(C). For this choice, the corresponding LP problem (18) are bounded by 1 and 0, and so are their duals (19). Let ŷ be an optimal solution for (19)_{l}. Thus ⟨ŷ, μ⟩ is a feasible solution of (20)_{l}. Similarly, let ẑ be an optimal solution for (19)_{r}. Thus ⟨ẑ, μ̃⟩ is a feasible solution of (20)_{r}. Hence μ̃ ∈ ms(C).

Moreover, even for the case of the linear fragment of CLP(ℚ) —and CLP(ℝ)— checking for the existence of a linear ranking function is a weakly polynomial problem.

Corollary 4.5. Let C be the binary CLP(ℚ) clause p(x̄) :− c[x̄, x̄′], p(x̄′), where c[x̄, x̄′] is a linear satisfiable constraint. The decision problem lrf(C) = ∅ is in P.

A space of ranking functions can be obtained (at a computational price that is no longer polynomial) by projecting the constraints of (20) onto μ̃. Any μ̃ satisfying all the projected constraints corresponds to one ranking function that, subject to c[x̄, x̄′], is bounded from below by 0 and that decreases by at least 1 at each iteration. From these “normalized” ranking functions, the opposite of the transformation outlined in Section 3 allows to recover all affine ranking functions: these are induced by the set of parameters

{ ⟨h, kμ⟩ | ⟨μ_{0}, μ⟩ ∈ lrf(C), h ∈ ℚ, k ∈ ℚ_{+} \ {0} }.  (21)

4.3. Application to the Analysis of Imperative While Loops

The generalization of Mesnard and Serebrenik can be used, almost unchanged, to analyze the termination behavior of imperative while loops with integer- or rational-valued variables. Consider a loop of the form (2), i.e., { I } while B do C where I is known to hold before any evaluation of B and C is known to always terminate in that loop. Termination analysis is conducted as follows:

1. Variables are duplicated: if x̄ are the n variables of the original loop, we introduce a new tuple of variables x̄′.

2. An analyzer based on convex polyhedra [29] is used to analyze the following program:

\begin{matrix}\{I\} \\ x′_{1} := x_{1}; …; x′_{n} := x_{n}; \\ if B[x̄′/x̄] then \\ C[x̄′/x̄] \\ ⋆ \end{matrix}  (22)

Let the invariant obtained for the program point marked with ‘⋆’ be c[x̄, x̄′]; this is a finite conjunction of linear constraints.

3. The method of Mesnard and Serebrenik is now applied to the CLP(ℚ) clause p(x̄) :− c[x̄, x̄′], p(x̄′): if termination can be established for that clause, then the while loop we started with will terminate.

Notice how the clause p(x̄) :− c[x̄, x̄′], p(x̄′) approximates the termination behavior of the loop: if we interpret the predicate p applied to x̄ as “the loop guard is evaluated on values x̄,” then the clause can be read as “if the loop guard is evaluated on values x̄, and c[x̄, x̄′] holds, then the loop guard will be evaluated again on values x̄′.”

We illustrate the overall methodology with an example.

Example 4.6. The following program, where x_{1} and y take values in ℤ, computes and stores in x_{2} the integer base-2 logarithm of x_{1} if x_{1} > 0, 0 otherwise:

x_{2} := 0;
{x_{1} ≥ 0 ∧ x_{2} ≥ 0}
while x_{1} ≥ 2 do
x_{1} := x_{1} div 2;
x_{2} := x_{2} + 1

where the loop invariant {x_{1} ≥ 0∧x_{2} ≥ 0} has been obtained by static analysis. After the duplication of variables, we submit to the analyzer the program

{x_{1} ≥ 0 ∧ x_{2} ≥ 0}
x′_{1} := x_{1}; x′_{2} := x_{2};
if x′_{1} ≥ 2 then
x′_{1} := x′_{1} div 2;
x′_{2} := x′_{2} + 1
⋆

and we obtain, for program point ‘⋆’, the invariant

x_{1} ≥ 2 ∧ 2x′_{1} + 1 ≥ x_{1} ∧ 2x′_{1} ≤ x_{1} ∧ x′_{2} = x_{2} + 1 ∧ x′_{2} ≥ 1.

Applying the method of Mesnard and Serebrenik we obtain that, for each μ_{0}, μ_{1}, μ_{2} ∈ ℚ such that μ_{1} − μ_{2} ≥ 1, μ_{2} ≥ 0, and μ_{0} + 2μ_{1} ≥ 0, the function f(x_{1}, x_{2}) := μ_{0} + μ_{1}x_{1} + μ_{2}x_{2} is a ranking function for the given while loop. It is interesting to observe that the first constraint guarantees strict decrease (at least 1), the addition of the second constraint guarantees boundedness from below, while the further addition of the third constraint ensures nonnegativity, i.e., that 0 is a lower bound.

4.4. Application to Conditional Termination Analysis

An important observation is that the method of Mesnard and Serebrenik is immediately applicable in conditional termination analysis. This is the problem of (automatically) inferring the preconditions under which code that does not universally terminate (i.e., there are inputs for which it does loop forever) is guaranteed to terminate. This problem has been recently studied in [30], where preconditions are inferred under which functions that are either decreasing or bounded become proper ranking functions. The two systems in (20), projected onto μ̃, exactly define the space of non-negative candidate ranking functions and the space of decreasing candidate ranking functions, respectively. While this is subject for future research, we believe that the availability of these two spaces allows to improve the techniques presented in [30].

5. The Approach of Podelski and Rybalchenko

Andreas Podelski and Andrey Rybalchenko [8] introduce a method for finding linear ranking functions for a particular class of unnested while loops that, with the help of a preliminary analysis phase, is indeed completely general.

Consider a while loop of the form

\begin{matrix}\{I\} \\ while B do \\ ▼ \\ C \\ ⋆ \end{matrix}  (23)

in which variables x_{1}, …, x_{n} occur. Suppose we have determined (e.g., by a data-flow analysis based on convex polyhedra) that the invariant

∑_{i=1}^{n} g_{k,i}x_{i} ≤ b_{k},  for k = 1, …, r,  (24)

holds at the program point marked with ‘▼’, while the invariant

∑_{i=1}^{n} a′_{k,i}x′_{i} ≤ ∑_{i=1}^{n} a_{k,i}x_{i} + b_{k},  for k = r + 1, …, r + s,  (25)

holds at the program point marked with ‘⋆’, where unprimed variables represent the values before the update and primed variables represent the values after the update, and all the coefficients and variables are assumed to take values in ℚ.^{11}

The inequalities in (24) can be expressed in the form (25) by just defining a′_{k,i} := 0 and a_{k,i} := −g_{k,i} for i = 1, …, n and k = 1, …, r. The conjunction of (24) and (25) can now be stated in matrix form as

(A A′) (\begin{matrix}x \\ x & ′ \end{matrix}) ≤ b,  (26)

where the matrix (A A′) is obtained by juxtaposition of the two (r + s) × n matrices A := (−a_{k,i}) and A′ := (a′_{k,i}), b := ⟨b_{1}, b_{2}, …, b_{r+s}⟩ and, as explained in Section 2, ⟨x, x′⟩ is obtained by juxtaposing the vectors x := ⟨x_{1}, x_{2}, …, x_{n}⟩ and x′ := ⟨x′_{1}, x′_{2}, …, x′_{n}⟩.

Podelski and Rybalchenko have proved that (23) is guaranteed to terminate on all possible inputs if there exist two (r+s)-dimensional non-negative rational vectors λ_{1} and λ_{2} such that:^{12}

λ1^{T}A′ = 0,  (27a)
(λ1^{T} − λ2^{T})A = 0,  (27b)
λ2^{T}(A + A′) = 0,  (27c)
λ2^{T}b < 0.  (27d)

Note that we have either zero or infinitely many solutions, since if the pair of vectors λ_{1} and λ_{2} satisfies the constraints, then the pair kλ_{1} and kλ_{2} satisfies them as well, for any k ∈ ℚ_{+} \ {0}. Podelski and Rybalchenko proved also the following completeness result: if the iterations of (23) are completely characterized by conditions (24) and (25) —in which case they call it a “simple linear loop”— then λ_{1}, λ_{2} ∈ ℚ_{+}^{r+s} satisfying conditions (27a)–(27d) exist if and only if the program terminates for all inputs.

5.1. Generation of Ranking Functions

For each pair of vectors λ_{1} and λ_{2} satisfying the conditions (27a)–(27d), a linear ranking function for the considered program can be obtained as

f(x) := λ2^{T}A′x.  (28)

11In [8] variables are said to have integer domain, but this restriction seems unnecessary and, in fact, it is not present in [31].

12For an informal justification of these equations, see Section 5.2; a more detailed explanation is available in [3].

In [8] a slightly more complex form is proposed, namely:

g(x) := \begin{cases}λ2^{T}A′x, & if there exists x′ such that (A A′) (\begin{matrix}x \\ x & ′ \end{matrix}) ≤ b , \\ (λ2^{T} − λ1^{T})b, & otherwise, \end{cases}(29)

but the extra provisions are actually necessary only if one is interested into an “extended ranking function” that is strictly decreasing also on the very last iteration of the loop, that is, when the effect of the command C is such that x would violate the loop guard B at the following iteration. As this more complex definition does not seem to provide any additional benefit, we disregard it and consider only the linear ranking function (28).

Example 5.1. Consider again the program of Example 4.6. The invariants in the forms dictated by (24) and (25) are given by the systems {−x_{1} ≤ −2, −x′_{2} ≤ −1} and {2x′_{1} ≤ x_{1}, −2x′_{1}−1 ≤ −x_{1}, −x′_{2} ≤ −x_{2}−1, x′_{2} ≤ x_{2}+1}, respectively. These can be expressed in the matrix form (26) by letting

A := (\begin{matrix}−1 & 0 \\ −1 & 0 \\ 1 & 0 \\ 0 & 1 \\ 0 & −1 \\ 0 & 0 \end{matrix}) , A′ := (\begin{matrix}0 & 0 \\ 2 & 0 \\ −2 & 0 \\ 0 & −1 \\ 0 & 1 \\ 0 & −1 \end{matrix}) , b = (\begin{matrix}−2 \\ 0 \\ 1 \\ −1 \\ 1 \\ −1 \end{matrix}) .

Two non-negative rational vectors solving the system (27) are, for instance, λ_{1} = ⟨2, 0, 0, 0, 0, 0⟩^{T} and λ_{2} = ⟨1, 1, 0, 0, 0, 0⟩^{T}.

5.2. Justification of the Approach

A reader of [8] wonders where the method of Podelski and Rybalchenko comes from. In fact, the paper does not give an intuition about why conditions (27a)–(27d) imply termination of (23). Those conditions can be mapped into a strengthening, tailored to the linear case, of the well known Floyd termination verification conditions,^{13} but such a higher level view needs to be extracted, with some effort, from the details of the proof of [8, Theorem 1]. The relative completeness of the approach is then proved in [8, Theorem 2] by exploiting the affine form of Farkas’ Lemma, showing that such a strengthening is unconsequential for the case of linear ranking functions and simple linear loops.

The intuitive reading hidden in the proof details is made explicit in [24, Section 6.2], where Patrick Cousot hints that the method by Podelski and Rybalchenko can be derived from the Floyd termination verification conditions by application of Lagrangian relaxation. The interested reader can find more details on this connection in the technical report version of this paper [3].

13A. Rybalchenko, personal communication, 2011.

5.3. An Alternative Implementation Approach

As long as the distinction between invariants (24) and (25) is retained, the method of Podelski and Rybalchenko can be implemented following an alternative approach. The linear invariants (26) are more precisely described by

( \begin{matrix}^{A_{B}} & 0 \\ ^{A_{C}} & ^{A′_{C}} \end{matrix} ) (\begin{matrix}x \\ x & ′ \end{matrix}) ≤ (\begin{matrix}b & B \\ b & C \end{matrix})  (30)

where A_{B} ∈ ℚ^{r×n}, A_{C} ∈ ℚ^{s×n}, A′_{C} ∈ ℚ^{s×n}, b_{B} ∈ ℚ^{r}, b_{C} ∈ ℚ^{s}. As shown in the technical report version of this paper [3, Section 5.3], the existence of a linear ranking function for the system (30) is equivalent to the existence of three vectors v_{1} ∈ ℚ_{+}^{r}, v_{2} ∈ ℚ_{+}^{r}, v_{3} ∈ ℚ_{+}^{s} such that

(v_{1} − v_{2})^{T}A_{B} − v3^{T}A_{C} = 0,  (31a)
v2^{T}A_{B} + v3^{T}(A_{C} + A′_{C}) = 0,  (31b)
v_{2}b_{B} + v_{3}b_{C} < 0.  (31c)

As already noted, the two vectors of the original Podelski and Rybalchenko method can be reconstructed as λ_{1} = ⟨v_{1}, 0⟩ and λ_{2} = ⟨v_{2}, v_{3}⟩.

Note that the same approach is still valid when starting from the single matrix form (26) in full generality, i.e., when we can’t assume that the distinction between invariants (24) and (25) has been retained or that invariants are listed in the order we used to build the matrix form (26): it is enough to apply a straightforward permutation to (26) to rearrange it in the form (30). In that case, due to the permutation involved, we would solve a different linear programming problem; however, we still obtain the same space of linear ranking functions we would have obtained by applying the original method starting from the matrix form (26), as we prove using the following lemma.

Lemma 5.2. Let S be the space of linear ranking functions obtained by applying the method of Podelski and Rybalchenko to (A A′) ⟨x, x′⟩ ≤ b, i.e.,

S := { ⟨λ2^{T}A′, λ1^{T}b⟩ ∈ ℚ^{n+1} | ⟨λ_{1}, λ_{2}⟩ is a solution of (27) },

and let P ∈ ℚ^{(r+s)×(r+s)} be a permutation matrix.^{14} Then the application of the method of Podelski and Rybalchenko to P (A A′) ⟨x, x′⟩ ≤ P b yields the same space of linear ranking functions S.

Proof. The system (27) corresponding to P (A A′) ⟨x, x′⟩ ≤ P b becomes

η1^{T}P A′ = 0,  (32a)
(η1^{T} − η2^{T})P A = 0,  (32b)
η2^{T}P (A + A′) = 0,  (32c)
η2^{T}P b < 0,  (32d)

14We recall that a k-dimensional permutation matrix is a square matrix obtained by a permutation of the rows or columns of the k-dimensional identity matrix.

to be solved for the two (r + s)-dimensional non-negative rational vectors η_{1} and η_{2}.

Now, ⟨λ_{1}, λ_{2}⟩ is a solution of (27) if and only if ⟨λ_{1}, λ_{2}⟩P^{−1} is a solution of (32): on one side, if ⟨λ_{1}, λ_{2}⟩ is a solution of (27) then ⟨η_{1}, η_{2}⟩ defined as ⟨η_{1}, η_{2}⟩ := ⟨λ_{1}, λ_{2}⟩P^{−1} is a solution of (32); on the other side, if ⟨η_{1}, η_{2}⟩ is a solution of (32) then ⟨λ_{1}, λ_{2}⟩ defined as ⟨λ_{1}, λ_{2}⟩ := ⟨η_{1}, η_{2}⟩P is a solution of (27) and the desired property can be verified by right-multiplying by P^{−1} both solutions.

The space of linear ranking functions for the permuted system is

S_{P} = { ⟨η2^{T}P A′, η1^{T}P b⟩ ∈ ℚ^{n+1} | ⟨η_{1}, η_{2}⟩ is a solution of (32) }

= { ⟨λ2^{T}P^{−1}P A′, λ1^{T}P^{−1}P b⟩ ∈ ℚ^{n+1} | ⟨λ_{1}, λ_{2}⟩ is a solution of (27) }

= S,

and thus it is unaltered with respect to the space of linear ranking functions S corresponding to the non-permuted system.

Since the system (30) is obtained by applying a suitable permutation to (26), a straightforward application of this lemma proves that the space of linear ranking functions obtained is the same in both cases.

Moreover, as λ2^{T}A′ = v3^{T}A′_{C} and λ_{1}b = v_{1}b_{B}, we can express the space of linear ranking functions as

S := { ⟨v3^{T}A′_{C}, v_{1}b_{B}⟩ ∈ ℚ^{n+1} | ⟨v_{1}, v_{2}, v_{3}⟩ is a solution of (31) }.

6. Comparison of the Two Methods

In this section we compare the method by Mesnard and Serebrenik with the method by Podelski and Rybalchenko: we first prove that they have the same “inferential power”, then we compare their worst-case complexities, then we experimentally evaluate them on a representative set of benchmarks.

6.1. Equivalence of the Two Methods

We will now show that the method proposed in [8] is equivalent to the one given in [32] on the class of simple linear loops, i.e., that if one of the two methods can prove termination of a given simple linear loop, then the other one can do the same. This is an expected result since both methods claim to be complete on the class of programs considered.

It is worth noting that a completeness result was already stated in [16, Theorem 5.1] for the case of single predicate CLP(ℚ_{+}) procedures, which can be seen to be a close variant of the binary, directly recursive CLP(ℚ_{+}) programs considered in Theorem 4.2 and Corollary 4.3. Probably due to the programming paradigm mismatch, Podelski and Rybalchenko [8] fail to recognize the actual strength and generality of the mentioned result, thereby claiming originality for their completeness result.

Theorem 6.1. Let C be the binary CLP(ℚ) clause p(x̄) :− c[x̄, x̄′], p(x̄′), where p is an n-ary predicate and c[x̄, x̄′] is a linear satisfiable constraint. Let pr(C) and m̂s(C) be the spaces of linear ranking functions for C obtained through the method of Podelski and Rybalchenko and through the method of Mesnard and Serebrenik, respectively, that is,

pr(C) := { ⟨λ2^{T}A′, λ1^{T}b⟩ ∈ ℚ^{n+1} | ⟨λ_{1}, λ_{2}⟩ is a solution of (27) },
m̂s(C) := \{ kμ̃ ∈ ℚ^{n+1} | \begin{matrix}⟨y, μ⟩ and ⟨z, μ̃⟩ are solutions of (20), \\ k ∈ ℚ_{+} \ \{0\} \end{matrix} \}.

where c[x̄, x̄′] is equivalent to (A A′) ⟨x, x′⟩ ≤ b or to A_{c}⟨x, x′⟩ ≥ b_{c}, respectively. Then pr(C) = m̂s(C).

Proof. We will, as customary, prove the two inclusions pr(C) ⊆ m̂s(C) and pr(C) ⊇ m̂s(C).

pr(C) ⊆ m̂s(C). Suppose that there exist two non-negative rational vectors λ_{1} and λ_{2} satisfying (27), i.e., λ1^{T}A′ = (λ1^{T} − λ2^{T})A = λ2^{T}(A + A′) = 0 and λ2^{T}b < 0. We need to show that ⟨λ2^{T}A′, λ1^{T}b⟩ ∈ m̂s(C), which is equivalent to proving that there exists a positive coefficient (that we can denote with \frac{1}{k} without loss of generality) \frac{1}{k} ∈ ℚ_{+} \ \{0\} such that ⟨\frac{1}{k} λ2^{T}A′, \frac{1}{k} λ1^{T}b⟩ ∈ ms(C), or, by Theorem 4.4, that ⟨\frac{1}{k} λ2^{T}A′, \frac{1}{k} λ1^{T}b⟩ ∈ lrf(C), which is in turn equivalent, by definition, to λ2^{T}A′x − λ2^{T}A′x′ ≥ k and λ1^{T}b + λ2^{T}A′x ≥ 0. We have

(A A′) (\begin{matrix}x \\ x & ′ \end{matrix}) ≤ b ⟹ Ax + A′x′ ≤ b
⟹ − Ax ≥ A′x′ − b
⟹ − λ2^{T}Ax ≥ λ2^{T}A′x′ − λ2^{T}b  [by λ_{2} ≥ 0]
⟹ λ2^{T}A′x ≥ λ2^{T}A′x′ − λ2^{T}b  [by (27c)]

and the former property is satisfied if we choose k = −λ2^{T}b, which is nonnegative by relation (27d). For the latter property, we have

Ax + A′x′ ≤ b ⟹ λ1^{T}Ax + λ1^{T}A′x′ ≤ λ1^{T}b as λ1^{T} is non-negative
⟹ λ1^{T}Ax ≤ λ1^{T}b  [by (27a)]
⟹ λ2^{T}Ax ≤ λ1^{T}b  [by (27b)]
⟹ − λ2^{T}A′x ≤ λ1^{T}b  [by (27c)]

and both properties are thus proved.

pr(C) ⊇ m̂s(C). In order to prove the inverse containment, we will need to recall the affine form of Farkas’ Lemma (see [28]).

Lemma 6.2 (Affine form of Farkas’ lemma). Let P be a nonempty polyhedron defined by the inequalities Cx + d ≥ 0. Then an affine function f(x) is non-negative everywhere in P if and only if it is a positive affine combination of the columns of Cx + d: f(x) = λ_{0} + λ^{T}(Cx + d) with λ_{0} ≥ 0, λ ≥ 0.

Let μ̃ ∈ m̂s(C). Then there exists h ∈ ℚ_{+} \ {0} such that hμ̃ ∈ lrf(C) describes a linear ranking function f for C.

The inequalities (A A′) ⟨x, x′⟩ ≤ b define a polyhedron; according to the affine form of Farkas’ lemma, a function is non-negative on this polyhedron, i.e., throughout the loop, if and only if it is a positive affine combination of the column vectors (A A′) ⟨x, x′⟩ ≤ b. In particular this holds for the ranking function f and its two properties: f(x) ≥ 0 and f(x) − f(x′) ≥ 1.

Hence there exist two non-negative rational vectors λ_{1} and λ_{2} and two nonnegative numbers λ_{0,1} and λ_{0,2} such that

f(x) = λ_{0,1} + λ1^{T} (− (A A′) ⟨x, x′⟩ + b)

and

f(x) − f(x′) − 1 = λ_{0,2} + λ2^{T} (− (A A′) ⟨x, x′⟩ + b) .

Replacing f(x) by hμx + hμ_{0}, we get two equalities —one for the part containing variables and one for the remaining part— for each expression. After simplification we obtain the following equalities:

−λ1^{T} (A A′) ⟨x, x′⟩ = hμx  (33a)
−λ2^{T} (A A′) ⟨x, x′⟩ = hμx − hμx′  (33b)
−λ2^{T}b = 1 + λ_{0,2}  (33c)

From (33a) and (33b) we obtain λ1^{T}A = −hμ^{T}, λ1^{T}A′ = 0, λ2^{T}A = −hμ^{T} and λ2^{T}A′ = hμ^{T}. We can rewrite it as 0 = λ1^{T}A′ = (λ1^{T} − λ2^{T})A = λ2^{T}(A + A′). From (33c) we deduce λ2^{T}b < 0.

The four conditions (27) to prove termination by [8] are thus satisfied.

The combination of Theorems 4.4 and 6.1 gives:

Theorem 6.3. Let C be the binary CLP(ℚ) clause p(x̄) :− c[x̄, x̄′], p(x̄′), where p is an n-ary predicate and c[x̄, x̄′] is a linear satisfiable constraint. Let l̂rf(C) be the set of (positive multiples of) linear ranking functions for C, m̂s(C) be the set of (positive multiples of) solutions of the Mesnard and Serebrenik system (20) projected onto μ and pr(C) be the set of the ranking function coefficients obtained through the method of Podelski and Rybalchenko, that is,

l̂rf(C) := \{ kμ̃ ∈ ℚ^{n+1} | \begin{matrix}∀x̄, x̄′ ∈ ℚ^{n} : c[x̄, x̄′] ⟹ \\ ∑_{i=1}^{n} μ_{i}x_{i} − ∑_{i=1}^{n} μ_{i}x′_{i} ≥ 1 \\ ∧ μ_{0} + ∑_{i=1}^{n} μ_{i}x_{i} ≥ 0, \\ k ∈ ℚ_{+} \ \{0\} \end{matrix} \},

m̂s(C) := \{ kμ̃ ∈ ℚ^{n+1} | \begin{matrix}⟨y, μ⟩ and ⟨z, μ̃⟩ are solutions of (20), \\ k ∈ ℚ_{+} \ \{0\} \end{matrix} \},
pr(C) := { ⟨λ2^{T}A′, λ1^{T}b⟩ ∈ ℚ^{n+1} | ⟨λ_{1}, λ_{2}⟩ is a solution of (27) },

where c[x̄, x̄′] is equivalent to (A A′) ⟨x, x′⟩ ≤ b or to A_{c}⟨x, x′⟩ ≥ b_{c}, respectively. Then l̂rf(C) = m̂s(C) = pr(C).

6.2. Worst-Case Complexity Using the Simplex Algorithm

The computationally most expensive component in both methods is the resolution of a linear optimization problem that can always be expressed in the standard form

minimize c^{T}x
subject to Ax = b
x ≥ 0

by applying well known transformations: inequalities and unconstrained (i.e., not subject to lower or upper bounds) variables can be replaced and the resulting equivalent problem in standard form has one more variable for each inequality or unconstrained variable appearing in the original problem.

The most common way to solve this linear optimization problems involves using the simplex algorithm [33], an iterative algorithm that requires \binom{e+u}{e} pivoting steps in the worst-case scenario, where e and u denote the number of equalities in A and unknowns in x respectively.

For a simple linear loop of m inequalities over n variables, Podelski and Rybalchenko require to solve a linear problem in standard form having 3n equalities over 2m variables (the opposite of the expression appearing in (27d) can be used as the quantity to be minimized); this gives a worst-case complexity of \binom{3n+2m}{3n} pivoting steps, corresponding, by Stirling’s formula, to an exponential complexity of exponent 3n + 2m approximately.

If the alternative formalization of the Podelski and Rybalchenko method is adopted for the same loop, then we will have the same m constraints as above for the ‘⋆’ invariant, while the ‘▼’ invariant will be described by other ℓ constraints. If redundant constraints are removed, we will have ℓ ≤ m. Hence, the alternative approach will result in a linear programming problem having 2n equalities over m + 2ℓ variables. Hence, the worst-case number of pivoting steps will be an exponential of exponent approximately 2n + m + 2ℓ.

For the same simple linear loop, Serebrenik and Mesnard require the resolution of two linear problems, that can be rewritten to contain 2n equalities over m + n variables (with n unconstrained variables) and 2n + 1 equalities over (m + 2) + (n + 1) variables (with n + 1 unconstrained variables), respectively. They can then be merged to generate a single linear problem of 4n + 1 equalities over m + (m + 2) + (n + 1) variables, n + 1 of which unconstrained, and an extra inequality replacing one of the two objective functions. In the end, we get a linear problem in standard form with 4n + 2 equalities over 2m + 2n + 5 variables. This means a worst-case complexity of \binom{6n+2m+7}{4n+2} pivoting steps and an exponential complexity of exponent 6n + 2m approximately.

So the method proposed by Podelski and Rybalchenko has, in general, a lower worst-case complexity than the one proposed by Mesnard and Serebrenik, if the single linear problem approach is chosen. The comparison of the two alternative implementation approaches for the Podelski and Rybalchenko method depends on the relations between quantities n, m and ℓ. On the one hand, if ℓ is significantly smaller than m, then the alternative approach could result in an efficiency improvement. On the other hand, if the number of constraints is much higher than the number of variables, then the original implementation approach should be preferred. Note that the need for two loop invariants instead of a single one should not be seen as a big practical problem: in fact, most analysis frameworks will provide the ‘▼’ invariant as the original input to the termination analysis tool, which will then use it to compute the ‘⋆’ invariant (via the abstract execution of a single iteration of the loop); that is, the computational cost for the ‘▼’ invariant is implicitly paid anyway.

It is well known, though, that the worst-case scenario for the simplex algorithm is extremely uncommon in practice. An average complexity analysis and, more recently, a smoothed complexity analysis [34] have been carried out on the simplex algorithm and showed why it usually takes polynomial time. Besides the theoretical studies, several experimental evaluations of implementations of the simplex algorithm reported that the average number of pivoting steps seems to grow linearly with the sum e + u of the number of equalities and unknowns of the problem. Therefore, for a more informative and meaningful comparison, the next section presents an experimental evaluation of the methods on a representative set of while loops.

7. Implementation and Experimental Evaluation

The Parma Polyhedra Library (PPL) is a free software, professional library for the handling of numeric approximations targeted at static analysis and computer-aided verification of hardware and software systems [6, 35]. The PPL, which features several unique innovations [36, 37, 38, 39, 40], is employed by numerous projects in this field, most notably by GCC, the GNU Compiler Collection, probably the most widely used suite of compilers.^{15}

As an integral part of the overall project to which the present paper belongs —whose aim is to make the technology of the automatic synthesis of linear ranking functions thoroughly explained and generally available—, we have extended the PPL with all the methods discussed in the present paper. Previously, only a rather limited demo version of RankFinder was available, only in x86/Linux binary format, implementing the method by Podelski and Rybalchenko.^{16} In contrast, the PPL implementation is completely general and available, both in source and binary formats, with high-level interfaces to C, C^{++}, Java, OCaml and six different Prolog systems.

For each of the methods —Mesnard and Serebrenik (MS) or Podelski and Rybalchenko (PR)—, for each of the two possibilities to encode the input — either the single ⋆ invariant of (22) in Section 4.3, or the two ▼ and ⋆ invariants of (23) in Section 5—, for each numerical abstractions supported by the PPL —including (not necessarily closed) convex polyhedra, bounded-difference shapes and octagonal shapes—, the PPL provides three distinct functionalities to investigate termination of the loop being analyzed:

15See "http://bugseng.com/products/ppl for more information.

16See http://www7.in.tum.de/~rybal/rankfinder/, last checked on August 18th, 2011.

1. a Boolean termination test;

2. a Boolean termination test that, in addition, returns the coefficients of one (not further specified) affine ranking function;

3. a function returning a convex polyhedron that encodes the space of all affine ranking functions.

In addition, using the MS method and for each input method, the PPL provides

4. a function returning two convex polyhedra that encode the space of all decreasing functions (also known as quasi-ranking functions) and all bounded functions, respectively, for use in conditional termination analysis.

We have evaluated the performance of the new algorithms implemented in the PPL using the termination analyzer built into Julia, a state-of-the-art analyzer for Java bytecode [41]. We have thus taken several Java programs in the Julia test suite and, using Julia, we have extracted the constraint systems that characterize the loops in the program that Julia cannot quickly resolve with syntax-based heuristics. This extraction phase allowed us to measure the performance of the methods described in the present paper, factoring out the time spent by Julia in all the analyses (nullness, sharing, path-length, unfolding, …) that allow to obtain such constraint systems.

We first tested the performance (and correctness) of the new PPL implementation with the implementation of the MS method, based on CLP(Q), previously used by Julia and with the implementation of PR, still based on CLP(Q), provided by the demo version of RankFinder. The reason we did this comparison is that, while we know that the infinite precision implementation of the simplex algorithm available in the PPL performs better than its direct competitors [6, Section 4, Table 3],^{17} we know there is much room for improvement: it could have been the case that the constraint solver employed in modern CLP systems made our implementation useless. The result was quite satisfactory: the PPL implementation is one to two orders of magnitude faster over the considered benchmark suite.

The benchmark programs are: CaffeineMark, from Pendragon Software Corporation, measures the speed of Java; JLex is a lexical analyzer generator developed by Elliot Berk and C. Scott Ananian; JavaCC is a parser generator from Sun Microsystems; Java_CUP is a parser generator developed by Scott Hudson, Frank Flannery and C. Scott Ananian; Jess is a rule engine written by Ernest Friedman-Hill; Kitten is a didactic compiler for a simple imperative object-oriented language written by Fausto Spoto; NQueens is a solver of the n-queens problem which includes a library for binary decision diagrams; Raytracer is a ray-tracing program; Termination is a JAR file containing all

17I.e., Cassowary (http://www.cs.washington.edu/research/constraints/cassowary/)

and Wallaroo (http://sourceforge.net/projects/wallaroo/). While GLPK, the GNU Linear Programming Toolkit (http://www.gnu.org/software/glpk/) includes a solver that is termed “exact,” it still depends critically on floating point computations; moreover, it has not yet been made available in the public interface.

Table 1: Benchmarks used in the experimental evaluation

benchmark  loops  n  \overline{n}  σ_{n}  m  \overline{m}  σ_{m}
CaffeineMark 151 [1,9]  6.0 1.3 [2,26] 17.  3.8
JLex  467 [1,14] 7.2 2.5 [2,45] 17.  6.7
JavaCC  136 [1,14] 8.6 4.1 [1,45] 22. 12.
Java_CUP  29 [2,14] 8.3 4.3 [5,45] 23. 13.
Jess  151 [1,9]  6.0 1.3 [2,26] 17.  3.8
Kitten  1484 [1,15] 11.  3.6 [2,45] 29. 10.
NQueens  359 [1,14] 6.3 3.6 [2,45] 17. 10.
Raytracer  8 [2,9]  4.5 2.7 [5,26] 11.  7.8
Termination  121 [1,9]  4.2 3.5 [2,27] 12.  9.9

Table 2: MS vs PR: CPU time in seconds

term. test  one r. f.  all r. f.
benchmark  MS  PR  MS  PR  MS  PR
CaffeineMark 0.42 0.26 0.43 0.25 0.31 0.34
JLex  1.62 0.83 1.64 0.84 1.17 1.14
JavaCC  0.86 0.43 0.87 0.45 0.67 0.65
Java_CUP  0.35 0.14 0.35 0.14 0.29 0.22
Jess  0.42 0.26 0.43 0.26 0.29 0.34
Kitten 11.8 6.87 11.9 6.84 8.41 10.2
NQueens  1.43 0.76 1.44 0.74 0.99 1.03
Raytracer  0.04 0.03 0.04 0.03 0.03 0.03
Termination  0.25 0.15 0.25 0.15 0.18 0.21

the programs of [41, Figure 16]. In Table 1 we report, for each benchmark, the number of loops for which termination was investigated, the interval, mean and standard deviations —with two significant figures— of the quantities n (number of variables) and m (number of constraints) that characterize those loops.

The results of the CPU-time comparison between the MS and PR methods are reported in Table 2. Measurements took place on a GNU/Linux system equipped with an Intel Core 2 Quad CPU Q9400 at 2.66 GHz and 8 Gbytes of main memory; a single core was used and the maximum resident set size over the entire set of tests was slightly above 53 Mbytes. From these we can conclude that the difference in performance between the two methods is rather limited. The PR method is more efficient on the problem of semi-deciding termination, with or without the computation of a witness ranking function, while the MS method is superior on the problem of computing the space of all affine ranking functions.

We also present, in Table 3, the precision results. For each benchmark, along with the total number of loops, we have the number of loops for which termination is decided positively, either with the MS or the PR method (column ‘term’); the remaining loops are divided, using the MS method, between those that admit a linear decreasing function (column ‘w/ d.f.’) and those who do

Table 3: Precision results and application to conditional termination

benchmark  loops term w/ d.f. w/o d.f.

            CaffeineMark 151 149  0  2

JLex  467 453  3  11
JavaCC  136 120  4  12

            Java_CUP  29  27  0  2

Jess  151 149  0  2
Kitten  1484 1454  3  27
NQueens  359 271  4  84

            Raytracer  8  6  0  2
            Termination  121 119  0  2

not (column ‘w/o d.f.’). It can be seen that the percentage of loops for which termination is decided positively ranges from 75% to 99%, depending on the benchmark. This means that we are conducting the experimental evaluation with a termination analyzer, Julia, whose analysis algorithms —though certainly improvable— very often provide enough information for termination analysis. This is crucial for the meaningfulness of the experimental evaluation presented in this section.

8. Conclusions

Linear ranking functions play a crucial role in termination analysis, as the termination of many programs can be decided by the existence on one such function. In this paper we have addressed the topic of the automatic synthesis of linear ranking functions with the aim of clarifying its origins, thoroughly explaining the underlying theory, and presenting new, efficient implementations that are being made available to the general public.

In particular, we have introduced, in general terms independent from any programming paradigm, the problem of automatic termination analysis of individual loops —to which more general control flows can be reconducted— and its solution technique based on the synthesis of ranking functions.

We have then presented and generalized a technique originally due to Sohn and Van Gelder, that was virtually unknown outside the logic programming field despite its general applicability and its relative completeness (given a linear constraint system approximating the behavior of a loop, if a linear ranking function exists for that system, then the method will find it). This method, due to its ability to characterize the spaces of all the linear decreasing functions and all the linear bounded functions, is also immediately applicable to conditional termination analysis; this theme is an excellent candidate for future work.

We have also presented and, for the first time, fully justified, a more recent technique by Podelski and Rybalchenko. For this we also present an alternative formulation that can lead to efficiency improvements.

We have compared the two methods, first proving their equivalence — thus obtaining an independent confirmation on their correctness and relative completeness— and then studying their worst-case complexity.

Finally, we have presented the implementation of all the techniques described in the paper recently included in the Parma Polyhedra Library, along with an experimental evaluation covering both the efficiency and the precision of the analysis.

Acknowledgments

We would like to express our gratitude to Amir Ben-Amram, David Merchat, Andrey Rybalchenko, Alessandro Zaccagnini and the anonymous reviewers for their useful comments and suggestions. The connection between the MS method and conditional termination analysis was indicated to us during a discussion with Samir Genaim.

References

[1] J. C. Lagarias, The 3x+1 problem and its generalizations, American Mathematical Monthly 92 (1985) 3–23.

[2] E. W. Dijkstra, A Discipline of Programming, Prentice Hall, 1976.

[3] R. Bagnara, F. Mesnard, A. Pescetti, E. Zaffanella, The automatic synthesis of linear ranking functions: The complete unabridged version,

arXiv:cs.PL/1004.0944v2, available from http://arxiv.org/. (2012).

[4] P. Cousot, R. Cousot, ‘A la Floyd’ induction principles for proving inevitability properties of programs, in: M. Nivat, J. C. Reynolds (Eds.), Algebraic Methods in Semantics, Cambridge University Press, 1985, pp. 277–312.

[5] R. W. Floyd, Assigning meanings to programs, in: J. T. Schwartz (Ed.), Mathematical Aspects of Computer Science, Vol. 19 of Proceedings of Symposia in Applied Mathematics, American Mathematical Society, Providence RI, New York City, USA, 1967, pp. 19–32.

[6] R. Bagnara, P. M. Hill, E. Zaffanella, The Parma Polyhedra Library: Toward a complete set of numerical abstractions for the analysis and verification of hardware and software systems, Science of Computer Programming 72 (1–2) (2008) 3–21.

[7] K. Sohn, A. Van Gelder, Termination detection in logic programs using argument sizes (extended abstract), in: Proceedings of the Tenth ACM SIGACT-SIGMOD-SIGART Symposium on Principles of Database Systems, ACM, Association for Computing Machinery, Denver, Colorado, United States, 1991, pp. 216–226.

[8] A. Podelski, A. Rybalchenko, A complete method for the synthesis of linear ranking functions, in: B. Steffen, G. Levi (Eds.), Verification, Model Checking and Abstract Interpretation: Proceedings of the 5th International Conference, VMCAI 2004, Vol. 2937 of Lecture Notes in Computer Science, Springer-Verlag, Berlin, Venice, Italy, 2004, pp. 239–251.

[9] M. A. Colón, H. B. Sipma, Practical methods for proving program termination, in: E. Brinksma, K. G. Larsen (Eds.), Computer Aided Verification: Proceedings of the 14th International Conference (CAV 2002), Vol. 2404 of Lecture Notes in Computer Science, Springer-Verlag, Berlin, Copenhagen, Denmark, 2002, pp. 442–454.

[10] K. Sohn, Automated termination analysis for logic programs, Ph.D. thesis, University of California Santa Cruz, Santa Cruz, CA, USA (1993).

[11] J. Fischer, Termination analysis for Mercury using convex constraints, Honours report, Department of Computer Science and Software Engineering, The University of Melbourne, Melbourne, Australia (Aug. 2002).

[12] C. Speirs, Z. Somogyi, H. Søndergaard, Termination analysis for Mercury, in: P. Van Hentenryck (Ed.), Static Analysis: Proceedings of the 4th International Symposium, Vol. 1302 of Lecture Notes in Computer Science, Springer-Verlag, Berlin, Paris, France, 1997, pp. 157–171.

[13] A. R. Bradley, Z. Manna, H. B. Sipma, Linear ranking with reachability, in: Computer Aided Verification: Proceedings of the 15th International Conference, Vol. 3576 of Lecture Notes in Computer Science, Springer-Verlag, Berlin, Edinburgh, Scotland, UK, 2005, pp. 491–504.

[14] K. Verschaetse, D. De Schreye, Deriving termination proofs for logic programs, using abstract procedures, in: K. Furukawa (Ed.), Proceedings of the 8th International Conference on Logic Programming, The MIT Press, Paris, France, 1991, pp. 301–315.

[15] M. Codish, C. Taboch, A semantic basis for the termination analysis of logic programs, Journal of Logic Programming 41 (1) (1999) 103–123.

[16] F. Mesnard, Inferring left-terminating classes of queries for constraint logic programs by means of approximations, in: M. J. Maher (Ed.), Logic Programming: Proceedings of the Joint International Conference and Symposium on Logic Programming, MIT Press Series in Logic Programming, The MIT Press, Bonn, Germany, 1996, pp. 7–21.

[17] F. Mesnard, A. Serebrenik, A polynomial-time decidable class of terminating binary constraint logic programs, Tech. Rep. 05-11, Université de la Réunion (2005).

[18] F. Mesnard, A. Serebrenik, Recurrence with affine level mappings is P-time decidable for CLP(R), Theory and Practice of Logic Programming 8 (1) (2008) 111–119.

[19] A. Podelski, A. Rybalchenko, Transition invariants, in: Logic in Computer Science, Proceedings of the 19th IEEE Symposium, LICS 2004, IEEE Computer Society, Turku, Finland, 2004, pp. 32–41.

[20] B. Cook, A. Podelski, A. Rybalchenko, Abstraction refinement for termination, in: C. Hankin, I. Siveroni (Eds.), Static Analysis: Proceedings of the 12th International Symposium, Vol. 3672 of Lecture Notes in Computer Science, Springer-Verlag, Berlin, London, UK, 2005, pp. 87–101.

[21] M. T. Nguyen, D. De Schreye, Polynomial interpretations as a basis for termination analysis of logic programs, in: M. Gabbrielli, G. Gupta (Eds.), Proceedings of the 21st International Conference on Logic Programming, no. 3668 in Lecture Notes in Computer Science, Springer-Verlag, Berlin, Sitges, Spain, 2005, pp. 311–325.

[22] D. S. Lankford, A finite termination algorithm, Internal memo, Southwestern University, Georgetown, TX, USA (1976).

[23] H. Zantema, Termination of term rewriting, Tech. Rep. UU-CS-2000-04, Department of Computer Science, Universiteit Utrecht, Utrecht, The Netherlands (2000).

[24] P. Cousot, Proving program invariance and termination by parametric abstraction, lagrangian relaxation and semidefinite programming, in: R. Cousot (Ed.), Verification, Model Checking and Abstract Interpretation: Proceedings of the 6th International Conference (VMCAI 2005), Vol. 3385 of Lecture Notes in Computer Science, Springer-Verlag, Berlin, Paris, France, 2005, pp. 1–24.

[25] M. A. Colón, H. B. Sipma, Synthesis of linear ranking functions, in: T. Margaria, W. Yi (Eds.), Tools and Algorithms for Construction and Analysis of Systems, 7th International Conference, TACAS 2001, Vol. 2031 of Lecture Notes in Computer Science, Springer-Verlag, Berlin, Genova, Italy, 2001, pp. 67–81.

[26] J. Jaffar, M. J. Maher, Constraint logic programming: A survey, Journal of Logic Programming 19&20 (1994) 503–582.

[27] M. R. Garey, D. S. Johnson, Computers and Intractability; A Guide to the Theory of NP-Completeness, W. H. Freeman & Co., New York, NY, 1990.

[28] A. Schrijver, Theory of linear and integer programming, Wiley, Chichester, New York, 1986.

[29] P. Cousot, N. Halbwachs, Automatic discovery of linear restraints among variables of a program, in: Conference Record of the Fifth Annual ACM Symposium on Principles of Programming Languages, ACM Press, Tucson, Arizona, 1978, pp. 84–96.

[30] B. Cook, S. Gulwani, T. Lev-Ami, A. Rybalchenko, M. Sagiv, Proving conditional termination, in: A. Gupta, S. Malik (Eds.), Computer Aided Verification: Proceedings of the 20th International Conference (CAV 2008), Vol. 5123 of Lecture Notes in Computer Science, Springer-Verlag, Berlin, Princeton, NJ, USA, 2008, pp. 328–340.

[31] A. Rybalchenko, Temporal verification with transition invariants, Ph.D. thesis, Universität des Saarlandes, Saarbrücken, Germany (2004).

[32] A. Serebrenik, F. Mesnard, On termination of binary CLP programs, in: S. Etalle (Ed.), Logic Based Program Synthesis and Transformation: 14th International Symposium, Revised Selected Papers, no. 3573 in Lecture Notes in Computer Science, Springer-Verlag, Berlin, Verona, Italy, 2005, pp. 231–244.

[33] C. H. Papadimitriou, K. Steiglitz, Combinatorial Optimization: Algorithms and Complexity, Prentice Hall, Englewood Cliffs, NJ, USA, 1982.

[34] D. A. Spielman, S.-H. Teng, Smoothed analysis: Why the simplex algorithm usually takes polynomial time, Journal of the ACM 51 (2004) 385– 463.

[35] R. Bagnara, P. M. Hill, E. Zaffanella, Applications of polyhedral computations to the analysis and verification of hardware and software systems, Theoretical Computer Science 410 (46) (2009) 4672–4691.

[36] R. Bagnara, P. M. Hill, E. Ricci, E. Zaffanella, Precise widening operators for convex polyhedra, Science of Computer Programming 58 (1–2) (2005) 28–56.

[37] R. Bagnara, P. M. Hill, E. Zaffanella, Not necessarily closed convex polyhedra and the double description method, Formal Aspects of Computing 17 (2) (2005) 222–257.

[38] R. Bagnara, P. M. Hill, E. Zaffanella, Widening operators for powerset domains, Software Tools for Technology Transfer 8 (4/5) (2006) 449–466.

[39] R. Bagnara, P. M. Hill, E. Zaffanella, Weakly-relational shapes for numeric abstractions: Improved algorithms and proofs of correctness, Formal Methods in System Design 35 (3) (2009) 279–323.

[40] R. Bagnara, P. M. Hill, E. Zaffanella, Exact join detection for convex polyhedra and other numerical abstractions, Computational Geometry: Theory and Applications 43 (5) (2010) 453–473.

[41] F. Spoto, F. Mesnard, É. Payet, A termination analyzer for Java bytecode based on path-length, ACM Transactions on Programming Languages and Systems 32 (3).