Ali Asadi
email: aasadi@ce.sharif.edu

Krishnendu Chatterjee
email: krishnendu.chatterjee@ist.ac.at

Amir Hongfei Fu

Kafshdar Goharshady
email: goharshady@cse.ust.hk

Mohammad Mahdavi
email: mimmahdavi@ce.sharif.edu

Mohammad Mahdavi Polynomial

Amir Kafshdar Goharshady

Polynomial Reachability Witnesses via Stellensätze

Keywords: CCS Concepts:, Theory of computation → Logic and verification, Automated reasoning, • Software and its engineering → Formal software verification Reachability, Inductive Reasoning, Stellensätze

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Reachability. Reachability analysis is a basic and fundamental problem in computer science, starting from the halting problem of Turing machines. It is a core problem in program verification, constitutes the most basic liveness property, and has been widely studied in program analysis and model checking [START_REF] Clarke | Handbook of model checking[END_REF][START_REF] Robert | Assigning meanings to programs[END_REF][START_REF] Antony | An Axiomatic Basis for Computer Programming[END_REF][START_REF] Manna | Temporal verification of reactive systems: safety[END_REF][START_REF] Pnueli | The Temporal Logic of Programs[END_REF]. The target states considered in reachablity analysis can be either desirable so that reachability to these states should be guaranteed, or undesirable so that the goal is to find an execution path leading to an unwanted behavior, hence proving incorrectness of the system. Reachability to desirable states encodes the most basic type of liveness property, but reachability to undesirable states is also ubiquitous in verification problems and useful when one needs to identify realistic bugs in software implementations (see e.g. [START_REF] O'hearn | Incorrectness logic[END_REF]). Indeed, in real-world software development, most bugs are identified by finding an execution path that leads to a specific error [START_REF] Distefano | Scaling static analyses at Facebook[END_REF][START_REF] Godefroid | Compositional dynamic test generation[END_REF][START_REF] Majumdar | Hybrid Concolic Testing[END_REF]. This idea led to developments such as incorrectness logic [START_REF] O'hearn | Incorrectness logic[END_REF].

Previous Works on Formal Models. A large body of research on reachability analysis is conducted over formal models [START_REF] Clarke | Handbook of model checking[END_REF], such as finite-state systems [START_REF] Baier | Principles of model checking[END_REF], Petri nets [START_REF] Faouzi | Approximating Petri Net Reachability Along Context-free Traces[END_REF][START_REF] Czerwinski | The reachability problem for Petri nets is not elementary[END_REF][START_REF] Darondeau | Petri Net Reachability Graphs: Decidability Status of First Order Properties[END_REF][START_REF] Mayr | An Algorithm for the General Petri Net Reachability Problem[END_REF] and timed automata [START_REF] Alur | Automata For Modeling Real-Time Systems[END_REF]. For these models, precise decidability and complexity results are attained. Although these models serve as an important abstraction mechanism for realistic systems, the techniques for reachability analysis over them cannot be applied directly to imperative programs, because the values taken by variables in a program typically come from an infinite, even uncountable, domain and the underlying program structure might be irregular. In many cases a given program cannot be translated into any of the formal models above.

Reachability in Software Model Checkers. Many model checkers rely heavily on reachability analysis [START_REF] Ball | The SLAM project: debugging system software via static analysis[END_REF][START_REF] Beyer | The software model checker Blast[END_REF][START_REF] Beyer | CPAchecker: A Tool for Configurable Software Verification[END_REF][START_REF] Holzmann | The Model Checker SPIN[END_REF]. Notably, the BLAST project [START_REF] Beyer | The software model checker Blast[END_REF] describes itself as "a verification tool for the C language that solves the reachability problem". Even when considering safety properties, all approaches based on Counterexample-Guided Abstraction Refinement (CEGAR) [START_REF] Alur | Timing Verification by Successive Approximation[END_REF][START_REF] Balarin | An Iterative Approach to Language Containment[END_REF][START_REF] Clarke | Counterexample-Guided Abstraction Refinement[END_REF][START_REF] Gurfinkel | YASM: A Software Model-Checker for Verification and Refutation[END_REF][START_REF] Hajdu | Efficient strategies for CEGARbased model checking[END_REF], including SLAM [START_REF] Ball | A decade of software model checking with SLAM[END_REF][START_REF] Ball | The SLAM project: debugging system software via static analysis[END_REF] and BLAST [START_REF] Beyer | The software model checker Blast[END_REF], need to constantly perform reachability analyses to obtain their counterexamples. These model checkers rely on predicate abstraction refinement and, assuming a termination requirement, guarantee completeness when the variables have finite domains [START_REF] Henzinger | Lazy abstraction[END_REF].

Previous Works on (Imperative) Programs. When considering imperative programs, the reachability problem, and in particular the special case of termination analysis, has been widely studied over the past decades. Previous works include symbolic execution [START_REF] Cadar | Symbolic execution for software testing: three decades later[END_REF][START_REF] Cavada | The nuXmv Symbolic Model Checker[END_REF][START_REF] Jaffar | TRACER: A Symbolic Execution Tool for Verification[END_REF], termination analysis [START_REF] Robert | Assigning meanings to programs[END_REF], abstract interpretation [START_REF] Cousot | Abstract Interpretation: A Unified Lattice Model for Static Analysis of Programs by Construction or Approximation of Fixpoints[END_REF] and recent results on incorrectness logic [START_REF] De | Reverse Hoare Logic[END_REF][START_REF] O'hearn | Incorrectness logic[END_REF].

• Symbolic execution runs program code statically in a symbolic fashion, and is thus effective for programs without general unbounded loops. For programs with loops, symbolic execution can only unfold the loop up to a bounded depth, and hence cannot handle general loops with an unbounded number of iterations. This is also applicable to other approaches that rely on loop unrolling, such as [START_REF] Albarghouthi | From Under-Approximations to Over-Approximations and Back[END_REF]. • Termination analysis is a special kind of reachability which is usually guaranteed by well-foundedness reasoning such as (lexicographic) ranking functions [START_REF] Alias | Multi-dimensional Rankings, Program Termination, and Complexity Bounds of Flowchart Programs[END_REF][START_REF] Amir | Complexity of Bradley-Manna-Sipma Lexicographic Ranking Functions[END_REF][START_REF] Amir | On Multiphase-Linear Ranking Functions[END_REF][START_REF] Aaron R Bradley | Linear ranking with reachability[END_REF][START_REF] Bradley | Termination of Polynomial Programs[END_REF][START_REF] Michael | Synthesis of linear ranking functions[END_REF][START_REF] Cousot | Proving Program Invariance and Termination by Parametric Abstraction, Lagrangian Relaxation and Semidefinite Programming[END_REF][START_REF] Podelski | A Complete Method for the Synthesis of Linear Ranking Functions[END_REF][START_REF] Podelski | Transition Invariants[END_REF]. It does not consider target states defined through numerical constraints over program variables. • Abstract interpretation is mainly used to generate overapproximations of reachable states (i.e. certain states may be reached), but there are also several abstraction-based approaches that compute under-approximations [START_REF] Albarghouthi | From Under-Approximations to Over-Approximations and Back[END_REF][START_REF] Albarghouthi | UFO: A Framework for Abstraction-and Interpolation-Based Software Verification[END_REF][START_REF] Giacobazzi | Making abstract interpretations complete[END_REF][START_REF] Ranzato | Complete Abstractions Everywhere[END_REF][START_REF] Rival | Understanding the Origin of Alarms in Astrée[END_REF][START_REF] Schmidt | A calculus of logical relations for over-and underapproximating static analyses[END_REF]. However, they cannot provide guarantees of completeness except in specific special cases [START_REF] Giacobazzi | Completeness in Abstract Interpretation: A Domain Perspective[END_REF]. • Finally, incorrectness logic [START_REF] O'hearn | Incorrectness logic[END_REF] is a sound and complete logic that is similar to Hoare logic but performs underapproximation for reachable program states. A disadvantage of incorrectness logic, much like Hoare logic, is that it requires a considerable amount of manual effort for writing assertions, and cannot be directly automated.

Previous Works on Invariants. An invariant is, in a sense, a dual notion of reachability, and invariant generation is also prominent in the PL literature. Informally, an invariant is an over-approximation of the set of reachable states that can be used to prove safety properties over programs. Invariant generation is widely-studied in program analysis and verification. Many approaches are considered, e.g. abstract interpretation [START_REF] Bagnara | Precise Widening Operators for Convex Polyhedra[END_REF][START_REF] Singh | Fast polyhedra abstract domain[END_REF], constraint solving [START_REF] Chatterjee | Polynomial Invariant Generation for Non-deterministic Recursive Programs[END_REF][START_REF] Michael A Colón | Linear invariant generation using non-linear constraint solving[END_REF][START_REF] Sriram Sankaranarayanan | Non-linear loop invariant generation using Gröbner bases[END_REF][START_REF] Sankaranarayanan | Constraint-Based Linear-Relations Analysis[END_REF] and abductive inference [START_REF] Dillig | Inductive invariant generation via abductive inference[END_REF].

Our Focus. We consider reachability analysis over imperative programs and study the problem of automatically verifying that a set of target states can be reached in program execution. While invariants provide an over-approximation of the set of reachable states, we consider their natural dual, i.e. under-approximations of the set of states that can reach a target. We consider programs with non-determinism and distinguish between existential and universal reachability.

Existential reachability is the more classical notion and requires that target states are reachable under some resolution of the nondeterminism. Universal reachability requires the program to reach the target states no matter how the nondeterminism is resolved. We focus on existential reachability, but our results generalize to the universal case, too. Our Approach. Our methods are based on classical theorems (stellensätze) in polyhedral and real algebraic geometry. We first extend ideas from both ranking functions and inductive invariant generation to cover the reachability problem. Informally, we use techniques from inductive invariant generation to capture a subset T ♦ of program states from which the execution steps of the program will either reach our target states or stay in T ♦ itself. Simultaneously, we use arguments similar to ranking functions to ensure that every state in T ♦ can reach a target state in finitely many steps. While using an invariant and a ranking function to prove termination is a classical approach, the key distinction is that our set T ♦ is an under-approximation of the set of states that can eventually reach a target state, whereas invariants are overapproximations of reachable states. Moreover, our inductive sets are closed existentially (for each state in the inductive set, some successor should be in the same set), whereas inductive invariants are closed universally (all successors should be in the set). A similar approach based on danger invariants has been proposed in [START_REF] David | Danger invariants[END_REF]. A detailed comparison will be presented in the sequel.

Our Contributions. We propose a novel approach for reachability analysis over programs with real variables * . Our contributions are:

• We propose the notion of Inductive Reachability Witnesses (IRWs) for existential reachability, which consists of a set T ♦ of program states and a ranking function 𝑓 over T ♦ . The state set T ♦ satisfies certain invariant-like conditions. The ranking function 𝑓 serves as a proof that every state in T ♦ can indeed reach a target. We also propose the notion of Universal Inductive Reachability Witnesses (UIRWs), the counterpart of IRWs for the universal case. • From a theoretical point-of-view, we show that IRWs and UIRWs are sound and complete for proving existential and universal reachability, respectively. • Similar to several previous template-based works [START_REF] Chatterjee | Polynomial Invariant Generation for Non-deterministic Recursive Programs[END_REF][START_REF] Michael A Colón | Linear invariant generation using non-linear constraint solving[END_REF][START_REF] Sriram Sankaranarayanan | Non-linear loop invariant generation using Gröbner bases[END_REF][START_REF] Sankaranarayanan | Constraint-Based Linear-Relations Analysis[END_REF], we use Farkas' Lemma, Putinar's Positivstellensatz, and Handelman's Theorem for automatically synthesizing linear and polynomial IRWs/UIRWs. However, we face new challenges regarding satisfiability in the polynomial case and address them with novel methods based on Hilbert's Strong Nullstellensatz. This is the main technical novelty of this work. To the best of our knowledge, such methods (Section 3.2 and Theorem 8) are a novel contribution to constraint-based analysis of polynomial programs, and have no parallels even in the context of invariant generation or termination analysis. Moreover, our synthesis algorithms are complete in the linear case and semi-complete in the polynomial case. Comparison with Danger Invariants. Structures similar to IRWs, considering inductive sets that are used as underapproximations for reachability, have been studied in various contexts, including in [START_REF] Bakhirkin | Finding recurrent sets with backward analysis and trace partitioning[END_REF][START_REF] Tewodros A Beyene | Solving existentially quantified horn clauses[END_REF][START_REF] Chen | Proving nontermination via safety[END_REF][START_REF] Cousot | On fixpoint/iteration/variant induction principles for proving total correctness of programs with denotational semantics[END_REF][START_REF] Urban | Abstract interpretation of CTL properties[END_REF]. The most similar previous work is [START_REF] David | Danger invariants[END_REF] in which danger invariants are introduced. Like our notion of IRW, a (partial) danger invariant consists of an inductive set and a ranking function. However, we differ from [START_REF] David | Danger invariants[END_REF] in a number of significant ways:

• Problem: We can handle both universal and existential reachability, whereas [START_REF] David | Danger invariants[END_REF] is limited to the existential case. • Algorithm: Our synthesis algorithms are based on a novel combination of null-and positivstellensätze, while [START_REF] David | Danger invariants[END_REF]'s approach is based on symbolic execution, evolutionary mutations, and second-order SAT solving. Our stellensätzebased approach is the central contribution of this work. • Completeness of the Synthesis Algorithms: We consider linear and polynomial programs and witnesses. Our pushbutton synthesis algorithms are complete for the linear (Theorem 4) and semi-complete for the polynomial case (Theorem 10), while [START_REF] David | Danger invariants[END_REF] does not provide completeness. • Conceptual: [START_REF] David | Danger invariants[END_REF] attempts to synthesize an explicit scheduler in order to handle non-determinism. In contrast, we directly encode the non-determinism in the inputs to our stellensätze and leave the scheduler implicit. • Practical: Our experimental results in Section 4 demonstrate that our approach can handle instances that are beyond the reach of [START_REF] David | Danger invariants[END_REF]. There are also programs which are not expressible in our framework, due to being nonpolynomial, but can be handled by [START_REF] David | Danger invariants[END_REF]. Deep Bugs. A major challenge in reachability analysis is posed by deep bugs. If undesirable states are only reachable after a very long execution, then finite model checking or symbolic execution will have a hard time identifying the bug. Works that address this problem include [START_REF] David | Danger invariants[END_REF][START_REF] Ganty | Underapproximation of procedure summaries for integer programs[END_REF][START_REF] Kroening | Underapproximating loops in C programs for fast counterexample detection[END_REF][START_REF] Kumar Madhukar | Compositional Safety Refutation Techniques[END_REF]. Our runtime does not depend on the depth of the bug and we can identify various deep bugs with ease. See Section 4.

Comparison with Approaches for Deep Bugs. Although the main focus of our approach is providing completeness results in the linear and polynomial cases, which no previous method can provide, our approach is also very effective in proving the existence of deep bugs, i.e. proving reachability properties that can only be attained through long execution paths. As such, we now compare with several other approaches for deep bugs:

• The work of [START_REF] Kroening | Underapproximating loops in C programs for fast counterexample detection[END_REF] focuses mostly on deep bugs in programs that manipulate arrays. It considers under-approximating the behavior of loops using the notion of "auxiliary paths". This is similar to loop acceleration/summarization but can handle more general loops, given that it under-approximates the behavior of the loop rather than trying to match it exactly. This method is sound and applicable to a wide family of programs. However, it does not provide completeness guarantees. In contrast, our approach is sound and complete for linear and polynomial programs. • Another work in this category is [START_REF] Ganty | Underapproximation of procedure summaries for integer programs[END_REF]. This work focuses on summarizing and under-approximating the behavior of recursive functions and identifying bugs that can only be reached using deep recursion (a long stack). Another major difference between our work and [START_REF] Ganty | Underapproximation of procedure summaries for integer programs[END_REF] is that [START_REF] Ganty | Underapproximation of procedure summaries for integer programs[END_REF] focuses on integer programs, whereas our transition systems are defined with real variables. Finally, we also provide completeness guarantees. • Finally, [START_REF] Kumar Madhukar | Compositional Safety Refutation Techniques[END_REF] explores compositional approaches to proving reachability (safety refutation) and provides three algorithms with various levels of completeness. However, these algorithms are only applicable to loop-free programs, whereas we can handle (nested) loops.

Inductive Reachability Witnesses

We now provide the basic definitions needed for reachability analysis and introduce the concept of Inductive Reachability Witnesses (IRWs/UIRWs). In the sequel, we use transition systems with real variables to model programs.

Valuations. Let V be a finite set of variables. A valuation over V is a function ν : V → R. We denote the set of all valuations over V by R V . Transition Systems. A transition system (or simply system) is a tuple 𝑆 = (V, L, ℓ 0 , 𝐼, Θ), in which V is a finite set of variables, L is a finite set of locations, ℓ 0 ∈ L is the initial or starting location, 𝐼 is an assertion over V which defines the set of possible initial valuations, and Θ is a finite set of transitions.

Each transition 𝜃 ∈ Θ is of the form 𝜃 = (ℓ, ℓ ′ , 𝜑, 𝜇) where ℓ, ℓ ′ ∈ L are the pre and post locations, 𝜑 is an assertion over V that serves as the transition condition, and 𝜇 : R V → R V is an update function. For brevity, in the sequel, we assume that we have fixed a system 𝑆 = (V, L, ℓ 0 , 𝐼, Θ) which is under study. For a location ℓ ∈ L, we write Θ ℓ to denote the set of transitions out of ℓ. We say that a system is 𝛽-branching if |Θ ℓ | ≤ 𝛽 for every location ℓ.

Example 1. Consider the program in Figure 1 (top), in which □ denotes non-deterministic choice between transitioning to 𝑏 or 𝑐. The transition system in Figure 1 (bottom) represents this program. Note that we have ℓ 0 = 𝑎 and assume the initial valuations satisfy 𝑥, 𝑦, 𝑧 ≥ 0.

States. A state in 𝑆 is a pair 𝜎 = (ℓ, ν) ∈ L × R V , consisting
of a location and a valuation. We denote the set of states by Σ. A subset Σ ′ ⊆ Σ is called bounded if the set of valuations that appear in the elements of Σ ′ is bounded.

Successors. A state 𝜎 ′ = (ℓ ′ , ν ′ ) is called a successor of a state 𝜎 = (ℓ, ν) if there exists a transition 𝜃 = (ℓ, ℓ ′ , 𝜑, 𝜇) ∈ Θ such that ν |= 𝜑 and ν ′ = 𝜇 (ν).
For theoretical elegance, we assume that every state has at least one successor, i.e. we assume that there is a transition from the final state to itself. Runs. A run of the system 𝑆 is an

𝐼 : 𝑥 ≥ 0 ∧ 𝑦 ≥ 0 ∧ 𝑧 ≥ 0 𝑎 : while 𝑥 ≥ 𝑦 : 𝑏 : □ (𝑥, 𝑦) := (𝑥 + 1, 𝑦 + 2) 𝑐 : □ (𝑥, 𝑦, 𝑧) := (𝑥 + 𝑧, 𝑦 + 𝑧, 𝑧 -1) 𝑑 : a d c b ( x , y , z ) ↦ ( x + 1 , y + 2 , z ) ( x , y , z ) ↦ ( x + z , y + z , z -1 ) x < y x ≥ y x ≥ y θ 1 θ 2 θ 4 θ 5 θ 3 θ 6 Figure 1. A Simple
infinite sequence r = {𝜎 𝑖 , 𝜃 𝑖 } ∞ 𝑖=0 = {(ℓ 𝑖 , ν 𝑖 ), 𝜃 𝑖 } ∞ 𝑖=0
, where each 𝜎 𝑖 ∈ Σ is a state consisting of ℓ 𝑖 ∈ L and ν 𝑖 ∈ R V , and each 𝜃 𝑖 = (ℓ 𝑖 , ℓ 𝑖+1 , 𝜑 𝑖 , 𝜇 𝑖 ) ∈ Θ is a transition from ℓ 𝑖 to ℓ 𝑖+1 , such that: • r starts in the initial location ℓ 0 ; • ν 0 |= 𝐼 , i.e. the initial valuation satisfies the assertion 𝐼 ; • For every 𝑖, we have ν 𝑖 |= 𝜑 𝑖 and ν 𝑖+1 = 𝜇 𝑖 (ν 𝑖 ), i.e. 𝜎 𝑖+1 is a successor of 𝜎 𝑖 through 𝜃 𝑖 . Semi-runs. A semi-run is defined similarly to a run, except that it does not have to start at ℓ 0 or satisfy 𝐼 . A path of length 𝑛 is a finite prefix 𝜋 = 𝜎 0 , 𝜃 0 , . . . , 𝜎 𝑛-1 , 𝜃 𝑛-1 , 𝜎 𝑛 of a run. A path must always end at a state. Similarly, a semi-path is a finite prefix of a semi-run that ends at a state. Non-determinism. The system 𝑆 is called deterministic if there is exactly one possible transition at every state. Formally, 𝑆 is deterministic if for every 𝜎 = (ℓ, ν) ∈ Σ, there exists exactly one 𝜃 ∈ Θ such that 𝜃 = (ℓ, ℓ ′ , 𝜑, 𝜇) and ν |= 𝜑. Otherwise, 𝑆 is non-deterministic. Existential Reachability. A set T ⊆ Σ is called existentially reachable or simply reachable if there exists an integer 𝑛 and a run r = {𝜎 𝑖 , 𝜃 𝑖 } ∞ 𝑖=0 such that 𝜎 𝑛 ∈ T. In other words, T is reachable if there exists a run that visits T. Universal Reachability. A set T ⊆ Σ is called universally reachable if there exists a valuation ν 0 ∈ R V and an integer 𝑛, such that (i) ν 0 |= 𝐼 , and (ii) every run r = {(ℓ 𝑖 , ν 𝑖 ), 𝜃 𝑖 } ∞ 𝑖=0 visits T in its first 𝑛 steps. In other words, for each such r, there exists an index 𝑖 ≤ 𝑛 such that (ℓ 𝑖 , ν 𝑖 ) ∈ T.

Intuitively, this requires that we can fix an initial valuation such that no matter how the non-determinism is resolved, the execution is forced to visit T after at most 𝑛 steps. Note that our notion of universal reachability is bounded by 𝑛 and different from having ♢T in every execution path.

Motivation for Universal Reachability. Universal reachability is useful for modeling angelic non-determinism. Consider that T is a set of undesirable (bug) states and the question is whether for all initial valuations, there exists an execution that avoids the bug. For example, consider a cyber-physical system modeled as a transition system and a controller who can resolve non-determinism. Deciding (the negation of) universal reachability corresponds to whether for all choices of initial state, a controller can avoid the bug. Example 3. Consider the system in Figure 1 (bottom), and let

T = {(𝑑, ν) | ν ∈ R V }.
Reaching T is equivalent to the termination of the program in Figure 1 (top). T is existentially reachable, i.e. there are runs of the system that reach label 𝑑, for example the following:

(𝑎, 0, 0, 0) 𝜃 1 -→ (𝑏, 0, 0, 0) 𝜃 4 -→ (𝑎, 1, 2, 0) 𝜃 3 -→ (𝑑, 1, 2, 0) → . . . .
It is also universally reachable, because every execution starting from (𝑎, 1, 2, 3) will reach 𝑑 in a single step. As another example, consider the target set T ′ = {(𝑑, ν) | ν(𝑥) < 0}. This corresponds to reaching 𝑑 (ending the program) with a negative value for 𝑥. This time, the set T ′ is existentially reachable, for example through the following run:

(𝑎, 0, 0, 0)

𝜃 2 -→ (𝑐, 0, 0, 0) 𝜃 5 -→ (𝑎, 0, 0, -1) 𝜃 2
-→ (𝑐, 0, 0, -1)

𝜃 5 -→ (𝑎, -1, -1, -2) 𝜃 2 -→ (𝑐, -1, -1, -2) 𝜃 5 -→ (𝑎, -3, -3, -3) 𝜃 1 -→ (𝑏, -3, -3, -3) 𝜃 4 -→ (𝑎, -2, -1, -3) 𝜃 3
-→ (𝑑, -2, -1, -3) → . . . , but it is not universally reachable. To see this, note that if an initial value satisfies 𝑥 < 𝑦, then it does not enter the while loop at all, and hence when it reaches 𝑑 it satisfies 𝑥 ≥ 0 (the initial condition). On the other hand, if an initial value satisfies 𝑥 ≥ 𝑦, there is a run that always chooses the transition 𝜃 2 when at 𝑎, and hence never reaches T ′ .

We now look into proof concepts for reachability. T-inductive Sets. Given a set T ⊆ Σ of target states, a set T ♦ ⊆ Σ is called T-inductive if for every 𝜎 ∈ T ♦ \ T, there exists a successor 𝜎 ′ of 𝜎 such that 𝜎 ′ ∈ T ♦ .

Intuitively, if T ♦ is T-inductive, then if we start the execution of the program from a state in T ♦ , there exists a way for resolving the non-determinism so that we either reach T or can inductively prove that we will never leave T ♦ . Example 4. Consider the system in Figure 1 and let T = {(𝑑, ν) | ν(𝑥) < 0}, i.e. the target is reaching 𝑑 with 𝑥 having a negative value. Let

T ♦ := {(ℓ, ν) | ℓ ∈ L, ν ∈ R V , ν |= 𝐴 ℓ }
be the set of states satisfying the following assertions:

ℓ 𝐴 ℓ 𝑎 𝑥, 𝑦, 𝑧 ≤ 0 ∧ (𝑥 -𝑦) • (𝑥 -𝑦 + 1) = 0 𝑏 𝑥 ≤ -2 ∧ 𝑦, 𝑧 ≤ 0 ∧ 𝑥 = 𝑦 𝑐 𝑥, 𝑦, 𝑧 ≤ 0 ∧ 𝑥 = 𝑦 𝑑 𝑥 < 0
Then, we can verify that T ♦ is a T-inductive set. Concretely, consider a state (𝑎, ν 𝑎 ) ∈ T ♦ . In other words, ν 𝑎 |= 𝐴 𝑎 . In such a state, we have (𝑥 -𝑦) • (𝑥 -𝑦 + 1) = 0. Therefore, either ν 𝑎 (𝑥) = ν 𝑎 (𝑦) or ν 𝑎 (𝑥) = ν 𝑎 (𝑦) -1. In the former case, we can take transition 𝜃 2 , and it is easy to verify that the new state satisfies 𝐴 𝑐 , hence there is a successor that is also in T ♦ . In the latter case, we can take 𝜃 3 and reach 𝑑 with a valuation that satisfies 𝑥 < 0, because ν 𝑎 |= (𝑦 ≤ 0 ∧ 𝑥 = 𝑦 -1). Similarly, if (𝑏, ν 𝑏 ) ∈ T ♦ , we know that ν 𝑏 |= (𝑥 ≤ -2 ∧𝑦, 𝑧 ≤ 0 ∧ 𝑥 = 𝑦). Therefore, taking the transition 𝜃 4 , corresponding to the update (𝑥, 𝑦) := (𝑥+1, 𝑦+2), leads to a state in 𝑎 that satisfies (𝑥, 𝑦, 𝑧 ≤ 0 ∧ 𝑥 = 𝑦 -1). Note that 𝑥 = 𝑦 -1 ⇒ (𝑥 -𝑦) • (𝑥 -𝑦 + 1) = 0, therefore 𝐴 𝑎 is satisfied and we have a successor in T ♦ . It is easy to verify the same property at 𝑐. Finally, if we have a state (𝑑, ν 𝑑 ) ∈ T ♦ , then it is also in T.

If we start at an initial state that satisfies 𝐴 𝑎 , we can find a run of the system that either reaches T or stays inside T ♦ . However, this is not enough for reachability to T. Such a run might stay inside T ♦ forever without visiting T. For example, we can keep taking the transition 𝜃 2 when at 𝑎, and hence never reach 𝑑. To avoid this, we need a T-ranking function.

T-ranking Functions. Given a T-inductive set T ♦ , a function 𝑓 : T ♦ → [0, ∞) is called a T-ranking function with parameter 𝜖 > 0, if for every 𝜎 ∈ T ♦ \ T, there exists a successor 𝜎 ′ ∈ T ♦ of 𝜎, for which we have 𝑓 (𝜎 ′ ) ≤ 𝑓 (𝜎) -𝜖. Inductive Reachability Witnesses (IRWs). Given a set T of target states in a system 𝑆 = (V, L, ℓ 0 , 𝐼, Θ), an Inductive Reachability Witness for T is a tuple (T ♦ , 𝑓 , 𝜖) such that:

• T ♦ is a T-inductive set; • 𝜖 ∈ (0, ∞); • 𝑓 : T ♦ → [0, ∞) is T-

ranking with parameter 𝜖;

• There exists ν ∈ R V such that (ℓ 0 , ν) ∈ T ♦ and ν |= 𝐼 .

Informally, an IRW serves as a proof of existential reachability for a target set T. The inductivity of T ♦ ensures that starting from the initial state (ℓ 0 , ν) ∈ T ♦ , we will never be forced to leave T ♦ unless we reach T, while the existence of the T-ranking function 𝑓 proves that we cannot avoid T forever. It is also noteworthy that the T-inductive set T ♦ is similar to an inductive invariant, but the main difference is that while an invariant is by definition a superset of all reachable states, a T-inductive set T ♦ is a subset of those states from which we can reach the target set T. An IRW (T ♦ , 𝑓 , 𝜖) is called bounded if T ♦ is bounded.

Example 5. Consider the system in Figure 1, with the same target set as in Example 4, i.e.

T = {(𝑑, ν) | ν(𝑥) < 0}. Let T ♦ := {(ℓ, ν) | ν |= 𝐴 ℓ } and 𝑓 (ℓ, ν) := 𝑓 ℓ (ν) be defined as: ℓ 𝐴 ℓ 𝑓 ℓ 𝑎 -10 ≤ 𝑥, 𝑦, 𝑧 ≤ 0 ∧ 𝑥 = 𝑦 -1 ∨ 𝑥 = 𝑦 = -𝑧•(𝑧+1) 2 100 + 𝑥 -𝑦 + 𝑧 𝑏 -10 ≤ 𝑥 ≤ -2 ∧ 𝑧 ≤ 0 ∧ 𝑥 = 𝑦 = -𝑧•(𝑧+1) 2 99.5 + 𝑧 𝑐 -2 ≤ 𝑥 ≤ 0 ∧ 𝑧 ≤ 0 ∧ 𝑥 = 𝑦 = -𝑧•(𝑧+1) 2 99.5 + 𝑧 𝑑 𝑥 ≤ -0.5 0
Note that the 𝐴 ℓ 's are more restrictive than in Example 4. We can verify that T ♦ is a T-inductive set in the same manner as in Example 4. We should also verify that 𝑓 is a valid T-ranking function. Whenever we take either transition 𝜃 1 or 𝜃 2 (from 𝑎 to 𝑏 or 𝑐), we are assured that 𝑥 = 𝑦, hence the value of 𝑓 goes from 100 + 𝑧 to 99.5 + 𝑧 and decreases by 0.5. Also, because in 𝐴 𝑎 we have -10 ≤ 𝑥, 𝑦, 𝑧 ≤ 0, the value of 𝑓 at 𝑎 is at least 80, and hence transition 𝜃 3 (from 𝑎 to 𝑑) decreases 𝑓 by more than 0.5. Now consider transition 𝜃 4 (from 𝑏 to 𝑎). This transition does not change the value of 𝑧, but makes it so that 𝑦 = 𝑥 + 1. So it changes the value of 𝑓 from 99.5 + 𝑧 to 99 + 𝑧. Note that transition 𝜃 5 (from 𝑐 to 𝑎), decreases 𝑧 by 1 while keeping 𝑥 = 𝑦. Hence, it decreases 𝑓 by 0.5. Also, 𝜃 6 (the self-transition from 𝑑 to 𝑑) is irrelevant in this case, because our 𝐴 𝑑 entails inclusion in T. Finally, (𝑎, 0, 0, 0) satisfies both 𝐼 and 𝐴 𝑎 .

UIRWs. A Universal Inductive Reachability Witness (UIRW) is defined similarly, except that existential quantifiers in T ♦ and 𝑓 are replaced by universal quantifiers. See Appendix A. Remark 1. In existential IRWs (such as Example 5), the set T ♦ is an under-approximation of the states from which there exists a way of resolving the non-determinism so that we eventually reach T. Similarly, in UIRWs (Appendix A), T ♦ underapproximates the set of states from which every execution of the program is forced to visit T. Hence, our T-inductive sets T ♦ are essentially natural duals of inductive invariants [START_REF] Chatterjee | Polynomial Invariant Generation for Non-deterministic Recursive Programs[END_REF][START_REF] Michael A Colón | Linear invariant generation using non-linear constraint solving[END_REF].

We prove existential (universal) reachability by synthesizing an IRW (a UIRW). This is both sound and complete. Theorem 1 (Soundness, Proof in Appendix B). Let T ⊆ Σ be a set of states in the system 𝑆.

(i) If there exists an IRW (T ♦ , 𝑓 , 𝜖) for T, then T is existentially reachable. (ii) If there exists a UIRW (T ♦ , 𝑓 , 𝜖) for T, then T is universally reachable. Theorem 2 (Completeness, Proof in Appendix B). Let T ⊆ Σ be a set of states in the system 𝑆.

(i) If T is existentially reachable, then there exists an IRW (T ♦ , 𝑓 , 𝜖) for T. (ii) If T is universally reachable, then there exists a UIRW (T ♦ , 𝑓 , 𝜖) for T.

Undecidability. Synthesis of IRWs (UIRWs) is equivalent to proving existential (universal) reachability, which are undecidable due to Rice's theorem [START_REF] Gordon | Classes of recursively enumerable sets and their decision problems[END_REF]. We consider the special case of linear or polynomial systems, with target sets that are defined by linear/polynomial inequalities, and focus on synthesizing linear/polynomial IRWs and UIRWs.

Polynomial Systems. A transition system 𝑆 = (V, L, ℓ 0 , 𝐼, Θ) is called (𝑑, 𝑘)-polynomial if • 𝐼 is a conjunction of at most 𝑘 polynomial inequalities of degree at most 𝑑 over V, and • for every 𝜃 = (ℓ, ℓ ′ , 𝜑, 𝜇) ∈ Θ, the transition condition 𝜑 is a conjunction of at most 𝑘 polynomial inequalities of degree at most 𝑑 over V, and

• for every 𝜃 = (ℓ, ℓ ′ , 𝜑, 𝜇) ∈ Θ and variable 𝑣 ∈ V, we have 𝜇 (𝑣) ∈ R[V] and deg(𝜇 (𝑣)) ≤ 𝑑, i.e. 𝜇 (𝑣) is a polynomial of degree at most 𝑑 over V. A (1, 𝑘)-polynomial system is also called 𝑘-linear. Linear IRWs/UIRWs. An IRW/UIRW (T ♦ , 𝑓 , 𝜖) is called 𝑘-linear if for every location ℓ ∈ L :

• The set T ♦ ℓ := T ♦ ∩ ({ℓ } ×R V ) is a closed polyhedron which is an intersection of at most 𝑘 half-spaces, i.e. there exists a set 𝐴 ℓ of at most 𝑘 non-strict linear inequalities over V such that a valuation ν satisfies 𝐴 ℓ iff (ℓ, ν) ∈ T ♦ .

• The function 𝑓 ℓ : Sat(𝐴 ℓ ) → [0, ∞), defined as 𝑓 ℓ (ν) = 𝑓 (ℓ, ν), is a linear function over V. Here, Sat(𝐴 ℓ ) is the set of all valuations that satisfy 𝐴 ℓ . Polynomial IRWs/UIRWs. An IRW/UIRW (T ♦ , 𝑓 , 𝜖) is called (𝑑, 𝑘)-polynomial if for every ℓ ∈ L :

• The set T ♦ ℓ := T ♦ ∩ ({ℓ } × R V
) is a closed semi-algebraic set defined by at most 𝑘 non-strict polynomial inequalities of degree 𝑑 or less. Equivalently, there exists a set 𝐴 ℓ of at most 𝑘 non-strict polynomial inequalities of degree at most 𝑑 over V such that ν |= 𝐴 ℓ iff (ℓ, ν) ∈ T ♦ .

• The function 𝑓 ℓ , defined as 𝑓 ℓ (ν) = 𝑓 (ℓ, ν), is a polynomial of degree at most 𝑑 over V. A (𝑑, 𝑘)-polynomial IRW/UIRW is explicitly bounded if each set 𝐴 ℓ contains a polynomial inequality 𝑔 ≥ 0 such that Sat(𝑔 ≥ 0) is bounded.

Automated Synthesis of IRWs/UIRWs

We now provide our algorithms for synthesizing linear or polynomial witnesses for linear or polynomial systems.

Synthesizing Linear IRWs/UIRWs

Problem Definition. Given a 𝑘-linear system 𝑆 = (V, L, ℓ 0 , 𝐼, Θ), together with a set 𝜏 ℓ of at most 𝑘 non-strict linear inequalities at every location ℓ ∈ L, synthesize a 𝑘-linear IRW/UIRW for the target set T := ∪ ℓ ∈L {ℓ } × Sat(𝜏 ℓ ) or report that no such IRW/UIRW exists. In the sequel, we assume V = {𝑣 1 , . . . , 𝑣 𝑟 }, and L = {ℓ 0 , . . . , ℓ 𝑛 }. Mathematical Tool. Our approach in this section is based on a well-known theorem in linear programming, called Farkas' Lemma. The presentation we use is similar to [START_REF] Chatterjee | Computational approaches for stochastic shortest path on succinct MDPs[END_REF][START_REF] Michael A Colón | Linear invariant generation using non-linear constraint solving[END_REF]. Lemma 1 (Farkas' Lemma [START_REF] Farkas | Theorie der einfachen Ungleichungen[END_REF][START_REF] Matousek | Understanding and using linear programming[END_REF]). Consider a set V = {𝑣 1 , . . . , 𝑣 𝑟 } of real-valued variables and the following system of 𝑚 linear inequalities over V:

Φ :          𝑎 1,0 + 𝑎 1,1 • 𝑣 1 + . . . + 𝑎 1,𝑟 • 𝑣 𝑟 ≥ 0 . . . 𝑎 𝑚,0 + 𝑎 𝑚,1 • 𝑣 1 + . . . + 𝑎 𝑚,𝑟 • 𝑣 𝑟 ≥ 0
When Φ is satisfiable, it entails a given linear inequality

𝜓 : 𝑐 0 + 𝑐 1 𝑣 1 + . . . + 𝑐 𝑟 𝑣 𝑟 ≥ 0
iff 𝜓 can be written as a non-negative linear combination of 1 ≥ 0 and the inequalities in Φ, i.e. iff there exist non-negative real numbers 𝑦 0 , 𝑦 1 , . . . , 𝑦 𝑚 , such that:

𝑐 0 = 𝑦 0 + 𝑚 𝑖=1 𝑦 𝑖 • 𝑎 𝑖,0 ; 𝑐 1 = 𝑚 𝑖=1 𝑦 𝑖 •𝑎 𝑖,1 ; . . . ; 𝑐 𝑟 = 𝑚 𝑖=1 𝑦 𝑖 •𝑎 𝑖,𝑟
. Moreover, Φ is unsatisfiable iff -1 ≥ 0 can be derived as above.

We often encounter Φ's with both strict and non-strict linear inequalities. So, we use the following corollary: Corollary 1 (Proof in Appendix D). Take a set V = {𝑣 1 , . . . , 𝑣 𝑟 } of real-valued variables and the following system of 𝑚 linear inequalities over V:

Φ :          𝑎 1,0 + 𝑎 1,1 • 𝑣 1 + . . . + 𝑎 1,𝑟 • 𝑣 𝑟 1 0 . . . 𝑎 𝑚,0 + 𝑎 𝑚,1 • 𝑣 1 + . . . + 𝑎 𝑚,𝑟 • 𝑣 𝑟 𝑚 0
in which 𝑖 ∈ {>, ≥}. When Φ is satisfiable, it entails a given non-strict linear inequality 𝜓 : 𝑐 0 + 𝑐 1 𝑣 1 + . . . + 𝑐 𝑟 𝑣 𝑟 ≥ 0 iff 𝜓 can be written as a non-negative linear combination of 1 ≥ 0 and the inequalities in Φ, i.e. iff there exist non-negative real numbers 𝑦 0 , 𝑦 1 , . . . , 𝑦 𝑚 , such that:

𝑐 0 = 𝑦 0 + 𝑚 𝑖=1 𝑦 𝑖 • 𝑎 𝑖,0 ; 𝑐 1 = 𝑚 𝑖=1 𝑦 𝑖 •𝑎 𝑖,1 ; . . . ; 𝑐 𝑟 = 𝑚 𝑖=1 𝑦 𝑖 •𝑎 𝑖,𝑟 .
Moreover, Φ is unsatisfiable iff either -1 ≥ 0 can be derived as above, or 0 > 0 can be derived as above while requiring 𝑖 ∈ {> } 𝑦 𝑖 > 0 (i.e. in order to derive a strict inequality, we should use at least one of the strict inequalities in Φ with non-zero coefficient).

Overview of the Algorithm. Before presenting our algorithm in detail, we provide a high-level overview of its steps:

• Step 1.
The algorithm creates a template for the desired IRW/UIRW. It considers every expression that should be synthesized as part of an IRW/UIRW, i.e. the descriptions of T ♦ and 𝑓 , and creates a template for it in which the coefficients are unknown (to be synthesized). • Step 2. The algorithm generates a series of so-called "constraint pairs". These constraint pairs are of a specific form that is amenable to Farkas' Lemma. They encode the requirements that T ♦ should be a T-inductive set and that 𝑓 should be a valid T-ranking function. • Step 3. The algorithm applies Farkas' lemma to the constraints generated in Step 2 and translates them to an equivalent system of quadratic (in)equalities over the unknown template variables. After this step, no program variable appears in the quadratic constraints. • Step 4. The algorithm adds a few additional constraints that ensure the existence of a suitable initial valuation. • Step 5. Finally, it solves the constraints by calling an off-the-shelf Quadratic Programming (QP) solver. It then plugs back the solution into the templates generated in

Step 1 to obtain the desired IRW/UIRW. We now provide the details of each step of our algorithm.

Step 1. Setting up a template. Take a 𝑘-linear IRW/UIRW for reaching T in 𝑆. It has a 𝑘-linear set T ♦ , defined by a set 𝐴 ℓ of 𝑘 linear inequalities at every location ℓ, and a function 𝑓 , also defined by a linear expression 𝑓 ℓ at every location ℓ. The algorithm sets up a symbolic template for each 𝐴 ℓ 𝐼 : 𝑥, 𝑦 ≥ 10 1 : i f 𝑥 < 𝑦 : 2 : □ 𝑥 := 𝑥 + 10 3 :

□ 𝑥 := 𝑥 + 5 4 :

1 4 3 2 ( x , y ) ↦ ( x + 1 0 , y ) ( x , y ) ↦ ( x + 5 , y ) x ≥ y x < y x < y θ 1 θ 2 θ 4 θ 5 θ 3 θ 6
𝐴 ℓ :          𝑐 ℓ,1,0 + 𝑐 ℓ,1,1 • 𝑣 1 + . . . + 𝑐 ℓ,1,𝑟 • 𝑣 𝑟 ≥ 0 . . . 𝑐 ℓ,𝑘,0 + 𝑐 ℓ,𝑘,1 • 𝑣 1 + . . . + 𝑐 ℓ,𝑘,𝑟 • 𝑣 𝑟 ≥ 0 𝑓 ℓ = 𝑑 ℓ,0 + 𝑑 ℓ,1 • 𝑣 1 + . . . + 𝑑 ℓ,𝑟 • 𝑣 𝑟
Intuitively, the goal of the algorithm is to find suitable real values for the unknown coefficients (i.e. 𝑐 ℓ,𝑖,𝑗 's and 𝑑 ℓ,𝑗 's) so that when we plug them into 𝑓 ℓ 's and 𝐴 ℓ 's, they yield a valid IRW/UIRW. The algorithm also defines a new unknown 𝜖 which serves as the decrease parameter for 𝑓 . Example 6. Consider the system in Figure 2. We will use this system as our running example and aim to synthesize a 2-linear IRW and a 2-linear UIRW for it. For the IRW case, suppose that the target set is

T = {(4, ν) | ν |= (𝑥 ≥ 𝑦 + 8)}.
For the UIRW case, we let

T ′ = {(4, ν) | ν |= (𝑥 ≥ 𝑦 + 4)}.
The algorithm generates a variable 𝜖 and the following templates:

𝐴 1 : 𝑐 0 + 𝑐 1 • 𝑥 + 𝑐 2 • 𝑦 ≥ 0 𝑐 3 + 𝑐 4 • 𝑥 + 𝑐 5 • 𝑦 ≥ 0 𝐴 2 : 𝑐 6 + 𝑐 7 • 𝑥 + 𝑐 8 • 𝑦 ≥ 0 𝑐 9 + 𝑐 10 • 𝑥 + 𝑐 11 • 𝑦 ≥ 0 𝐴 3 : 𝑐 12 + 𝑐 13 • 𝑥 + 𝑐 14 • 𝑦 ≥ 0 𝑐 15 + 𝑐 16 • 𝑥 + 𝑐 17 • 𝑦 ≥ 0 𝐴 4 : 𝑐 18 + 𝑐 19 • 𝑥 + 𝑐 20 • 𝑦 ≥ 0 𝑐 21 + 𝑐 22 • 𝑥 + 𝑐 23 • 𝑦 ≥ 0 𝑓 1 = 𝑑 0 + 𝑑 1 • 𝑥 + 𝑑 2 • 𝑦 𝑓 2 = 𝑑 3 + 𝑑 4 • 𝑥 + 𝑑 5 • 𝑦 𝑓 3 = 𝑑 6 + 𝑑 7 • 𝑥 + 𝑑 8 • 𝑦 𝑓 4 = 𝑑 9 + 𝑑 10 • 𝑥 + 𝑑 11 • 𝑦
The goal is to synthesize values for 𝜖, 𝑐 0 , . . . , 𝑐 23 and 𝑑 0 , . . . , 𝑑 11 , so that the templates above become an IRW/UIRW.

Step 2a. Computing IRW Constraint Pairs. This step is only performed when we want to synthesize an IRW. In an IRW, the existential inductive set T ♦ should satisfy the condition that for every state 𝜎 ∈ T ♦ \ T, there exists a successor 𝜎 ′ ∈ T ♦ of 𝜎 with 𝑓 (𝜎 ′ ) ≤ 𝑓 (𝜎) -𝜖. Let ℓ ∈ L be a location and Θ ℓ be the set of transitions out of ℓ, † We use the notation • to denote variables/expressions whose values should be synthesized by the algorithm.

i.e. transitions whose pre-location is ℓ. The IRW properties at ℓ are equivalent to:

∀ν ∈ R V , ν |= 𝐴 ℓ ⇒ ν |= 𝜏 ℓ ∨ 𝜃 =(ℓ,ℓ ′ ,𝜑,𝜇) ∈Θ ℓ 𝜉 (𝜃 ) (1)
where 𝜉 (𝜃 ) = 𝜉 (ℓ, ℓ ′ , 𝜑, 𝜇) is defined as:

𝜉 (𝜃 ) := ν |= 𝜑 ∧ 𝜇 (ν) |= 𝐴 ℓ ′ ∧ 𝑓 ℓ ′ (𝜇 (ν)) ≤ 𝑓 ℓ (ν) -𝜖 (2)
Intuitively, the constraint in (1) says that if ν |= 𝐴 ℓ or equivalently (ℓ, ν) ∈ T ♦ , then either (ℓ, ν) ∈ T which is equivalent to ν |= 𝜏 ℓ , or there exists a transition 𝜃 ∈ Θ ℓ , using which we can obtain a successor

(ℓ ′ , 𝜇 (ν)) ∈ T ♦ such that 𝑓 (ℓ ′ , 𝜇 (ν)) ≤ 𝑓 (ℓ, ν) -𝜖.
The latter is formalized by 𝜉 (𝜃 ).

In this step, the algorithm symbolically computes (1) and writes it in the following equivalent format:

∀ν ∈ R V , ν |= 𝐴 ℓ ∧ 𝜃 =(ℓ,ℓ ′ ,𝜑,𝜇) ¬𝜉 (𝜃 ) ⇒ ν |= 𝜏 ℓ (3)
Let 𝑃 ℓ be the LHS assertion in (3) above. Then 𝑃 ℓ is constructed from logical operations and atomic strict/non-strict linear inequalities over V. The coefficients in these linear inequalities contain the unknown 𝑐 ℓ,𝑖,𝑗 's and 𝑑 ℓ,𝑗 's.

The algorithm writes 𝑃 ℓ in disjunctive normal form, obtaining 𝑃 ℓ = 𝑃 ℓ,1 ∨𝑃 ℓ,2 ∨. . . ∨𝑃 ℓ,𝑝 , where each 𝑃 ℓ,𝑖 is a conjunction of linear inequalities over V. It then symbolically computes the following "constraint pair" for every 𝑃 ℓ,𝑖 :

𝛾 ℓ,𝑖 := 𝑃 ℓ,𝑖 , 𝜏 ℓ (4)
The algorithm computes these constraint pairs for every ℓ ∈ L and stores them in a set Γ. Every constraint pair 𝛾 = (𝜆, 𝜚 ) ∈ Γ consists of two parts. 𝜆 is a set of strict/nonstrict linear inequalities, while 𝜚 is a set of only non-strict linear inequalities. Informally, 𝛾 encodes the requirement that every inequality in 𝜚 be entailed by inequalities in 𝜆.

Example 7. Consider the system in Figure 2 together with the templates generated in Example 6. In this step, the algorithm considers location 1 ∈ L, and writes the constraint in (3):

𝐴 1 ∧ ¬𝜉 (𝜃 1 ) ∧ ¬𝜉 (𝜃 2 ) ∧ ¬𝜉 (𝜃 3 ) ⇒ 𝜏 1 (5)
Intuitively, this constraint says if we are at a T ♦ state in location 1, and cannot transition to another T ♦ with smaller 𝑓 -value, in other words ¬𝜉 (𝜃 1 ) ∧ ¬𝜉 (𝜃 2 ) ∧ ¬𝜉 (𝜃 3 ), then we must already be in a target state (satisfy 𝜏 1 ). There is no target state at location 1, so we can assume 𝜏 1 ≡ (-1 ≥ 0). The algorithm computes (5) symbolically:

𝑐 0 + 𝑐 1 • 𝑥 + 𝑐 2 • 𝑦 ≥ 0 ∧ 𝑐 3 + 𝑐 4 • 𝑥 + 𝑐 5 • 𝑦 ≥ 0 ∧ ¬(𝑥 < 𝑦 ∧ 𝑐 6 + 𝑐 7 • 𝑥 + 𝑐 8 • 𝑦 ≥ 0 ∧ 𝑐 9 + 𝑐 10 • 𝑥 + 𝑐 11 • 𝑦 ≥ 0 ∧ 𝑑 3 + 𝑑 4 • 𝑥 + 𝑑 5 • 𝑦 ≤ 𝑑 0 + 𝑑 1 • 𝑥 + 𝑑 2 • 𝑦 -𝜖) ∧ ¬(𝑥 < 𝑦 ∧ 𝑐 12 + 𝑐 13 •𝑥 + 𝑐 14 •𝑦 ≥ 0 ∧ 𝑐 15 + 𝑐 16 •𝑥 + 𝑐 17 •𝑦 ≥ 0 ∧ 𝑑 6 + 𝑑 7 • 𝑥 + 𝑑 8 • 𝑦 ≤ 𝑑 0 + 𝑑 1 • 𝑥 + 𝑑 2 • 𝑦 -𝜖) ∧ ¬(𝑥 ≥ 𝑦 ∧ 𝑐 18 + 𝑐 19 •𝑥 + 𝑐 20 •𝑦 ≥ 0 ∧ 𝑐 21 + 𝑐 22 •𝑥 + 𝑐 23 •𝑦 ≥ 0 ∧ 𝑑 9 + 𝑑 10 • 𝑥 + 𝑑 11 • 𝑦 ≤ 𝑑 0 + 𝑑 1 • 𝑥 + 𝑑 2 • 𝑦 -𝜖) ⇒ (-1 ≥ 0)
Intuitively, the first part of the constraint above models a state in T ♦ at location 1. The second part models the fact that it is not possible to take transition 𝜃 1 and reach another state in T ♦ at location 2 such that the 𝑓 -value decreases by at least 𝜖. The next two parts model similar constraints for 𝜃 2 and 𝜃 3 . Finally, the last part says that if no suitable transition is possible, then the current state must itself be a target, which is impossible in this case because there are no target states at location 1. Next, the algorithm writes the constraint above in DNF:

𝑃 1,1 ∨ 𝑃 1,2 ∨ . . . ∨ 𝑃 1,𝑝 ⇒ (-1 ≥ 0)
Just as before, the algorithm computes each of 𝑃 1,1 , . . . , 𝑃 1,𝑝 concretely in terms of 𝑥, 𝑦, 𝜖, 𝑐 𝑖 's and 𝑑 𝑖 's, but to save space, we omit the full expansion here, e.g. we can assume 𝑃 1,1 is:

𝑐 0 + 𝑐 1 • 𝑥 + 𝑐 2 • 𝑦 ≥ 0 ∧ 𝑐 3 + 𝑐 4 • 𝑥 + 𝑐 5 • 𝑦 ≥ 0 ∧ 𝑥 ≥ 𝑦 ∧ 𝑑 9 + 𝑑 10 • 𝑥 + 𝑑 11 • 𝑦 > 𝑑 0 + 𝑑 1 • 𝑥 + 𝑑 2 •
𝑦 -𝜖 This corresponds to the case where we cannot use either transition 𝜃 1 or 𝜃 2 because 𝑥 ≥ 𝑦, and also taking transition 𝜃 3 will lead to a state whose 𝑓 -value is not small enough. For each such 𝑃 1,𝑖 it generates a constraint pair (𝑃 1,𝑖 , 𝜏 1 ) = (𝑃 1,𝑖 , -1 ≥ 0).The algorithm handles other locations similarly.

Step 2b. Computing UIRW Constraint Pairs. This step is only performed when synthesizing a UIRW and is similar to its IRW variant in Step 2a above. Due to space constraints, we have relegated the details of this step to Appendix C.

Step 2c. Computing Non-negativity Constraints. In an IRW/UIRW, the ranking function 𝑓 should have non-negative value over T ♦ . Let ℓ ∈ L. The non-negativity condition at ℓ is equivalent to:

∀ν ∈ R V , ν |= 𝐴 ℓ ⇒ 𝑓 ℓ (ν) ≥ 0
To ensure this constraint, for every ℓ ∈ L, the algorithm adds the constraint pair ( 𝐴 ℓ , 𝑓 ℓ ≥ 0) to Γ. Example 8. Based the templates of Example 6, the algorithm creates the following non-negativity constraint pair 𝛾 = (𝜆, 𝜚 ), encoding 𝜆 ⇒ 𝜚, at location 1 ∈ L:

𝜆 : 𝑐 0 + 𝑐 1 • 𝑥 + 𝑐 2 • 𝑦 ≥ 0 𝑐 3 + 𝑐 4 • 𝑥 + 𝑐 5 • 𝑦 ≥ 0 𝜚 : ( 𝑑 0 + 𝑑 1 • 𝑥 + 𝑑 2 • 𝑦 ≥ 0)
Step 3. Applying Farkas' Lemma. The algorithm applies Corollary 1 to every constraint pair generated in the previous step to obtain a non-linear constraint system based on the template variables (i.e. 𝑐 ℓ,𝑖,𝑗 's and 𝑑 ℓ,𝑗 's), the ranking parameter 𝜖, and new variables defined in this step. Crucially, this non-linear constraint system does not include any of the variables in V. We now explain this step more concretely.

For every constraint pair 𝛾 = (𝜆, 𝜚 ) ∈ Γ, we know that 𝜆 is a set of strict/non-strict linear inequalities

{𝜆 𝑖,0 + #» 𝜆 𝑖 • #» V 𝑖 0} 𝑚
𝑖=1 , in which 𝑖 ∈ {>, ≥}. Moreover, 𝜚 is a set of non-strict inequalities and every inequality in 𝜚 should be entailed by

𝜆. Let 𝛼 0 + 𝛼 1 • 𝑣 1 + . . . + 𝛼 𝑟 • 𝑣 𝑟 ≥ 0 ≡ 𝛼 0 + #» 𝛼 • #» V ≥ 0 be
an inequality in 𝜚 . According to Corollary 1, there are three cases in which {𝜆 𝑖,0

+ 𝜆 𝑖 • V 𝑖 0} 𝑚 𝑖=1 entails 𝛼 0 + 𝛼 • V ≥ 0 : (i) 𝛼 0 + 𝛼 • V ≥ 0 is a non-negative combination of 1 ≥ 0 and {𝜆 𝑖,0 + 𝜆 𝑖 • V 𝑖 0} 𝑚 𝑖=1
, or (ii) -1 ≥ 0 can be derived as above, or (iii) 0 > 0 can be derived as above. The algorithm writes constraints that model each of the three cases above and then combines them disjunctively. Given their similarity, we only explain (i). The algorithm creates new variables 𝑦 0 , . . . , 𝑦 𝑚 ≥ 0 and computes the equality

𝛼 0 + 𝛼 • V = 𝑦 0 + 𝑚 𝑖=1 𝑦 𝑖 • (𝜆 𝑖,0 + 𝜆 𝑖 • V) (6)
The two sides above are linear expressions over V. As such, they are equal iff they agree on the coefficient of every term.

The algorithm equates the corresponding coefficients, and adds these equalities to the constraint system:

𝛼 0 = 𝑦 0 + 𝑚 𝑖=1 𝑦 𝑖 • 𝜆 𝑖,0
i.e. the constant factor should be equal on both sides, and

∀𝑗 ≠ 0 𝛼 𝑗 = 𝑚 𝑖=1 𝑦 𝑖 • 𝜆 𝑖,𝑗
i.e. the coefficient of every variable 𝑣 𝑗 ∈ V should be equal. The algorithm handles (ii) and (iii) similarly, except that in (iii) we should ensure that at least one strict inequality is used when trying to obtain 0 > 0. Hence, in this case, the algorithm also adds the extra constraint 𝑖 ∈ {> } 𝑦 𝑖 > 0.The algorithm performs the same operations for every constraint pair 𝛾 = (𝜆, 𝜚 ) and every linear inequality in 𝜚 and combines the resulting non-linear constraint systems conjunctively. Example 9. Consider the constraint pair 𝛾 = (𝜆, 𝜚 ) below, which was obtained in Example 7:

𝜆 :            𝑐 0 + 𝑐 1 • 𝑥 + 𝑐 2 • 𝑦 ≥ 0 𝑐 3 + 𝑐 4 • 𝑥 + 𝑐 5 • 𝑦 ≥ 0 𝑥 -𝑦 ≥ 0 𝑑 9 + 𝑑 10 • 𝑥 + 𝑑 11 • 𝑦 -𝑑 0 -𝑑 1 • 𝑥 -𝑑 2 • 𝑦 + 𝜖 > 0 𝜚 : (-1 ≥ 0)
We want to make sure that 𝜆 entails 𝜚 . By Corollary 1, either 𝜚 or -1 ≥ 0 or 0 > 0 should be a non-negative combination of inequalities in 𝜆. Here, 𝜚 is itself -1 ≥ 0, so we take two cases:

• -1 ≥ 0 is obtainable from 𝜆: The algorithm creates unknowns 𝑦 0 , 𝑦 1 , . . . , 𝑦 4 ≥ 0 and computes this equality:

𝑦 0 + 𝑦 1 • ( 𝑐 0 + 𝑐 1 •𝑥 + 𝑐 2 •𝑦) + 𝑦 2 • ( 𝑐 3 + 𝑐 4 •𝑥 + 𝑐 5 •𝑦) + 𝑦 3 • (𝑥 -𝑦)+ 𝑦 4 • ( 𝑑 9 + 𝑑 10 • 𝑥 + 𝑑 11 • 𝑦 -𝑑 0 -𝑑 1 • 𝑥 -𝑑 2 • 𝑦 + 𝜖) = -1.
Our program variables are 𝑥 and 𝑦. All other variables are created by the algorithm and we need to synthesize a value for them. The above is an equality between two polynomials in R[𝑥, 𝑦] that has to hold for all values of 𝑥 and 𝑦. Hence, the algorithm equates its corresponding coefficients:

• 𝑦 1 • 𝑐 1 + 𝑦 2 • 𝑐 4 + 𝑦 3 + 𝑦 4 • 𝑑 10 -𝑦 4 • 𝑑 1 = 0 (for 𝑥), • 𝑦 1 • 𝑐 2 + 𝑦 2 • 𝑐 5 -𝑦 3 + 𝑦 4 • 𝑑 11 -𝑦 4 • 𝑑 2 = 0 (for 𝑦), • 𝑦 0 + 𝑦 1 • 𝑐 0 + 𝑦 2 • 𝑐 3 + 𝑦 4 • 𝑑 9 -𝑦 4 • 𝑑 0 + 𝑦 4 • 𝜖 = -1.
• 0 > 0 is obtainable from 𝜆: The algorithm creates 5 new variables 𝑦 5 , . . . , 𝑦 9 and proceeds to obtain equalities over non-program variables in the exact same manner as in the previous case, except that it also adds the condition 𝑦 9 > 0.

Step 4. Computing Initial Constraints. An IRW/UIRW should have an initial state (ℓ 0 , ν) with ν |= 𝐼 . Equivalently,

∃ν 0 = (ν 0,1 , . . . , ν 0,𝑟 ) ∈ R V , ν |= 𝐴 ℓ 0 ∧ 𝐼 . (7)
By 𝑘-linearity of the system 𝑆, we know that the initial assertion 𝐼 is a conjunction of at most 𝑘 linear inequalities. Thus, the assertion above is a conjunction of at most 2𝑘 linear inequalities, and is equivalent to Sat( 𝐴 ℓ 0 ∧ 𝐼 ) ≠ ∅. In this step, the algorithm creates 𝑟 new variables ν 0,1 , . . . , ν 0,𝑟 , symbolically computes the linear inequalities in [START_REF] Alur | Timing Verification by Successive Approximation[END_REF], and adds them (conjunctively) to the non-linear constraint system. Example 10. For our example (Figure 2), the algorithm creates new variables ν 0,𝑥 and ν 0,𝑦 and computes the following:

𝑐 0 + 𝑐 1 • ν 0,𝑥 + 𝑐 2 • ν 0,𝑦 ≥ 0 ν 0,𝑥 ≥ 10 𝑐 3 + 𝑐 4 • ν 0,𝑥 + 𝑐 5 • ν 0,𝑦 ≥ 0 ν 0,𝑦 ≥ 10
The left constraints ensure that the valuation ν 0 = ( ν 0,𝑥 , ν 0,𝑦 ) satisfies 𝐴 1 and the right constraints enforce 𝐼 .

Step 5. Solving the Resulting Constraint System. Finally, the algorithm uses an off-the-shelf solver to solve the resulting non-linear constraint system. If the system is unsatisfiable, it reports that no 𝑘-linear IRW/UIRW exists. Otherwise, it obtains a solution 𝔰 of the non-linear constraint system. Let 𝔰( 𝑥) denote the value assigned by 𝔰 to variable 𝑥, and extend this definition in the natural way so to any expression 𝑒. The algorithm outputs 𝐴 ℓ := 𝔰( 𝐴 ℓ ) and 𝑓 ℓ := 𝔰( 𝑓 ℓ ), for all ℓ ∈ L, as the IRW/UIRW. Moreover, 𝔰( ν 0,1 , . . . , ν 0,𝑟 ) is the initial state, and 𝔰( 𝜖) is the decrease parameter for 𝑓 . Example 11. When the algorithm solves the non-linear constraints obtainted in the previous steps, it successfully synthesizes the following IRW (left) for

T = {(4, ν) | ν |= (𝑥 ≥ 𝑦 +8)} and UIRW (right) for T ′ = {(4, ν) | ν |= (𝑥 ≥ 𝑦 + 4)}: ℓ 𝐴 ℓ 𝑓 ℓ 1 𝑦 -2 ≤ 𝑥 ≤ 𝑦 -1 2 2 𝑦 -2 ≤ 𝑥 ≤ 𝑦 -1 1 3 -1 ≥ 0 -1 4 𝑥 ≥ 𝑦 + 8 0 𝜖 = 1, ν 0 = (11, 12) ℓ 𝐴 ℓ 𝑓 ℓ 1 𝑦 -0.6 ≤ 𝑥 ≤ 𝑦 -0.5 2 2 𝑦 -0.6 ≤ 𝑥 ≤ 𝑦 -0.5 1 3 𝑦 -0.6 ≤ 𝑥 ≤ 𝑦 -0.5 1 4 𝑥 ≥ 𝑦 + 4.4 0 𝜖 = 1, ν 0 = (11, 11.55)
Theorem 3 (Soundness, Proof in Appendix D). Given a 𝑘-linear system 𝑆 = (V, L, ℓ 0 , 𝐼, Θ), and a 𝑘-linear set T of target states, every solution of the non-linear constraint system solved in Step 5 of the algorithm above produces a valid 𝑘-linear IRW/UIRW for T in 𝑆. Theorem 4 (Completeness, Proof in Appendix D). Given a 𝑘-linear system 𝑆 = (V, L, ℓ 0 , 𝐼, Θ), and a 𝑘-linear set T of target states, every 𝑘-linear IRW/UIRW for T in 𝑆 is produced by some solution to the non-linear constraint system solved in Step 5 of the algorithm above.

Theorem 5 (Complexity, Proof in Appendix D). For fixed constants 𝑘 and 𝛽, given a 𝑘-linear 𝛽-branching system 𝑆 = (V, L, ℓ 0 , 𝐼, Θ), and a 𝑘-linear set T of target states, Steps 1-4 of the algorithm above lead to a polynomial-time reduction from the problem of generating a 𝑘-linear IRW/UIRW to solving a Quadratic Programming (QP) instance.

Synthesizing Polynomial IRWs/UIRWs

We now extend our algorithm to the case where the system, the target set, and the IRW/UIRW are all polynomial. Outline. Our algorithm is similar to the linear case. The main difference is that we should now handle polynomial constraint pairs. In other words, we have entailments of the form 𝑃 → 𝑄 in which both 𝑃 and 𝑄 are conjunctions of polynomial inequalities with unknown coefficients. Our goal is to rewrite such entailments as a system of quadratic equations over the unknown coefficients, so that we can use off-the-shelf QP-solvers. To do this, we will need polynomial counterparts of Farkas' Lemma and Corollary 1. Note that Farkas' lemma also provided a direct way of checking unsatisfiability of systems of linear inequalities. We will similarly need an unsatisfiability criterion for the polynomial case, since 𝑃 → 𝑄 holds vacuously if 𝑃 is unsatisfiable. We first introduce and develop these required mathematical tools and then present the algorithm. Strong Positivity. Let 𝑋 ⊆ R V be a set of valuations and 𝑔 ∈ R[V] a polynomial over V. We say that 𝑔 is strongly positive over 𝑋 , and write 𝑋 |= 𝑔 ≫ 0 (or simply 𝑔 ≫ 0 when 𝑋 is clear from context), if inf 𝑥 ∈𝑋 𝑔(𝑥) > 0. The real value 𝛿 := inf 𝑥 ∈𝑋 𝑔(𝑥) is called the positivity gap of 𝑔 over 𝑋 . 𝑔 is strongly greater than ℎ, denoted 𝑔 ≫ ℎ, iff 𝑔ℎ ≫ 0. Problem Definition. Given technical constants Υ 1 , . . . , Υ 4 ∈ N, a (𝑑, 𝑘)-polynomial system 𝑆 = (V, L, ℓ 0 , 𝐼, Θ), together with a set 𝜏 ℓ of at most 𝑘 strong polynomial inequalities of degree at most 𝑑 at every location ℓ ∈ L, synthesize a (𝑑, 𝑘)polynomial IRW/UIRW for a target set T that satisfies 𝜏 ℓ at every ℓ ∈ L, i.e. T ∩ {ℓ } × R V |= 𝜏 ℓ , or report that no such IRW/UIRW exists. The technical constants Υ 𝑖 are bounds on the degrees of various polynomials we construct as part of our algorithm. We will soon discuss them more concretely. Theorem 6 (Putinar's Positivstellensatz [START_REF] Putinar | Positive polynomials on compact semi-algebraic sets[END_REF]). Consider a set V = {𝑣 1 , . . . , 𝑣 𝑟 } of real-valued variables and the following system of 𝑚 polynomial inequalities over V:

Φ : 𝑔 1 (𝑣

∈ R[V].
This theorem automatically provides a criterion for satisfiability of Φ, given that the inequality -1 > 0 is entailed by Φ iff Φ is unsatisfiable. As before, we need a variant of this theorem that can handle strict inequalities in Φ. Corollary 2. Take the following system of 𝑚 polynomial inequalities over V = {𝑣 1 , . . . , 𝑣 𝑟 }: Φ : 𝑔 1 (𝑣 1 , . . . , 𝑣 𝑟 ) 1 0, . . . , 𝑔 𝑚 (𝑣 1 , . . . , 𝑣 𝑟 ) 𝑚 0 in which every 𝑔 𝑖 ∈ R[V] is a polynomial and 𝑖 ∈ {>, ≥}. Also, assume that there is some 𝑖 such that the set Sat(𝑔 𝑖 ≥ 0) is compact, or equivalently Sat(𝑔 𝑖 𝑖 0) is bounded. If Φ is satisfiable, then it entails a strong polynomial inequality 𝑔(𝑣 1 , . . . , 𝑣 𝑟 ) ≫ 0, iff there exist a constant 𝑦 0 ∈ (0, ∞) and polynomials ℎ 0 , . . . , ℎ 𝑚 ∈ R[V] such that

𝑔 = 𝑦 0 + ℎ 0 + 𝑚 𝑖=1 ℎ 𝑖 • 𝑔 𝑖 (8)
and every ℎ 𝑖 is a sum of squares, i.e.

ℎ 𝑖 = ℎ 2 𝑖,𝑗 for some polynomials ℎ 𝑖,𝑗 ∈ R[V].
Proof. It is obvious that any polynomial 𝑔 that can be represented as in Equation ( 8) is strongly positive over Sat(Φ) and has a positivity gap of at least 𝑦 0 > 0. Suppose Φ is satisfiable and entails 𝑔 ≫ 0 with positivity gap 𝛿, and choose 0

< 𝑦 0 < 𝛿 ′ < 𝛿. We have Sat(Φ) ⊆ Sat(𝑔 > 𝛿 ′ ) so Sat(Φ) = Sat(Φ) ⊆ Sat(𝑔 > 𝛿 ′ ) = Sat(𝑔 ≥ 𝛿 ′ ) ⊆ Sat(𝑔 > 𝑦 0 ). Hence, Φ entails 𝑔 -𝑦 0 > 0.
Applying Theorem 6 to Φ and 𝑔 -𝑦 0 leads to the desired result. □

This corollary characterizes strongly positive polynomials over Sat(Φ). However, it is only applicable when Φ is satisfiable. We are using this assumption in the proof above when we write Sat(Φ) = Sat(Φ). If Φ is unsatisfiable, then it vacuously entails any polynomial inequality, but the corollary above is inapplicable. Therefore, we also need a criterion for unsatisfiability of Φ. We use Hilbert's Strong Nullstellensatz to obtain a suitable criterion (Theorem 8). Theorem 7 (Strong Nullstellensatz [START_REF] Francis | Introduction to Commutative Algebra[END_REF]). Take V = {𝑣 1 , . . . , 𝑣 𝑟 } and let 𝑔 1 , . . . , 𝑔 𝑚 , 𝑔 ∈ R[V] be polynomials over V. Then exactly one of the following statements holds:

• There exists a valuation ν ∈ R V , such that 𝑔 1 (ν) = 𝑔 2 (ν) = . . . = 𝑔 𝑚 (ν) = 0, but 𝑔(ν) ≠ 0. • There exist 𝛼 ∈ N∪ {0} and polynomials ℎ 1 , . . . , ℎ 𝑚 ∈ R[V] such that 𝑚 𝑖=1 ℎ 𝑖 • 𝑔 𝑖 = 𝑔 𝛼 .
The following theorem is the main theoretical basis of our synthesis algorithm and characterizes unsatisfiability of Φ when it contains both strict and non-strict inequalities: Theorem 8. Take a set V = {𝑣 1 , . . . , 𝑣 𝑟 } of real-valued variables and a system of 𝑚 polynomial inequalities over V: Φ : 𝑔 1 (𝑣 1 , . . . , 𝑣 𝑟 ) 1 0, . . . , 𝑔 𝑚 (𝑣 1 , . . . , 𝑣 𝑟 ) 𝑚 0 where each 𝑔 𝑖 ∈ R[V] is a polynomial and 𝑖 ∈ {>, ≥}. Φ is unsatisfiable iff at least one of the following statements holds: (i) There exist a constant 𝑦 0 ∈ (0, ∞) and sum-of-square polynomials ℎ 0 , . . . , ℎ 𝑚 ∈ R[V] such that

-1 = 𝑦 0 + ℎ 0 + 𝑚 𝑖=1 ℎ 𝑖 • 𝑔 𝑖 .
(ii) There exist 𝛼 ∈ N ∪ {0} and ℎ 1 , . . . , ℎ 𝑚 ∈ R[V * ] for V * = V ∪ {𝑤 1 , . . . , 𝑤 𝑚 }, such that for some 1 ≤ 𝑗 ≤ 𝑚 with 𝑗 ∈ {>}, we have

𝑤 2•𝛼 𝑗 = 𝑚 𝑖=1 ℎ 𝑖 • (𝑔 𝑖 -𝑤 2 𝑖 ).
Proof. If Φ is satisfiable, then it cannot entail -1 > 0, so (i) is impossible. We now show that (ii) implies unsatisfiability of Φ as well. Define g𝑖 (𝑣

g𝑖 𝑗 • 𝑤 2•(𝛼-𝑖) 𝑗 = 𝑤 2•𝛼 𝑗 + ℎ ′ 𝑗 • g𝑗 for some ℎ ′ 𝑗 ∈ R[V * ]. So, letting ℎ ′′ 𝑖 = ℎ 𝑖 for 𝑖 ≠ 𝑗 and ℎ ′′ 𝑗 = ℎ 𝑗 + ℎ ′ 𝑗 , we have 𝑔 𝛼 𝑗 = 𝑚 𝑖=1 ℎ ′′ 𝑖 • (𝑔 𝑖 -𝑤 2 𝑖 ) Let ν ∈ R V ∩Sat(Φ).
We extend ν to ν * ∈ R V * as follows: for every 𝑤 𝑖 , let ν * (𝑤 𝑖 ) = ν(𝑔 𝑖 ). So, we have ν * (𝑔 𝑖 -𝑤 2 𝑖 ) = 0, and hence the RHS of the equation above is 0 at ν * . On the other hand, we have ν

* (𝑔 𝛼 𝑗 ) = ν(𝑔 𝛼 𝑗 ) = ν(𝑔 𝑗 ) 𝛼 > 0.
This contradiction shows that Φ is unsatisfiable.

We now prove the other side. Suppose that Φ is unsatisfiable. If Φ is unsatisfiable, then it entails -1.5 > 0 and hence we can apply Theorem 6 to write -1.5 = ℎ 0 + 𝑚 𝑖=1 ℎ 𝑖 • 𝑔 𝑖 for some sum-of-squares polynomials ℎ 𝑖 , which is equivalent to -1 = 0.5 + ℎ 0 + 𝑚 𝑖=1 ℎ 𝑖 • 𝑔 𝑖 , hence leading to case (i) above. The only remaining case is if Φ is satisfiable but Φ is not. Reorder the inequalities in Φ so that the non-strict inequalities appear first. Let 𝑗 be the smallest index for which Φ[1 . . . 𝑗], i.e. the set of first 𝑗 inequalities in Φ, is unsatisfiable. By definition, Φ[1 . . . 𝑗 -1] is satisfiable and hence

Sat(Φ[1 . . . 𝑗 -1]) = Sat(Φ[1 . . . 𝑗 -1]). Moreover, since Φ[1 . . . 𝑗] = Φ[1 . . . 𝑗 -1] ∧ 𝑔 𝑗 > 0 is unsatisfiable, we know that Φ[1 . . . 𝑗 -1] entails 𝑔 𝑗 ≤ 0. In other words, Sat(Φ[1 . . . 𝑗 -1]) ⊆ Sat(𝑔 𝑗 ≤ 0). Taking closures from both sides shows that Φ[1 . . . 𝑗 -1] entails 𝑔 𝑗 ≤ 0. So, Φ[1 . . . 𝑗] en- tails 𝑔 𝑗 = 0. Define g𝑖 (𝑣 1 , . . . , 𝑣 𝑟 , 𝑤 1 , . . . , 𝑤 𝑚 ) := 𝑔 𝑖 (𝑣 1 , . . . , 𝑣 𝑟 )- 𝑤 2
𝑖 . We claim there is no valuation ν * ∈ R V * such that for all 1 ≤ 𝑖 ≤ 𝑗, g𝑖 (ν * ) = 0, but 𝑔 𝑗 (ν * ) ≠ 0. To prove this, suppose that such a valuation exists, and let ν be its restriction to V. For each 1 ≤ 𝑖 ≤ 𝑗, since g𝑖 (ν * ) = 0, we have 𝑔 𝑖 (ν) ≥ 0. Moreover, 𝑔 𝑗 (ν) = 𝑔 𝑗 (ν * ) ≠ 0. This is a contradiction with the previously proven fact that Φ[1 . . . 𝑗] entails 𝑔 𝑗 = 0. Applying the Strong Nullstellensatz (Theorem 7) to the g𝑖 's and 𝑔 𝑗 , we conclude that there exist a non-negative integer 𝛼 and polynomials h1 , . . . ,

h𝑗 ∈ R[V * ] such that 𝑔 𝛼 𝑗 = 𝑗 𝑖=1 h𝑖 • g𝑖 . Note that 𝑔 𝛼 𝑗 = g𝑗 + 𝑤 2 𝑗 𝛼 = 𝛼 𝑖=0 𝛼 𝑖 g𝑖 𝑗 • 𝑤 2•(𝛼-𝑖) 𝑗 = 𝑤 2•𝛼 𝑗 + ℎ ′ 𝑗 • g𝑗 for some ℎ ′ 𝑗 ∈ R[V * ].
Defining ℎ 𝑖 = h𝑖 for all 𝑖 ≠ 𝑗, and ℎ 𝑗 = h𝑗ℎ ′ 𝑗 , we get

𝑤 2•𝛼 𝑗 = 𝑗 𝑖=1 ℎ 𝑖 • g𝑖 = 𝑗 𝑖=1 ℎ 𝑖 • (𝑔 𝑖 -𝑤 2 𝑖
). The Synthesis Algorithm. As in the linear case, our algorithm consists of 5 steps. The main differences are in Steps 1 and 3. In Step 1, our algorithm should now generate a polynomial template. Moreover, in Step 3, it employs Corollary 2 and Theorem 8 for characterizing entailment. The other steps are exactly like our previous algorithm.

Step 1. Setting up a template. The algorithm symbolically computes the set of monomials of degree at most 𝑑 over V:

𝑀 𝑑 (V) := {𝔪 1 , 𝔪 2 , . . . , 𝔪 𝑢 } := {𝑣 𝛼 1 1 • 𝑣 𝛼 2 2 • . . . • 𝑣 𝛼 𝑟 𝑟 | 𝛼 1 , . . . , 𝛼 𝑟 ∈ N ∪ {0} ∧ 𝛼 1 + . . . + 𝛼 𝑟 ≤ 𝑑 }.
It sets up templates for 𝐴 ℓ and 𝑓 ℓ at every ℓ ∈ L :

𝐴 ℓ :          𝑐 ℓ,1,1 • 𝔪 1 + . . . + 𝑐 ℓ,1,𝑢 • 𝔪 𝑢 ≥ 0 . . . 𝑐 ℓ,𝑘,1 • 𝔪 1 + . . . + 𝑐 ℓ,𝑘,𝑢 • 𝔪 𝑢 ≥ 0 𝑓 ℓ = 𝑑 ℓ,1 • 𝔪 1 + . . . + 𝑑 ℓ,𝑢 • 𝔪 𝑢
The 𝑐 ℓ,𝑖,𝑗 's and 𝑑 ℓ,𝑗 's are unknown variables for which we should synthesize a value such that the 𝐴 ℓ 's and 𝑓 ℓ 's form an IRW or a UIRW. We do not need to add a separate constant factor to our templates because 1 ∈ 𝑀 𝑑 (V).

Step 2. Computing Constraint Pairs. Steps 2a-2c are the same as in Section 3.1. However, note that the resulting constraint pairs 𝛾 = (𝜆, 𝜚 ) ∈ Γ are now polynomial. Concretely, 𝜆 is a set of strict or non-strict polynomial inequalities over V and 𝜚 is a set of strong polynomial inequalities over V.

Step

𝑔 = 𝑦 0 + ℎ 0 + 𝑚 𝑖=1 ℎ 𝑖 • 𝑔 𝑖 .
Finally, the algorithm translates this equality to quadratic equalities over template variables. The systems of quadratic equalities generated in (i)-(iii) above are combined disjunctively. Steps 4 and 5. These steps are the same as in Section 3.1.

Example 12. Suppose that Υ 1 = Υ 2 = Υ 3 = Υ 4 = 1, and the algorithm is in Step 3, handling the following constraint pair:

𝜆 : 𝑐 1 • 𝑥 > 0 𝑐 2 • 𝑦 ≥ 0 𝜚 : ( 𝑐 3 • 𝑥 • 𝑦 + 𝑐 4 ≫ 0)
The algorithm considers the following cases: (i) It generates three new template polynomials

ℎ 0 = 𝜂 0,1 + 𝜂 0,2 • 𝑥 + 𝜂 0,3 • 𝑦 ℎ 1 = 𝜂 1,1 + 𝜂 1,2 • 𝑥 + 𝜂 1,3 • 𝑦 ℎ 2 = 𝜂 2,1 + 𝜂 2,2 • 𝑥 + 𝜂 2,3 • 𝑦
and computes a quadratic system of (in)equalities over the 𝜂 𝑖,𝑗 's that ensures every ℎ 𝑖 is a sum-of-squares (See Appendix F for details). The algorithm then computes the following equality symbolically (with 𝑦 0 > 0):

-1 = 𝑦 0 + ℎ 0 + ℎ 1 • 𝑐 1 • 𝑥 + ℎ 2 • 𝑐 2 • 𝑦
and rewrites it as quadratic equalities between the unknown variables by equating the coefficients of corresponding terms on the two sides of the polynomial equality. Intuitively, if there is a valuation for the unknown variables that satisfies these constraints, then -1 is a combination of 𝑐 1 • 𝑥, 𝑐 2 • 𝑦 and sum-of-square polynomials. Hence, 𝜆 is unsatisfiable. (ii) The algorithm creates two new program variables 𝑤 1 , 𝑤 2 and sets up the following templates:

ℎ 3 = 𝜂 3,1 + 𝜂 3,2 • 𝑥 + 𝜂 3,3 • 𝑦 + 𝜂 3,4 • 𝑤 1 + 𝜂 3,5 • 𝑤 2 ℎ 4 = 𝜂 4,1 + 𝜂 4,2 • 𝑥 + 𝜂 4,3 • 𝑦 + 𝜂 4,4 • 𝑤 1 + 𝜂 4,5 • 𝑤 2 .
Unlike the previous case, ℎ 3 and ℎ 4 need not be sum-ofsquares. It then writes the equality:

𝑤 2 1 = ℎ 3 • ( 𝑐 1 • 𝑥 -𝑤 2 1 ) + ℎ 4 • ( 𝑐 2 • 𝑦 -𝑤 2
2 ), and converts this polynomial equality to quadratic equalities by equating corresponding coefficients. Note that the LHS and RHS of the polynomial equality above are in R[𝑥, 𝑦, 𝑤 1 , 𝑤 2 ]. By Theorem 8, any solution to these constraints serves as a proof for unsatisfiability of 𝜆. (iii) The algorithm generates the following templates:

ℎ 5 = 𝜂 5,1 + 𝜂 5,2 • 𝑥 + 𝜂 5,3 • 𝑦, ℎ 6 = 𝜂 6,1 + 𝜂 6,2 • 𝑥 + 𝜂 6,3 • 𝑦, ℎ 7 = 𝜂 7,1 + 𝜂 7,2 • 𝑥 + 𝜂 7,3 • 𝑦;
enforces them to be sum-of-squares just as in case (i) above (Appendix F) and writes the polynomial equality:

𝑐 3 • 𝑥 • 𝑦 + 𝑐 4 = 𝑦 1 + ℎ 5 + ℎ 6 • 𝑐 1 • 𝑥 + ℎ 7 • 𝑐 2 • 𝑦 in which 𝑦 1 > 0.
It handles it similarly to the previous cases. This is again a polynomial equality in R[𝑥, 𝑦]. The algorithm combines the systems of quadratic inequality in (i)-(iii) above disjunctively.

It is now easy to obtain the following theorems, whose proofs are similar to the linear case: Theorem 9 (Soundness). Given a (𝑑, 𝑘)-polynomial system 𝑆 = (V, L, ℓ 0 , 𝐼, Θ), and a set 𝜏 ℓ of at most 𝑘 polynomial inequalities of degree 𝑑 or less at every ℓ ∈ L, every solution of the non-linear constraint system solved in Step 5 of the algorithm above produces a valid (𝑑, 𝑘)-polynomial IRW/UIRW for a target set T that satisfies 𝜏 ℓ at every ℓ ∈ L. Theorem 10 (Semi-completeness). Take a (𝑑, 𝑘)-polynomial system 𝑆 = (V, L, ℓ 0 , 𝐼, Θ) and a set 𝜏 ℓ of at most 𝑘 strong polynomial inequalities of degree 𝑑 or less at every ℓ ∈ L. Let 𝑊 = (T ♦ , 𝑓 , 𝜖) be an explicitly bounded (𝑑, 𝑘)-polynomial IR-W/UIRW for a target set T that satisfies 𝜏 ℓ at every ℓ ∈ L.

If large enough values are assigned to technical constants Υ 1 , . . . , Υ 4 , the witness 𝑊 is produced by some solution of the non-linear constraint system solved in Step 5 of the algorithm. Theorem 11 (Complexity). For fixed constants 𝑘, 𝑑 and 𝛽, and technical constants Υ 1 , . . . , Υ 4 , given a (𝑘, 𝑑)-polynomial 𝛽-branching system 𝑆 = (V, L, ℓ 0 , 𝐼, Θ), and a set 𝜏 ℓ of at most 𝑘 polynomial inequalities of degree 𝑑 or less at every ℓ ∈ L, Steps 1-4 of the algorithm above lead to a polynomial-time reduction from the problem of generating a (𝑘, 𝑑)-polynomial IRW/UIRW to solving a QP instance. Remark 2. Theorem 10 provides semi-completeness, i.e. completeness when the chosen technical constants are large enough, because the stellensätze have no bound on the degree of polynomials in their characterizations, but we have to fix a degree in our algorithm. We use Υ 1 , . . . , Υ 4 for this purpose. Such results arise routinely in constraint-based termination analysis [START_REF] Chatterjee | Termination analysis of probabilistic programs through Positivstellensatz's[END_REF] and invariant generation [START_REF] Feng | Finding polynomial loop invariants for probabilistic programs[END_REF]. In practice, solutions are often found with small technical constants (see Section 4). Remark 3. We presented our approach for conjunctive IRWs, i.e. at each label, the set of states that are in T ♦ was defined by a conjunction of linear/polynomial inequalities. However, the approach can be extended to handle disjunction, as well. In the linear case, this can be achieved using Motzkin's Transposition Theorem [START_REF] Ben-Israel | Motzkin transposition theorem[END_REF][START_REF] Ts Motzkin | Beitrage zur Theorie der linearen Ungleichungen[END_REF]. In the polynomial case, we can apply Stengle's Positivstellensatz [START_REF] Stengle | A nullstellensatz and a positivstellensatz in semialgebraic geometry[END_REF] instead of Putinar's. In both cases, the completeness will be unaffected but the runtime will increase.

Experimental Results

Implementation. We implemented our algorithms for IRW synthesis in Python using SymPy [START_REF] Meurer | SymPy: symbolic computing in Python[END_REF] for symbolic computations. We also added several heuristics for improving performance. Notably, we used Z3 [START_REF] De | Z3: An efficient SMT solver[END_REF] to identify and discard unsatisfiable or tautological constraint pairs. The QPs were then solved by LOQO [START_REF] Robert | LOQO: An interior point code for quadratic programming[END_REF]. All results were obtained on an Intel Core i5-2540M (2.6 GHz) machine with 8 GB RAM running Ubuntu, with a limit of 1800 seconds per task. Linear Benchmarks. We used C benchmarks from SV-COMP 2020 [START_REF] Beyer | Advances in automatic software verification: SV-COMP 2020[END_REF]. We considered tasks in the "Reachability/Safety" category of the competition in which the error states were reachable. Of these, we removed benchmarks that cannot be soundly modeled as transition systems with real variables, i.e. benchmarks containing pointers, arrays, etc. This left us with 30 benchmarks. The only source of non-determinism in these benchmarks is their input. We handle this by synthesizing an initial valuation as part of the IRW. In absence of integer division/mod, our approach can soundly handle integer variables by requiring that the synthesized initial valuation (in Step 4) is integral, i.e. ν 0,𝑖 's are integer unknowns. This leads to a mixed QP instance, but there are very few integral variables, and the problem remains tractable. Polynomial Benchmarks. For the polynomial case, all standard benchmarks focused on safety. Hence, we created 6 simple programs with complex reachability structure to showcase the strengths of our approach. Due to space restrictions, these examples are put in Appendix G. Specifically, it is noteworthy that these benchmarks demonstrate the fact that our algorithm's success is not dependent on the length or proportion of paths that reach the target set T. Previous Tools. We compare our approach against the two best-performing tools in the Reachability/Safety category of SV-COMP 2020, namely VeriAbs [START_REF] Afzal | Veri-Abs : Verification by Abstraction and Test Generation[END_REF] and CPAchecker [START_REF] Beyer | CPAchecker: A Tool for Configurable Software Verification[END_REF], as well as the most related previous work, Dangerzone [START_REF] David | Danger invariants[END_REF]. Linear Results. The results over linear benchmarks are summarized in Table 1. Our approach could handle every linear reachability benchmark in SV-COMP 2020. It is noteworthy that according to the SV-COMP results, none of the participating model checkers could handle all benchmarks of Table 1. CPAchecker times out on 9 of the instances, whereas VeriAbs fails on only 1 instance. Additionally, Dangerzone is not applicable to 4 of the instances and times out on 14.

By manual inspection, we realized that CPAChecker and Ve-riAbs are faster than our approach when reachability can be attained using liberal abstractions and a relatively short path (benchmarks towards the top of Table 1). This is not surprising, given that in these situations, abstract interpretation and symbolic execution are considerably faster than quadratic constraint solving. However, as the paths become longer and sparser (towards the the bottom of Table 1), the advantages of our approach begin to show. When the paths are long, i.e. thousands of steps of program execution, CPAchecker always fails to verify the instance. VeriAbs manages to handle these instances by a combination of ideas from loop pruning, loop summarization, abstract interpretation and bounded model checking. However, this comes with a considerable overhead, leading to a much worse performance in comparison with our approach. Finally, Dangerzone, while being able to efficiently find some of the deep bugs, fails on other instances. This is also not surprising, given that Dangerzone does not provide completeness guarantees and relies on evolutionary methods in its synthesis process.

Our runtime is mainly dependent on the size of the QP instance, which in turn depends on the size of the program and the degree of the IRW. In contrast, previous approaches' runtimes depend on the depth of the bug. Some previous approaches, e.g. VeriAbs, also perform abstraction, so their runtime becomes dependent on the number of refinements needed to detect the bug. Informally speaking, other approaches are faster in detecting shallow and wide bugs, whereas ours performs better on deep and narrow bugs. Of course, our performance also depends on the QP solver.

Nested Loop Benchmark. The program below illustrates the core of deep-nested, the only linear benchmark that could be handled by neither VeriAbs nor CPAchecker: In this program, 𝑀 is the largest value that fits into a 32-bit integer. We also ran these tools over this benchmark with an extended time limit of 12 hours, but they timed out. Moreover, according to SV-COMP results, no other participating model checker could handle this example, either. We believe this is because the target state can only be reached after an enormously-long path. Moreover, the target set is quite thin and even the smallest loss of precision in abstraction can cause a failure to prove reachability. However, this particular benchmark is not at all challenging for our method. The runtime of our method does not depend on the length of the paths, and we do not perform abstraction. Moreover, our approach is complete for linear IRWs. As such, it can easily prove reachability in this program.

Polynomial Results. Table 2 shows our experimental results over 6 polynomial instances. Informally, sqrt1 is a simple program that given an input integer 𝑛 ≥ 1 computes 𝑠 = ⌊ √ 𝑛⌋ by trying every possible integer starting from 1. The goal is to (choose a value for 𝑛 so as to) reach a state with 𝑛 -𝑠 > 10 5 . sqrt2 is a more clever variant of the same program that doubles the current value in a single step when the doubled value does not exceed ⌊ √ 𝑛⌋. sum is a program that sums up all the integers from 1 to 𝑛. The goal is to synthesize a value for 𝑛 such that the sum falls in a specific interval. sum2 is a similar benchmark in which the program sums squares of all integers from 1 to 𝑛. In robot1, two robots are put in the same position in a 2d plane. At each step, each robot non-deterministically chooses to move one unit either upwards or to the right. The goal is to reach a state where the square of the distance between the robots is more than 10 5 . In robot2, the same two robots are placed on the lower-right and upper-left corners of a square of side length 10 4 . The goal is to show that they can reach a distance of less than 10 from each other. See Appendix G for details.

As in the linear case, we observe that CPAchecker and VeriAbs can handle cases where the path reaching the targets is short, and when there is no combinatorial explosion in the number of paths due to repeated nondeterministic choice. Notably, CPAchecker can handle sqrt2 but not sqrt1. The only difference between these two programs is that sqrt2 is more efficient and hence the path to targets is shorter. We also observe that the various other techniques used by VeriAbs, which made it more successful in the linear case, do not extend well to the polynomial case. In contrast, our approach and Dangerzone are able to handle reachability with long paths. Our approach is the only one that can handle all of the polynomial benchmarks. This is due to its semicompleteness over polynomial IRWs.

Conclusion and Future Work

We proposed IRWs for reachability analysis over imperative programs. Our approach synthesizes an under-approximation of the set of program states that can reach the target, then uses a ranking argument to ensure eventual reachability. We proved that our approach is sound and complete when there is no restriction over the form of inductive reachability witnesses, and presented automated sound and semi-complete algorithms for synthesizing linear and polynomial inductive reachability witnesses. In practice, our experimental results show that our automated approaches can solve instances beyond the reach of previous methods.

An interesting future direction is to incorporate more advanced ranking-function synthesis methods such as lexicographic ranking functions [START_REF] Alias | Multi-dimensional Rankings, Program Termination, and Complexity Bounds of Flowchart Programs[END_REF][START_REF] Amir | On Multiphase-Linear Ranking Functions[END_REF][START_REF] Aaron R Bradley | Linear ranking with reachability[END_REF] into reachability analysis. Another direction is to consider how our approach can be extended to automate the search for proofs in incorrectness logic [START_REF] O'hearn | Incorrectness logic[END_REF]. Invariant generation has been successfully used in automating aspects of Hoare logic and termination analysis [START_REF] Chatterjee | Termination analysis of probabilistic programs through Positivstellensatz's[END_REF][START_REF] Chatterjee | Non-polynomial worst-case analysis of recursive programs[END_REF][START_REF] Michael | Synthesis of linear ranking functions[END_REF][START_REF] Podelski | A Complete Method for the Synthesis of Linear Ranking Functions[END_REF]. Given that our witnesses are natural duals of inductive invariants, we expect that this direction will be fruitful.

A Universal Reachability Witnesses

𝐼 : 𝑖 = 𝑠 = 0 ∧ 𝑛 ≥ 0 𝑎 : while 𝑖 ≤ 𝑛 : 𝑏 :

(𝑠, 𝑖) := (𝑠 + 1, 𝑖 + 1)

𝑐 : □ (𝑠, 𝑖) := (𝑠 + 2, 𝑖 + 1) 𝑑 : a d c b ( s , i ) ↦ ( s + 1 , i + 1 ) ( s , i ) ↦ ( s + 2 , i + 1 )
i In this section, we provide detailed definitions and an example of the Universal Inductive Reachability Witnesses (UIRWs).

> n i ≤ n i ≤ n θ 1 θ 2 θ 4 θ 5 θ 3 θ 6
Universal T-inductive Sets. Given a set T ⊆ Σ of target states, a set T ♦ ⊆ Σ is called universally T-inductive if for every 𝜎 ∈ T ♦ \ T and every successor 𝜎 ′ of 𝜎, we also have

𝜎 ′ ∈ T ♦ .
The idea behind universal T-inductive sets is that any execution of the program that starts in such a set T ♦ will either reach T or one can prove using induction that it will never leave T ♦ , no matter how the non-determinism is resolved.

Universal T-ranking Functions. Given a universal T-inductive set T ♦ , a function 𝑓 : T ♦ → [0, ∞) is called a universal Tranking function with parameter 𝜖 > 0, if for every 𝜎 ∈ T ♦ \T and every successor 𝜎 ′ of 𝜎, we have 𝑓 (𝜎 ′ ) ≤ 𝑓 (𝜎) -𝜖.

Universal Inductive Reachability Witnesses (UIRWs). Given a set T of target states in a system 𝑆 = (V, L, ℓ 0 , 𝐼, Θ), a Universal Inductive Reachability Witness for T is a tuple (T ♦ , 𝑓 , 𝜖) such that:

• T ♦ is a universal T-inductive set; • 𝜖 ∈ (0, ∞); • 𝑓 : T ♦ → [0, ∞) is a universal T-ranking function with
parameter 𝜖; • There exists a valuation ν ∈ R V such that (ℓ 0 , ν) ∈ T ♦ and ν |= 𝐼 .

Example 13. Figure 3 shows a simple program together with its representation as a transition system. Let T = {(𝑑, ν) | ν(𝑠) ≥ 20}, i.e. the target is reaching point 𝑑 with an 𝑠 value of more than 20. Let T ♦ := {(ℓ, ν) | ν |= 𝐴 ℓ } and 𝑓 (ℓ, ν) := 𝑓 ℓ (ν) be defined as follows:

ℓ 𝐴 ℓ 𝑓 ℓ 𝑎 𝑛 ≥ 50 ∧ 𝑠 ≥ 𝑖 ≥ 0 ∧ 𝑛 + 1 ≥ 𝑖 𝑛 + 1.5 -𝑖 𝑏 𝑛 ≥ 50 ∧ 𝑠, 𝑛 ≥ 𝑖 ≥ 0 𝑛 + 1 -𝑖 𝑐 𝑛 ≥ 50 ∧ 𝑠, 𝑛 ≥ 𝑖 ≥ 0 𝑛 + 1 -𝑖 𝑑 𝑠 ≥ 50 0
It is easy to check that (T ♦ , 𝑓 , 0.5) is a UIRW for T. Intuitively, this guarantees that if a run starts with an initial valuation that satisfies 𝐴 𝑎 , it will definitely reach a target state.

(i) Given that T is reachable, by definition there exists a path 𝜋 = (ℓ 0 , ν 0 ), 𝜃 0 , . . . , (ℓ 𝑛 , ν 𝑛 ) such that (ℓ 𝑛 , ν 𝑛 ) ∈ T and ν 0 |= 𝐼 . Without loss of generality, we choose such a 𝜋 that is prefix-minimal, i.e. that no prefix of 𝜋 has the same properties. Let

T ♦ = {(ℓ 𝑖 , ν 𝑖 )|0 ≤ 𝑖 ≤ 𝑛}, then T ♦ is T-inductive, because (ℓ 𝑛 , ν 𝑛 ) ∈
T and for every 𝑖 ≠ 𝑛, the state (ℓ 𝑖 , ν 𝑖 ) can be succeeded by (ℓ 𝑖+1 , ν 𝑖+1 ). Let 𝑓 : T ♦ → [0, ∞) be defined as follows: 𝑓 (ℓ 𝑖 , ν 𝑖 ) := 𝑛 -𝑖.

It is easy to verify that (T ♦ , 𝑓 , 1) is an IRW for T. (ii) We define Σ 𝑘 ⊆ Σ as the set of all states such that every semi-path of length 𝑘 starting in these states is guaranteed to visit T. Note that Σ 0 = T and if 𝜎 ∈ Σ 𝑘 \ T, then by definition every successor 𝜎 ′ of 𝜎 must be in Σ 𝑘-1 . Let T ♦ = ∞ 𝑖=0 Σ 𝑘 , and for every 𝜎 ∈ T ♦ , define 𝑓 (𝜎) := min{𝑘 | 𝜎 ∈ Σ 𝑘 }. It is easy to prove by definition-chasing that (T ♦ , 𝑓 , 1) is a UIRW. □

C Computing Universal Constraint Pairs

In this section, we provide a detailed description of Step 2b of our algorithm, which aims to generate constraint pairs for universal inductive reachability witnesses.

Step 2b. Computing UIRW Constraint Pairs. This step is only performed when synthesizing a UIRW and is similar to its IRW variant in Step 2a. In a UIRW, the universal T-inductive set T ♦ should satisfy the condition that for every state 𝜎 ∈ T ♦ \ T, every successor 𝜎 ′ of 𝜎 is also in T ♦ . Moreover, given that 𝑓 is a universal T-ranking function, we must have 𝑓 (𝜎 ′ ) ≤ 𝑓 (𝜎) -𝜖 for every such 𝜎 ′ . Let ℓ ∈ L be a location. The UIRW properties at ℓ are equivalent to:

∀ν ∈ R V , ν |= 𝐴 ℓ ⇒ ν |= 𝜏 ℓ ∨ 𝜃 =(ℓ,ℓ ′ ,𝜑,𝜇) 𝜁 (𝜃 ) (9)
where 𝜁 (𝜃 ) = 𝜁 (ℓ, ℓ ′ , 𝜑, 𝜇) is defined as:

𝜁 (𝜃 ) := ν |= 𝜑 ⇒ 𝜇 (ν) |= 𝐴 ℓ ′ ∧ 𝑓 ℓ ′ (𝜇 (ν)) ≤ 𝑓 ℓ (ν) -𝜖 (10)
Informally, the constraint in [START_REF] Francis | Introduction to Commutative Algebra[END_REF] says that if ν |= 𝐴 ℓ or equivalently (ℓ, ν) ∈ T ♦ , then either (ℓ, ν) ∈ T, i.e. ν |= 𝜏 ℓ , or for every transition 𝜃 from ℓ the assertion 𝜉 (𝜃 ) holds, i.e. if the transition is possible (ν |= 𝜑), then the successor state (ℓ ′ , 𝜇 (ν)) is also in T ♦ , and the 𝑓 value decreases by at least 𝜖 when going to this successor. As in the previous case, the algorithm computes (9) symbolically and writes it in the following equivalent format:

∀ν ∈ R V , ν |= 𝐴 ℓ ∧ 𝜃 =(ℓ,ℓ ′ ,𝜑,𝜇) ¬𝜁 (𝜃 ) ⇒ ν |= 𝜏 ℓ (11)
Let 𝑄 ℓ be the LHS assertion above. Similar to Step 2a, 𝑄 ℓ is constructed form logical operations and atomic strict/nonstrict linear inequalities over V, and its coefficients include the unknown template variables 𝑐 ℓ,𝑖,𝑗 's and 𝑑 ℓ,𝑗 's defined in Step 1. The algorithm writes 𝑄 ℓ in disjunctive normal form, hence obtaining 𝑄 ℓ = 𝑄 ℓ,1 ∨𝑄 ℓ,2 ∨. . .∨𝑄 ℓ,𝑞 in which each 𝑄 ℓ,𝑖 is a conjunction of strict/non-strict linear inequalities over V. It then computes the following constraint pair symbolically: 𝛾 ′ ℓ,𝑖 := (𝑄 ℓ,𝑖 , 𝜏 ℓ ) The algorithm performs these operations for every location ℓ ∈ L and stores all the resulting 𝛾 ′ ℓ,𝑖 constraint pairs in a set Γ.

Example 14. In our running example (Figure 2), we are looking for a linear UIRW for the target set T ′ = {(4, ν) | ν |= (𝑥 ≥ 𝑦 + 4)}. In this step, the algorithm creates constraints at every location. We now demonstrate how the process works for location 3. At location 3, the algorithm considers

𝐴 3 ∧ ¬𝜁 (𝜃 5 ) ⇒ 𝜏 3
and symbolically computes it as:

𝑐 12 + 𝑐 13 • 𝑥 + 𝑐 14 • 𝑦 ≥ 0 ∧ 𝑐 15 + 𝑐 16 • 𝑥 + 𝑐 17 • 𝑦 ≥ 0 ∧ ¬(1 ≥ 0 ⇒ ( 𝑐 18 +5• 𝑐 19 + 𝑐 19 •𝑥+ 𝑐 20 •𝑦 ≥ 0 ∧ 𝑐 21 +5• 𝑐 22 + 𝑐 22 •𝑥+ 𝑐 23 •𝑦 ≥ 0 ∧ 𝑑 9 + 5 • 𝑑 10 + 𝑑 10 • 𝑥 + 𝑑 11 • 𝑦 ≤ 𝑑 6 + 𝑑 7 • 𝑥 + 𝑑 8 • 𝑦 -𝜖)) ⇒ (-1 ≥ 0)
Note that the transition 𝜃 5 is unconditional, as such we can assume that its condition is simply 1 ≥ 0. Similarly, because there is no target state at location 3, we assume 𝜏 3 ≡ (-1 ≥ 0). Moreover, the transition 𝜃 5 updates the value of 𝑥 to 𝑥 + 5. This is taken into account when generating the constraint above.

The algorithm writes the LHS of the constraint in DNF and handles it exactly as in Example 7.

D Proofs of Theorems Presented in Section 3.1

Notation. Given a set 𝑋 ⊆ R V , we write 𝑋 to denote the closure of 𝑋, i.e. the smallest closed subset of R V that contains 𝑋 . Similarly, for Φ defined as below, we use the notation Φ to denote the system of linear inequalities obtained by replacing every 𝑖 in Φ with ≥ .

Corollary 1. Consider a set V = {𝑣 1 , . . . , 𝑣 𝑟 } of real-valued variables and the following system of 𝑚 linear inequalities over V:

Φ :          𝑎 1,0 + 𝑎 1,1 • 𝑣 1 + . . . + 𝑎 1,𝑟 • 𝑣 𝑟 1 0 . . . 𝑎 𝑚,0 + 𝑎 𝑚,1 • 𝑣 1 + . . . + 𝑎 𝑚,𝑟
• 𝑣 𝑟 𝑚 0 in which 𝑖 ∈ {>, ≥}. When Φ is satisfiable, it entails a given non-strict linear inequality 𝜓 : 𝑐 0 + 𝑐 1 𝑣 1 + . . . + 𝑐 𝑟 𝑣 𝑟 ≥ 0 iff 𝜓 can be written as a non-negative linear combination of 1 ≥ 0 and the inequalities in Φ, i.e. iff there exist non-negative real numbers 𝑦 0 , 𝑦 1 , . . . , 𝑦 𝑚 , such that: As was the case with Farkas' Lemma, it is useful to have a variant of Handelman's theorem that can handle strict inequalities in Φ. We present such a variant, which is a direct corollary of Theorem 12 and characterizes strongly positive polynomials over bounded, but not necessarily closed, polyhedra: Corollary 3. Consider a set V = {𝑣 1 , . . . , 𝑣 𝑟 } of real-valued variables and the following system of 𝑚 linear inequalities over V:

𝑐 0 = 𝑦 0 + 𝑚 𝑖=1 𝑦 𝑖 • 𝑎 𝑖,0 , 𝑐 1 = 𝑚 𝑖=1 𝑦 𝑖 • 𝑎 𝑖,
Φ :          𝑎 1,0 + 𝑎 1,1 • 𝑣 1 + . . . + 𝑎
Φ :          𝑎 1,0 + 𝑎 1,1 • 𝑣 1 + . . . + 𝑎 1,𝑟
𝑔 = 𝑦 0 + 𝑢 𝑖=1 𝑦 𝑖 • ℎ 𝑖 . ( 12
)
Proof. It is obvious that every 𝑔 in the form of ( 12) is strongly positive over Sat(Φ), given that Φ trivially entails 𝑔 ≥ 𝑦 0 > 0. We now prove the other side. Suppose Φ entails 𝑔 ≫ 0.

Let 𝛿 > 0 be the positivity gap of 𝑔 over Sat(Φ) and choose 𝛿 ′ , 𝑦 0 such that 0 < 𝑦 0 < 𝛿 ′ < 𝛿. So, Sat(Φ) ⊆ Sat(𝑔 > 𝛿 ′ ) and hence Sat(Φ) = Sat(Φ) ⊆ Sat(𝑔 > 𝛿 ′ ) = Sat(𝑔 ≥ 𝛿 ′ ). Therefore, Φ entails 𝑔 -𝛿 ′ ≥ 0. So, it also entails 𝑔 -𝑦 0 > 0.

Applying Theorem 12 to Φ and 𝑔 -𝑦 0 , we have:

𝑔 -𝑦 0 = 𝑢 𝑖=1 𝑦 𝑖 • ℎ 𝑖
which is equivalent to Equation [START_REF] Baier | Principles of model checking[END_REF]. □

The Synthesis Algorithm. Our synthesis algorithm is similar to the one in Section 3.1 and consists of five steps. The main difference is in Step 3, in which constraint pairs are translated to non-linear constraints over template variables.

In the previous section, our main tool for this translation was Farkas' Lemma. In this section, due to the more complicated nature of our target sets, we now supplement Farkas' Lemma with Handelman's theorem. For brevity, we do not repeat the presentation of other steps, which are the same as our previous algorithm. Recall that Step 2 (either Steps 2a and 2c for IRWs, or Steps 2b and 2c for UIRWs) has already generated a set Γ of constraint pairs. Each constraint pair 𝛾 ∈ Γ is of the form 𝛾 = (𝜆, 𝜚 ) and encodes the requirement that every inequality in 𝜚 should be entailed by 𝜆. Moreover, 𝜆 is a set of strict/nonstrict linear inequalities over V, whereas 𝜚 is a set of strong polynomial inequalities of degree at most 𝑑. Let 𝑔 ≫ 0 be a strong inequality in 𝜚 . Either 𝜆 is satisfiable and 𝑔 should be represented in the form of Equation 12(cf. Corollary 3) or 𝜆 is unsatisfiable, in which case -1 ≥ 0 or 0 > 0 can be derived as non-negative combinations of inequalities in 𝜆 and 1 ≥ 0 (cf. Corollary 1).

Step 3. Applying Handelman's Theorem and Farkas' Lemma. For every 𝛾 = (𝜆, 𝜚 ) ∈ Γ and strong polynomial inequality 𝑔 ≫ 0 in 𝜚, the algorithm performs the following operations:

• Let Monoid 𝑑 (𝜆) = {ℎ 1 , ℎ 2 , . . . , ℎ 𝑢 } be the set of all polynomials in Monoid(𝜆) whose degree is at most 𝑑. The algorithm symbolically computes Monoid 𝑑 (𝜆) and all of it elements. • The algorithm considers the following three cases, writes constraints that model each of them, and then combines them disjunctively: (i) Writing 𝑔 as in Equation 12. The algorithm creates 𝑢 + 1 new variables 𝑦 0 , 𝑦 1 , . . . , 𝑦 𝑢 with the constraints 𝑦 0 > 0 and 𝑦 1 , . . . , 𝑦 𝑢 ≥ 0, and symbolically computes the equation

𝑔 = 𝑦 0 + 𝑢 𝑖=1 𝑦 𝑖 • ℎ 𝑖 .
Note that both sides of this equation are polynomials of degree 𝑑 over V. Hence, they are equal iff they agree on the coefficient of every monomial. The algorithm equates the coefficients of corresponding monomials in the LHS and RHS of the equation above, hence obtaining a set of equalities over template variables. (ii) Obtaining -1 ≥ 0 as a non-negative combination of 𝜆 and 1 ≥ 0. (iii) Obtaining 0 > 0 as a non-negative combination of 𝜆 and 1 ≥ 0. Cases (ii) and (iii) are handled using Farkas' Lemma in the exact same manner as in our previous algorithm (Section 3.1).

• The algorithm adds the resulting constraints to the nonlinear constraint system Example 15. Consider our running example (Figure 2) together with the templates generated in Example 6. Assume that we aim to synthesize an IRW for 𝜏 3 := 𝑥 2 -𝑥 -100 ≫ 0 , and no target sets in other locations. When Step 2 of the algorithm is applied to location 3 (in exactly the same manner as in Section 3.1) it creates several constraint pairs, including the following:

𝜆 :        𝑐 12 + 𝑐 13 • 𝑥 + 𝑐 14 • 𝑦 ≥ 0 𝑐 15 + 𝑐 16 • 𝑥 + 𝑐 17 • 𝑦 ≥ 0 -𝑐 3 -5 • 𝑐 4 -𝑐 4 • 𝑥 -𝑐 5 • 𝑦 > 0 𝜚 : (𝑥 2 -𝑥 -100 ≫ 0)
In Step 3 of the algorithm, the constraint pair 𝛾 = (𝜆, 𝜚 ) is handled as follows:

• The algorithm computes Monoid 2 (𝜆) which consists of all products of polynomials in 𝜆 up to degree 2. Explicitly, it computes an expanded version of the following polynomials:

ℎ 1 := 1 ℎ 2 := 𝑐 12 + 𝑐 13 • 𝑥 + 𝑐 14 • 𝑦 ℎ 3 := 𝑐 15 + 𝑐 16 • 𝑥 + 𝑐 17 • 𝑦 ℎ 4 := -𝑐 3 -5 • 𝑐 4 -𝑐 4 • 𝑥 -𝑐 5 • 𝑦 ℎ 5 := ℎ 2 2 ℎ 6 := ℎ 2 • ℎ 3 ℎ 7 := ℎ 2 • ℎ 4 ℎ 8 := ℎ 2 3 ℎ 9 := ℎ 3 • ℎ 4 ℎ 10 := ℎ 2 4
• The algorithm considers cases (i)-(iii) as above. Cases (ii) and (iii) are similar to Section 3.1, so we focus on (i). The algorithm introduces 11 new variables 𝑦 0 , . . . , 𝑦 10 , adds the constraints 𝑦 0 > 0 and 𝑦 1 . . . 𝑦 10 ≥ 0 and symbolically computes the following equality:

𝑥 2 -𝑥 -100 = 𝑦 0 + 10 𝑖=1 𝑦 𝑖 • ℎ 𝑖
As before, this is a polynomial equality in R[𝑥, 𝑦], and must hold for all values of 𝑥, 𝑦. So, the corresponding coefficients of the two sides should be equal. The algorithm generates these equalities. For example, given that the constant factor must be the same in the LHS and RHS, the algorithm generates this equality:

-100 = 𝑦 0 + 𝑦 1 + 𝑦 2 • 𝑐 12 + 𝑦 3 • 𝑐 15 -𝑦 4 • 𝑐 3 -5 • 𝑦 4 • 𝑐 4 + 𝑦 5 • 𝑐 12 2 + 𝑦 6 • 𝑐 12 • 𝑐 15 -𝑦 7 • 𝑐 12 • 𝑐 3 -5 • 𝑦 7 • 𝑐 12 • 𝑐 4 + 𝑦 8 • 𝑐 15 2 -𝑦 9 • 𝑐 15 • 𝑐 3 -5 • 𝑦 9 • 𝑐 15 • 𝑐 4 + 𝑦 10 • 𝑐 3 2 + 10 • 𝑦 10 • 𝑐 3 • 𝑐 4 + 25 • 𝑦 10 • 𝑐 4 2 .
The algorithm generates similar equalities for the coefficients of 𝑥, 𝑦, 𝑥 2 , 𝑥 • 𝑦, and 𝑦 2 .

Note that Steps 4 and 5 are also exactly the same as in our previous algorithm and are omitted here. This being said, we have the following theorems, whose proofs are similar to Section 3.1: Theorem 13 (Soundness). Given an input 𝑘-linear system 𝑆 = (V, L, ℓ 0 , 𝐼, Θ), and a set 𝜏 ℓ of at most 𝑘 polynomial inequalities of degree 𝑑 or less at every ℓ ∈ L, every solution of the non-linear constraint system solved in Step 5 of the algorithm above produces a valid 𝑘-linear IRW/UIRW for a target set T that satisfies 𝜏 ℓ at every ℓ ∈ L. Theorem 14 (Completeness). Given a 𝑘-linear system 𝑆 = (V, L, ℓ 0 , 𝐼, Θ), and a set 𝜏 ℓ of at most 𝑘 strong polynomial inequalities of degree 𝑑 or less at every ℓ ∈ L, every bounded 𝑘-linear IRW/UIRW for a target set T that satisfies 𝜏 ℓ at every ℓ ∈ L, is produced by some solution of the non-linear constraint system solved in Step 5 of the algorithm above. Theorem 15 (Complexity). For fixed constants 𝑘, 𝑑 and 𝛽, given a 𝑘-linear 𝛽-branching system 𝑆 = (V, L, ℓ 0 , 𝐼, Θ), and a set 𝜏 ℓ of at most 𝑘 polynomial inequalities of degree 𝑑 or less at every ℓ ∈ L, Steps 1-4 of the algorithm above lead to a polynomial-time reduction from the problem of generating a 𝑘-linear IRW/UIRW to solving a QP instance. Remark 4. Unlike the linear case, our completeness result in Theorem 14 requires strong inequalities and boundedness. This is because Handelman's theorem is only applicable when Sat(Φ) is compact, and hence Corollary 3 can only handle strong inequalities over bounded polyhedra. These requirements do not apply to our soundness result, and although they are theoretically necessary, they have very little impact in practice. If there is an IRW/UIRW for a target set T that ensures reachability within 𝑛 steps, it is easy to verify that there is also a bounded IRW/UIRW with the same property, i.e. the semi-runs starting at ν 0 and taking 𝑛 transitions cannot visit an unbounded set of valuations. Moreover, if the target set contains a non-strong inequality such as 𝑔 ≥ 0 or 𝑔 > 0, one can replace this inequality with 𝑔 + 𝜖 ≫ 0 for a new variable 𝜖 ≥ 0 and solve a quadratic programming instance with the goal of minimizing 𝜖. This trick will slightly change the problem, but it rarely has practical significance.

F Enforcing a Polynomial to be a Sum of Squares

In several places in our algorithm, we have a sum-of-squares polynomial ℎ defined by a template

ℎ := 𝜂 1 • 𝔪 1 + . . . + 𝜂 𝔫 • 𝔪 𝔫
in which {𝔪 1 , . . . , 𝔪 𝔫 } are monomials over a set V of variables, the 𝜂 𝑖 's are unknown reals, and the algorithm depends on ensuring that ℎ is indeed a sum-of-squares polynomial. Given that our algorithm reduces the problem of generating an IRW/UIRW to quadratic programming, we would like to similarly reduce the problem of ℎ being a sum-of-squares to quadratic programming over the 𝜂 𝑖 's. In this section, we show how such a reduction can be obtained. This is a standard procedure and has previously been used in many other constraint-based program analysis algorithms. The presentation we use is taken from [START_REF] Chatterjee | Polynomial Invariant Generation for Non-deterministic Recursive Programs[END_REF]. Our main tools are two well-known theorems: Theorem 16 ([72, Chapter 7]). A polynomial ℎ ∈ R[V] of even degree 𝔡 is a sum-of-squares iff there exists an 𝔯-dimensional symmetric positive semi-definite matrix P such that ℎ = 𝑦 𝑇 P𝑦, where 𝔯 is the number of monomials of degree no greater than 𝔡/2 and 𝑦 is a column vector consisting of every such monomial.

Theorem 17 ([69]

). A symmetric square matrix P is positive semi-definite iff it has a Cholesky decomposition of the form P = LL 𝑇 where L is a lower-triangular matrix with nonnegative diagonal entries.

Given the two theorems above, we use the following standard process for generating quadratic equations that are equivalent to ℎ being a sum-of-squares: Generating Sum-of-Squares Constraints. The algorithm generates the set 𝑀 𝔡 2 consisting of all monomials of degree at most 𝔡 2 over V and creates a vector 𝑦 of these monomials. It then symbolically computes the following equality:

ℎ = 𝑦 𝑇 L L 𝑇 𝑦.
Here, L is a lower-triangular matrix. Every entry of L is a new unknown variable, and every diagonal entry is constrained to be non-negative. As usual, the algorithm equates the corresponding terms on both sides of this polynomial equality to obtain quadratic equations over the unknown variables. It follows directly from the two theorems above that this reduction is both sound and complete. Example 16 (Taken from [START_REF] Chatterjee | Polynomial Invariant Generation for Non-deterministic Recursive Programs[END_REF]). Let V = {𝑎, 𝑏} be the set of variables and ℎ ∈ R[V] a quadratic polynomial, i.e. ℎ(𝑎, 𝑏)

= 𝜂 1 + 𝜂 2 • 𝑎 + 𝜂 3 • 𝑏 + 𝜂 4 • 𝑎 2 + 𝜂 5 • 𝑎 • 𝑏 + 𝜂 6 • 𝑏 2 .
We aim to encode the property that ℎ is a sum-of-squares as a system of quadratic equalities and inequalities. To do so, we first generate all monomials of degree at most ⌊𝔡/2⌋ = 1, which are 1, 𝑎 and 𝑏. Hence, we let 𝑦 = 1 𝑎 𝑏 𝑇 . We then generate a template for a lower-triangular matrix L whose every non-zero entry is a new variable:

L =        𝑙 1 0 0 𝑙 2 𝑙 3 0 𝑙 4 𝑙 5 𝑙 6        .
We also add the inequalities 𝑙 1 ≥ 0, 𝑙 3 ≥ 0 and 𝑙 6 ≥ 0 to our system. Now, we write the equation ℎ = 𝑦 𝑇 L L 𝑇 𝑦 and compute it symbolically:

ℎ = 1 𝑎 𝑏        𝑙 1 0 0 𝑙 2 𝑙 3 0 𝑙 4 𝑙 5 𝑙 6               𝑙 1 𝑙 2 𝑙 4 0 𝑙 3 𝑙 5 0 0 𝑙 6              1 𝑎 𝑏      
, which leads to:

𝜂 1 + 𝜂 2 • 𝑎 + 𝜂 3 • 𝑏 + 𝜂 4 • 𝑎 2 + 𝜂 5 • 𝑎 • 𝑏 + 𝜂 6 • 𝑏 2 = 𝑙 1 2 + 2 • 𝑙 1 • 𝑙 2 • 𝑎 + 2 • 𝑙 1 • 𝑙 4 • 𝑏 + ( 𝑙 2 2 + 𝑙 3 2 ) • 𝑎 2 +(2 • 𝑙 2 • 𝑙 4 + 2 • 𝑙 3 • 𝑙 5 ) • 𝑎 • 𝑏 + ( 𝑙 4 2 + 𝑙 5 2 + 𝑙 6 2 ) • 𝑏 2 .
Note that both sides of the equation above are polynomials over {𝑎, 𝑏}, hence they are equal iff their corresponding coefficients are equal. So, we get the following quadratic equalities over the 𝜂 𝑖 's and 𝑙 𝑖 's:

𝜂 1 = 𝑙 1 2 , 𝜂 2 = 2 • 𝑙 1 • 𝑙 2 , . . . , 𝜂 6 = 𝑙 4 2 + 𝑙 5 2 + 𝑙 6 2
. This concludes the construction of our quadratic system.

G Polynomial Programs used in the Experimental Results

In this section, we provide details of the polynomial programs that were used in our experiments in Section 4. Figures 4-9 show the programs. We now discuss each benchmark in more detail:

• sqrt1: This program gets a value 𝑛 as input and computes ⌊ √ 𝑛⌋ by simply iterating through every integer starting from 1. The goal is to reach the end of the program with 𝑛 -𝑠 > 10 5 . Therefore, to solve this task, a verifier has to assign a value to 𝑛 such that 𝑛 -⌊ √ 𝑛⌋ > 10 5 . It is easy to see that any 𝑛 > 10 5 + 316 works. However, this example is interesting because the shortest path to a target state needs to go through 316 iterations of the while loop. Moreover, the loop has a quadratic guard. As such, a verifier that is based on abstract interpretation needs to obtain a relatively fine abstraction, whereas approaches based on loop-unrolling and symbolic execution need to unroll this non-linear loop 316 times. As mentioned in Table 2, CPAchecker times out on this instance. However, VeriAbs succeeds in proving reachability in 207.3 seconds. In contrast, our approach synthesizes an IRW in just 40.8 seconds.

𝐼 : 𝑛 ≥ 1 ∧ 𝑠 = 1 while (𝑠 + 1) 2 ≤ 𝑛 : • sqrt2: This is a variant of sqrt1 in which the value of 𝑠 is doubled in every step if 2 • 𝑠 ≤ √ 𝑛. This simple change means that there are now many short paths that reach the target. For example, by setting 𝑛 = 2 18 , one can reach the target in just 9 iterations. Unsurprisingly, both CPAchecker and VeriAbs can handle this example (Table 2). This being said, note that the complexity of our approach does not depend on the length of paths. As such, while sqrt1 is much harder than sqrt2 for other approaches, our runtimes on these two benchmarks are very close. Indeed, our method solves sqrt1 a bit faster than sqrt2 (40.8𝑠 vs 46.4𝑠). This is because sqrt1 is a smaller program. • sum: This program simply sums up all the integers from 1 to a given value 𝑛. Note that the program itself is linear (the loop guard and the updates are linear). However, we need polynomial arguments given that for every integer 𝑛, at the end of this program we will have 𝑠 = 𝑛 •(𝑛+1)

𝑠 := 𝑠 + 1

2

. As in the previous examples, the target set can only be reached after many iterations. To reach the target, it suffices to choose 10000 ≤ 𝑛 ≤ 11000. Our algorithm is exact and can handle tight inequalities. We chose this liberal interval in order to make the instance solvable for abstract interpretation approaches with good precision. Nevertheless, CPAchecker timed out and VeriAbs terminated with no result, i.e. returned "unknown".

𝐼 : 𝑠 = 0 ∧ 𝑖 = 1 while 𝑖 ≤ 𝑛 :

(𝑠, 𝑖) := (𝑠 + 𝑖, 𝑖 + 1)

Figure 6. The program sum. Our target is to reach the end of this program with 50005000 ≤ 𝑠 ≤ 60505500.

• sum2: This program is similar to sum but it adds the squares of integers from 1 to 𝑛. Because this program has nonlinear assignments, it is intuitively harder to verify in comparison with sum.

𝐼 : 𝑠 = 0 ∧ 𝑖 = 1 while 𝑖 ≤ 𝑛 : (𝑠, 𝑖) := (𝑠 + 𝑖 2 , 𝑖 + 1)

Figure 7. The program sum2. Our target is to reach the end of this program with 333383335000 ≤ 𝑠 ≤ 443727168500.

• robot1: This program models the behavior of two robots on a 2d plane. One robot is located at (𝑥 1 , 𝑦 1 ) and the other at (𝑥 2 , 𝑦 2 ). Initially, we have (𝑥 1 , 𝑦 1 ) = (𝑥 2 , 𝑦 2 ). At each iteration, each robot moves one unit upwards or to the right. The direction is chosen non-deterministically. The goal is to prove reachability to the endpoint of the program. This is equivalent to proving that it is possible for the robots to move in such a way that makes their distance from each other more than √ 10 5 . The main difficulty in this program is the combinatorial explosion in the number of paths. Nevertheless, note that a relatively large proportion of the paths lead to the desired target. As such, it was surprising for us to see that our approach was the only one that succeeded in handling this example.

𝐼 : 𝑥 1 = 𝑥 2 ∧ 𝑦 1 = 𝑦 2 while (𝑥 1 -𝑥 2 ) 2 + (𝑦 1 -𝑦 2 ) 2 ≤ 10 5 : • robot2: This is a variant of robot1 which intuitively seems to be a bit harder. The same two robots are put on opposite sides of a square with side-length 10 4 and the goal is to prove that they can move in a way that decreases their distance to less than 10. This example creates the same combinatorial explosion in the number of possible paths as in robot1, but this time, only a very small proportion of these paths reach the target. Nevertheless, our approach can handle this example in virtually the same amount of time as robot1. This is because our approach is (semi-)complete and finds a polynomial IRW if one exists. It does not depend on the proportion of paths that lead to a target state.

𝑥 1 := 𝑥 1 + 1 □ 𝑦 1 := 𝑦 1 + 1 𝑥 2 := 𝑥 2 + 1 □ 𝑦 2 := 𝑦 2 + 1

H Related Works

Below, we compare our approach with several families of previous results. CEGAR-based Model-Checkers. Counterexample-guided abstraction refinement [START_REF] Parosh | Counter-Example Guided Program Verification[END_REF][START_REF] Alur | Timing Verification by Successive Approximation[END_REF][START_REF] Balarin | An Iterative Approach to Language Containment[END_REF][START_REF] Clarke | Counterexample-Guided Abstraction Refinement[END_REF] is one of the most successful ideas in software verification and has been implemented in many model-checkers, including the well-known tools SLAM [START_REF] Ball | The SLAM project: debugging system software via static analysis[END_REF] and BLAST [START_REF] Beyer | The software model checker Blast[END_REF], which handle not only safety properties, but also problems such as test-case generation [START_REF] Beyer | Generating Tests from Counterexamples[END_REF]. These model-checkers repeatedly run reachability analyses in order to obtain the required counterexamples for refining their abstractions. A significant challenge is that when variable domains are infinite or uncountable, e.g. R, they cannot guarantee both termination and completeness. They either provide a complete approach that might not terminate, or a sound terminating approach with no completeness guarantee. Another significant challenge arises when there are many spurious counterexamples. Approaches for mitigating this problem also rely on reachability (e.g. see [START_REF] Berdine | Resourceful Reachability as HORN-LA[END_REF]). Invariant Generation. Invariant generation aims to automatically generate over-approximations of reachable sets, while our approach targets reachability analysis that aims to check whether certain undesirable states can be reached (whether a bug is present in the program). Note that although we have the same inductiveness idea as in inductive invariant generation, the idea leads to under-approximation (rather than over-approximation) of the set of states that can reach some target state. Invariant generation is very well-studied and many approaches are developed for automating it, including recurrence analysis [START_REF] Farzan | Compositional Recurrence Analysis[END_REF][START_REF] Kincaid | Compositional recurrence analysis revisited[END_REF][START_REF] Kincaid | Non-linear reasoning for invariant synthesis[END_REF], abstract interpretation [START_REF] Bagnara | Generation of Basic Semi-algebraic Invariants Using Convex Polyhedra[END_REF][START_REF] Chakarov | Expectation Invariants for Probabilistic Program Loops as Fixed Points[END_REF][START_REF] Feautrier | Accelerated Invariant Generation for C Programs with Aspic and C2fsm[END_REF][START_REF] Martin | Invariant Generation through Strategy Iteration in Succinctly Represented Control Flow Graphs[END_REF], constraint-solving [START_REF] Chatterjee | Polynomial Invariant Generation for Non-deterministic Recursive Programs[END_REF][START_REF] Michael A Colón | Linear invariant generation using non-linear constraint solving[END_REF], inference [START_REF] Sharma | From Invariant Checking to Invariant Inference Using Randomized Search[END_REF] and symbolic execution [START_REF] Csallner | DySy: dynamic symbolic execution for invariant inference[END_REF]. Symbolic Execution. Symbolic execution [START_REF] Burch | Symbolic Model Checking: 10 20 States and Beyond[END_REF][START_REF] Cadar | Symbolic execution for software testing: three decades later[END_REF] runs program codes in a static and symbolic manner, thus it is very effective for analyzing programs without loops or with bounded loops. When tackling loops, symbolic execution can only unfold the loop a bounded number of steps, hence it is unsuitable for loops with an unbounded number of iterations.

In contrast, our approach can handle loops with unbounded iterations, and provides soundness and completeness. Abstract Interpretation. Abstract interpretation mainly focuses on over-approximation of reachable states through the widening operator, which often leads to a loss of precision [START_REF] Gonnord | Abstract acceleration in linear relation analysis[END_REF]. There are also several abstraction-based results on under-approximation [START_REF] Giacobazzi | Making abstract interpretations complete[END_REF][START_REF] Ranzato | Complete Abstractions Everywhere[END_REF][START_REF] Rival | Understanding the Origin of Alarms in Astrée[END_REF][START_REF] Schmidt | A calculus of logical relations for over-and underapproximating static analyses[END_REF]. However, a theory with completeness guarantees for generating underapproximations such as our T-inductive sets through abstract interpretation is still lacking. Termination Analysis. Termination analysis only considers whether a program terminates in a finite number of steps, i.e. whether the program can reach the terminal program counter or not, without constraints over program variables.

In our approach, we consider reachability to program states defined by numerical constraints over program variables, which is a considerable extension of the termination property. The primary method of proving termination is to synthesize a ranking function, and there are template-based algorithms for the synthesis of linear/polynomial ranking functions [START_REF] Chatterjee | Termination analysis of probabilistic programs through Positivstellensatz's[END_REF][START_REF] Michael | Synthesis of linear ranking functions[END_REF][START_REF] Leike | Ranking Templates for Linear Loops[END_REF][START_REF] Podelski | A Complete Method for the Synthesis of Linear Ranking Functions[END_REF]. Termination and reachability properties have also been extensively studied in the context of probabilistic programs, especially through martingale-based approaches (e.g. see [START_REF] Chakarov | Probabilistic Program Analysis with Martingales[END_REF][START_REF] Huang | Modular verification for almost-sure termination of probabilistic programs[END_REF][START_REF] Takisaka | Ranking and Repulsing Supermartingales for Reachability in Probabilistic Programs[END_REF]). Incorrectness Logic. Incorrectness logic [START_REF] O'hearn | Incorrectness logic[END_REF] and reverse Hoare logic [START_REF] De | Reverse Hoare Logic[END_REF] are similar to Hoare logic but target underapproximations. The logical background behind our approach is a bit different from incorrectness logic. Incorrectness logic obtains under-approximations of the set of reachable states. Hence, a bug can be found by taking the intersection of the under-approximation obtained by incorrectness logic and the set T of undesirable states. In contrast, we find underapproximations of the sets of states from which reachability to an undesirable state (or bug) in T is guaranteed. Intuitively, the relationship between our approach and incorrectness logic is similar to the relationship between inductive invariants and Hoare logic. It is also noteworthy that incorrectness logic needs manual effort to write assertions, while our approach is entirely automated when we consider linear/polynomial IRWs/UIRWs.

Figure 2 .

 2 Figure 2. Our Running Example

  f o r ( 𝑎 := 0 ; 𝑎 < 𝑀 -1 ; 𝑎 := 𝑎 + 1 ) : f o r (𝑏 := 0 ; 𝑏 < 𝑀 -1 ; 𝑏 := 𝑏 + 1 ) : f o r ( 𝑐 := 0 ; 𝑐 < 𝑀 -1 ; 𝑐 := 𝑐 + 1 ) : f o r ( 𝑑 := 0 ; 𝑑 < 𝑀 -1 ; 𝑑 := 𝑑 + 1 ) : f o r ( 𝑒 := 0 ; 𝑒 < 𝑀 -1 ; 𝑒 := 𝑒 + 1 ) : i f 𝑀 -2 ≤ 𝑎, 𝑏, 𝑐, 𝑑, 𝑒 : print ( "target reached" )

Figure 3 .

 3 Figure 3. A Non-deterministic Program and its Representation as a Transition System

Figure 4 .

 4 Figure 4. The program sqrt1. Our target is to reach the end of this program with 𝑛 -𝑠 > 10 5 .

𝐼 : 𝑛 ≥ 1 ∧

 1 𝑠 = 1 while (𝑠 + 1) 2 ≤ 𝑛 : i f 4 • 𝑠 2 ≤ 𝑛 : 𝑠 := 2 • 𝑠 e l s e : 𝑠 := 𝑠 + 1

Figure 5 .

 5 Figure 5. The program sqrt2. Our target is to reach the end of this program with 𝑛 -𝑠 > 10 5 .

Figure 8 .

 8 Figure 8. The program robot1. Our target is to reach the end of this program.

𝐼 : 𝑦 1 = 1 Figure 9 .

 119 Figure 9. The program robot2. Our target is to reach the end of this program.

  1 , . . . , 𝑣 𝑟 ) ≥ 0, . . . , 𝑔 𝑚 (𝑣 1 , . . . , 𝑣 𝑟 ) ≥ 0 where 𝑔 1 , . . . , 𝑔 𝑚 ∈ R[V] are polynomials. If there exists a 𝑔 𝑖 such that the set Sat(𝑔 𝑖 ≥ 0) is compact, and Φ entails a given polynomial inequality 𝑔(𝑣 1 , . . . , 𝑣 𝑟 ) > 0 then there exist polynomials ℎ 𝑚 𝑖=1 ℎ 𝑖 • 𝑔 𝑖 and every ℎ 𝑖 is a sum of squares, i.e. ℎ 𝑖 = ℎ 2 𝑖,𝑗 for some polynomials ℎ 𝑖,𝑗

0 , ℎ 1 , . . . , ℎ 𝑚 ∈ R[V] such that 𝑔 = ℎ 0 +

  1 , . . . , 𝑣 𝑟 , 𝑤 1 , . . . , 𝑤 𝑚 ) := 𝑔 𝑖 (𝑣 1 , . . . , 𝑣 𝑟 ) -𝑤 2 𝑖 . So, we have 𝑤 2•𝛼

			𝑗	= 𝑚 𝑖=1 ℎ 𝑖 • g𝑖 . More-
	𝛼
	over, 𝑔 𝛼 𝑗 = g𝑗 + 𝑤 2 𝑗	= 𝛼 𝑖=0	𝛼 𝑖

  𝜆 is a set of polynomial inequalities of the form {𝑔 𝑖 𝑖 0} 𝑚 𝑖=1 . Let 𝑔 ≫ 0 be a strong polynomial inequality in 𝜚 . 𝜆 must entail 𝑔 ≫ 0. The algorithm considers three cases:(i) 𝜆 is unsatisfiable due to case (i) in Theorem 8: The algorithm considers the set 𝑀 Υ 1 (V) := {𝔪 1 , 𝔪 2 , . . . , 𝔪 𝔫 } of all monomials of degree at most Υ 1 over V. Recall that Υ 1 is the first technical parameter given in input. It then generates the following templates ℎ 𝑖 for 0 ≤ 𝑖 ≤ 𝑚:ℎ 𝑖 := 𝜂 𝑖,1 • 𝔪 1 + . . . + 𝜂 𝑖,𝔫 • 𝔪 𝔫by introducing new variables 𝜂 𝑖,𝑗 . It also adds certain constraints on 𝜂 𝑖,𝑗 's that ensure every ℎ 𝑖 is a sum-ofsquares. See Appendix F for more details. Then, the algorithm introduces 𝑦 0 > 0 and symbolically computes-1 = 𝑦 0 + ℎ 0 + 𝑚 𝑖=1 ℎ 𝑖 • 𝑔 𝑖 Finally,the algorithm equates the corresponding coefficients on the two sides of the equality above, and obtains quadratic equalities over the unknown variables. No program variable appears in these equalities. (ii) 𝜆 is unsatisfiable due to case (ii) in Theorem 8: The algorithm considers the set 𝑀 * Υ 2 := {𝔪 * 1 , . . . , 𝔪 * 𝔫 * } of all monomials of degree at most Υ 2 (our second technical parameter) over the extended variable set V * = V ∪ {𝑤 1 , . . . , 𝑤 𝑚 }. It generates templates ℎ 𝑖 for 1 ≤ 𝑖 ≤ 𝑚 : ℎ 𝑖 := 𝜂 𝑖,1 • 𝔪 * 1 + . . . + 𝜂 𝑖,𝔫 • 𝔪 * 𝑖 • (𝑔 𝑖 -𝑤 2 𝑖 ). Here Υ 3 is our third technical parameter and both sides are polynomials in R[V * ]. The algorithm equates the corresponding coefficients on the LHS and RHS and obtains quadratic equalities over unknown variables. The systems of quadratic equalities generated for each index 𝑗 are combined disjunctively. (iii) 𝑔 is a combination of 𝑔 𝑖 's as in Corollary 2: The algorithm considers the set 𝑀 Υ

		𝔫 *
	and symbolically computes the equality below for every
	index 𝑗 with a strict inequality 𝑔 𝑗 > 0 in 𝜆:
	𝑤 2•Υ 3 𝑗	= 𝑚 𝑖=1 ℎ

[START_REF] Albarghouthi | From Under-Approximations to Over-Approximations and Back[END_REF]

. Applying the Stellensätze. The algorithm applies Corollary 2 and Theorem 8 to every constraint pair generated in the previous step to obtain a non-linear constraint system based on the template variables, 𝜖, and new variables defined in this step. Let 𝛾 = (𝜆, 𝜚 ) ∈ Γ be a constraint pair. 4 := {𝔪 1 , . . . , 𝔪 𝔫 } of monomials of degree at most Υ 4 over V, and generates the following templates ℎ 𝑖 for 0 ≤ 𝑖 ≤ 𝑚:

ℎ 𝑖 := ℎ 𝑖 := 𝜂 𝑖,1 • 𝔪 1 + . . . + 𝜂 𝑖,𝔫 • 𝔪 𝔫

by introducing new variables 𝜂 𝑖,𝑗 and adding constraints that ensure every ℎ 𝑖 is a sum-of-squares (Appendix F). It then introduces 𝑦 0 > 0 and symbolically computes

Table 1 .

 1 Experimental

	Benchmark		|L |	|Θ |	|V | 𝑘	|QP |	Gen Solve	Ours CPAchecker VeriAbs Dangerzone	Depth
	gcnr2008		8	14	4	2	1838	14.8	81.4	96.2	1.8	17.6	3.3	4
	trex02-2		5	7	1	2	260	1.8	3.4	5.3	4.3	16.6	2.4	3
	nec11		4	8	3	2	2871	13.2	45.7	58.9	4.2	10.8	2.7	6
	terminator_02-1		5	8	3	3	1962	12.0	19.4	31.4	4.2	17.2	2.9	5
	simple_2-2		3	5	1	2	267	1.7	1.5	3.2	4.3	10.4	T/O	3
	count_up_down-2		3	4	3	2	244	1.6	1.9	3.5	4.4	5.9	T/O	4
	nested_1b		4	7	1	2	819	3.3	11.0	14.3	4.7	3.5	T/O	9
	while_infinite_loop_4	10	14	1	2	1223	5.1	7.9	13.0	4.1	15.3	N/A	6
	multivar_1-2		3	6	2	2	900	5.8	19.8	25.6	4.4	9.0	T/O	4
	trex01-1		14	27	6	3	9491	69.7	228.2	297.8	4.5	17.3	4.5	7
	for_bounded_loop1		10	13	5	2	1579	9.8	30.1	39.9	5.6	16.8	4.6	12
	underapprox_1-1		3	6	2	2	901	3.7	5.6	9.3	4.7	14.0	T/O	21
	sum03-1		9	14	2	2	20963	77.9	413.1	491.0	6.1	16.3	4.8	57
	underapprox_2-1		3	6	2	2	901	3.7	14.5	18.2	4.9	10.6	T/O	21
	sum04-1		6	10	2	2	1082	5.4	8.0	13.4	5.1	17.0	T/O	22
	terminator_03-1		6	11	2	3	1740	10.0	25.7	35.6	5.1	9.7	T/O	5
	trex03-1		4	12	6	2	8500	49.2	197.9	247.0	5.2	9.1	5.3	5
	Mono1_1-1		3	5	1	2	262	1.3	4.4	5.7	T/O	377.2	T/O	5.50e+7
	sum01_bug02_base		7	13	3	2	7972	38.0	133.4	171.4	6.0	17.3	24.1	18
	id_trans		5	11	5	2	11192	68.7	171.2	239.8	6.4	19.8	10.8	11
	sum01_bug02		7	12	4	3	17632	60.0	218.6	278.6	6.5	17.3	71.8	26
	sum01-1		7	12	3	2	7316	36.9	55.1	92.0	7.6	16.7	T/O	32
	nested_1-2		4	6	2	2	329	2.9	8.0	10.9	T/O	86.0	N/A	6.17e+9
	Mono5_1		5	7	3	4	1048	8.1	31.8	39.8	T/O	332.4	19.6	3.50e+7
	const_1-2		3	6	2	2	901	4.6	17.6	22.2	T/O	49.6	T/O	3076
	Mono3_1		6	8	2	3	660	4.0	20.0	24.0	T/O	369.9	T/O	4.00e+6
	Mono4_1		5	7	2	4	949	5.3	22.1	27.3	T/O	635.8	T/O	3.50e+6
	Mono6_1		5	7	3	5	1502	11.6	48.5	60.1	T/O	382.2	T/O	3.50e+7
	deep_nested		7	17	5	5	3686	28.6	69.6	98.2	T/O	T/O	N/A	1.46e+48
	nested_delay_notd2	6	9	6	3	12858	55.3	396.7	451.9	T/O	247.4	N/A	1385
	Benchmark	|L |	|Θ |	|V | 𝑘 𝑑	|QP|	Gen Solve	Ours	CPAchecker VeriAbs Dangerzone	Depth
	sqrt2	5	7	2	5	2	2494	24.1	22.4	46.4	10.5	19.7	2.7	207
	sqrt1	3	4	2	4	2	920	10.7	30.1	40.8	T/O	207.3	2.7	633
	robot2	5	8	4	5	2	5537	71.1	681.7	752.9	T/O	F	3.1	29982
	robot1	5	8	4	5	2	5537	69.8	724.2	794.0	T/O	F	5.97	2127
	sum	3	4	3	5	2	1826	20.4	59.7	80.1	T/O	F	T/O	30003
	sum2	3	4	3	5	3	2476	36.8	167.5	204.2	T/O	T/O	T/O	30003

Results over SV-COMP Benchmarks. Times are in seconds. "Gen"/"Solve" are the times spent generating/solving the QP. "Ours" is our total runtime. "T/O" denotes timeout. "Depth" is the length of the shortest path that reaches the target states (bugs). Since our approach does not output a path, but rather an IRW, we calculated the depths manually.

Table 2 .

 2 Experimental Results over Polynomial Programs. 'F' denotes failure. We set our Υ variables equal to 𝑑.

  1 , . . . , 𝑐 𝑟 = Let 𝑋 ⊆ R V be a set of valuations and 𝑔 ∈ R[V] a polynomial over V. We say that 𝑔 is strongly positive over 𝑋 , and write 𝑋 |= 𝑔 ≫ 0 (or simply 𝑔 ≫ 0 when 𝑋 is clear from context), if inf 𝑥 ∈𝑋 𝑔(𝑥) > 0. The real value 𝛿 := inf 𝑥 ∈𝑋 𝑔(𝑥) is called the positivity gap of 𝑔 over 𝑋 . Moreover, 𝑔 is strongly greater than ℎ, denoted 𝑔 ≫ ℎ, iff 𝑔ℎ ≫ 0. Problem Definition. In this section, we consider the following problem: Given a 𝑘-linear system 𝑆 = (V, L, ℓ 0 , 𝐼, Θ) together with a set 𝜏 ℓ of at most 𝑘 strong polynomial inequalities of degree at most 𝑑 at every location ℓ ∈ L, synthesize a 𝑘-linear IRW/UIRW for a target set T that satisfies 𝜏 ℓ at every ℓ ∈ L, or report that no such IRW/UIRW exists. Mathematical Tool. Our main mathematical tool in this section is a theorem, due to Handelman, that characterizes positive polynomials over compact polyhedra. Before presenting this theorem, it is useful to define the notion of monoid. Monoid. Consider a set V = {𝑣 1 , . . . , 𝑣 𝑟 } of real-valued variables and the following system of 𝑚 linear inequalities over V:+ 𝑎 1,1 • 𝑣 1 + . . . + 𝑎 1,𝑟• 𝑣 𝑟 1 0 . . . 𝑎 𝑚,0 + 𝑎 𝑚,1 • 𝑣 1 + . . . + 𝑎 𝑚,𝑟 • 𝑣 𝑟 𝑚 0 in which 𝑖 ∈ {>, ≥}. Let 𝑔 𝑖 be the LHS of the 𝑖-th inequality, i.e. 𝑔 𝑖 (𝑣 1 , . . . , 𝑣 𝑟 ) := 𝑎 𝑖,0 + 𝑎 𝑖,1 • 𝑣 1 + . . . + 𝑎 𝑖,𝑟 • 𝑣 𝑟 . The monoid of Φ is defined as: Consider a set V = {𝑣 1 , . . . , 𝑣 𝑟 } of realvalued variables and the following system of 𝑚 linear inequalities over V:

	E Synthesizing Linear IRWs/UIRWs for
	Linear Systems with Polynomial Target
	Sets
	For technical reasons, we need the concept of strong positiv-
	ity.
	Strong Positivity. Φ :          𝑎 1,0 Monoid(Φ) :=	𝑚 𝑖=1 𝑔 𝜅 𝑖 𝑖	∀1 ≤ 𝑖 ≤ 𝑚, 𝜅 𝑖 ∈ N ∪ {0} .
	Basically, Monoid(Φ) is the set of all polynomials that can
	be obtained by multiplying the linear expressions on the LHS
	of Φ together. Note that each such expression can appear
	zero or multiple times in the multiplication. Specifically, it is
	noteworthy that 1 ∈ Monoid(Φ).
	Theorem 12 ([67]).
			𝑚
			𝑦 𝑖 • 𝑎 𝑖,𝑟 .
			𝑖=1

  1,𝑟 • 𝑣 𝑟 ≥ 0 . . . 𝑎 𝑚,0 + 𝑎 𝑚,1 • 𝑣 1 + . . . + 𝑎 𝑚,𝑟 • 𝑣 𝑟 ≥ 0 If Φ is satisfiable, Sat(Φ) is compact, and Φ entails a given polynomial inequality 𝑔(𝑣 1 , . . . , 𝑣 𝑟 ) > 0 then there exist reals 𝑦 1 , 𝑦 2 , . . . , 𝑦 𝑢 ∈ [0, ∞) and ℎ 1 , ℎ 2 , . . . , ℎ 𝑢 ∈ Monoid(Φ) such that: 𝑔 = 𝑢 𝑖=1 𝑦 𝑖 • ℎ 𝑖 .

  • 𝑣 𝑟 1 0 . . . 𝑎 𝑚,0 + 𝑎 𝑚,1 • 𝑣 1 + . . . + 𝑎 𝑚,𝑟 • 𝑣 𝑟 𝑚 0 in which 𝑖 ∈ {>, ≥}. If Φ is satisfiable and Sat(Φ) is bounded, then Φ entails a given strong polynomial inequality 𝑔(𝑣 1 , . . . , 𝑣 𝑟 ) ≫ 0, or in other words Sat(Φ) |= 𝑔(𝑣 1 , . . . , 𝑣 𝑟 ) ≫ 0, if and only if there exist constants 𝑦 0 ∈ (0, ∞) and 𝑦 1 , 𝑦 2 , . . . , 𝑦 𝑢 ∈ [0, ∞), and polynomials ℎ 1 , ℎ 2 , . . . , ℎ 𝑢 ∈ Monoid(Φ) such that:

Acknowledgments

This research was partially supported by the ERC CoG 863818 (ForM-SMArt), the National Natural Science Foundation of China (NSFC) Grant No. 61802254, the Huawei Innovation Research Program, the Facebook PhD Fellowship Program, and DOC Fellowship No. 24956 of the Austrian Academy of Sciences (ÖAW).

B Proofs of Theorems Presented in Section 2

In this section we provide proofs of our basic soundness and completeness theorems. Theorem 1 (Soundness). Let T ⊆ Σ be a set of states in the system 𝑆 = (V, L, ℓ 0 , 𝐼, Θ).

(i) If there exists an IRW (T ♦ , 𝑓 , 𝜖) for T, then T is existentially reachable. (ii) If there exists a UIRW (T ♦ , 𝑓 , 𝜖) for T, then T is universally reachable.

Proof. We handle each case separately.

(i) We construct a run of 𝑆 that visits T. By definition of IRW, there exists a state 𝜎 0 = (ℓ 0 , ν 0 ) ∈ T ♦ such that ν 0 |= 𝐼 . We start our run with 𝜎 0 and inductively find the next transitions and states as follows: when we are in a state 𝜎 𝑖 ∈ T ♦ , either (a) 𝜎 𝑖 ∈ T in which case the path until this point has already reached T and we can extend it to an arbitrary run, or (b) 𝜎 𝑖 ∈ T ♦ \ T, in which case there exists a successor 𝜎 𝑖+1 ∈ T ♦ of 𝜎 𝑖 such that 𝑓 (𝜎 𝑖+1 ) ≤ 𝑓 (𝜎 𝑖 ) -𝜖, and we transition to 𝜎 𝑖+1 . Using this procedure, it is not possible to avoid case (a) forever, because each application of (b) decreases the value of 𝑓 by at least 𝜖 and 𝑓 is bounded from below. Hence, the constructed run will reach T. (ii) We choose 𝜎 0 = (ℓ 0 , ν 0 ) as in the previous case. We now prove that every path of length 𝑛 := 1 + ⌈𝑓 (𝜎 0 )/𝜖⌉ starting from 𝜎 0 will reach T. Let r = 𝜎 0 , 𝜃 0 , 𝜎 Theorem 2 (Completeness). Let T ⊆ Σ be a set of states in the system 𝑆 = (V, L, ℓ 0 , 𝐼, Θ).

(i) If T is existentially reachable, then there exists an IRW (T ♦ , 𝑓 , 𝜖) for T. (ii) If T is universally reachable, then there exists a UIRW (T ♦ , 𝑓 , 𝜖) for T.

Proof. In each case, we construct the required IRW/UIRW. Moreover, Φ is unsatisfiable iff either -1 ≥ 0 can be derived as above, or 0 > 0 can be derived as above with the extra requirement that 𝑖 ∈ {> } 𝑦 𝑖 > 0 (i.e. in order to derive a strict inequality, we should use at least one of the strict inequalities in Φ with non-zero coefficient).

Proof. For the first part, suppose that𝜓 is entailed by Φ, hence

The former is a closed set, hence it also includes the closure of the latter, which is the set of points that satisfy Φ. Hence, we can apply Lemma 1 to Φ and 𝜓 to obtain the desired result.

For the second part, if Φ is satisfiable, then obviously no non-negative combination of inequalities in Φ can sum up to a contradiction such as 0 > 0 or -1 ≥ 0. If Φ is not satisfiable, then by Lemma 1, we can obtain -1 ≥ 0. The only remaining case is if Φ is satisfiable but Φ is not. Reorder the inequalities in Φ so that the non-strict inequalities appear first. Then, consider the smallest 𝑖 for which the first 𝑖 inequalities in Φ are unsatisfiable. Let Φ[1 . . . 𝑖] denote the first 𝑖 inequalities. Based on our ordering, we know that the 𝑖-th inequality is strict and of the form

Therefore, by the first part above, we can write

as a non-negative combination of the first 𝑖 -1 inequalities. Moreover, the 𝑖-th inequality is:

Summing up these two, we get 0 > 0. □ Theorem 3 (Soundness). Given an input 𝑘-linear system 𝑆 = (V, L, ℓ 0 , 𝐼, Θ), and a 𝑘-linear set T of target states, every solution of the non-linear constraint system solved in Step 5 of the algorithm in Section 3.1 produces a valid 𝑘-linear IRW/UIRW for T in 𝑆.

Proof. Every solution 𝔰 satisfies the constraints generated in Step 3. Therefore, for every constraint pair 𝛾 = (𝜆, 𝜚 ) ∈ Γ generated in Step 2 and inequality 𝛼 0 + 𝛼 • V ≥ 0 in 𝜚 , either 𝔰(𝜆) is unsatisfiable, i.e. a non-negative linear combination of its inequalities sums up to 0 ≥ 1 or 0 > 0, or there is such a linear combination that sums up to 𝛼 0 +𝛼 •V ≥ 0. In each case, the coefficients of the combination are given by 𝔰( 𝑦 𝑖 ) for the corresponding 𝑦 𝑖 variables. Moreover, no matter which case happens, the inequalities in 𝜚 are entailed by 𝜆. By definition, the constraint pairs generated in Step 2 modeled inductivity, non-negativity and ranking conditions and hence 𝔰 satisfies these properties. Finally, 𝔰 satisfies the constraints generated in Step 4. Therefore, we have 𝔰( ν 0,1 , . . . , ν 0,𝑟 ) |= 𝔰( 𝐴 ℓ 0 ) ∧ 𝐼 . So, all the requirements for IRW/UIRW are met. □ Theorem 4 (Completeness). Given a 𝑘-linear system 𝑆 = (V, L, ℓ 0 , 𝐼, Θ), and a 𝑘-linear set T of target states, every 𝑘-linear IRW/UIRW for T in 𝑆 is produced by some solution to the non-linear constraint system solved in Step 5 of the algorithm in Section 3.1.

Proof. We construct the required solution. Let (T ♦ , 𝑓 , 𝜖) be a 𝑘-linear IRW/UIRW for T in 𝑆. Let 𝐴 ℓ be the set of inequalities defining T ♦ ∩ ℓ × R V , and 𝑓 ℓ the linear expression defining 𝑓 at ℓ. We use the coefficients in 𝐴 ℓ 's and 𝑓 ℓ 's as the corresponding values for 𝔰( 𝑐 ℓ,𝑖,𝑗 )'s and 𝔰( 𝑑 ℓ,𝑗 )'s. Moreover, we let 𝔰( 𝜖) = 𝜖. By definition, T ♦ is an existential/universal T-inductive set, and 𝑓 is an existential/universal T-ranking function with parameter 𝜖. Therefore, 𝐴 ℓ 's and 𝑓 ℓ 's satisfy the constraint pairs generated at Step 2 of the algorithm. By Corollary 1, there are suitable values for each variable 𝑦 𝑖 such that the constraints in Step 3 are satisfied. We use these values as 𝔰( 𝑦 𝑖 ). Finally, by definition of IRW/UIRW, there exists a valuation ν ∈ R V such that ν |= 𝐴 ℓ 0 ∧ 𝐼 = 𝔰( 𝐴 ℓ 0 ) ∧ 𝐼 . We let 𝔰( ν 0,𝑖 ) = ν 𝑖 . It is easy to verify that 𝔰 is a solution to the system of non-linear constraints solved in Step 5. □ Theorem 5 (Complexity). For fixed constants 𝑘 and 𝛽, given a 𝑘-linear 𝛽-branching system 𝑆 = (V, L, ℓ 0 , 𝐼, Θ), and a 𝑘-linear set T of target states, Steps 1-4 of the algorithm in Section 3.1 lead to a polynomial-time reduction from the problem of generating a 𝑘-linear IRW/UIRW to solving a Quadratic Programming (QP) instance.

Proof. It is easy to verify that all steps of the algorithm run in polynomial time ‡ , and that all the generated (in)equalities over non-program variables are quadratic. However, these (in)equalities are not always combined conjunctively. Specifically, in Step 3, the constraints corresponding to cases (i)-(iii) are combined disjunctively. This being said, we can perform the following actions to obtain a QP instance in polynomial time:

• We first convert every inequality of the form 𝔢 0 to 𝔢 -𝑥 𝔢 = 0 by introducing a new variable 𝑥 𝔢 0.

• We rewrite every disjunction 𝔢 1 = 0 ∨ 𝔢 2 = 0 as 𝔢 1 •𝔢 2 = 0.

Note that this might create polynomial equalities of higher degree. • We eliminate terms of degree more than 2 by defining new variables that are equal to their proper divisors, e.g. we rewrite 𝑐 1 • 𝑐 2 • 𝑐 3 2 as 𝜐 1 • 𝜐 2 where 𝜐 1 , 𝜐 2 are new variables, and add the equalities 𝜐 1 = 𝑐 1 • 𝑐 2 and 𝜐 2 = 𝑐 3 2 .

The steps above lead to a polynomial blow-up in the size of the system, given that in Step 3 of the algorithm we have disjunctions of at most 3 different boolean formulas. □ ‡ The reason for fixing 𝑘 and 𝛽 is to avoid exponential blow-up when rewriting boolean expressions in DNF.