Porous Invariants

We introduce the notion of porous invariants for multipath (or branching/nondeterministic) affine loops over the integers; these invariants are not necessarily convex, and can in fact contain infinitely many 'holes'. Nevertheless, we show that in many cases such invariants can be automatically synthesised, and moreover can be used to settle (non-)reachability questions for various interesting classes of affine loops and target sets.


Introduction
We consider the reachability problem for multipath (or branching) affine loops over the integers, or equivalently for nondeterministic integer linear dynamical systems. A (deterministic) integer linear dynamical system consists of an update matrix M ∈ Z d×d together with an initial point x (0) ∈ Z d . We associate to such a system its infinite orbit (x (i) ) consisting of the sequence of reachable points defined by the rule x (i+1) = M x (i) . The reachability question then asks, given a target set Y , whether the orbit ever meets Y , i.e., whether there exists some time i such that x (i) ∈ Y . The nondeterministic reachability question allows the linear update map to be chosen at each step from a fixed finite collection of matrices.
When the orbit does eventually hit the target, one can easily substantiate this by exhibiting the relevant finite prefix. However, establishing non-reachability is intrinsically more difficult, since the orbit consists of an infinite sequence of points. One requires some sort of finitary certificate, which must be a relatively simple object that can be inspected and which provides a proof that the set Y is indeed unreachable. Typically, such a certificate will consist of an overapproximation I of the set R of reachable points, in such a manner that one can check both that Y ∩ I = ∅ and R ⊆ I; such a set I is called an invariant.
Formally we study the following problem for inductive invariants: Meta Problem 1. Consider a system with update functions f 1 , . . . , f n . A set I is an inductive invariant if f i (I) ⊆ I for all i. Given a reachability query (x, Y ) we search for a separating inductive invariant I such that x ∈ I and Y ∩ I = ∅.
Meta Problem 1 is parametrised by the type of invariants and targets that are considered; that is, what are the classes of allowable invariant sets I and target sets Y , or equivalently how are such sets allowed to be expressed.
Fixing a particular invariant and target domain, a reachability query has three possible scenarios: (1) the instance is reachable, (2) the instance is unreachable and a separating invariant from the domain exists, or (3) the instance is unreachable but no separating invariant exists. Ideally, one would wish to provide a sufficiently expressive invariant domain so that the latter case does not occur, whilst keeping the resulting invariants as simple as possible and computable. For some classes of systems, it is known that distinguishing reachability (1) from unreachability (2,3) is undecidable; it can also happen that determining whether a separating invariant exists (i.e., distinguishing (2) from (3)) is undecidable.
We note that the existence of strongest inductive invariants 3 is a desirable property for an invariant domain-when strongest invariants exist (and can be computed), separating (2) from (1,3) is easy: compute the strongest invariant, and check whether it excludes the target state or not; if so, then you are done, and if not, no other invariant (from that class) can possibly do the trick either. However, unless (3) is excluded, computing the strongest invariant does not necessarily imply that reachability is decidable. Unfortunately, strongest invariants are not always guaranteed to exist for a particular invariant domain, although some separating inductive invariant may still exist for every target (or indeed may not).
In prior work from the literature, typical classes of invariants are usually convex, or finite unions of convex sets. In this paper we consider certain classes of invariants that can have infinitely many 'holes' (albeit in a structured and regular way); we call such sets porous invariants. These invariants can be represented via Presburger arithmetic 4 . We shall work instead with the equivalent formulation of semi-linear sets, generalising ultimately periodic sets to higher dimensions, as finite unions of linear sets of the form {b + p 1 N + · · · + p m N} (by which we mean {b + a 1 p 1 + · · · + a m p m | a 1 , . . . , a m ∈ N}, see Definition 2).
Let us first consider a motivating example: Example 1 (Hofstadter's MU Puzzle [7]). Consider the following term-rewriting puzzle over alphabet {M, U, I}. Start with the word M I, and by applying the following grammar rules (where y and z stand for arbitrary words over our alphabet), we ask whether the word M U can ever be reached.
The answer is no. One way to establish this is to keep track of the number of occurrences of the letter 'I' in the words that can be produced, and observe that this number (call it x) will always be congruent to either 1 or 2 modulo 3. In other words, it is not possible to reach the set {x | x ≡ 0 mod 3}. Indeed, Rules 2 and 3 are the only rules that affect the number of I's, and can be described by the system dynamics x → 2x and x → x − 3. Hence the MU Puzzle can be viewed as a one-dimensional system with two affine updates, 5 or a twodimensional system with two linear updates. 6 The set {1 + 3Z} ∪ {2 + 3Z} is an inductive invariant, and we wish to synthesise this. The problem can be rephrased as a safety property of the following multipath loop, verifying that the 'bad' state x = 0 is never reached, or equivalently that the above loop can never halt, regardless of the nondeterministic choices made.
The MU Puzzle was presented as a challenge for algorithmic verification in [4]; the tools considered in that paper (and elsewhere, to the best of our knowledge) rely upon the manual provision of an abstract invariant template. Our approach is to find the invariant fully automatically (although one must still abstract from the MU Puzzle the correct formulation as the program x → 2x || x → x − 3).
Main Contributions. Our focus is on the automatic generation of porous invariants for multipath affine loops over the integers, or equivalently nondeterministic integer linear dynamical systems.
-We first consider targets consisting of a single vector (or 'point targets'), and present the classes of invariants and systems for which invariants can and cannot be automatically computed for the reachability question. A summary of the results for linear and semi-linear invariants for these targets is given in Table 1. For completeness we also consider R, R + -(semi)-linear sets, where we complete the picture from prior work by showing that strongest R-semilinear invariants are computable.
• We establish the existence of strongest Z-linear invariants, and show that they can be found algorithmically (Theorem 2). These invariants may or may not separate the target under consideration. • If a Z-linear invariant is not separating, we may instead look for an Nsemi-linear invariant (which generalises both Z-semi-linear and N-linear invariants), and we show that such an invariant can always be found for any unreachable point target when dealing with deterministic integer linear dynamical systems (Theorem 4). 5 One-dimensional affine updates are functions of the form f (x) = ax + b. 6 a b 0 1 • However, for nondeterministic integer linear dynamical systems, computing an N-semi-linear invariants is an undecidable problem in arbitrary dimension (Theorem 5). Nevertheless we show how such invariants can be constructed in a low-dimensional setting, in particular for affine updates in one dimension (Theorem 6). As an immediate consequence, this establishes that the multipath loop associated with the MU Puzzle belongs to a class of programs for which we can automatically synthesise N-semi-linear invariants. -For full-dimensional 7 Z-linear targets we show that reachability is decidable, and, in the case of unreachability that a Z-semi-linear invariant can always be exhibited as a certificate (Theorem 3). If the target is not fulldimensional then the reachability problem is Skolem-hard and undecidable for deterministic and nondeterministic systems respectively. -In Section 6 we present our tool porous which handles one-dimensional affine systems for both point and Z-linear targets, solving both the reachability problem and producing invariants. Inter alia, this allows one to handle the multipath loop derived from the MU Puzzle in fully automated manner.

Related Work
The reachability problem (in arbitrary dimension) for loops with a single affine update, or equivalently for deterministic linear dynamical systems, is decidable in polynomial time for point targets (that is Y = {y}), as shown by Kannan and Lipton [16]. However for nondeterministic systems (where the update matrix is chosen nondeterministically from a finite set at each time step), reachability is undecidable, by reduction from the matrix semigroup membership problem [22]. In particular this entails that for unreachable nondeterministic instances we cannot hope always to be able to compute a separating invariant. In some cases we may compute the strongest invariant (which may suffice if this invariant happens to be separating for the given reachability query), or we may compute an invariant in sub-cases for which reachability is decidable (for example in low dimensions). For some classes of invariants, it is also undecidable whether an invariant exists (e.g., polyhedral invariants [8]).
Various types of invariants have been studied for linear dynamical systems, including polyhedra [23,8], algebraic [15], and o-minimal [1] invariants. For certain classes of invariants (e.g., algebraic [15]), it is decidable whether a separating invariant exists, notwithstanding the reachability problem being undecidable. Other works (e.g., [5]) use heuristic approaches to generate invariants, without aiming for any sort of completeness.
Kincaid, Breck, Cyphert and Reps [18] study loops with linear updates, studying the closed forms for the variables to prove safety and termination properties. Such closed forms, when expressible in certain arithmetic theories, can be interpreted as another type of invariant and can be used to over-approximate the reachable sets. The work is restricted to a single update function (deterministic loops) and places additional constraints on the updates to bring the closed forms into appropriate theories.
Bozga, Iosif and Konecný's FLATA tool [2] considers affine functions in arbitrary dimension. However, it is restricted to affine functions with finite monoids; in our one-dimensional case this would correspond to limiting oneself to counterlike functions of the form f (x) = x + b.
Finkel, Göller and Haase [9], extending Fremont [10], show that reachability in a single dimension is PSPACE-complete for polynomial update functions (and allowing states can be used to control the sequences of updates which can be applied). The affine functions (and single-state restriction) we consider are a special case, but we focus on producing invariants to disprove reachability.
Other tools, e.g., AProVE [11] and Büchi Automizer [14] may (dis-)prove termination/reachability on all branches, but may not be able to prove termination/reachability on some branch.
Inductive invariants specified in Presburger arithmetic have been used to disprove reachability in vector addition systems [20]. A generalisation, 'almost semi-linear sets' [21] are also non-convex and can capture exactly the reachable points of vector addition systems. Our nondeterministic linear dynamical systems can be seen as vector addition systems over Z extended with affine updates (rather than only additive updates).

Preliminaries
We denote by Z the integers and N the non-negative integers. We say that  A point y is reachable if there exists m ∈ N and B 1 , . . . , B m such that B 1 · · · B m x (0) = y and B i ∈ {M 1 , . . . , M k } for all 1 ≤ i ≤ m.
The reachability set O ⊆ Z d of an LDS is the set of reachable points.
Definition 2 (K-(semi)-linear sets). A linear set L is defined by a base vector b ∈ Z d and period vectors p 1 , . . . , p d ∈ Z d such that For convenience we often write {b + p 1 K + · · · + p d K} for L. A set is semi-linear if it is the finite union of linear sets.
N-semi-linear sets are precisely those definable in Presburger arithmetic (FO(Z, +, ≤)) [12]. However, we can also consider Z-semi-linear sets (corresponding to FO(Z, +) without order), and the real counterparts (R and R + ). Note that even if K = N we still allow p i ∈ Z d .

Definition 3. Given an integer linear dynamical system
Note in particular that every inductive invariant contains the reachability set (O ⊆ I). We are interested in the following problem: In our setting, we are interested in classes D of invariants that are linear, or semi-linear. When a separating inductive invariant I exists, we also wish to compute it. Since (semi)-linear invariants are enumerable, the decision problem is, in theory, sufficient-although all of our proofs are constructive.

R Invariants: R-linear and R-semi-linear
Before delving into porous invariants, let us consider invariants over the real numbers, i.e., described as R-(semi)-linear sets.
Strongest R-linear invariants are given precisely by the affine hull of the reachability set, and can be computed using Karr's algorithm [17]. Moreover, we will show that strongest R-semi-linear invariants also exist and can be computed by combining techniques for algebraic invariants [15] and R-linear invariants.

R-linear.
Recall that a set L is R-linear if L = {v 0 + v 1 R + · · · + v t R} for some v 0 , . . . , v t ∈ Z d that can be assumed to be linearly-independent 8 without loss of generality (and thus t ≤ d). Given two distinct points of L, every point on the infinite line connecting them must also be in L. Generalising this idea to higher dimensions, given a set S ⊆ R d , let the affine hull be Fix an LDS (x (0) , {M 1 , . . . , M k }) and consider its reachability Then O a is precisely the strongest R-linear invariant. Karr's algorithm [17,26] can be used to compute this strongest invariant in polynomial time. The next lemma follows from Theorem 3.1 of [26]. Let R 0 = x (0) , r 1 , . . . , r d be obtained as per Lemma 1, with d ≤ d. The R-linear invariant of the LDS is the affine span R 0 a , which can be written as the R-semi-linear. Let us now generalise this approach to R-semi-linear sets. The collection of R-semi-linear sets, { m i=1 L i | m ∈ N, L 1 , . . . , L m are R-linear sets}, is closed under finite unions and arbitrary intersections 9 . Thus for any given set X, the smallest R-semi-linear set containing X is simply the intersection of all R-semi-linear sets containing X. Let us denote by X R this smallest R-semi-linear set. We are interested in O R .
Algebraic sets are those that are definable by finite unions and intersections of zeros of polynomials. For example, {(x, y) | xy = 0} describes the lines x = 0 and y = 0. The (real) Zariski closure X z of a set X is the smallest algebraic subset of R d containing the set X. The Zariski closure of the set of reachable points, O z , can be computed algorithmically [15].
An algebraic set A is irreducible if whenever A ⊆ B ∪ C, where B and C are algebraic sets, then we have A ⊆ B or A ⊆ C. Any algebraic set (and in particular a Zariski closure) can be written effectively as a finite union of irreducible sets [3].
Proof. Since A i ⊆ X R = ∪ j L j , and A i is irreducible, we have A i ⊆ L j for some j (as the L j 's are algebraic sets). Since L j is R-linear, and A i a is the smallest , where y ∈ A \ W . Such a point y can always be found using quantifier elimination in the theory of the reals. Each step necessarily increases the dimension, which can occur at most d times, ensuring termination, at which point one has A a = W .

Strongest Z-linear Invariants
Recall that a Z-linear set {q + p 1 Z + · · · + p n Z} is defined by a base vector q ∈ Z d and period vectors p 1 , . . . , p n ∈ Z d . Equivalently, a Z-linear set describes a lattice, i.e., {p 1 Z + · · · + p n Z}, in d-dimensional space, translated to start from q rather than 0. The image of a Z-linear set L = {q + p 1 Z + · · · + p n Z} by a matrix M is the Z-linear set: M (L) = {M q + (M p 1 )Z + · · · + (M p n )Z}. The following lemma asserts that when two points are in a Z-linear set, the direction between these two points can be applied from any reachable point, and hence this direction can be included as a period without altering the set. Proposition 2. Let L = {q + a 1 p 1 + · · · + a n p n | a 1 , . . . , a n ∈ Z} be a Z-linear set. If x, y ∈ L then for all z ∈ L and all a ∈ Z we have z + (y − x)a ∈ L. In particular, we have L = {q + a 1 p 1 + · · · + a n p n + a (y − x) | a 1 , . . . , a n , a ∈ Z}.
Next we show minimality as a straightforward consequence of Proposition 2.
Clearly the vectors p 1 , . . . , p n can be added by Proposition 2 because any two points of L 1 differing by p i guarantees that adding p i does not alter the resulting set. Similarly, t 1 , . . . , t m can also be included. Finally, by Proposition 2, the vector s − q can be included because q and s both belong to L 1 ∪ L 2 .
A d-dimensional lattice can always be defined by at most d vectors; and thus if d is the dimension of the matrices, no more than d period vectors are needed in total. However, Proposition 3 induces a representation which may over-specify the lattice by producing more than d vectors to define the lattice. The Hermite normal form can be used to obtain a basis of the vectors that define the lattice. Consider a lattice L i = {p 1 Z + · · · + p d Z}. The lattice remains the same if p i is swapped with p j , if p i is replaced by −p i , or if p i is replaced by p i + αp j where α is any fixed integer 10 .
These are the unimodular operations. The Hermite normal form of a matrix M is a matrix H such that M = U H, where U is a unimodular matrix (formed by unimodular column operations) and H is lower triangular, non-negative and each row has a unique maximum entry which is on the main diagonal. Such a form always exists, and the columns of H form a basis of the same lattice as the columns of M , because they differ up to unimodular (lattice-preserving) operations. There are many texts on the subject; we refer the reader to the lecture notes of Shmonin [25] for more detailed explanations.
The columns of a matrix in Hermite normal form constitute a unique basis for the lattice (up to additional redundant zero columns). Hence a basis of minimal dimension can be obtained by computing the Hermite normal form of the matrix formed by placing the period vectors into columns.
We now prove the main theorem: Proof (Proof of Theorem 2). We claim that Algorithm 1 returns the strongest Z-linear invariant I. Algorithm 1 proceeds in two phases: -First find a necessary subset L 0 ⊆ I of the invariant having already the same dimension as I. -Then compute a growing sequence L 0 L 1 · · · L m−1 = L m = I, where at each step the algorithm merely increases the density of the attendant sets in order to 'fill in' missing points of the invariant.
Recall the set . Applying M 1 , . . . , M k can only increase the density, but not the dimension.
As each r i and x (0) are in O, by Proposition 2 we can assume that each of the directions (r i − x (0) ) must be represented in any Z-linear set containing O, and we therefore have that L 0 ⊆ I.
In the second phase, we 'fill in' the lattice as required to cover the whole of O. To do this we repeatedly apply the covering procedure of Proposition 3. That is, To keep the number of vectors small, we keep the period vectors of the Z-linear set in Hermite normal form.
The vectors ) form a parallelepiped (hyper-parallelogram) that repeats regularly. There are a finite number of integral points inside this parallelepiped. If new points are added in some step, they are added to every parallelepiped. Thus we can add new points finitely many times before saturating or becoming fixed. The volume of the parallelepiped is bounded above by |p 1 | · · · |p d |.
At each step, the volume of the parallelepiped must at least halve, thus the volume at step t is vol t ≤ |p 1 | · · · |p d |/2 t . The procedure must saturate at or before the volume becomes 1, which occurs after at most log(|p 1 | · · · |p d |) = i log(|p i |) steps. At each step, for efficiency considerations, we convert the Z-linear set into Hermite normal form to retain exactly d period vectors.
Claim (I is the strongest invariant). For every invariant J, we have I ⊆ J.
By induction, let us prove that every invariant J must contain L i . Clearly this is the case for L 0 because all points of R 0 ⊆ O must be in J and every period vectors in L 0 can be present, without loss of generality, thanks to Proposition 2. Assume L i ⊆ J. Then it must be the case that J contains every M j (L i ), as otherwise it would not be an invariant. It therefore follows that J must contain L i+1 , since the latter is the minimal Z-linear set containing L i and M j (L i ) for all j ≤ k. Finally, since I is itself one of the L i 's, we have I ⊆ J as required.
Remark 1. Note that a Z-linear set is not sufficient for the MU puzzle: both 1 and 2 are in the reachability set, thus {1 + 1Z} = Z is the strongest Z-linear invariant.

Extensions of Z-linear sets without strongest invariants
In this section we show that several generalisations of Z-linear domains fail to admit strongest invariants.
Z-semi-linear sets are unions of Z-linear sets, and therefore can include singletons. Consider the deterministic dynamical system starting from point 1 and doubling at each step M = (1, (x → 2x)). This system has reachability set O = 2 k | k ∈ N , which is not even N-semi-linear (our most general class). For this LDS we can construct the invariant 2, 4, 8, ..., 2 k ∪ 2 k+1 p 1 | p 1 ∈ Z for each k. For any proposed strongest Z-semi-linear invariant, one can find a k for which the corresponding invariant is an improvement.
N-linear sets generalise Z-linear sets (observe that Z-linear sets are a proper subclass, since {x + p i Z} can be expressed as {x + (−p i )N + p i N}, but {x + p i N} is clearly not Z-linear). Consider the LDS ((x 1 , x 2 ), ( 0 1 1 0 )), with a reachability set consisting of just two points x = (x 1 , x 2 ) and y = (x 2 , x 1 ). There are two incomparable candidates for the minimal N-linear invariant: {x + (y − x)N} and {y + (x − y)N}. Similarly for R + -linear invariants, the sets {y + (x − y)R + } and {x + (y − x)R + } are incomparable half-lines.

Z-linear targets
We have so far only considered invariants for point targets. We now turn to lattice-like targets, in particular targets specified as full-dimensional Z-linear sets.
Furthermore, for unreachable instances, a Z-semi-linear inductive invariant can be provided.
Theorem 3 requires the targets to be full-dimensional. For nondeterministic systems reachability is undecidable for non-full-dimensional targets (in particular point targets) [22]. However, even for deterministic systems, when Z-linear targets fail to be full-dimensional the reachability problem becomes as hard as the Skolem problem (see, e.g. [24]), for example by choosing as target the set Towards proving Theorem 3, we first show that full-dimensional linear sets can be expressed as 'square' hybrid-linear sets. Hybrid-linear sets are semi-linear sets in which all the components share the same period vectors, and thus differ only in starting position (whereas semi-linear sets allow each component to have distinct period vectors). By square, we mean that all period vectors are the same multiple of standard basis vectors.
Then there is an integral combination of p 1 , . . . , p d such that m i e i is an admissible direction in Y .
By the presence of And therefore Y can be written as b∈B {b We now prove Theorem 3.
Proof (Proof of Theorem 3). Choose m and B as in Lemma 2, so that Y is of the form b∈B {b + me 1 Z + · · · + me d Z}. We build an invariant I of the form If there exists y ∈ B ∩ I then return Reachable. This is because the same sequence of matrices applied to x (0) to produce y ∈ I would, thanks to the modulo step, wind up inside the set {y + me 1 Z + · · · + me d Z}, which is a part of the target.
Otherwise, return Unreachable and I as invariant. By construction, I is indeed an inductive invariant disjoint from the target set.
Remark 2. By the same argument, Theorem 3 extends to a restricted class of Z-semi-linear targets: the finite union of full-dimensional Z-linear sets.

N-semi-linear Invariants
We now consider N-semi-linear invariants, our most general class. N-semi-linear invariants gain expressivity thanks to the 'directions' provided by the period vectors. For example, the only possible Z-semi-linear invariant for the LDS (0, (x → x + 1)) is Z, yet the reachability set, N, is captured exactly by an N-linear invariant. We show that a separating N-semi-linear invariant can always be found for unreachable instances of deterministic integer LDS, although the computed invariant will depend on the target. However, finding invariants is undecidable for nondeterministic systems, at least in high dimension. Nevertheless, we show decidability for the low-dimensional setting of the MU Puzzle-one dimension with affine updates.

Existence of sufficient (but non-minimal) N-semi-linear invariants for point reachability in deterministic LDS
Kannan and Lipton showed decidability of reachability of a point target for deterministic LDS [16]. In this subsection, we establish the following result to provide a separating invariant in unreachability instances. To do so, we will invoke the results from [8] to compute an R + -semi-linear inductive invariant, and then extract from it an N-semi-linear inductive invariant. More precisely, the authors of [8] show how to build polytopic inductive invariants for certain deterministic LDS. Such polytopes are either bounded or are R + -semi-linear sets. In the first case, the polytope contains only finitely many integral points, which can directly be represented via an N-semi-linear set. In the second case, we build an N-semi-linear set containing exactly the set of integral points included in the R + -semi-linear invariant, thanks to the following lemma. Proof (Proof of Theorem 4). We note that every invariant produced in [8] has rational period vectors, as the vectors are given by the difference of successive point in the orbit of the system, and thus Lemma 3 can be applied. The authors of [8] build an inductive invariant in all cases except those for which every eigenvalue of the matrix governing the evolution of the LDS is either 0 or of modulus 1 and at least one of the latter is not a root of unity. This situation however cannot occur in our setting. Indeed, the eigenvalues of an integer matrix are algebraic integers, and an old result of Kronecker [19] asserts that unless all of the eigenvalues are roots of unity, one of them must have modulus strictly greater than 1 (the case in which all eigenvalues are 0 being of course trivial).
This concludes the proof of Theorem 4.

Undecidability of N-semi-linear invariants for nondeterministic LDS
If the enhanced expressivity of N-semi-linear sets allows us always to find an invariant for deterministic LDS, it contributes in turn to making the invariantsynthesis problem undecidable when the LDS is not deterministic. We establish this through a reduction from the infinite Post correspondence problem (ω-PCP) that can be defined in the following way: given m pairs of non-empty words This problem is known to be undecidable when m is at least 8 [13,6].
Theorem 5. The invariant synthesis problem for N-semi-linear sets and linear dynamical systems with at least two matrices of size 91 is undecidable.
Proof (Sketch). We first establish the result in the case of several matrices in low dimension; this can then be transformed in a standard way to two larger matrices (of size 91). The proof is by reduction from the infinite Post correspondence problem. Given an instance of this problem the pair of words corresponding to each sequence of tiles has an integer representation, using base-4 encoding. An important property of our encoding is that the operation of appending a new tile to an existing pair of words can be encoded by matrix multiplication.
Recall that if the instance of ω-PCP is negative, then every generated pair of words will differ at some point. Our encoding is such that this difference of letters creates a difference in their numerical encodings that can be identified with an N-semi-linear invariant. On the other hand, when there is a positive answer to the ω-PCP instance, there can be no N-semi-linear invariant.

Nondeterministic one-dimensional affine updates
The previous section shows that point reachability for nondeterministic LDS is undecidable once there sufficiently many dimensions, motivating an analysis at lower dimensions. The MU Puzzle requires a single dimension with affine updates (or equivalently two dimensions in matrix representation, with the coordinate along the second dimension kept constant). We consider this one-dimensional affine-update case, and therefore, rather than taking matrices as input, we directly work with affine functions of the form f i (x) = a i x + b i . Theorem 6. Given x (0) , y ∈ Z, along with a finite set of functions Moreover, when y is unreachable, an N-semi-linear separating inductive invariant can be algorithmically computed.
We note that decidability of reachability is already known [9,10]. We refine this result by exhibiting an invariant which can be used to disprove reachability. In fact our procedure will produce an N-semi-linear set which can be used to decide reachability, and which, in instances of non-reachability, will be a separating inductive invariant. We have implemented this algorithm into our tool porous, enabling us to efficiently tackle the MU Puzzle as well as its generalisation to arbitrary collections of one-dimensional affine functions. We report on our experiments in Section 6.
We build a case distinction depending on the type of functions that appear: Simplifying assumptions Lemma 4. Without loss of generality, redundant functions are redundant; more precisely, we can reduce the computation of an invariant for a system having redundant functions to finitely many invariant computations for systems devoid of such functions.
Proof. Clearly the identity function has no impact on the reachability set, and so can be removed outright. For any other redundant function, its impact on the reachability set does not depend on when the function is used, and we may therefore assume that it was used in the first step, or equivalently, using an alternative starting point. Hence the invariant-computation problem can be reduced to finitely many instances of the problem over different starting points, with redundant functions removed. Finally, taking the union of the resulting invariants yields an invariant for the original system. Proof and b = c (as f = g) these two functions are opposing.
Two opposing counters. Let us first observe that when there are two opposing counters, we essentially move in either direction by some fixed amount. This will entail that only Z-(semi)-linear invariants can be produced, rather than proper N-(semi)-linear invariants. Therefore, starting with x (0) + dZ ∈ I we can 'saturate' the invariant under construction using the following lemma: Consider the function f (x) = ax + b. If x = y + dk ∈ I, then f (x) = ax + b = ay + adk + b = f (y) + adk ∈ I. Now thanks to the presence of counter h(x) = x + d, by choosing the initial k ∈ Z appropriately and applying h(x) sufficiently many times (say m ∈ N times), one can reach f (x) + adk + dm = f (x) + dn for any desired n ∈ Z.
Without loss of generality if {x + dZ} is in the invariant, then 0 ≤ x < d. We then repeatedly use Lemma 8 to find the required elements of the invariant. Since there are only finitely many residue classes (modulo d), every reachable residue class {c 1 , . . . , c n } can be found by saturation (in at most d steps), yielding invariant {c 1 + dZ} ∪ · · · ∪ {c n + dZ}.
Thanks to Lemma 6, in all remaining cases there is without loss of generality at most one pure inverter. No Counters. If we are not in the preceding case and there are no counters, then there must be growing functions and by Lemma 6, without loss of generality at most one pure inverter. We show that all growing functions increase the modulus outside of some bounded region.
Proof. By the triangle inequality we have: This is the only situation in which the invariant is not exactly the reachability set, and requires us to take an overapproximation. If there is one pure inverter g(x) = −x + d then observe that −C is mapped to C + d and C + d is mapped to −C. Thus intuitively we want to use the interval (−C, C + d). However two problems may occur: (a) since d could be less than 0 then C + d may no longer be growing (under the application of the growing functions), and (b) an inverting growing function only ensures that −C is mapped to a value greater than or equal to C, rather than C+d. Hence, we choose C to ensure that C ± d is still growing by at least |d| (under the application of our growing functions). Let C = max C Non-opposing counters. The only remaining possibility (if there do not exist two opposing counters, and not all functions are growing or pure inverters), is that there are counter-like functions, but they are all counting in the same direction. There may also be a single pure inverter, and possibly some growing functions.
Pick a counter h(x) = x+d to be the reference counter; the choice is arbitrary, but it is convenient to pick a counter with minimal |d|. As a starting point, we have x (0) + dN ⊆ I. Proof. Let r = g(x) + dm for m ∈ Z. We show r ∈ I. Consider x + dn for n ∈ N, then g(x + dn) = −a(x + dn) + b = −ax + b − adn = g(x) − adn. Hence g(x) − adn + dk, n, k ∈ N, is reachable by applying k times the function h(x). Hence for any m ∈ Z there exists k, n ∈ N such that k − na = m, so that r is indeed reachable.
Similarly to the situation with two opposing counters, whenever the invariant contains some Z-linear set, Lemma 8 allows us to saturate amongst the finitely many reachable residue classes.
However, the invariant may contain subsets that are not Z-linear. Consider {x + dN} ⊆ I, which is not yet invariant. We repeatedly apply non-inverting functions to {x + dN} to obtain new N-linear sets (not Z-linear sets). When the function applied 'moves' in the direction of the counters this will ultimately saturate (in particular when applying other counter functions). However, in the opposite direction, we may generate infinitely many such classes. Clearly we can examine all reachable residue classes defined by our reference counter. Any residue class reachable after an inverting function induces a Z-linear set. So it remains to consider those N-linear sets reachable without inverting functions. The remaining case to handle occurs when we repeatedly induce Nlinear sets until they repeat a residue class in the direction opposite to that of the reference counter.
We consider the case for h(x) = x + d with d ≥ 0. The case with h(x) = x − d is symmetric. It remains to detect when a set {x + dN} leads to {y + dN} by a sequence of non-inverting functions with x ≡ y mod d. Then by repeated application of these functions one can reach sets {z + dN} with z arbitrarily small, hence we can replace {x + dN} by {x + dZ}. We give further details in the appendix.
Reachability. The above procedure is sufficient to decide reachability. In all cases apart from that in which there are no counters, the invariants produced co-incide precisely with the reachability sets. A reachability query therefore reduces to asking whether the target belongs to the invariant.
In the remaining case, the invariant obtained is parametrised by the target via the bound C . The target lies within the region (−C , C +d), within which we can compute all reachable points. Thus once again, the target is reachable precisely if it belongs to the invariant. However, for a new target of larger modulus, a different invariant would need to be built.
Complexity. Lemma 11. Assume that all functions, starting point, and target point are given in unary. Then the invariant can be computed in polynomial time.
Without the unary assumption, the invariant could have exponential size, and hence require at least exponential time to compute. That is because the invariant we construct could include every value in an interval, for example, (−C, C), where C is of size polynomial in the largest value.
As shown in [10], the reachability problem is at least NP-hard in binary, because one can encode the integer Knapsack problem (which allows an object to be picked multiple times rather at most once). Moreover the Knapsack problem is efficiently solvable in pseudo-polynomial time via dynamic programming; that is, polynomial time assuming the input is in unary, matching the complexity of our procedure.

The POROUS Tool
Our invariant-synthesis tool porous 11 computes N-semi-linear invariants for point and Z-linear targets on systems defined by one-dimensional affine functions. porous includes implementations of the procedures of Theorem 3 (restricted to one-dimensional affine systems) and Theorem 6. porous is built in Python and can be used by command-line file input, a web interface, or by directly invoking the Python packages.
porous takes as input an instance (a start point, a target, and a collection of functions) and returns the generated invariant. Additionally it provides a proof that this set is indeed an inductive invariant: the invariant is a union of N-linear sets, so for each linear set and each function, porous illustrates the application of that function to the linear set and shows for which other linear set in the invariant this is a subset. Using this invariant, porous can decide reachability; if the specific target is reachable the invariant is not in itself a proof of reachability (since the invariant will often be an overapproximation of the global reachability set). Rather, equipped with the guarantee of reachability, porous searches for a direct proof of reachability: a sequence of functions from start to target (a process which would not otherwise be guaranteed to terminate).  Table 2. Results varying by size parameter (last row includes all instances tested). Times are given in seconds, with the average and maximum shown (except reachability proof time, which are all approximately 30s due to instances that terminate just before the timeout).
Experimentation. porous was tested on all 2 7 − 1 possible combinations of the following function types, with a ≥ 2, b ≥ 1: positive counters (x → x + b), negative counters (x → x − b), growing (x → ax ± b), inverting and growing (x → −ax ± b), inverters with positive counters (x → −x + b), inverters with negative counters (x → −x − b) and the pure inverter (x → −x). For each such combination a random instance was generated, with a size parameter to control the maximum modulus of a and b, ranging between 8 and 1024. The starting point was between 1 and the size parameter and the target was between 1 and 4 times the size parameter. Ten instances were tested for each size parameter and each of the 2 7 − 1 combinations, with between 1 and 9 functions of each type (with a bias for one of each function type).
Our analysis, summarised in Table 2, illustrates the effect of the size parameter. The time to produce the proof of invariant is separated from the process of building the invariant, since producing the proof of invariant can become slower as |I| becomes larger; it requires finding L k ∈ I such that f i (L j ) ⊆ L k for every linear set L j ∈ I and every affine function f i . In every case porous successfully built the invariant, and hence decided reachability very quickly (on average well below 1 second) and also produced the proof of invariance in around half a second on average. To demonstrate correctness in instances for which the target is reachable porous also attempts to produce a proof of reachability (a sequence of functions from start to target). Since our paper is focused on invariants as certificates of non-reachability, our proof-of-reachability procedure was implemented crudely as a simple breadth-first search without any heuristics, and hence a timeout of 30 seconds was used for this part of the experiment only.
Our experimental methodology was partially limited due to the high prevalence of reachable instances. A random instance will likely exhibit a large (often universal) reachability set. When two random counters are included, the chance that gcd(b 1 , b 2 ) = 1 (whence the whole space is covered) is around 60.8% and higher if more counters are chosen.
Overall around 86% of instances were reachable (of which 84% produced a proof within 30 seconds). Of the 14% of unreachable instances, all produced a proof, with the invariant taking around 0.2 seconds to build and 0.6 seconds to produce the proof. The 30-second timeout when demonstrating reachability directly is several orders of magnitudes longer than answering the reachability query via our invariant-building method.
A typical academic/consumer laptop was used to conduct the timing and analysis (a four-year-old, four-core MacBook Pro).

Conclusions and Open Directions
We introduced the notion of porous invariants, which are not necessarily convex and can in fact exhibit infinitely many 'holes', and studied these in the context of multipath (or branching/nondeterministic) affine loops over the integers, or equivalently nondeterministic integer linear dynamical systems. We have in particular focused on reachability questions. Clearly, the potential applicability of porous invariants to larger classes of systems (such as programs involving nested loops) or more complex specifications remains largely unexplored.
Our focus is on the boundary between decidability and undecidability, leaving precise complexity questions open. Indeed, the complexity of synthesising invariants could conceivably be quite high, except where we have highlighted polynomial-time results. On the other hand, the invariants produced should be easy to understand and manipulate, from both a human and machine perspective.
On a more technical level, in our setting the most general class of invariants that we consider are N-semi-linear. There remains at present a large gap between decidability for one-dimensional affine functions, and undecidability for linear updates in dimension 91 and above. It would be interesting to investigate whether decidability can be extended further, for example to dimensions 2 and 3.  [26] refers to linear independence, this can be converted to affine independence by increasing the dimension by one.
The procedure works via a pruned version breadth-first search, with nodes only expanded if its children are linearly independent of the current set. Hence, the first point found in the tree is the initial point x (0) , and therefore this point is included. The maximum depth of the tree that needs to be explored is d, and so every point included is reached with at most d applications of matrices to x (0) . Hence, if the largest absolute value of a point or matrix entry is µ, after d iterations, the largest absolute value is d d−1 µ d . This is by induction on the largest possible value µ for every entry: Base case: The result of [26] is in polynomial time in the number of arithmetic operations, we observe that this is also polynomial time in the bit-size. The independence checking in the algorithm involves checking linear independence of at most d vectors all having bit size at most log((dµ) d ) = d log(d) + d log(µ), which can be done in polynomial time in the bit-size (for example by Bareiss algorithm for calculating the determinant). Proof. We will prove the result for m + 5 matrices of size 7. This can then be transformed in a usual way to two matrices of size 7m + 35 (See Theorem 9 of [8] for instance).

B Proof of Lemma 3
In order to simplify the main part of the proof, let us first show that one can enforce an order between the matrices using affine transformations on one dimension. Let us denote p this dimension, it is initially equal to 1 and its target value is 0. Consider the three following affine transformation: f 1 (p) = 2p − 1, f 2 (p) = 2p − 2 and f 3 (p) = 2p, then the only sequences of transformation allowing to reach the target are of the form f * 3 f 2 f * 1 . Indeed, let I = {p | p ≥ 2 ∨ p ≤ −1}, we have (1) if p ∈ I, then for all i ∈ {1, 2, 3}, f i (p) ∈ I, (2) f 1 (1) = 1 and f 1 (0) ∈ I, (3) f 2 (1) = 0 and f 2 (0) ∈ I and (4) f 3 (1) ∈ I and f 3 (0) = 0. As a consequence, the inductive invariant I ensure that any sequence of transformation that do not have the desired order cannot reach the target. In the following, we will call type 1, 2 or 3 the transformations we define, depending on whether they implictly contain the function f 1 , f 2 or f 3 .
We reduce an instance {(u 1 , v 1 ), . . . , (u m , v m )} of the ω-PCP problem to the invariant synthesis problem. In order to simplify future notations, given a finite or infinite word w, we denote by |w| the length of the word w and given an integer i ≤ |w|, we write w i for the i'th letter of w. Given a finite or infinite word w on alphabet {1, . . . , m} we denote by u w and v w the words on the alphabet {0, 2} such that u w = u w1 u w2 . . . and v w = v w1 v w2 . . We work with 5 dimension, (s, c, d, n, k), and define the following transformations: -For i ≤ m, the type 1 transformation Simulate i on (s, c, d, n, k) encode the action of reading the pair (u i , v i ) and increases the counters n and k: it These m+5 transformations need 7 dimensions in total: the five above, (s, c, d, n, k), the one used for ordering the transformations,p, and one dimension constantly equal to 1, required to use affine transformations.
We now show that there is a solution to the given instance of the ω-PCP problem iff there does not exist a N-semi-linear invariant for the system with initial point x = (0, 1, 1, 0, 0, 1, 1), target y = (0, 0, 0, 1, 0, 0, 1) and using the matrices inducing the transformations defined above.
Assume first that there is a solution w to the ω-PCP instance. Consider the sequence of points (x n ) obtained as follows: for all j ∈ N, denoting w ≤j the prefix of w of length, x j = (s j , c j , 0, n j , k j , 0, 1) = Transfer Simulate w ≤j x where Simulate w ≤j represents the transformation Simulate w j . . . Simulate w2 Simulate w1 . We have that s j and c j are negative. Indeed, let (s, c, d) be the three first components of Simulate w ≤j x, we have that s = c[u wi ] − d[v wi ]. As w ≤j is a prefix of a solution to the ω-PCP instance, assuming |u w i| ≤ |v w i| this implies that Due to the above, by applying to the points x j a number of time the transformations Inc s and Inc c , we obtain the sequence of points (y j ) where y j = (0, 0, 0, n j , k j , 0, 1). We claim that any semi-linear invariant containing all the points y j also contains a point of the shape (0, 0, 0, 0, n j + d, k j , 0, 1) where d is a positive integer. This will imply the result as from such a point, one can reach the target by d − 1 applications of Dec k and k j applications of Dec and thus there is no semi-linear invariant of the system that does not intersect the target.
Let us now prove the above claim. Let I be a semi-linear set containing every point (y j ) (which we will see as two-dimensional objects by projecting on the 4th and 5th dimension). Then there exists a linear set I ⊆ I that contains infinitely many vectors of (y j ). This set I is defined by an initial vector, and a set of period vectors. As I contains infinitely many vectors of (y j ) where the ratios between the first and second component is increasing, one of the period vectors is of the form (d, 0) where d is a strictly positive integer. Let j be such that y j ∈ I , then (n j + d, k j ) ∈ I which implies the claim.
As a consequence, every N-semi-linear set over-approximating the system intersects with the target.
Conversely, assume that there is no solution to the ω-PCP instance. There exists n 0 ∈ N such that for every infinite word w on alphabet {0, . . . , m} there exists n ≤ n 0 such that u w n = v w n . Indeed, consider the tree which root is labelled by (ε, ε) and, given a node (u, v) of the tree, if for all n ≤ min(|u|, |v|) we have u n = v n , then this node has m children: the nodes (uu i , vv i ) for i = 1 . . . m. This tree is finitely branching and does not contain any infinite path (which would induce a solution to the ω-PCP instance). Thus, according to König's lemma, it is finite. We can therefore choose the height of this tree as our n 0 .
We define the invariant I = I 1 ∪ I 2 ∪ I 3 where . . , m} * ∧ |w| ≤ n 0 + 1 , ∧ w ∈ {1, . . . , m} * ∧ |w| ≤ n 0 + 1 ∧ s, t, n, k ∈ N and I 3 = (s, c, d, n, k, p, By definition, I is an N-semi-linear set, contains x and does not contain y. The difficulty is to show stability under the transformations. • Let z = Simulate w (x) ∈ I 1 , for some w ∈ {1, . . . , m} * with |w| ≤ n 0 + 1. By ordering if we apply a transformation outside Transfer or a Simulate i for some i, we reach I 3 .
As c ≥ d, this shows that Simulate wi z ∈ I 3 . -Transferz ∈ I 2 .
• Let z ∈ I 2 and f be one of the transformations, then f (z) ∈ I 2 if f increased (resp. decreased) a negative (resp. positive) component. Otherwise f (z) ∈ I 3 .
There is three possibilities (1) p = 2 and thus f (z) ∈ I 3 , (2) f = Transfer then p = 0 and either s ≥ 1 or c ≥ 1 and thus f (z) ∈ I 3 or (3) f = Simulate i for i ≤ m. In the latter case without loss of generality, assume that d c (this is completely symmetric in c and d ). We have that by assumption on |s| since max i c ≥ c , max i d/3 ≥ d (as m i ≥ 4) and max i 4. This shows that f (z) ∈ I 3 .
Therefore I is inductive and thus a N-semi-linear invariant of the system. This concludes the reduction. Lemma 12. For , k coprime, the sequence a n = (n mod k) for n ∈ N cycles through every modulo class {0, . . . , k − 1}.
Proof. Any path longer than k visits some class twice, and if the shortest cycle is k, then it visits every class.
Suppose there is a cycle of length less than k; then n = c+mk and (n+i) = c + m k and hence i = (m − m)k, with i < k. Since is an integer i divides (m − m)k then i = pr for p, r ∈ N such that m −m p is integer and k r is integer. Observe that since r ≤ i < k we have k r > 1. But this implies that k r divides k and , contradicting gcd(k, ) = 1. Proof. Let b = kd, c = d, where k, are co-prime.
We show there exists m, n ≥ 0 such that mb − cn = d. We have mb − cn = d ⇐⇒ mkd − n d = d ⇐⇒ mk − n = 1. Then choose m = 1+n k . By Lemma 12 there exists n such that n is in any modulo class modulo k, and thus too for 1 + n and so k divides 1 + n for some n.
Hence the set {x + dN} is included in the reachability set: we obtain x + jd by applying function f mj times and applying function g nj times. Similarly, we can find m , n ≥ 0 such that m b − cn = −d and thus {x + dZ} is within the reachability set.

D.2 Extended argument for non opposing counters
The following shows that if {x + dN} does lead to {y + dN}, with y < x and y ≡ x mod d, then indeed we can reach {z + dN} for any z ≡ x mod d by reapplying the same set of functions which lead from x to y. . Consider x (0) ∈ I and a path x (0) , f i1 , x (1) , f i2 , . . . , f im , x (m) such that x (j) = f ij (x (j−1) ), x (j) ≤ −B, x (m) < x (0) and x (0) ≡ x (m) mod d.
Remark 3. By symmetry, Lemma 13 also holds for the opposite direction. That is when h(x) = x − d, d > 0, inequalities are inverted and C is used in place of −B.
We now consider inductively applying non-inverting functions to sets {x + dN} ∈ I. Then add {f i (x) + dN} provided it is not already a subset of some set already in I. If {f i (x) + dN} is new and a new modulo class we can again apply Lemma 10, from whence we may also need to apply Lemma 8.
However, when this procedure does not saturate there eventually exists be a sequence of actions in which {x + dN} leads to {y + dN} with x ≡ y mod d according to a path in Lemma 13. In particular y < x < −B since if x < y then {y + dN} ⊆ {x + dN}, some modulo class must repeat after at most d steps, and eventually the procedure must stay < −B for at least d steps. Then, according to Lemma 13, a new Z-linear set can be added ({x + dZ}) (which again can be saturated using Lemma 8). We repeat this process until all N-linear sets are invariant. This process terminates, as each application of Lemma 13 adds a new Z-linear set with period d, of which there are at most d.

D.3 Proof of Lemma 11
Lemma 11. Assume that all functions, starting point, and target point are given in unary. Then the invariant can be computed in polynomial time.
Proof. In the no-counter case, by Lemma 9, there is an interval [−C, C] of size approximately |b|+|M | |a|−1 , where |b|, |M |, |a| are all numbers represented in the input, and thus is of polynomial of size. This means the gap is of polynomial size, and thus the saturation algorithm, which must in each step add a point or terminate, is of polynomial time.
In each counter-case there is a reference counter period d arising directly from the input as the counter part of some function, or in the case of two opposing counters, possibly the sum of two counter parts. For this period d there are at most 3d possible types of non-singleton invariant ({x + dN} or {x − dN} for some x and x + dZ for x ∈ {0, . . . , d} ). Singletons only arise in the interval [−C, C] if they exist. Hence, there are at most O(2C + 3d) steps which change the invariant.
In the case of two opposing counters, immediately all invariants are of the form x + dZ for x ∈ {0, . . . , d}, and the reachable modulo classes can be found in O(dk) (recall k is the number of functions), by breadth first search.
In the case of all counters in the same direction, there are two phases, each has a bounded number of steps. First we consider updates which move in the direction of the counters and secondly we consider updates which move against the counters.
In the case of moving with the counters, outside of [−C, C] all functions are growing. Hence, by conducting breadth first search on a priority queue that always expands the minimal element we can find the sets of the form x + dN for x ∈ {0, . . . , d} in polynomial time. Only inside [−C, C] does the search result in smaller elements (which there are at most 2C such steps), and in the remaining case we either expand to find an element already covered, or we find the smallest element in that modulo class. Thus this step takes O(dk + 2C) time.
Secondly we search for cycles in the direction opposing the counters, to see if we can turn any x + dN sets into x + dZ sets, that is invariants induced by Lemma 13. There can be a path of length at most d steps outside of [−C, C] before a cycle is found, so the running time is O(2Cd).

E Tool
The tool's output, when, applied to the MU Puzzle can be seen to produce the invariant {1 + 3Z} ∪ {2 + 3Z} of Example 1: The web-interface can be found at http://invariants.davidpurser.net.