# Automated Expected Value Analysis of Recursive Programs

Martin Avanzini, Georg Moser, Michael Schaper

► To cite this version:

Martin Avanzini, Georg Moser, Michael Schaper. Automated Expected Value Analysis of Recursive Programs. Proceedings of the ACM on Programming Languages, 2023, 7 (PLDI), pp.1050-1072. 10.1145/3591263. hal-04345663

HAL Id: hal-04345663

https://inria.hal.science/hal-04345663v1

Submitted on 14 Dec 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à 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.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Automated Expected Value Analysis of Recursive Programs

MARTIN AVANZINI, INRIA Sophia Antipolis Méditerranée, France

GEORG MOSER, Universität Innsbruck, Austria

MICHAEL SCHAPER, Build Informed, Austria

In this work, we study the fully automated inference of expected result values of probabilistic programs in the presence of natural programming constructs such as procedures, local variables and recursion. While crucial, capturing these constructs becomes highly non-trivial. The key contribution is the definition of a term representation, denoted as infer[·], translating a pre-expectation semantics into first-order constraints, susceptible to automation via standard methods. A crucial step is the use of logical variables, inspired by previous work on Hoare logics for recursive programs. Noteworthy, our methodology is not restricted to tail-recursion, which could unarguably be replaced by iteration and wouldn't need additional insights. We have implemented this analysis in our prototype ev-imp. We provide ample experimental evidence of the prototype's algorithmic expressibility.

CCS Concepts: • **Theory of computation** → **Program analysis**; Automated reasoning.

Additional Key Words and Phrases: probabilistic programming, expected value analysis, weakest pre-expectation semantics, automation

**ACM Reference Format:**

Martin Avanzini, Georg Moser, and Michael Schaper. 2023. Automated Expected Value Analysis of Recursive Programs. *Proc. ACM Program. Lang.* 7, PLDI, Article 149 (June 2023), 27 pages. https://doi.org/10.1145/3591263

# 1 INTRODUCTION

The verification of the quantitative behaviour of probabilistic programs is a highly active field of research, motivated partly by the recent success of machine learning methodologies (see eg. [Batz et al. 2019; Eberl et al. 2020; Vasilenko et al. 2022]). Without verification, bugs may hide in correctness arguments and subsequent implementations, in particular, as reasoning about probabilistic programs (or data structures for that matter) is highly non-trivial and error prone.

Instead of verifying quantitative program behaviour semi-automatically, it would be desirable to fully automatically infer such estimates. For example, the goal of inference could be an approximate but still precise computation of *expected (amortised) costs* (eg. [Avanzini et al. 2020; Leutgeb et al. 2022; Wang et al. 2020]) or *quantitative invariants* (eg. [Bao et al. 2022; Wang et al. 2018]).<fnref n="1" /> The goal of inference is to compute approximations that are as precise as possible, which requires the use of optimisation techniques. As there are now powerful optimising constraint solvers available, such as Z3 [de Moura and Bjørner 2008] or OptiMathSAT [Sebastiani and Trentin 2020], we may

Authors' addresses: Martin Avanzini, INRIA Sophia Antipolis Méditerranée, Route des Lucioles - BP 93, France, martin.avanzini@inria.fr; Georg Moser, Universität Innsbruck, Innsbruck, Austria, georg.moser@uibk.ac.at; Michael Schaper, Build Informed, Innsbruck, Austria, mschaper@posteo.net.

Permission to make digital or hard copies of part or all of this work for personal or classroom use is granted without fee provided that copies are not made or distributed for profit or commercial advantage and that copies bear this notice and the full citation on the first page. Copyrights for third-party components of this work must be honored. For all other uses, contact the owner/author(s).

© 2023 Copyright held by the owner/author(s).

2475-1421/2023/6-ART149

https://doi.org/10.1145/3591263

<footnote n="1">¹The need for automated techniques that analyse eg. the computational cost of code has also been recognised in large software companies. For example at Facebook, one routinely runs a cost analysis on the start-up routines in order to ensure a fast loading of the Facebook web page [Distefano et al. 2019].</footnote>
even hope to automatically conduct optimisations that earlier required intricate analysis with pen and paper, cf. [Leutgeb et al. 2021, 2022].

The central contribution of our work is the (fully) automated inference of expected result values of probabilistic programs in the presence of natural programming constructs such as procedures, local variables and recursion. While crucial, capturing these constructs becomes highly non-trivial. To analyse recursive procedures successfully, it is required to properly model the call-stack. In particular, for each (recursive) procedure call the program context may change, that is, the analysis of (general) recursion requires some form of parametricity.<fnref n="2" /> To the best of our knowledge, prior work fails in general to deal with these natural constructs.

Upper invariants (on the expected value) constitute a liberalisation of exact quantitative invariants. This approximated rendering of invariants allows for easier (and thus more powerful) automation. Such an analysis may also act as a stepping stone towards the incorporation of support for automation into ITPs. To wit, a recent and partly motivating work is Vasilenko et al. [2022]. We take up one of their motivating examples—see Listing 1(a) below—and show how a variant of this example becomes susceptible to full automation.

Our analysis is based on a weakest pre-expectation semantics [McIver and Morgan 2005] (see also [Gretz et al. 2014; Kaminski et al. 2018])—an axiomatic semantic in Dijkstra's spirit—for a simple imperative language PWhile endowed with sampling instructions, non-deterministic choice, lexically scoped local variables, and, crucially, recursive procedures. On a conceptual level, this semantics—denoted as et[[P]] for a program P—capture the expected behaviour of probabilistic programs.

Automation of a weakest pre-expectation semantics in the context of recursive procedures is highly non-trivial.<fnref n="3" /> Technically, this is due to the fact that well-definedness of the semantics requires the existence of *higher-order* fixed-points. Full automation, however, requires the inference of (upper) bounds on closed-forms of such fixed-points. We overcome this challenge through the definition of a suitable *term representation*—denoted as infer[P]—of the aforementioned pre-expectation semantics et[[P]]. This syntactic representation translates the pre-expectation semantics into first-order constraints, susceptible to automation via standard methods. To provide for the aforementioned parametricity in this analysis, the use of logical variables—inspired by the work on Hoare logics for recursive programs [Kleymann 1999]—is essential. We make use of a template approach and employ Z3 as suitable optimising SMT solver. Perhaps the closest comparison is to work on amortised cost analysis of functional languages (eg. [Leutgeb et al. 2022; Wang et al. 2020]).

*Contributions.* In sum, we present a novel methodology for the *automated expected value* analysis of non-deterministic, probabilistic and recursive programs. Our starting point is a natural weakest pre-expectation semantics for PWhile in the form of an *expectation transformer* et[[P]], capturing also natural program constructs such as lexically scoped local variables, procedure parameters, unrestricted return statements, etc. Our main contribution lies in a novel *term representation* of the aforementioned expectation transformer, denoted as infer[P]. Its definition follows the pattern of et[[P]], but notably differs in the definition of procedure calls and loops, where we replace the underlying fixed-point constructions via suitably constrained first-order templates. Through this step, we manage to tightly over-approximate the precise, but higher-order semantics via first-order constraint generation susceptible to automation. Second, we establish a new and original automation of this quantitative analysis in our prototype implementation ev-imp.

<footnote n="2">²In the context of automated cost analysis this problem has been partly addressed through the notion of *resource parametricity* [Hoffmann 2011].</footnote>

<footnote n="3">³Automation is here understood here as “push-button” automation; no user-interaction is required.</footnote>
<figure><img src="image_2.png" /><figcaption>Fig. 1. Textbook Examples on Expected Value Analysis</figcaption></figure>

*Outline.* In Section 2 we provide a bird’s eye view on the contributions of this work. In Section 3 we detail the syntactic structure of our language PWhile and provide the definition of the aforementioned weakest pre-expectation semantics. Section 4 constitutes the main technical part of the work, detailing the conception and definition of a term representation of the weakest pre-expectation semantics, susceptible to automation. Conclusively in Section 5 we discuss implementation choices for our prototype implementation, the chosen benchmark suites and the experimental evaluation. Finally, we conclude with related works and future work in Sections 6 and 7, respectively. Omitted proofs can be found in the extended version [Avanzini et al. 2023].

## 2 EXPECTED VALUE ANALYSIS, AUTOMATED

We motivate the central contribution of our work: an *automated* expected result value analysis of non-deterministic, probabilistic programs, featuring *recursion*. First, we present three textbook examples (and their encoding in PWhile) motivating the interest and importance of an *expected value* analysis, in contrast to eg. an *expected cost analysis* (see eg. [Avanzini et al. 2020; Kaminski et al. 2018; Leutgeb et al. 2022; Ngo et al. 2018; Wang et al. 2020, 2019]) or the inference of quantitative invariants (see eg. [Bao et al. 2022; Katoen et al. 2010; Wang et al. 2018]). Second, we emphasise the need for formal verification techniques going beyond pen-and-paper proofs and provide the intuition behind the thus chosen methodology of so-called *expectation transformers*. Third, we detail challenges posed by *recursive procedures* for the formal definition of our expectation transformers. Finally, we highlight the challenges of automated verification techniques in this context.

*Textbook examples.* To motivate expected value analysis in general, we study corresponding textbook examples—taken from Cormen et al. [Cormen et al. 2009]—and suitable code representations thereof, depicted in Figure 1. Consider an unspecified number of bins and suppose we throw balls towards the bins. Let us attempt to answer the following three questions: (i) *How many balls will fall in a given bin?* (ii) *How many balls must we toss, on average, until a given bin contains a ball?* and finally (iii) *How many balls must we toss until every bin contains at least one ball?*

The first two questions can be easily answered given some mild background in probability theory. If we throw say $n$ balls, the number of balls per bin follows a binomial distribution. Assuming we will always hit at least one bin, then the success probability is $1/\text{bins}$, where $\text{bins}$ denotes the number of bins. Thus, the expected value of the number of balls in a single bin is given by $1/\text{bins} \cdot n$. Success in the second problem follows a geometric distribution. Thus the answer is $1/\text{bins} = \text{bins}$. But the last one is more tricky and involves a more intricate argumentation. Note that this problem and thus the procedure every is equivalent to the *Coupon Collector* problem, cf. [Avanzini et al. 2020; Cormen et al. 2009; Mitzenmacher and Upfal 2005]. Following the argument in [Cormen et al. 2009, Chapter 5.4], we say that we have a “hit”, if we have successfully hit a specific bin. Using the notion of hits, we can split the task into stages, each corresponding to a completed hit. If all stages are complete, we are done. To complete stage $k$, we need to hit the $k^{\text{th}}$ bin. Thus, the success probability
to complete this stage is given as $^{bins-k+1/bins}$. Now let $n_k$ denote the number of throws in the $k^{th}$ stage. Thus, the total number required to fill the bins is $n = \sum_{k=1}^{bins} n_k$. As $\mathbb{E}[n_k] = \frac{bins}{bins-k+1}$, we obtain from linearity of expectations

$$ \mathbb{E}[n] = \mathbb{E}[\sum_{k=1}^{bins} n_k] = \sum_{k=1}^{bins} \mathbb{E}[n_k] = \sum_{k=1}^{bins} \frac{bins}{bins-k+1} = bins \cdot \sum_{k=1}^{bins} \frac{1}{k}, $$

Further, we obtain $bins \cdot \sum_{k=1}^{bins} \frac{1}{k} \in \Theta(bins \cdot \log(bins))$, utilising that the harmonic number $H_{bins} = \sum_{k=1}^{bins} \frac{1}{k}$ is asymptotically bounded by $\log(bins)$, cf. [Graham et al. 1994].

*Probabilistic, recursive procedures.* As depicted in Listings 1(a-c), it is easy to encode the above questions as (probabilistic, recursive) procedures. The encoding fixes the number of bins to five. The command *Bernoulli(p)* draws a value from a Bernoulli distribution, returning 1 with probability $p$ and 0 with probability $(1-p)$.<fnref n="4" /> Our simple imperative language PWhile follows the spirit of Dijkstra's *Guarded Command Language*, see Section 3 for the formal details. Apart from randomisation and recursion, the encoding makes—we believe natural—use of local variables, formal parameters and return statements.

The textbook questions above now become (possibly intricate) questions on program behaviours, that is, on the *expected value* of program variables wrt. the memory distributions in the final program state. Expected value analysis encompasses expected cost analysis, as we can often, that is, for almost surely terminating programs, represent program costs through a dedicated counter. Similarly, upper bounds to expected values as considered here, form approximations of quantitative invariants. In the experimental validation of our prototype implementation, we will see that wrt. the latter often our bounds are in fact exact, that is, constitute *invariants* (see Section 5).

It is well-known that in general pen-and-paper analyses of the expected value of programs are not an easy matter—as we have for example seen in the answer to the last question above—and more principled methodologies are essential, in particular if our focus is on (full) automation. To this avail, we build upon earlier work on *weakest pre-expectation semantics* of (probabilistic) programs and develop an *expectation transformer* for PWhile.

*Expectation Transformers.* Predicate transformer semantics—introduced in the seminal works of Dijkstra [1975]—map each program statement to a function between two predicates on the state-space of the program. Their semantics can be viewed as a reformulation of Floyd-Hoare logic [Hoare 1969]. Subsequently, this methodology has been extended to randomised programs by replacing predicates with so-called expectations—real valued functions on the program’s state-space. *Mem V*—leading to the notion of *expectation transformers* (see [Kaminski et al. 2018; McIver and Morgan 2005]) and the development of *weakest pre-expectation semantics* [Gretz et al. 2014].

In the following, we impose a pre-expectation semantics on programs in PWhile, thereby providing a formal definition of the expected value of program values. Adequacy of such semantics is well understood in the literature, for example wrt. operational semantics based on Markov Decision Processes [Kaminski et al. 2018] or probabilistic abstract rewrite systems [Avanzini et al. 2020]. The difference to our work seems incremental to us. Thus, our starting point are pre-expectation transformer semantics, permitting us to focus on inference.

For a given command C, its *expectation transformer* et$\llbracket C \rrbracket$ (detailed in Figure 2 in Section 3) maps a post-expectation $f$ to a pre-expectation, measuring the expected value that $f$ takes on the distribution of states resulting from running C, as a function in the initial state of C. Slightly simplifying, the expectation transformer assigns semantics to *commands*

$$ \text{et} \llbracket C \rrbracket : (\text{Mem } V \to \mathbb{R}^{+\infty}) \to (\text{Mem } V \to \mathbb{R}^{+\infty}). $$

<footnote n="4">⁴For ease of encoding, Listing 1(c) employs the sampling from *Bernoulli*($i/\epsilon$). Setting $i = bins - k + 1$ and $bins = 5$, we recover the textbook argumentation.</footnote>
In the special case where $f$ is a predicate, ie. a $\{0, 1\}$ valued function on memories, et $[[C]]$ $f \sigma$ yields the probability that $f$ holds after completion of C for any initial state $\sigma$; thereby generalising Djikstra's and Nielson's weakest-precondition transformer, cf. [Kaminski et al. 2018; Nielson 1987]. This definitions is probably most easily understood as a denotational semantics in continuation passing style, with post-expectation $f$ interpreted as continuation, cf. [Avanzini et al. 2021; Friedman and Wand 2008].

In this reading, $f: \text{Mem} V \to \mathbb{R}^{+\infty}$ measures a quantity on the continuations outcome given its inputs—in expectation— while et $[[C]]$ $f$ performs a backward reasoning lifting $f$ to a function on initial states of C. (In the following, we use “continuation” and “expectation” interchangingly.)

Let us consider the procedure balls depicted in Figure 1(a). As mentioned, the final value of $b$, the number of balls hitting a dedicated bin, follows a binomial distribution. In each recursive call the success probability is $1/5$. Thus, the expected value of the number of balls in a single bin is $1/5 \cdot n$. To provide a bird's view on the working of expectation transformers in this context, we verify this simple equality in the sequel. Calculating the pre-expectation from the continuation $f$ is straightforward, as long as the program under consideration contains neither recursive calls nor loops. To wit, let $\langle b \rangle$ be the real-valued function that measures (the positive value of) variable $b$ in a given memory. Then

$$ \text{et}[[\text{if (Bernoulli}(1/5)) \{b=b+1\}]] \langle b \rangle = 1/5 \cdot \text{et}[[b=b+1]] \langle b \rangle + 4/5 \cdot \langle b \rangle = 1/5 \cdot \langle b+1 \rangle + 4/5 \cdot \langle b \rangle. $$

Here, addition and scaling by a constant on expectations should be understood point-wise. Notice how this expression captures the fact $b$ is incremented only with probability $1/5$. Composability of the transformer extends this kind of reasoning to sequences of command, in the sense that command composition is interpreted as the composition of the corresponding expectation transformers. For straight-line programs C, et $[[C]]$ $f$ can be computed in a bottom up fashion.

*Loops and Recursive Procedures.* When dealing with loops or recursive procedures though, the definition of et $[[C]]$ becomes more involved. As usual in giving denotational program semantics, the definition of the expectation transformer of these self-referential program constructs is based on a fixed-point construction, cf. [Avanzini et al. 2020; Bao et al. 2022; Kaminski et al. 2018; Katoen et al. 2010; McIver and Morgan 2005; Wang et al. 2018]. This permits attributing precise semantics to programs. In our setting, procedures $p$ are interpreted via an expectation transformer

$$ \text{et}[[p]] : (\mathbb{Z} \times \text{GMem} \to \mathbb{R}^{+\infty}) \to (\mathbb{Z}^{\text{ar}(p)} \times \text{GMem} \to \mathbb{R}^{+\infty}), $$

where GMem the set of global memories. Following the definition of $p$, et $[[p]]$ turns a continuation, parameterised in the return value and global memory, into a pre-expectation in the formal parameters of $p$ and the global memory before execution of $p$. The definition (formalised in Section 3) is slightly technical to permit recursive definitions and to ensure proper passing of arguments and lexical scoping. For a procedure declared as def $p(\vec{x}) \{C\}$ and a post-expectation $f$, et $[[p]]$ $f$ is given by the transformer et $[[C]]$ associated to its body C, initialising the parameters $\vec{x}$ accordingly and with the post-expectation $f$ applied whenever C leaves its scope through a return-statement.

Let us re-consider the procedure balls from Listing 1(a). We use Iverson brackets $[\cdot]$ to lift predicates on memories to expectations, that is, $[B](\sigma) := 1$, if B holds in memory $\sigma$, and $[B](\sigma) := 0$, otherwise. Since the program does not make use of global variables, we elide representation of the empty global memory for brievity. Thus, for an arbitrary continuation $f: \mathbb{Z} \to \mathbb{R}^{+\infty}$, the expectation transformer of balls is given by the least functional satisfying:<fnref n="5" />

$$ \text{et}[[\text{balls}]] f = \lambda n. [n > 0] \cdot \text{et}[[\text{balls}]] (\lambda b. 1/5 \cdot f(b+1) + 4/5 \cdot f(b)) (n-1) + [n \le 0] \cdot f(0), (\dagger) $$

<footnote n="5">⁵Note that the functional is also ordered point-wise. Well-definedness of the employed fixed-point construction rests on the observations that expectation transformers form $\omega$-CPOs [Winskel 1993].</footnote>
where the continuation $f$ passed to the recursive call of $\text{et}[[\text{balls}]]$ is computed as above. Thus, the continuation to the call to balls probabilistically updates the returned value. Arguing inductively, we see that $\text{et}[[\text{balls}]]$ $f_r = \lambda n. \frac{1}{5} \cdot n + r$ for any continuation of the form $f_r := \lambda b. b + r$, where $r$ denotes a non-negative real. Since $f_0$ measures the return value, we conclude that on argument $n$, $\text{et}[[\text{balls}]]$ returns a value of $\frac{1}{5} \cdot n$ in expectation. Conclusively, we re-obtain that a fifth of the thrown balls will fall in the considered bin.

In the same vein, the analysis of the return value of throws and every—in expectation—is performed for Listings 1(b) and 1(c), respectively. In this fashion, the translation of the above pen-and-paper analysis for the textbook examples into the formalism of an expectation transformer results into a rigorous and principled methodology. This builds a sound (and suitable) foundation for subsequent automation.

*Automation.* The above calculation of the expected return value of procedure balls gives evidence on the advantages of a formal methodology over an ad-hoc pen-and-paper analysis. The crux of turning such a calculus into a fully automated analysis lies in automatically deriving closed forms for expected values of loops and recursive procedures. Related problems have been extensively studied in the literature, eg. [Contejean et al. 2005; Fuhs et al. 2007; Podelski and Rybalchenko 2004; Sinn et al. 2014]. One prominent approach lies in assigning templates to expected value functions, by which the definition of the expectation transformer can be reduced to a set of constraints treatable with off-the-shelf SMT solvers, like Z3 [de Moura and Bjørner 2008].

Following this approach, we express values as linear combination $\sum_i c_i \cdot b_i$ of base functions (or norms) $b_i$ with variable coefficients $c_i$, mapping program valuations to (non-negative) real numbers. Norms serve as a numerical abstraction of memories and encompass a variety of common abstractions, for example the absolute value of a variable, the difference between two variables and more generally arbitrary polynomial combinations thereof. Often, expected values can be computed symbolically on such *value expressions*. Concerning recursive procedures or loops, the definition of the expectation transform is in essence recursive itself. Rather than computing this fixed-point directly, we make use of Park's theorem [Kaminski et al. 2018; Wechler 1992] and seek an upper bound in closed-form. In the context of recursive procedure, we also have to model the call-stack appropriately.

Concerning recursive procedures, our solution is to represent every higher-order functional $\text{et}[[\text{p}]]$ through a pair of first-order templates $\langle H_p, K_p \rangle$, representing pre- and post-expectations respectively. Concretely, such a template signifies $\text{et}[[\text{p}]] K_p \le H_p$. Inspired by Hoare-calculi for recursive programs (cf. [Kleymann 1999]), we index templates by *logical variables* $x$, kept implicit in the sequel, but emphasised here for clarity of exposition. Thus indeed, a template represents a *families of pre-/post-expectations*. To wit, let us revisit procedure balls in Listing 1(a) once more. We take the templates:

$$H_{\text{balls}}^{x} := \lambda n. c_{0} + c_{1} \cdot \langle n \rangle + c_{2} \cdot \mathbf{x}$$

$$K_{\text{balls}}^{x} := \lambda b. \langle b \rangle + \mathbf{x},$$

Notice that the logical variable $x$ should only be instantiated non-negatively, so as to ensure that these templates indeed capture (non-negative) expectations. Thus, the indended meaning of this template for balls is captured by the constraint

$$\forall x \ge 0. \forall n. \text{et}[[\text{balls}]] K_{\text{balls}}^{x} n \le H_{\text{balls}}^{x} n, \quad (\ddagger)$$

Symbolic evaluation of the left-hand side—making use of the template for calls to balls—can now be used to sufficiently constrain the undetermined coefficients. Concretely $(\ddagger)$ is representable via the following three constraints, where the variables $b, n$ and $x$ are universally quantified, while the
unknowns $c_0$, $c_1$, $d_0$ and $d_1$ are existentially quantified.

$$0 \le \mathbf{x} \models \frac{1}{5} \cdot \langle b + 1 \rangle + \frac{4}{5} \cdot \langle b \rangle + \mathbf{x} \le \langle b \rangle + (d_0 + d_1 \cdot \mathbf{x}) \quad (c1)$$

$$0 \le \mathbf{x} \models 0 \le d_0 + d_1 \cdot \mathbf{x} \quad (c2)$$

$$0 \le \mathbf{x} \models [n > 0] \cdot (c_0 + c_1 \cdot \langle n - 1 \rangle + c_2 \cdot (d_0 + d_1 \cdot \mathbf{x})) + [n \le 0] \cdot \mathbf{x} \le c_0 + c_1 \langle n \rangle + c_2 \cdot \mathbf{x}. \quad (c3)$$

To see how these constraints have been formed, recall the fixed-point equation for et$_{[balls]}$ in (†), instantiating $f = K_{balls}^x$. Constraint (c1) expresses that the post-expectation

$$\lambda b. \frac{1}{5} \cdot K_{balls}^x (b + 1) + \frac{4}{5} \cdot K_{balls}^x (b) = \lambda b. \frac{1}{5} \cdot \langle b + 1 \rangle + \frac{4}{5} \cdot \langle b \rangle + \mathbf{x},$$

passed to the call of et$_{[balls]}$ in (†) is over-approximated by instance $K_{balls}^{d_1+d_1 \cdot \mathbf{x}}$, taking a (to be further determined) instantiation $\mathbf{x} \mapsto d_1 + d_1 \cdot \mathbf{x}$ of the logical variable $\mathbf{x}$. This substitution is guaranteed non-negative (and thus well-defined) by constraint (c2). Thus taking (c1) into account, $y$ (‡) e have

$$\mathrm{et}_{[\mathrm{balls}]} (\lambda b. \frac{1}{5} \cdot K_{\mathrm{balls}}^{x} (b+1) + \frac{4}{5} \cdot K_{\mathrm{balls}}^{x} (b)) (n-1) \leq H_{\mathrm{balls}}^{d_0+d_1 \cdot x} (n-1).$$

Making use of this over-approximation in (†), it is now evident that the final constraint (c3) witnesses satisfiability of the main constraint (‡). These three constraints can be met by taking $c_1, d_0 = \frac{1}{5}, c_2, d_1 = 1$, and $c_0 = 0$.

We re-obtain (unsurprisingly) that the expected return value of procedure balls is given as $\frac{1}{5} \cdot n$. We emphasise that parameterising templates through logical variables is crucial, even for relatively simple examples such as balls. Since the program is not tail-recursive, the continuation changes after each recursive call. Parametricity allows us to vary templates across recursive calls, as we have done above through instantiating $\mathbf{x}$ by $d_0 + d_1 \cdot \mathbf{x}$.

We have successfully automated this approach in our prototype implementation ev-imp, see Section 5. Automation is based on our first contribution, the development of a term representation of the expectation transformer inducing an inference method that describes the generation of constraints. In the context of procedure balls, this representation reduces the task of finding a functional satisfying (†) to the definition of a set of constraints like (‡), whose solution yields over-approximations of the function graph of et$_{[balls]}$. (See Section 4 for the details.) Based on this formal development, our tool ev-imp proceeds with an analysis as outlined, and integrates a dedicated constraint solver for solving constraints of the form (c1)–(c3). Apart from handling recursive procedures, we have incorporated the constraints for handling loop programs. Here, we take inspiration from [Avanzini et al. 2020] to guarantee a modular (and thus scalable) inference of upper bounding functions.

Wrt. the three motivating examples in Figure 1, our tool derives (upper-bounds to) the corresponding expected return values in milliseconds. The bounds for the procedures in Listings 1(a) and (b) are precise, respectively. For the procedure every we, however, only manage to derive a (sound) upper bound. The latter is not surprising. As shown in the textbook proof the expected number of throws is given as $\Theta(\mathrm{bins} \cdot \log(\mathrm{bins}))$ in general. Our template approach does not (yet) support logarithmic functions. Hence, we cannot hope to derive the precise bound fully automatically. To the best of our knowledge, our prototype ev-imp is the only existing tool able to fully automatically analyse the expected result value (or expected cost for that matter) of probabilistic, recursive programs, if formulated in an imperative language. The complete evaluation of our prototype implementation is given in Section 5.

## 3 AN IMPERATIVE PROBABILISTIC LANGUAGE

We consider a small imperative language in the spirit of Dijkstra's Guarded Command Language. In particular, the language features (i) dynamic sampling instructions; (ii) non-deterministic choice, formalised via a non-deterministic choice operator (<<); (iii) nested loops; and (iv) crucially
recursive procedure declarations. In this section, we first formalise the syntax and then endow it
with axiomatic, pre-expectation semantics.

Syntax. Let Var = {x, y, ...} be a set of (integer-valued) program variable. We fix three syntactic
categories of Boolean valued expressions BExpr V, (integer valued) expressions Expr V, and sampling
instructions SExpr V over (finitely many) variables V ⊆ Var, respectively. In the following, B
will range over Boolean, E, F over integer valued expressions and G over sampling instructions.
Furthermore, let Proc = {p, q, ...} be a set of procedure (symbols). The set of commands Cmd V
over program variables x ∈ V is given by

$$
\begin{array}{l}
C, D ::= \text{skip} | x \approx G | x \approx p(E_1, \dots, E_{ar(p)}) | \text{return}(E) | \text{var } x \leftarrow E \text{ in } \{C\} \\
      | C;D | \text{if } (B) \{C\} \text{else } \{C\} | \text{while } (B) \{C\} | C \Rightarrow D
\end{array}
$$

The interpretation of these commands is fairly standard. The command skip acts as a no-op. In an assignment x := R, the right-hand side R evaluates to a distribution, from which a value is sampled and assigned to x. As right-hand side R, we permit either built-in sampling expressions G ∈ SExpr such as unif(lo, hi) for sampling an integer uniformly between constants lo and hi, or calls to user-defined procedures p ∈ Proc. A command var x ← E in {C} declares a local variable x, initialised by E, within command C. Further, commands can be defined by composition, via a conditional conditional, via a while-loop construct or via the non-deterministic choice operator C <> D, interpreted in a demonic way. For ease of presentation, we elide a probabilistic choice command and probabilistic guards (see eg [Kaminski et al. 2016, 2018]), since they do not add to the expressiveness of the language. In examples, however, we make use of mild syntactic sugaring to simplify readability, as we did already above.

A *program P* is given as a finite sequence of procedure definitions. Each procedure *p* expects a number of integer-valued arguments and returns an integer upon termination. The number of arguments of a procedure *p* is called its *arity* and denoted as *ar(p)*. The body of *p*, denoted as *Bdy*<sub>*p*</sub>, consists of a command *C* ∈ *Cmd V* that may refer both to formal parameters, locally and globally declared variables. Note that since *p* may trigger a sampling instruction, its evaluation evolves probabilistically, yielding a final value and modifying the global state with a certain probability. Formally, a program is a tuple *P* = (*GVar*, *Decl*) where *GVar* ⊆ *Var* is a finite set of *global variables*, and where *Decl* is a finite sequence of *procedure declarations* of the form

$$
\text{def } p(x_1, \dots, x_{ar(p)}) \{C\}.
$$

To avoid notational overhead due to variable shadowing, we assume that the *formal parameters*
*Args<sub>p</sub>* := *x*₁, ..., *x*<sub>ar(*p*)</sub> and global variable are all pairwise different, and distinct from variables
locally bound within *C*. Using *α*-renaming, this can always been guaranteed. Throughout the
following, we keep the program *P* = (*GVar*, *Decl*) fixed.

Weakest Pre-Expectations Semantics. A memory (or state) over finite variables V ⊆ Var is a mapping σ ∈ Mem V := V → Z from variables to integers. We write GMem := Mem GVar for the set of global memories, and σ|<sub>V</sub> for its restriction to variables in V. Let x ∈ V be a variable and let v ∈ Z. Then we write σ[x ↦ v] for the memory that is as σ except that x is mapped to v. As short-forms, we write σ<sub>g</sub> for the global memory σ|<sub>GVar</sub>, and dual, σ<sub>l</sub> for the local memory σ|<sub>Var \ GVar</sub>. Let R<sup>+∞</sup> denote the set of non-negative real numbers extended with ∞, ie. R<sup>+∞</sup> := R<sub>≥0</sub> ∪ {∞}. A (discrete) subdistribution over A is a function δ: A → R<sub>≥0</sub> so that Σ<sub>a∈A</sub> (δ a) ≤ 1, and a distribution, if Σ<sub>a∈A</sub> (δ a) = 1. We may write (sub)distributions δ as {{δ a : a}}<sub>a∈A</sub>. The set of all subdistributions over A is denoted by D A. We restrict to distributions over countable sets A. The expectation of a function f: A → R<sup>+∞</sup> wrt. a distribution δ is given by E<sub>δ</sub> f := Σ<sub>a∈A</sub> (δ a) · (f a). We suppose that expressions E ∈ Expr V, Boolean expressions B ∈ BExpr V and sampling expressions G ∈ SExpr V
<table><caption>Fig. 2. Expectation Transformer Semantics of PWhile.</caption><tr><td>et[<b>p</b>]<sup>η</sup> : (Z → GMem → ℝ<sup>+∞</sup>) → (Z<sup>ar(p)</sup> → GMem → ℝ<sup>+∞</sup>)</td></tr><tr><td>et[<b>p</b>]<sup>η</sup> f</td><td>:= λ<sup>v̄</sup>σ. et[<b>Bdy<sub>p</sub></b>]<sup>η</sup><sub>f</sub> (λτ. f 0 τ<sub>g</sub>) (σ ⊎ {Args<sub>p</sub> ↦ v̄})</td></tr><tr><td>et[<b>C</b>]<sup>η</sup><sub>f</sub>: (Mem V → ℝ<sup>+∞</sup>) → (Mem V → ℝ<sup>+∞</sup>)</td></tr><tr><td>et[<b>skip</b>]<sup>η</sup> g</td><td>:= g</td></tr><tr><td>et[<b>x</b> ≈ <b>G</b>]<sup>η</sup> g</td><td>:= λσ. E<sub>[<b>G</b>]</sub> σ (λv. f[<b>x</b> ↦ <b>v</b>])</td></tr><tr><td>et[<b>x</b> ≈ <b>p</b>(<b>E</b>)]<sup>η</sup> g</td><td>:= λσ. η p (λv. τ. g (σ<sub>l</sub> ⊎ τ) [<b>x</b> ↦ <b>v</b>]) (<b>E</b><sup>v̄</sup> σ) σ<sub>g</sub></td></tr><tr><td>et[<b>return</b>(<b>E</b>)]<sup>η</sup> g</td><td>:= λσ. f (<b>E</b>] σ) σ<sub>g</sub></td></tr><tr><td>et[<b>var</b> <b>x</b> ← <b>E</b> in {<b>C</b>}]<sup>η</sup> g</td><td>:= λσ. et[<b>C</b>]<sup>η</sup> g σ [<b>x</b> ↦ (<b>E</b>] σ)</td></tr><tr><td>et[<b>C</b>;<b>D</b>]<sup>η</sup> g</td><td>:= λσ. et[<b>C</b>]<sup>η</sup> (et[<b>D</b>]<sup>η</sup> g) σ</td></tr><tr><td>et[<b>if</b> (B) {<b>C</b>} else {<b>D</b>}]<sup>η</sup> g</td><td>:= λσ. [<b>B</b>] σ) · et[<b>C</b>]<sup>η</sup> g σ + [<b>¬B</b>] σ) · et[<b>D</b>]<sup>η</sup> g σ</td></tr><tr><td>et[<b>while</b> (B) {<b>C</b>}]<sup>η</sup> g</td><td>:= λσ. lfp (λG. λτ. [<b>B</b>] τ) · et[<b>C</b>]<sup>η</sup> G τ + [<b>¬B</b>] τ) · g τ) σ</td></tr><tr><td>et[<b>C</b> <> <b>D</b>]<sup>η</sup> g</td><td>:= λσ. max et[<b>C</b>]<sup>η</sup> g, et[<b>D</b>]<sup>η</sup> g) σ</td></tr></table>

are equipped with interpretations [E] : Mem V → Z, [B] : Mem V → B and [G] : Mem V → DZ, respectively. Functions in A → ℝ<sup>+∞</sup> are called expectation (over a set A) in the literature and are usually denoted by f, g, h etc. We equip expectations and transformers with the point-wise ordering, that is, f ≤ g if fa ≤ ga for all a ∈ A. We also extend functions over A point-wise to expectations and denote these extensions in typewriter font, eg., f + g := λa. fa + ga etc. In particular, θ = λa.0 and ∞ = λa.∞. Equipped with this ordering, expectations form an ω-CPO [Winskel 1993], whose least element is the constant zero function θ.

Expectation transformer have the shape $F: (A \to \mathbb{R}^{+\infty}) \to (B \to \mathbb{R}^{+\infty})$ (for some sets A and B respectively). Ordered point-wise these again form an ω-CPO and have thus enough structure to give a denotational model to programs. In particular, when $A = B$ the least fixed-point lfp (F) of F is well-defined and given by $\sup_n (F^n \perp_A)$, for $F^n$ the n-fold composition of F [Winskel 1993]. Following Dijkstra [Dijkstra 1976], expectation transformers can be seen as giving rise to a (denotational) semantics.<fnref n="6" />

The transformer et[<b>P</b>] is defined in terms of an expectation transformer et[<b>p</b>], p ∈ P, for procedures, which in turn is mutual recursively defined via an expectation transformer et[<b>C</b>] on commands, cf. Figure 2. These transformers are parameterised in a procedure environment η, of type

$$p: \text{Proc} \to (\mathbb{Z} \times \text{GMem} \to \mathbb{R}^{+\infty}) \to (\mathbb{Z}^{\text{ar}(p)} \times \text{GMem} \to \mathbb{R}^{+\infty}),$$

associating an expectation semantic to each p ∈ Proc, which is used in the case of procedure calls x ≈ p(E). As already alluded to, the definitions given in Figure 2 are best understood as a denotational semantics in continuation passing style, with post-expectations f : Z × GMem → ℝ<sup>+∞</sup>

<footnote n="6">⁶Such a transformer constitutes a standard denotational semantic, where effects are interpreted in the continuation monad $\text{Cont}_{\mathbb{R}^{+\infty}}(-) = ((-)) \to \mathbb{R}^{+\infty}) \to \mathbb{R}^{+\infty}$, cf. [Avanzini et al. 2021]. To observe this, flip the arguments f and σ in Figure 2.</footnote>
<table><caption>Fig. 3. Expectation Transformer Laws.</caption><tr><td>monotonicity</td><td>η ≤ χ ∧ f<sub>1</sub> ≤ f<sub>2</sub> ∧ g<sub>1</sub> ≤ g<sub>2</sub> ⇒ et[C]<sub>f<sub>1</sub></sub><sup>η</sup> g<sub>1</sub> ≤ et[C]<sub>f<sub>2</sub></sub><sup>χ</sup> g<sub>2</sub></td></tr><tr><td>linearity</td><td>et[C]<sub>Σ<sub>i</sub> R<sub>i</sub> · f<sub>i</sub></sub><sup>η</sup> (Σ<sub>i</sub> R<sub>i</sub> · g<sub>i</sub>) = Σ<sub>i</sub> R<sub>i</sub> · et[C]<sub>f<sub>i</sub></sub><sup>η</sup> g<sub>i</sub></td></tr><tr><td>loop-invariant</td><td>[¬B] · g<sub>1</sub> ≤ g<sub>2</sub> ∧ [B] · et[C]<sub>f</sub><sup>η</sup> g<sub>2</sub> ≤ g<sub>2</sub> ⇒ et[while (B) {C}]<sub>f</sub><sup>η</sup> g<sub>1</sub> ≤ g<sub>2</sub></td></tr><tr><td>procedure-invariant</td><td>∀p ∈ Proc. et[p]<sub>f</sub><sup>χ</sup> ≤ χp ⇒ et[P] ≤ χ</td></tr></table>

and $g: \text{Mem } V \to \mathbb{R}^{+\infty}$ being interpreted as continuations, respectively. In short, the functional

$$ \text{et}[p] : (\mathbb{Z} \times \text{GMem} \to \mathbb{R}^{+\infty}) \to (\mathbb{Z}^{\text{ar}(p)} \times \text{GMem} \to \mathbb{R}^{+\infty}), $$

turns a continuation in the return-value and (possibly modified within the body of p) global memory, and lifts it to a function in the formal parameters of p and the initial global memory. To this end, et[p] links formal parameters $\text{Arg}_p$ to the input values $\vec{v}$, and yields to the transformer et[Bdy$_p$]$_f$ associated to its body. As denoted, this auxiliary transformer on commands is parameterised by a procedure environment $\eta$, and continuation f of the procedure p. The latter is necessary to model evaluation when Bdy$_p$ pre-maturely leaves scope through a return statement. To capture the sitation where evaluation completes without encountering a return, the continuation to et[Bdy$_p$]$_f$ supplies to f a default return value of zero.

To give some intuitions on the transformer of commands, it is again advisable to think of et[C]<sub>f</sub><sup>η</sup> g σ as running the command C on memory σ, and then proceeding with continuation g. The only exception to this reading lies in the treatment of sampling instructions.

For a sampling instructions $x \approx G$, et[C]<sub>f</sub><sup>η</sup> g σ computes the expected value of the continuation g on the distribution of stores obtained by updating x with elements sampled from G. (As a simplified and rough intuition, think of the assignment rule in Hoare logics.) For a procedure call $x \approx p(\vec{E})$, we make use of the semantics $\eta p$ of p, which is applied to the evaluated arguments and the current global memory $\sigma_g$. The continuation passed to $\eta p$ runs the continuation g on the combination of global memory $\tau$ yielded by p and the local pre-memory $\sigma_l$, with x updated by the return value yielded by p. When C is a return statement, the transformer skips continuation g and jumps directly to the continuation f defined by the enscoping procedure, supplying the returned value. For a local variable declarations, et[\varpi x \leftarrow E \text{ in } \{C\}]<sub>f</sub><sup>η</sup> g implement lexical scope, updating variable x in $\sigma$ by E. Due to the variable convention—as emphasised on page 8—the local variable x is fresh and thus cannot conflict with any variable in V.

The transformer for composed commands is given by the composition of the corresponding transformers. In the case of conditionals, we use Iverson bracket [·] to interpret Boolean values false and true as integers 0 and 1, respectively. Thus the transformer of conditionals effectively recurses on one of the two branches of the conditional, depending on the condition B. The expectation transformer for loops can be seen as (the least fixed-point) satisfying

$$ \text{et[while (B) \{C\}]}_f^{\eta} g \sigma = \begin{cases} \text{et[C;while (B) \{C\}]}_f^{\eta} g \sigma & \text{if } [B] \sigma = \text{true,} \\ g \sigma & \text{if } [B] \sigma = \text{false,} \end{cases} $$

running C once followed by the while loop in case the guard holds; and calling the continuation g otherwise. Finally, non-determinism is modelled as the maximum of the pre-expectations obtained from the alternatives. This is motivated by the fact that we are interested in worst-case bounds (on
expected values and costs) and follows the treatment of non-determinism in the context of program analysis, cf. [Nielson et al. 1999]. We note that our treatment of non-determinism constitutes a conceptual difference to the *weakest pre-expectation calculus* of [McIver and Morgan 2005]. There the focus is on (quantitative) program behaviours and thus a lower-bound to the pre-expectation is sought which results in choosing the *minimum* of the pre-expectations to handle non-deterministic choice.

What remains is to set up the procedure environment $\eta$ according to the declarations in P. To this end, we associate the semantics et$\llbracket P \rrbracket$ with the procedure environment $\eta$ that makes each $p \in \text{Proc}$ adhere to the semantics et$\llbracket p \rrbracket$' $f$, that is, that satisfies the (least) fixed point $\eta p f = \text{et} \llbracket p \rrbracket$' $f$. More precisely, et$\llbracket P \rrbracket := \text{lfp} (\lambda \Xi. \lambda p. \text{et} \llbracket p \rrbracket^\Xi)$. Again, this least (higher-order) fixed-point exists as expectation transformers form an $\omega$-CPO.

We emphasise that all the individual transformers are defined mutually, thus permitting mutual recursion on procedure declarations. In the following, we write et$\llbracket \cdot \rrbracket$ instead of et$\llbracket \cdot \rrbracket_f$ when $\eta = \text{et} \llbracket P \rrbracket$ and caller expectation $f$ is clear from context. Table 3 lists laws, tacitly employed above, for the expectation transformers et$\llbracket p \rrbracket$ and et$\llbracket C \rrbracket$, respectively.

# 4 INFERENCE

In this section we present the key contribution of this work, the development of a *term representation* of the expectation transformer et$\llbracket P \rrbracket$, see Figure 4. This forms the crucial basis of the automated inference of upper bounds to et$\llbracket P \rrbracket$, as implemented in our prototype ev-imp.

To a great extent, the definition of infer[p], p a procedure, follows the pattern of the definition of et$\llbracket p \rrbracket$. Notably, however, it differs in the definition of procedure calls, where we employ the templates $\langle H_p, K_p \rangle$ outlined in Section 2 and the definition of loops, where we replace the fixed-point construction via a suitably constrained template. Theorem 4.1 verifies that this approximation is sound; in proof we make essential use of the invariant laws depicted in Figure 3. We represent expectations syntactically as terms, denoting linear combinations of norms:

$$ \text{Norm } V Z \ni N ::= [B] \cdot E \quad (\text{norms}) \qquad \text{Term } V Z \ni T, S, U ::= \sum_i R_i \cdot N_i \quad (\text{terms}) $$

$E \in \text{Expr}(V \Join Z)$ denotes an arbitrary expression over *program variables* $V$ and *logical variables* $Z = \{x, y, \dots\}$. $B \in \text{BExpr}(V \Join Z)$ denotes a Boolean expression over program and logical variables. Coefficients $R$ denote terms yielding non-negative real numbers.

We require that for any norm $N$, $E$ is non-negative, whenever $B$ holds. A norm abstracts an expression over a program variable as a non-negative real number. For instance, $\max(x, 0) = [x \ge 0] \cdot x$, or $[x \ge y] \cdot (x - y)$ gives the distance from $x$ to $y$. For brevity, we set $\langle x \rangle := [x \ge 0] \cdot x$. Following the semantics, we set $[B] \cdot (\sum_i R_i \cdot ([B_i] \cdot E_i)) = \sum_i R_i \cdot ([B \wedge B_i] \cdot E_i)$. Note that, since norms are non-negative, a term $T \in \text{Term} V \emptyset$ can be interpreted as an expectation $\llbracket T \rrbracket := \text{Mem} V \to \mathbb{R}^{+\infty}$ over the state-space. Let $G \in \text{SExpr}$ denote a sampling instruction and let $T \in \text{Term} V \emptyset$, then $E_{x \leftarrow G} T$ denotes the term representation of the expectation $E_{\llbracket G \rrbracket} \sigma \llbracket T[x \mapsto v] \rrbracket$ of the continuation $\llbracket T[x \mapsto v] \rrbracket$ wrt. the distribution $\llbracket G \rrbracket$ applied to the current memory $\sigma$.

The definition of infer$\llbracket \cdot \rrbracket$ (see Figure 4) is best understood as a syntactic representation of the expectation transformer et$\llbracket \cdot \rrbracket$. It translates this denotational semantics into first-order constraints, susceptible to automation. Apart from returning a term, representing a pre-expectation, it generates a set of *side-conditions* of the form $\Gamma \vdash S \leq T$. Such a constraint is valid, if for all logical variables $\vec{x}$ occurring in the expressions $B, S$ and $T$, respectively, we have $\llbracket B[\vec{x} := \vec{v}] \rrbracket \models \llbracket S[\vec{x} := \vec{v}] \rrbracket \models \llbracket T[\vec{x} := \vec{v}] \rrbracket$ for all $\vec{v}$. In this reformulation, the procedure environment $\eta$ is kept implicit. The semantics of each $p \in \text{Proc}$ is thus representable as a pair of terms $K_p \in \text{Term} G \text{Var} \Join \{r\} Z$, $\langle H_p, K_p \rangle$, with
<table><caption>Fig. 4. Term Representation of Expectation Transformer Semantics.</caption><tr><td>infer[p] : Term GVar {r} → Term GVar {ā<sub>p</sub>}</td></tr><tr><td>infer[p] S</td><td>:= (infer[Bdy<sub>p</sub>]<sub>S</sub> S[r ↦→ 0])[Args<sub>p</sub> ↦→ ā<sub>p</sub>]</td></tr><tr><td>infer[C]<sub>S</sub> : Term V {z̄<sub>p</sub>} → Term V {z̄<sub>p</sub>}</td></tr><tr><td>infer[skip]<sub>S</sub> T</td><td>:= T</td></tr><tr><td>infer[x :≈ G]<sub>S</sub> T</td><td>:= E<sub>x←G</sub> T</td></tr><tr><td>infer[x :≈ q(E)]<sub>S</sub> T</td><td>:= H<sub>q</sub>[ā<sub>q</sub> ↦→ E, z̄<sub>q</sub> ↦→ U] where U ∈ Term (V \ GVar) {z̄<sub>p</sub>}<br>Γ<sub>q</sub> ⊢ Γ<sub>q</sub>[z̄<sub>q</sub> ↦→ U]; Γ<sub>q</sub> ⊢ T[x ↦→ r] ≤ K<sub>q</sub>[z̄<sub>q</sub> ↦→ U]</td></tr><tr><td>infer return(E)<sub>S</sub> T</td><td>:= S[r ↦→ E]</td></tr><tr><td>infer[var x ← E in {C}]<sub>S</sub> T</td><td>:= (infer[C]<sub>S</sub> T)[x ↦→ E]</td></tr><tr><td>infer[C;D]<sub>S</sub> T</td><td>:= infer[C]<sub>S</sub> (infer[D]<sub>S</sub> T)</td></tr><tr><td>infer[if (B) {C} else {D}]<sub>S</sub> T</td><td>:= [B] · infer[C]<sub>S</sub> T + [¬B] · infer[D]<sub>S</sub> T</td></tr><tr><td>infer[while (B) {C}]<sub>S</sub> T</td><td>:= U B ⊢ infer[C]<sub>S</sub> U ≤ U; ¬B ⊢ T ≤ U</td></tr><tr><td>infer[C <> D]<sub>S</sub> T</td><td>:= U ⊢ infer[C]<sub>S</sub> T ≤ U; ⊢ infer[D]<sub>S</sub> T ≤ U</td></tr></table>

$H_p \in \text{Term GVar} \Leftrightarrow \{\bar{a}_p\} Z$ and where $\bar{a}_p = a_1, \dots, a_{ar(p)}$ and $r$ are dedicated variables, referring to the formal parameters and the return value of $p$.

For each $p$, the terms $H_p$ and $K_p$ can be understood as families of terms, parameterised in the substitution of logical variables. To wit, let $\theta$ denote an arbitrary substitution of logical variables for values $\vec{v}$, then $[H_p]\theta$ and $[K_p]\theta$ denote pre- and post-expectations $h_{\vec{v}}: \mathbb{Z}^n \times \text{GMem} \to \mathbb{R}^{+\infty}$ and $k_{\vec{v}}: \mathbb{Z} \times \text{GMem} \to \mathbb{R}^{+\infty}$, respectively. Conceptually, the logical variables model resource parametricity [Hoffmann 2011], which is required to model the call-stack in recursive calls suitably. To improve upon the expressiveness of templates, we require that the terms $H_p$ and $K_p$ are non-negative only under an associated constraint $\Gamma_p$ on logical variables.

For instance, in Section 2, we implicitly used the logical context $\Gamma_{\text{balls}} = (0 \le x)$ to ensure that templates $H_{\text{balls}} = c_0 + c_1 \cdot \langle a \rangle + c_2 \cdot x$ and $K_{\text{balls}} = r + x$ are non-negative. To ensure that all the pairs $\langle H_p, K_p \rangle_{p \in \text{Prog}}$ adhere to the semantics of P, wrt. to their associated contexts $\Gamma_p$, we finally require (for all $p \in P$)

$$ \Gamma_p \vdash (\text{infer}[p] K_p)[\text{Args}_p \mapsto \bar{a}] \leq H_p. \quad (1) $$

The left-hand side of the inequality in (1) may reference the pair $\langle H_p, K_p \rangle$, namely when $p$ calls itself recursively.

As we have already seen in the informal description in Section 2 this recourse may be performed for an *instantiation* $[\vec{z}_p \mapsto \vec{U}]$ of logical variables. In Section 5 we detail how this instance is chosen. Under the intended meaning of $\langle H_p, K_p \rangle$, any application of $\langle H_p, K_p \rangle$ to an expectation $T$ (an alternation of $K_p$) can safely be replaced by $H_p$, provided the passed expectation $T$ is bounded again from above by $K_p$, viz, the corresponding constraint $\Gamma_q \vdash T[x \mapsto r] \le K_q[\vec{z}_q \mapsto \vec{U}]$ in Figure 4.

**THEOREM 4.1 (SOUNDNESS THEOREM).** If for all $p \in P$ the constraint (1) is fulfilled and all side-conditions in Figure 4 are met, then for all $p \in P$, et $[[p]] \le [[\text{infer}[p]]]$, that is, the inference algorithm is sound.
<figure><img src="image_3.png" /><figcaption>Fig. 5. Studies in Expected Value Analysis</figcaption></figure>

(a) Generating a biased coin [Wang et al. 2018].

(b) Binomial update [Katoen et al. 2010].

(c) Hire a new assistant [Cormen et al. 2009].

**Templating approach.** In the inference of upper invariants, we follow the *templating approach* standard in the literature (see eg. [Avanzini et al. 2020; Leutgeb et al. 2022; Wang et al. 2020]) in which the functions to be synthesised—in our case a closed form of infer[P]—are given as linear combination $\sum_i c_i \cdot b_i$ of pre-determined *base functions* $b_i$, with variable coefficient $c_i$. We emphasise that base functions can be *linear* or *non-linear*. Concretely, for straight-line commands C this is captured by the definition of the term representation infer[C] in the context of the continuation T. In Section 2 we have seen an informal account of this recipe. As a slightly more involved example, consider procedure *biasedcoin* depicted in Listing 5(a).<fnref n="7" /> The procedure incorporates a *non-deterministic choice*, denoted by the conditional with unspecified guard (*).

We focus on the *non-deterministic choice* and the *sampling instruction* incorporated. Wrt. C <> D the term representation eludes an explicit representation of the maximum function max. Instead non-determinism is resolved by asserting the constraints (i) $\vdash \text{infer[C]}_S T \le U$ and (ii) $\vdash \text{infer[D]}_S T \le U$, respectively. Wrt. sampling instructions, $\text{infer}[x \approx G]_S T$ is defined as the term $E_{x \leftarrow G} T$, representing the computation of the expected value of the continuation T symbolically on the distribution of the memory obtained by sampling elements according to the instruction G.

To illustrate, let C denote the body of the procedure *biasedcoin*. By default, we choose the return value $x_1$ as continuation, which we abstracted by the norm $\langle x_1 \rangle$. Restricting to the non-trivial branch of the non-deterministic choice, we obtain by symbolic calculation that $\text{infer}[C] \langle x_1 \rangle = \langle x_1 \rangle + [x_1 \ge \frac{1}{2} + x_2]$. Solely by symbolic execution, we obtain the constraint

$$ \vdash \langle x_1 \rangle + [x_1 \ge \frac{1}{2} + x_2] \le c_0 \cdot [x_1 \ge \frac{1}{2} + x_2] \cdot \langle x_1 \rangle + c_1 \cdot [x_2 \ge x_1] \cdot \langle x_1 \rangle + \frac{1}{2} \cdot c_2 \cdot [x_1 \ge \frac{1}{2} + x_2], $$
solvable with $c_0 = c_1 = c_2 = 1$, yielding $\langle x_1 \rangle + [x_1 \ge \frac{1}{2} + x_2]$ as the desired expected return value.

**Loop programs.** Considering loops, we employ the *loop-invariant law* to derive a closed form, cf. Figure 3. Conclusively, for a loop statement *while* (B) {C}, $\text{infer[while (B) {C}]}_S T$ asserts the constraints (i) $B \vdash \text{infer[C]}_S U \le U$ and (ii) $\neg B \vdash T \le U$, respectively. Ie. U represents an upper bound $I_T$, parameterised in the post-expectation T. Due to its reminiscence with a loops invariant, the function $I_T$ is called an *upper invariant* in the literature [Kaminski et al. 2018].

To illustrate, consider the procedure *binomial_update* from Figure 5(b).<fnref n="8" /> Again we approximate the return value by the norm $\langle x \rangle$. The above constraints on the term U induce the following constraints on an upper invariant $I_{\langle x \rangle}$ (i) $N \le n \vdash \langle x \rangle \le I_{\langle x \rangle}$ and (ii) $n < N \vdash \frac{1}{2} \cdot \text{et} \llbracket x \approx x+1; n \approx n+1 \rrbracket I_{\langle x \rangle} \le$

<footnote n="7">⁷Listing 5(a) is taken from Wang et al. [Wang et al. 2018] and—making use of a non-deterministic abstraction—constitutes a variant of an example considered by Katoen et al. [Katoen et al. 2010]. The latter example uses a stream of fair coin flips to generate a (single) biased coin.</footnote>

<footnote n="8">⁸The code sets variable $x$ to a value between 0 and $N$, following a binomial distribution, cf. [Katoen et al. 2010].</footnote>
$I_{\langle x \rangle}$, respectively. Making use of linear template based on base functions $1, \langle x \rangle$ and $\langle N - n \rangle$ this can be instantiated as the following constraints.

$$ 
\begin{array}{l}
N \le n \models \langle x \rangle \le c_0 + c_1 \cdot \langle N - n \rangle + c_2 \cdot \langle x \rangle \\
n < N \models c_0 + c_1 \cdot \langle N - (n+1) \rangle + \frac{1}{2} \cdot c_2 \cdot (\langle x \rangle + \langle x+1 \rangle) \le c_0 + c_1 \cdot \langle N - n \rangle + c_2 \cdot \langle x \rangle,
\end{array}
$$

solvable with $c_0 = c_1 = 0$ and $c_2 = \frac{1}{2}$, yielding the quantitative invariant $\frac{1}{2} \cdot x$, that is, the derived upper invariant is optimal, cf. [Katoen et al. 2010].

**Recursive procedures.** Since we represent et $[[p]]$ through the term representation infer $[p]$, a similar approach can be suited to recursively defined procedures $\text{def } p(\bar{x}) \{\text{Bdy}_p\}$. Instead of the loop-invariant law, we implicitly employ the law on procedure-invariants to derive a closed form, though (see Figure 3). Thus, the definition of infer $[x \approx q(\bar{E})]_S T$ asserts the constraints (i) $\Gamma_q \vdash \Gamma_q[\bar{z}_q \mapsto \bar{U}]$ and (ii) $\Gamma_q \vdash T[x \mapsto r] \le K_q[\bar{z}_q \mapsto \bar{U}]$, respectively. Here, the latter constraint guarantees that the continuation $T[x \mapsto r]$ of the procedure call is bounded by an instance of the bounding function $K_q$. On the other hand the first constraint guarantees that the substitution $[\bar{z}_q \mapsto \bar{U}]$ of logical variables employed is properly represented in the context information. Consider procedure hire depicted in Listing 5(c).<fnref n="9" /> Following the recipe employed for the procedure balls, we generate the templates (i) $K_{\text{hire}} := \langle \mathbf{r} \rangle + \mathbf{x}$ and (ii) $H_{\text{hire}} := c_0 + c_1 \cdot \langle \mathbf{n} \rangle + c_2 \cdot \mathbf{x}$ where $\mathbf{n}$ refers to the parameter taken by hire, and $\mathbf{x}$ is a logical variable subject to the constraint $\Gamma_{\text{hire}} := (0 \le \mathbf{x})$. This results in the following three constraints

$$ 
\begin{array}{l}
0 < \mathbf{n}, \Gamma_{\text{hire}} \models (1+r)/n + (1-1/n)\langle\mathbf{r}\rangle + \mathbf{x} \le \langle\mathbf{r}\rangle + (d_0 + d_1 \cdot \mathbf{x}) \\
\Gamma_{\text{hire}} \models 0 \le d_0 + d_1 \cdot \mathbf{x} \\
\Gamma_{\text{hire}} \models [0 < \mathbf{n}] \cdot (c_0 + c_1 \cdot \langle \mathbf{n} - 1 \rangle + c_2 \cdot (d_0 + d_1 \cdot \mathbf{x})) + [\mathbf{n} \le 0] \cdot \mathbf{x} \le c_0 + c_1 \cdot \langle \mathbf{n} \rangle + c_2 \cdot \mathbf{x},
\end{array}
$$

solvable with $c_0 = 0$, $c_1 = c_2 = d_0 = d_1 = 1$. Here, the additional premise $0 < \mathbf{n}$ in the first constraint stems from a (forward) analysis of the conditional governing the call to hire. This yields $\langle \mathbf{n} \rangle$ as (sub-optimal) upper bound to expected number of hires in hire. Note that the expected value for the number of expected hires is given as harmonic number $H_k \in \Theta(\log(n))$. So linear is the best we can do with the templates provided in our prototype implementation ev-imp.

# 5 AUTOMATION

We have implemented the outlined procedures, proven sound in Theorem 4.1, in a prototype implementation, dubbed ev-imp.<fnref n="10" /> Our tool ev-imp estimates upper-bounds on the expected (normalised) return value of procedures p, as a function of the inputs, that is, ev-imp computes an upper-bound to et $[[p]]$ $(\lambda v_{\_\_} \langle v \rangle)$. We note that the restriction to measurements of the expected return value does not constitute a real restriction but rather constitutes a slight design simplification, as long as the estimated continuation f is representable as a return expression E. The term representation from Figure 4 forms the basis of our prototype's implementation ev-imp. Below, we highlight the main design choices that lead us from the term representation to a concrete algorithm and provide ample experimental evidence of the algorithmic expressivity of our prototype implementation.

**Assignments.** Our implementation supports sampling from finite, discrete distributions G, assigning probabilities $p_i$ to values $e_i$ ($0 \le i \le k$ for constant k), where probabilities $p_i$ are expressed as rationals and values $e_i$ as integer expressions. Both, probabilities and expressions can possibly

<footnote n="9">⁹The procedure represents a hiring process, encoded as probabilistic program, cf. [Cormen et al. 2009].</footnote>

<footnote n="10">¹⁰Our prototype implementation employs libraries of the open-source tool, freely available at https://gitlab.inria.fr/mavanzin/ecoimp. In particular, we rely on its auxiliary functionality, such as program parsers etc., but also on its underlying constraint solver.</footnote>
depend on the current memory, thus our tool natively supports *dynamic sampling instructions* such as in *Bernoulli(1/n)*, used for example in our rendering of the textbook example Listing 5(c). Note that the probabilities depend on the value of *n* taken in the memory under which the distribution is evaluated. Such dynamic distributions are also required to represent our variant of the Coupon Collector's problem—procedure *every*, cf. Listing 1(c). Overall, *ev-imp* natively implements in this way a variety of standard distributions, in particular it supports sampling *uniform distributions with bounded support*, *Bernoulli*, *binomial* and *hypergeometric* distributions. Note, that through recursion more involved distributions can be built, with unbounded and dynamic support. Practical extensions, eg. further support for dynamic sampling are (in our opinion) engineering tasks that do not require novel theoretical insights.

*Template selection and instantiation.* The overall approach rests on templating to over-approximate the behaviour of procedures and loops. As indicated earlier, we describe these templates via linear combinations $\sum_i c_i \cdot N_i$ of pre-determined *norms* $N_i$ and undetermined coefficients $c_i$. To determine these templates, our implementation selects a set of candidate *base functions* from post-expectations, program invariants (determined through a simple forward-analysis) and loop-guards, loosely following the heuristics of Sinn et al. [2016] and Avanzini et al. [2020]. Specifically, an established invariant or guard $x \le y$ gives rise to the base function $\langle x - y \rangle$, modelling the (non-negative) distance between $x$ and $y$. This heuristic is extended from variables to expressions, and to arbitrary Boolean formulas, and turns out to work well in practice. These base functions are then combined to an overall *linear*, or *simple-mixed* template, cf. Contejean et al. [2005], of base functions. For instance, $c_1\langle x \rangle + c_2\langle y \rangle + c_3x^2 + c_4y^2 + c_5[x \ge 0 \land y \ge 0](x \cdot y) + c_6$ is a simple-mixed template over base functions $\langle x \rangle$ and $\langle y \rangle$. In a similar spirit, we make use of templates for the instantiations $\vec{U} \in \text{Term}(V \GVar) \{\bar{z}_p\}$, as linear functions in the local and logical variables, in the same vain as we have already done when presenting examples. Our prototype implements caching and backtracking to test different templates when operating in a modular setting (see the paragraph on modularity below). In particular, non-linear base functions are employed only if the linear ones fail. Here, we follow similar approaches in the literature, cf. [Avanzini et al. 2020; Leutgeb et al. 2021, 2022; Wang et al. 2018, 2020].

*Constraint solving.* Evaluation of $\inf[p]$ for a procedure *p* results in a set of constraints, whose solution is then used to assess bounding functions. In effect, these constraints are inequalities over real-valued polynomial expressions over unknown coefficients, enriched with conditionals (through Iverson's bracket). We emphasise that via our definition of $\inf[\cdot]$, the computational very complicated problem of computing the pre-expectation $\langle[p] f$ for a given continuation *f* (finding solutions to second-order fixed point equations to deal with general recursion) has been transformed into a problem suitable to be expressed in constraints, susceptible to automation. In particular this permit us to use the constraint solver implemented within Avanzini et al. [2020] tool eco-imp to reason about such constraints. In brief, the solver resolves conditionals through case analysis, resulting in a set of equivalent constraints over unconditional polynomials, and then makes essential use of Handelman's theorem [Handelman 1988] to turn these into constraints over undetermined coefficients. This in turn enables the use of off-the-shelf SMT solver supporting QF_NRA, in our case Z3.

*Improving upon modularity of the analysis.* As in many denotational semantics, the expectation transformer $\langle \cdot \rangle$ is compositional, which facilitates reasoning. Unfortunately, this compositionality does not directly give rise to modularity of inference. Indeed, in the context of our expected result value analysis, such a modular inference is unsound in general. Recursive procedures and loops cause cyclic dependencies that hinder modularity. In our setting, these cyclic dependencies are
reflected within the constraints generated by infer[·]. Templates assigned to nested loops, or potentially mutually recursive procedures, are thus defined through cyclic constraints. As a result, constraints cannot be solved in isolation, effectively rendering the approach described so far a whole program analysis.

In some cases, however, stratification present in the program can be exploited to run our machinery in an iterative way, thereby greatly improving upon modularity, and in consequence improving upon the performance of the overall implementation. To illustrate, eg. one case where our implementation departs from the definition of infer[·] lies within the analysis of nested loops. Here our implementation follows ideas originally proposed by Avanzini et al. [2020]. In particular, pre-expectations of nested loops can be determined in isolation. In essence, this modular treatment of loops depends on the linear shape of templates, and the linearity of expectations law (see Figure 3). While the approach can in general not be lifted to our setting with procedure calls, for a sizeable class of loops—for example those that do not contain recursive calls—the approach extends to our setting. In this specialised case the constraints imposed by the treatment of the loop are free of unknowns expect those mentioned in the template U over-approximating the expectation of the loop. Thus, the added two constraints can be solved independently, and consequently, a concrete upper-bound rather than a template term U can be substituted for infer[while (B) {C}]s T. Soundness of this specialisation follows by a straightforward adaption of [Avanzini et al. 2020, Theorem 7.5]. These efforts yield that our prototype implementation can analyse the entirety of our 53 benchmarks examples in around 5 minutes on a standard desktop. Accepting a slight deterioration of strength (loosing one example), this benchmark can be handled in less than 5 seconds.

**Evaluation.** To the best our knowledge there is currently no tool providing a fully automated expected value analysis of programs available, regardless whether the considered programming language is imperative or not, admits recursive programs or not. On the other hand, there is ample work on (automated) generation of *quantitative invariants*, see eg. [Bao et al. 2022; Chakarov and Sankaranarayanan 2014; Katoen et al. 2010; McIver and Morgan 2005; Wang et al. 2018], employing a variety of techniques. Furthermore, there is a large body work on *expected cost analysis*, see eg. [Avanzini et al. 2020; Kaminski et al. 2018; Ngo et al. 2018; Wang et al. 2019]. Thus we have chosen examples from these seminal works as basis of the developed benchmark suite. In addition, we have added a number of examples of our own, detailed below. In sum this amounts to a test-suite of 53 examples. The results of these evaluations are given in Tables 1 and 2. In short, we can handle ⁴⁹/₅³ of the benchmark suite, without any recourse to user interaction.

**Selection and evaluation results on examples on invariant generation.** We have chosen challenging examples from [Chakarov and Sankaranarayanan 2014; Katoen et al. 2010], detailed in benchmarks (a) and (b) in Table 1. Our tool can handle all of these examples, with the exception of uniform-dist, where we can only handle a restrictive instance (denoted as uniform-dist-100 in the benchmark). We also note that the expectations employed in [Katoen et al. 2010] are more sophisticated than ours. On the other hand only semi-automation is achieved. Further, we consider examples from [Bao et al. 2022; Wang et al. 2018], establishing fully automated methodologies. Wang et al. employ a templating approach similar to ours, while the very recent work by Bao et al. employs a conceptually highly interesting learning approach. The results are described in Table 1 (c), (d) as well as in Table 2 (d), respectively. To suit to the expressivity of ev-imp, we instantiate the probabilities by constant once in the majority of the benchmarks in (d). Variable probabilities are not (yet) expressible as upper invariants in our prototype.

We can handle all except one example from these benchmarks. In the one example eg (benchmark (c)) that we cannot handle the templates chosen are not precise enough, which is explained by the more liberal construction we use to handle non-determinism in the definition of infer[·].
Wrt. precision, the bounds generated by ev-imp are often as precise as those generated by the tools from Wang et al. [2018] and Bao et al. [2022]. Wrt. to speed, our prototype implementation ev-imp is typically able to handle the benchmark examples in milliseconds. As already mentioned the whole testsuite can be handled in around 5 minutes on a standard desktop.

*Selection and evaluation results on examples on expected costs.* As argued by Kaminski et al. [2018] an expected value analysis is in general unsound for expected cost analysis. However, if the program under consideration is *almost-surely terminating* (AST for short), then an expected cost analysis can be recovered from an expected value analysis by counter instrumentation. In this context, an expected value analysis constitutes a strict extension of an expected cost analysis. Note that all benchmarks considered in Table 2 (f) and (g) are AST (even *positive almost-sure terminating*).

We have incorporated recursive examples from Kaminski et al. [2018], suited to our possibilities, cf. Table 2 (e). While procedure geo can be handled instantly with ev-imp, the growth rate of the (expected) value of faulty_factorial is non-polynomial. We thus suited a variant—dubbed faulty_sum—to our benchmarks replacing the multiplication by a sum, thus featuring polynomial growth but similar algorithmic complexity. In addition, we considered challenging—newly introduced examples from Avanzini et al. [2020]—if expressible in our prototype. Wrt. both benchmarks, we adapted the original expected cost analysis to an expected return value analysis, as proposed above. The evaluation results are given in Table 2 (f) and (g). Remarkably, we can handle both recursive examples from Kaminski et al. and also make significant inroad into challenging non-linear example from [Avanzini et al. 2020]. Unsurprisingly, we cannot handle all the selected benchmarks due to the greater generality of our prototype ev-imp.

For example, the crucial motivating example by Avanzini et al.—*Coupon Collector*—cannot (yet) be handled by our prototype implementation ev-imp, as support for dynamic uniform distributions is lacking. This is left for future work. We can however express (and handle) concrete instances of this benchmark. (We included coupon-10, coupon-50 and coupon-100 as examples.) Wrt. precision, the bounds generated by ev-imp are on the same order of magnitude as those generated by the tools from eco-imp, while wrt. speed the generality of our tool takes its toll.

*Programs considered in this work.* We have already discussed the examples balls, throws, every-5 and benchmark example hire in Section 2 and 4, respectively. Benchmark every constitutes the general case to procedure every, where we do not restrict the number of bins to five, while every-while constitutes an encoding of this problem as loop program. For the moment, ev-imp cannot handle the general encoding of the question how many balls have to be thrown in average until all bins are filled with at least one ball (compare Section 2) due to complexity of the example. The remaining examples constitute recursive variants of standard examples and are given in full in [Avanzini et al. 2023].

# 6 RELATED WORK

Very briefly, we refer to the extensive literature of analysis methods for (non-deterministic) probabilistic programs introduced in the last years. These have been provided in the form of *abstract interpretations* [Chakarov and Sankaranarayanan 2014]; *martingales*, eg., ranking supermartingales [Agrawal et al. 2018; Takisaka et al. 2018; Wang et al. 2019]; or equivalently *Lyapunov ranking functions* [Bournez and Garnier 2005]; *model checking* [Katoen 2016]; *program logics* [Bao et al. 2022; Kaminski and Katoen 2017; Kaminski et al. 2018; Kaminski and Katoen 2015; McIver and Morgan 2005; McIver et al. 2018; Ngo et al. 2018; Wang et al. 2018]; *proof assistants* [Barthe et al. 2009]; *recurrence relations* [Sedgewick and Flajolet 1996]; methods based on *program analysis* [Celiku and McIver 2005; Katoen et al. 2010; Kozen 1985]; or *symbolic inference* [Gehr et al.
<table><caption>Table 1. Automatically Derived Bounds on the Expected Value via our Prototype ev-imp.</caption><thead><tr><th>Program</th><th>Return Value</th><th>Inferred Bound</th><th>Time (sec)</th></tr></thead><tbody><tr><td colspan="4">(a) examples from Katoen et al. [2010]</td></tr><tr><td>biased-coin</td><td>bool</td><td>1/2</td><td>0.029</td></tr><tr><td>binom-update</td><td>x</td><td>1/2 · &lt;N&gt;</td><td>0.012</td></tr><tr><td>uniform-dist</td><td>g</td><td>—</td><td>2.173</td></tr><tr><td>uniform-dist-100</td><td>g</td><td>99</td><td>0.09</td></tr><tr><td colspan="4">(b) examples from Chakarov and Sankaranarayanan [2014]</td></tr><tr><td>mot-ex</td><td>count</td><td>44/3</td><td>0.02</td></tr><tr><td colspan="4">(c) benchmark from Wang et al. [2018]</td></tr><tr><td>2d-walk</td><td>count</td><td>&lt;1 + count&gt;</td><td>0.012</td></tr><tr><td>aggregate-rv</td><td>x</td><td>1/2 · [x ≥ −1 ∧ 499 ≥ i] · (1 + x) + 1/2 · [499 ≥ i] · &lt;x&gt; + [i ≥ 500] · &lt;x&gt;</td><td>0.009</td></tr><tr><td>biased-coin</td><td>x<sub>1</sub></td><td>1/2 · &lt;x<sub>1</sub>&gt; + 1/2 · [x<sub>1</sub> &gt; x<sub>2</sub>]</td><td>0.014</td></tr><tr><td>binom-update</td><td>x</td><td>1/4 · [x ≥ −1 ∧ 99 ≥ n] · (1 + x) + 3/4 · [99 ≥ n] · &lt;x&gt; + [n ≥ 100] · x</td><td>0.009</td></tr><tr><td>coupon5</td><td>count</td><td>[count ≥ −1 ∧ 4 ≥ i] · &lt;1 + count&gt; + [i ≥ 5] · &lt;count&gt;</td><td>0.009</td></tr><tr><td>eg-tail</td><td>x</td><td>&lt;x&gt; + 19/16 · &lt;z&gt; + 7/4</td><td>0.056</td></tr><tr><td>eg</td><td>x</td><td>—</td><td>0.18</td></tr><tr><td>hare-turtle</td><td>h</td><td>[h ≥ −1 ∧ t ≥ h] · 1/22 · &sum;<sup>10</sup><sub>i=0</sub> (h + i) + [h &gt; t] · &lt;h&gt; + [h ≤ t]<sup>6/11</sup> · &lt;h&gt;</td><td>0.010</td></tr><tr><td>hawk-dove</td><td>count</td><td>1</td><td>0.013</td></tr><tr><td>mot-ex</td><td>count</td><td>&lt;1 + count&gt;</td><td>0.009</td></tr><tr><td>recursive</td><td>x</td><td>&lt;x&gt; + 10</td><td>0.03</td></tr><tr><td>uniform-dist</td><td>g</td><td>1/2 · ([g ≥ −1/2 ∧ 9 ≥ n] · (1 + 2 · y) + [9 ≥ n] · 2 · &lt;y&gt;) + [n ≥ 10] · &lt;g&gt;</td><td>0.009</td></tr><tr><td colspan="4">(d) benchmark from Bao et al. [2022]</td></tr><tr><td>biasdir</td><td>x</td><td>N1 − x + &lt;x&gt;</td><td>0.022</td></tr><tr><td>bin0</td><td>x</td><td>1/2 · &lt;n&gt; · &lt;y&gt; + &lt;x&gt;</td><td>0.075</td></tr><tr><td>bin1</td><td>x</td><td>1/2 · &lt;M − n&gt; + &lt;x&gt;</td><td>0.012</td></tr><tr><td>bin2</td><td>x</td><td>1/4 · (&lt;n&gt; + n<sup>2</sup>) + 1/2 · &lt;n&gt; · &lt;y&gt; + &lt;x&gt;</td><td>0.141</td></tr><tr><td>deprv</td><td>z</td><td>&lt;z&gt;</td><td>0.011</td></tr><tr><td>detm</td><td>count</td><td>&lt;11 − x&gt; + &lt;count&gt;</td><td>0.012</td></tr><tr><td>duel</td><td>turn</td><td>&lt;turn&gt; + 2 · &lt;continuing&gt;</td><td>0.036</td></tr><tr><td>fair</td><td>count</td><td>&lt;count&gt; + 2 · &lt;1 − c<sub>2</sub>&gt;</td><td>0.86</td></tr><tr><td>gambler0</td><td>z</td><td>[x ≥ 0 ∧ y ≥ x] · (x · y − x<sup>2</sup>) + &lt;z&gt;</td><td>0.063</td></tr><tr><td>geo0</td><td>z</td><td>&lt;1 − flip&gt; + &lt;z&gt;</td><td>0.014</td></tr></tbody></table>

2016]; and finally type systems [Avanzini et al. 2019; Leutgeb et al. 2022; Vasilenko et al. 2022; Wang et al. 2020].

**Invariant generation.** Generally, invariant generation is more challenging than the expected result value analysis that we study. Still, often ev-imp derives exact bounds, thus establishing invariants. On the other hand, our methodology is applicable in a more general framework, for example to expected cost analysis. Further, our methodology encompasses recursive (imperative) programs, which is—to the best of our knowledge—not the case for any of the approaches on invariant generation (or expected cost analysis, for that matter). Katoen et al. [2010] provide constraint-based methods for the semi-automated generation of linear quantitative invariants, based on sophisticated proof-based methods. The studied examples are highly interesting and have been integrated into our benchmarks (see Section 5). The form of expectations considered can be very expressive. We emphasise, however, that our method is fully automated, while the approach in [Katoen et al. 2010] is only partly automated, in particular nested loops require user-interaction. Related results have been reported by Chakarov and Sankaranarayanan [2014], suitable adapting an
<table><caption>Table 2. Automatically Derived Bounds on the Expected Value via our Prototype ev-imp.</caption><thead><tr><th>Program</th><th>Return Value</th><th>Inferred Bound</th><th>Time (sec)</th></tr></thead><tbody><tr><td colspan="4">(d) benchmark from Bao et al. [2022] (cont'd)</td></tr><tr><td>geo1</td><td>z &lt;1 - flip&gt; + &lt;z&gt;</td><td></td><td>0.014</td></tr><tr><td>geo2</td><td>z &lt;1 - flip&gt; + &lt;z&gt;</td><td></td><td>0.015</td></tr><tr><td>geoar0</td><td>x <sup>3</sup>&frasl;<sub>2</sub>(1 + &lt;1 + y&gt; + &lt;x&gt;)</td><td></td><td>0.09</td></tr><tr><td>linexp</td><td>z &lt;z&gt; + <sup>21</sup>&frasl;<sub>8</sub> &middot; &lt;n&gt;</td><td></td><td>0.014</td></tr><tr><td>mart</td><td>rounds &lt;rounds&gt; + 3 &middot; &lt;b&gt;</td><td></td><td>0.012</td></tr><tr><td>prinsys</td><td>bool 1</td><td></td><td>0.013</td></tr><tr><td>revbin</td><td>z &lt;z&gt; + 2 &middot; &lt;x&gt;</td><td></td><td>0.011</td></tr><tr><td>sum0</td><td>x <sup>1</sup>&frasl;<sub>4</sub> &middot; (&lt;n&gt; + n<sup>2</sup>) + &lt;x&gt;</td><td></td><td>0.034</td></tr><tr><td colspan="4">(e) examples from [Kaminski et al. 2018]</td></tr><tr><td>faulty_sum</td><td>x <sup>3</sup>&frasl;<sub>7</sub> &middot; &lt;x&gt;<sup>2</sup> + <sup>4</sup>&frasl;<sub>7</sub> &middot; &lt;x&gt; + 1</td><td></td><td>0.04</td></tr><tr><td>geo</td><td>0 0</td><td></td><td>0.009</td></tr><tr><td colspan="4">(f) benchmark from [Avanzini et al. 2020]</td></tr><tr><td>bridge</td><td>count [b &ge; x &and; x &ge; a] &middot; (-a &middot; b + a &middot; x + b &middot; x - x<sup>2</sup>)</td><td></td><td>0.067</td></tr><tr><td>coupon-10</td><td>count 110</td><td></td><td>0.145</td></tr><tr><td>coupon-50</td><td>count 2550</td><td></td><td>2.691</td></tr><tr><td>coupon-100</td><td>count 10100</td><td></td><td>10.63</td></tr><tr><td>nest-1</td><td>count 4 &middot; &lt;n&gt;</td><td></td><td>0.014</td></tr><tr><td>nest-2</td><td>count —</td><td></td><td>100.00</td></tr><tr><td>trader-5</td><td>&lt;price&gt; <sup>20</sup>&frasl;<sub>3</sub> &middot; [price &ge; -1] &middot; (1 + 2 &middot; price + price<sup>2</sup>) + <sup>35</sup>&frasl;<sub>6</sub> &middot; &lt;1 + price&gt;</td><td></td><td></td></tr><tr><td colspan="4">(g) examples from this work</td></tr><tr><td>balls</td><td>b <sup>1</sup>&frasl;<sub>5</sub> &middot; &lt;n&gt;</td><td></td><td>0.012</td></tr><tr><td>throws</td><td>throws 5</td><td></td><td>0.012</td></tr><tr><td>every</td><td>number —</td><td></td><td>100.0</td></tr><tr><td>every-5</td><td>number 20</td><td></td><td>0.017</td></tr><tr><td>every-while</td><td>number 25</td><td></td><td>0.973</td></tr><tr><td>double_recursive</td><td>y 0</td><td></td><td>0.011</td></tr><tr><td>hire</td><td>hire &lt;n&gt;</td><td></td><td>0.504</td></tr><tr><td>rdwalk</td><td>n 2 &middot; &lt;n&gt;</td><td></td><td>0.016</td></tr><tr><td>rec1</td><td>n <sup>1</sup>&frasl;<sub>2</sub> &middot; (1 + &lt;n&gt;)</td><td></td><td>0.014</td></tr></tbody></table>

abstract interpretation framework to the notion of invariant generation. Their motivating example can be handled by ev-imp fully automatically, establishing a slightly worse constant bound than the exact bound (see Section 5). In contrast to [Chakarov and Sankaranarayanan 2014; Katoen et al. 2010], Wang et al. [2018] and the very recent Bao et al. [2022] provide fully automated (linear) invariant generation methodologies. Wang et al. [2018] provide a compelling algebraic framework for the program analysis of probabilistic programs. Apart from the here relevant linear invariant generation, interprocedural Bayesian inference analysis and the Markov decision problem are conducted. Wrt. linear invariant generation, we have considered all provided examples in our benchmarks and obtained comparable results, while improving the speed of the analysis by a magnitude (in comparison to the analysis times reported in [Wang et al. 2018]), cf. Section 5. Automation of the method developed in [Wang et al. 2018] is based, like ours, on a template approach. To overcome the dependency on templates, Bao et al. [2022] have developed a highly interesting learning approach. Their approach is data-driven. Apart from invariants, Bao et al. also consider *sub-invariants* which are dual to our upper invariants, establishing lower bounds on the pre-expectations. In both cases, however, analysis times are high. We have incorporated
the benchmark examples from [Bao et al. 2022] into our test suites. In all cases our tool ev-imp provides precise invariants, cf. Section 5. This is remarkable, as reported in [Bao et al. 2022], their prototype implementation Exist handles only 12/18 of their benchmark suite fully automatically.

*Expected cost analysis.* Kaminski et al. [2018] establishes an expected cost analysis of recursive programs and we have suited the corresponding two example to our benchmark suite. Both examples can be handled fully automatically. Technically our development of recursive programs constitutes an extension as our language (and methodology) admits local variables, formal parameters and (unrestricted) return values. This allows a more natural representation of programs (see Section 3). Still the definition of our expectation transformer $\text{et}[[p]]$ is closely related to the corresponding definition in [Kaminski et al. 2018, Chapter 7]. We emphasise, however, that in [Kaminski et al. 2018] automation is discussed only superficially. Avanzini et al. [2020], on the other hand, take automation very seriously. As mentioned above, we have taken inspiration from their work in the modular analysis of recursion-free programs. Despite our efforts, however, our prototype implementation ev-imp lacks scalability and speed in comparison to their tool eco-imp. Wrt. precision, however, we often derive the same bounds on expected costs. Further, and as argued, our methodology focusing on expected value analysis is more general and our implementation incorporated the highly non-trivial handling of recursive programs. Technically the closest comparison is to work on amortised cost analysis of functional languages (eg. [Leutgeb et al. 2022; Wang et al. 2020]), as resource parametricity [Hoffmann 2011] has been studied in this context. In comparison to Wang et al. [2020], our tool ev-imp cannot infer higher-order moments but concerning expectations seems to have better support for recursion. For instance, we were not able to reproduce an expected cost analysis of the balls procedure in their RaML prototype. Similarly, ev-imp is no match to Leutgeb et al. [2022]'s ATLAS, when it comes to the precise analysis of (probabilistic) data structures. Remarkably, however, the balls benchmark can only be expressed convolutedly in their language. (Due to the lack of support for general, inductive, data structures.) Further, their automated analysis does not derive an optimal bound.

*Expected value analysis.* Finally, we briefly remark on very recent and partly motivating work by Vasilenko et al. [2022]. In [Vasilenko et al. 2022] a refinement type system—Liquid Haskell, cf. [Handley et al. 2020; Vazou 2016]—is updated, to reason about relational properties of probabilistic computations. Eg. Vasilenko et al. [2022] provide a semi-automated proof that for (a variant of) procedure balls, cf. Listing 1(a), the expected return value of b is given as $p \cdot n$.

# 7 CONCLUSION

We established a fully automated expected result value analysis for probabilistic programs in the presence of natural programming constructs, in particular recursion. Our analysis is in the form of an *expectation transformer* $\text{et}[[\cdot]]$ and its syntactic representation $\text{infer}[\cdot]$. As argued, automated inference of upper invariants is challenging for PWhile, due to the presence of recursion. We have overcome these challenges and implemented the established methodology in our novel prototype implementation ev-imp. In future work, we aim to incorporate (i) *more program features* and, like eg. support for *dynamic uniform distributions*; (ii) improve the *constraint solving* capabilities of our prototype implementation, to handle the analysis of further natural probabilistic programs and data structures fully automatically.

# ACKNOWLEDGMENTS

We would like to thank the anonymous reviewers for their work and invaluable suggestions, which greatly improved our presentation. This work is partly supported by the INRIA Associate Team TC(Pro)³ and by the ANR Project PPS: "Probabilistic Program Semantics".
# 8 AVAILABILITY OF DATA AND SOFTWARE

Our prototype implementation ev-imp is publicly available via the following Zenodo link: https://doi.org/10.5281/zenodo.7706691. The artifact is given as a Docker image running a minimal Linux distribution, containing the sources in directory /evimp, as well as pre-compiled executables. A second, larger, image contains the source distribution as well as the complete toolchain to compile ev-imp. We have successfully tested the images under Linux, MacOS and Windows, as long as these were non-ARM architectures. For ARM64 architectures, the artifact additionally provides instructions to install from source. Detailed instructions on the use of ev-imp are provided as well.

# REFERENCES

S. Agrawal, K. Chatterjee, and P. Novotný. 2018. Lexicographic Ranking Supermartingales: An Efficient Approach to Termination of Probabilistic Programs. PACMPL 2, POPL (2018), 34:1–34:32. https://doi.org/10.1145/3385412.3386002

M. Avanzini, G. Barthe, and U. Dal Lago. 2021. On Continuation-Passing Transformations and Expected Cost Analysis. *PACM on Programming Languages* 5, ICFP (2021), 1–30. https://doi.org/10.1145/3473592

M. Avanzini, U. Dal Lago, and A. Ghyselen. 2019. Type-Based Complexity Analysis of Probabilistic Functional Programs. In *Proc. of 34th LICS*. IEEE, 1–13. https://doi.org/10.1109/LICS.2019.8785725

Martin Avanzini, Georg Moser, and Michael Schaper. 2020. A modular cost analysis for probabilistic programs. *Proc. ACM Program. Lang.* 4, OOPSLA (2020), 172:1–172:30. https://doi.org/10.1145/3428240

Martin Avanzini, Georg Moser, and Michael Schaper. 2023. Automated Expected Value Analysis of Recursive Programs. arXiv:2304.01284 [cs.PL]

Jialu Bao, Nitesh Trivedi, Drashti Pathak, Justin Hsu, and Subhajit Roy. 2022. Data-Driven Invariant Learning for Probabilistic Programs. In *Proc. of 34th CAV (LNCS, Vol. 13371)*, 33–54. https://doi.org/10.1007/978-3-031-13185-1_3

G. Barthe, B. Grégoire, and S. Z. Béguelin. 2009. Formal Certification of Code-based Cryptographic Proofs. In *Proc. of 36th POPL*. ACM, 90–101. https://doi.org/10.1145/1480881.1480894

Kevin Batz, Benjamin Lucien Kaminski, Joost-Pieter Katoen, Christoph Matheja, and Thomas Noll. 2019. Quantitative separation logic: a logic for reasoning about probabilistic pointer programs. *PACM on Programming Languages* 3, POPL (2019), 34:1–34:29. https://doi.org/10.1145/3290347

O. Bournez and F. Garnier. 2005. Proving Positive Almost-Sure Termination. In *Proc. of 16th RTA (LNCS, Vol. 3467)*. Springer, 323–337. https://doi.org/10.1142/S0129054112400588

O. Celiku and A. McIver. 2005. Compositional Specification and Analysis of Cost-Based Properties in Probabilistic Programs. In *Proc. of FM 2005 (LNCS, Vol. 3582)*. Springer, 107–122. https://doi.org/10.1007/11526841_9

A. Chakarov and S. Sankaranarayanan. 2014. Expectation Invariants for Probabilistic Program Loops as Fixed Points. In *Proc. of 21st SAS (LNCS)*. Springer, 85–100. https://doi.org/10.1007/978-3-319-10936-7_6

E. Contejean, C. Marché, A.-P. Tomás, and X. Urbain. 2005. Mechanically Proving Termination Using Polynomial Interpretations. *JAR* 34, 4 (2005), 325–363. https://doi.org/10.1007/s10817-005-9022-x

Thomas H. Cormen, Charles E. Leiserson, Ronald L. Rivest, and Clifford Stein. 2009. *Introduction To Algorithms* (3rd ed.). MIT Press. http://mitpress.mit.edu/books/introduction-algorithms

Leonardo Mendonça de Moura and Nikolaj Bjørner. 2008. Z3: An Efficient SMT Solver. In *Proc. 14th TACAS (LNCS, Vol. 4963)*. 337–340.

E. W. Dijkstra. 1975. Guarded Commands, Nondeterminacy and Formal Derivation of Programs. ACM 18, 8 (1975), 453–457.

E. W. Dijkstra. 1976. *A Discipline of Programming*. Prentice-Hall.

Dino Distefano, Manuel Fähndrich, Francesco Logozzo, and Peter W. O'Hearn. 2019. Scaling static analyses at Facebook. *Commun. ACM* 62, 8 (2019), 62–70. https://doi.org/10.1145/3338112

Manuel Eberl, Max W. Haslbeck, and Tobias Nipkow. 2020. Verified Analysis of Random Binary Tree Structures. *J. Autom. Reason*. 64, 5 (2020), 879–910. https://doi.org/10.1007/s10817-020-09545-0

Daniel P. Friedman and Mitchell Wand. 2008. *Essentials of Programming Languages* (3rd ed.). The MIT Press.

C. Fuhs, J. Giesl, A. Middeldorp, P. Schneider-Kamp, R. Thiemann, and H. Zankl. 2007. SAT Solving for Termination Analysis with Polynomial Interpretations. In *Proc. of 10th SAT (LNCS, Vol. 4501)*. Springer, 340–354. https://doi.org/10.1007/978-3-540-72788-0_33

T. Gehr, S. Misailovic, and M. Vechev. 2016. PSI: Exact Symbolic Inference for Probabilistic Programs. In *Proc. of 28th CAV (LNCS, Vol. 9779)*. Springer, 62–83. https://doi.org/10.1007/978-3-319-41528-4_4

Ronald L. Graham, Donald E. Knuth, and Oren Patashnik. 1994. *Concrete Mathematics: A Foundation for Computer Science* (2nd ed.). Addison-Wesley. https://www-cs-faculty.stanford.edu/%7Eknuth/gkp.html

Friedrich Gretz, Joost-Pieter Katoen, and Annabelle McIver. 2014. Operational versus weakest pre-expectation semantics for the probabilistic guarded command language. *Performance Evaluation* 73 (2014), 110–132.
D. Handelman. 1988. Representing Polynomials by Positive Linear Functions on Compact Convex Polyhedra. *PJM* 132, 1 (1988), 35–62. https://doi.org/10.2140/pjm.1988.132.35

Martin A. T. Handley, Niki Vazou, and Graham Hutton. 2020. Liquidate your assets: reasoning about resource usage in liquid Haskell. *PACMPL* 4, POPL (2020), 24:1–24:27. https://doi.org/10.1145/3371092

C. A. R. Hoare. 1969. An Axiomatic Basis for Computer Programming. *Commun. ACM* 12, 10 (1969), 576–580. https://doi.org/10.1145/363235.363259

Jan Hoffmann. 2011. *Types with potential: polynomial resource bounds via automatic amortized analysis*. Ph. D. Dissertation. Ludwig Maximilians University Munich. http://edoc.ub.uni-muenchen.de/13955/

B. L. Kaminski and J.-P. Katoen. 2017. A Weakest Pre-expectation Semantics for Mixed-sign Expectations. In *Proc. of 32nd LICS*. IEEE, 1–12. https://doi.org/10.1109/LICS.2017.8005153

B. Lucien Kaminski, J.-P. Katoen, C. Matheja, and F. Olmedo. 2016. Weakest Precondition Reasoning for Expected Run-Times of Probabilistic Programs. In *Proc. of 25th ESOP (LNCS, Vol. 9632)*. Springer, 364–389. https://doi.org/10.1007/978-3-662-49498-1_15

B. L. Kaminski, J.-P. Katoen, C. Matheja, and F. Olmedo. 2018. Weakest Precondition Reasoning for Expected Runtimes of Randomized Algorithms. *JACM* 65, 5 (2018), 30:1–30:68. https://doi.org/10.1145/3208102

B. L. Kaminski and J.-P. Katoen. 2015. On the Hardness of Almost-Sure Termination. In *MFCS 2015, Part I (LNCS)*. Springer, 307–318. https://doi.org/10.1007/978-3-662-48057-1_24

J.-P. Katoen. 2016. The Probabilistic Model Checking Landscape. In *Proc. of 31st LICS*. ACM, 31–45. https://doi.org/10.1145/2933575.2934574

J.-P. Katoen, A. McIver, L. Meinicke, and C.C. Morgan. 2010. Linear-Invariant Generation for Probabilistic Programs: Automated Support for Proof-Based Methods. In *Proc. of 17th SAS (LNCS, Vol. 6337)*. Springer, 390–406. https://doi.org/10.1007/978-3-642-15769-1_24

Thomas Kleymann. 1999. Hoare Logic and Auxiliary Variables. *Formal Aspects Comput.* 11, 5 (1999), 541–566. https://doi.org/10.1007/s001650050057

D. Kozen. 1985. A Probabilistic PDL. *JCSC* 30, 2 (1985), 162–178. https://doi.org/10.1016/0022-0000(85)90012-1

Lorenz Leutgeb, Georg Moser, and Florian Zuleger. 2021. ATLAS: Automated Amortised Complexity Analysis of Self-adjusting Data Structures. In *Proc. of 33rd CAV (LNCS, Vol. 12760)*. 99–122. https://doi.org/10.1007/978-3-030-81688-9_5

Lorenz Leutgeb, Georg Moser, and Florian Zuleger. 2022. Automated Expected Amortised Cost Analysis of Probabilistic Data Structures. In *Proc. of 34th CAV (LNCS, Vol. 13372)*. 70–91. https://doi.org/10.1007/978-3-031-13188-2_4

Annabelle McIver and Carroll Morgan. 2005. *Abstraction, refinement and proof for probabilistic systems*. Springer Science & Business Media.

A. McIver, C. Morgan, B. L. Kaminski, and J-P Katoen. 2018. A New Proof Rule for Almost-sure Termination. *PACMPL* 2, POPL (2018), 33:1–33:28. https://doi.org/10.1145/3158121

M. Mitzenmacher and E. Upfal. 2005. *Probability and Computing: Randomized Algorithms and Probabilistic Analysis*. Cambridge University Press. https://doi.org/10.1017/CBO9780511813603

N. C. Ngo, Q. Carbonneaux, and J. Hoffmann. 2018. Bounded Expectations: Resource Analysis for Probabilistic Programs. In *Proc. of 39th PLDI*. ACM, 496–512. https://doi.org/10.1145/3296979.3192394

Flemming Nielson, Hanne Riis Nielson, and Chris Hankin. 1999. *Principles of program analysis*. Springer. https://doi.org/10.1007/978-3-662-03811-6

Hanne Riis Nielson. 1987. A Hoare-Like Proof System for Analysing the Computation Time of Programs. *Sci. Comput. Program.* 9, 2 (1987), 107–136. https://doi.org/10.1016/0167-6423(87)90029-3

A. Podelski and A. Rybalchenko. 2004. A Complete Method for the Synthesis of Linear Ranking Functions. In *Proc. of 5th VMCAI (LNCS, Vol. 2937)*. Springer, 239–251. https://doi.org/10.1007/978-3-540-24622-0_20

E. Schlechter. 1996. *Handbook of Analysis and Its Foundations*. Elsevier.

Roberto Sebastiani and Patrick Trentin. 2020. OptiMathSAT: A Tool for Optimization Modulo Theories. *J. Autom. Reason.* 64, 3 (2020), 423–460. https://doi.org/10.1007/s10817-018-09508-6

R. Sedgewick and P. Flajolet. 1996. *An Introduction to the Analysis of Algorithms*. Addison-Wesley-Longman. https://doi.org/10.1142/10875

M. Sinn, F. Zuleger, and H. Veith. 2014. A Simple and Scalable Static Analysis for Bound Analysis and Amortized Complexity Analysis. In *Proc. of 26th CAV (LNCS, Vol. 8559)*. Springer, Heidelberg, DE, 745–761.

M. Sinn, F. Zuleger, and H. Veith. 2016. A Simple and Scalable Static Analysis for Bound Analysis and Amortized Complexity Analysis. In *Proc. of 26th CAV (LNCS, Vol. 8559)*. Springer, 745–761. https://doi.org/10.1007/978-3-319-08867-9_50

T. Takisaka, Y. Oyabu, N. Urabe, and I. Hasuo. 2018. Ranking and Repulsing Smartartingales for Reachability in Probabilistic Programs. In *Proc. of 16th ATVA (LNCS, Vol. 11138)*. Springer, 476–493. https://doi.org/10.1007/978-3-030-01090-4_28

Elizaveta Vasilenko, Niki Vazou, and Gilles Barthe. 2022. Safe couplings: coupled refinement types. *Proc. ACM Program. Lang.* 6, ICFP (2022), 596–624. https://doi.org/10.1145/3547643
Niki Vazou. 2016. *Liquid Haskell: Haskell as a Theorem Prover*. Ph. D. Dissertation. University of California, San Diego, USA. http://www.escholarship.org/uc/item/8dm057ws

D. Wang, J. Hoffmann, and T. W. Reps. 2018. PMAF: An Algebraic Framework for Static Analysis of Probabilistic Programs. In *Proc. of 39th PLDI*. ACM, 513–528. https://doi.org/10.1145/3192366.3192408

Di Wang, David M. Kahn, and Jan Hoffmann. 2020. Raising expectations: automating expected cost analysis with types. *PACM on Programming Languages* 4, ICFP (2020), 110:1–110:31. https://doi.org/10.1145/3408992

P. Wang, H. Fu, A. K. Goharshady, K. Chatterjee, X. Qin, and W. Shi. 2019. Cost Analysis of Nondeterministic Probabilistic Programs. In *Proc. of 40th PLDI*. ACM, 204–220. https://doi.org/10.1145/3314221.3314581

Wolfgang Wechler. 1992. *Universal Algebra for Computer Scientists*. EATCS Monographs on Theoretical Computer Science, Vol. 25. Springer. https://doi.org/10.1007/978-3-642-76771-5

G. Winskel. 1993. *The Formal Semantics of Programming Languages*. MIT Press. https://doi.org/10.7551/mitpress/3054.003.0004
# A MATHEMATICAL BACKGROUND

**PROPOSITION A.1 (FUNCTION LIFTING OF ω-CPOs [WINSKEL 1993, SECTION 8.3.3])**. Let $(D, \sqsubseteq)$ be an ω-CPO. Then $(A \to D, \sqsubseteq)$ where $\sqsubseteq$ extends $\sqsubseteq$ point-wise forms an ω-CPO, with the supremum on $A \to D$ given point-wise. If $\sqsubseteq$ has least and greatest elements $\perp$ and $\top$, then $\perp$ and $\top$ are the least and greatest elements of $\sqsubseteq$, respectively.

**THEOREM A.2 (KLEENE'S FIXED-POINT THEOREM FOR ω-CPOs, [WINSKEL 1993, THEOREM 5.11])**. Let $(D, \sqsubseteq)$ be a ω-CPO with least element $\perp$. Let $\chi: D \to D$ be continuous (thus monotone). Then $\chi$ has a least fixed-point given by

$$ \text{lfp}(\chi) = \sup_{n \in \mathbb{N}} \chi^n(\perp). $$

**LEMMA A.3 (CONTINUITY OF EXPECTATION)**. $\mathbb{E}_\mu \sup_{n \in \mathbb{N}} f_n = \sup_{n \in \mathbb{N}} \mathbb{E}_\mu f_n$ for all ω-chains $(f_n)_{n \in \mathbb{N}}$.

**PROOF.** This is the discrete version of Lebesgue's Monotone Convergence Theorem [Schlechter 1996, Theorem 21.38]. □

# B OMITTED PROOFS

**PROPOSITION B.1 (EXPECTATION TRANFORMER LAWS).** For any program P, any procedure environment $\eta$, any command C and any expectations $f, f_1, f_2, g, g_1, g_2, \dots$ the laws in Figure 3 hold.

**PROOF.** The proofs follow the pattern of the proof of continuity of the expected cost transformer in [Avanzini et al. 2020, Lemma 6.2]. □

We define finite approximations of procedure environments et $[[P]]^{(i)}$ inductively, so that et $[[P]]^{(0)} := \lambda p. \lambda f \vec{v} \sigma. 0$ and et $[[P]]^{(i+1)} := \lambda p. \text{et}[[p]]^{(i)}$, where $p \in P$. In this way et $[[P]]^{(0)}$ represents the poorest approximation by the constant zero function, while approximations are refined iteratively.

**PROPOSITION B.2 (FINITE APPROXIMATIONS OF PROCEDURE ENVIRONMENTS).** For any program P, we have et $[[P]] = \sup_{i \ge 0} \text{et}[[P]]^{(i)}$.

**PROOF.** Direct consequence of the continuity of et $[[p]]$ and Knaster-Tarski fixed-point theorem. □

**THEOREM 4.1 (SOUNDNESS THEOREM).** If for all $p \in P$ the constraint (1) is fulfilled and all side-conditions in Figure 4 are met, then for all $p \in P$, et $[[p]] \le [[\text{infer}[p]]]$, that is, the inference algorithm is sound.

**PROOF.** Let $S, T \in \text{Term} V Z$ and $\theta: Z \to \mathbb{Z}$. Then we prove for all commands C that

$$ \text{et}[[C]]_{[S\theta]}^{\eta} [T\theta] \le \text{[}(\text{infer}[C]_S T)\theta\text{]}. \quad (2) $$

Let $s := [[S\theta]]$ and $t := [[T\theta]]$. In order to prove (2), we prove the following, for all $i \in \mathbb{N}$:

$$ \text{et}[[C]]_s^{\eta_i} t \le \text{[}(\text{infer}[C]_S T)\theta\text{]}. $$

where (i) $\eta_0 p := \lambda p. 0$; (ii) $\eta_{i+1} p := \text{et}[[p]]^{\eta_i}$; and $\eta := \sup_{i \ge 0} \eta_i$.

In proof, we elide the logical context $\Gamma_p$, employed in the definition of $\text{infer}[p]$ for notational convenience. As the essence of this context information is that logical variables are always instantiated non-negatively, this can be guaranteed globally.

We proceed by main induction on $i$ and side induction on C, where we focus on the (main) step case, as the case for $i = 0$ is similar, but simpler. Thus let $i > 0$ and we proceed by case induction on C, considering the most interesting cases, only.
- **CASE skip.** By unfolding of definitions, we easily obtain

$$ \text{et}[[\text{skip}]_s^{\eta_{i+1}} t = [T \theta] = [(\text{infer}[\text{skip}]_s T) \theta], $$

which concludes the case.

- **CASE x :≈ p(\vec{E}).** Suppose $p \in \text{Prog}$. By unfolding of definitions, we obtain for $\sigma \in \text{MemV}$

$$
\begin{align*}
\text{et}[[x :\approx p(\vec{E})]_s^{\eta_{i+1}} t \sigma &= \eta_{i+1} p \underbrace{(\lambda v \tau. t (\sigma_l \triangleright \tau_g)[x \mapsto v])}_{:=f} ([\vec{E}] \sigma) \sigma_g \\
&= \text{et}[[p]]^{\eta_i} f ([\vec{E}] \sigma) \sigma_g \\
&= \text{et}[[Bdy_p]]_f^{\eta_i} (\lambda \tau. t (\sigma_l \triangleright \tau_g)[x \mapsto 0] \sigma_g \triangleright \{\text{Args}_p \mapsto [\vec{E}] \sigma\}) \\
&\le \text{et}[[Bdy_p]]_{K_p \theta}^{\eta_i} [K_p [\mathbf{r} \mapsto 0] \theta] \sigma_g \triangleright \{\text{Args}_p \mapsto [\vec{E}] \sigma\} \\
&\le [(\text{infer}[Bdy_p]_{K_p} K_p [\mathbf{r} \mapsto 0]) \theta] \sigma_g \triangleright \{\text{Args}_p \mapsto [\vec{E}] \sigma\} \\
&= [(\text{infer}[Bdy_p]_{K_p} K_p [\mathbf{r} \mapsto 0]) [\text{Args}_p \mapsto \vec{a}_p] \theta] \sigma_g \triangleright \{\vec{a}_p \mapsto [\vec{E}] \sigma\} \\
&= [(\text{infer}[p] K_p) \theta] \sigma_g \triangleright \{\vec{a}_p \mapsto [\vec{E}] \sigma\} \\
&\le [H_p \theta] \sigma_g \triangleright \{\vec{a}_p \mapsto [\vec{E}] \sigma\} \\
&= [H_p [\vec{a}_p \mapsto \vec{E}] \theta] \sigma = [(\text{infer}[x :\approx p(\vec{E})]_s T \theta) \sigma.
\end{align*}
$$

Here, we have used in conjunction with reduction to definitions (i) two instances of the assumed side-condition $\vdash T[x \mapsto \mathbf{r}] \le K_q \theta$ in line three; (ii) together with monotonicity of the expectation transformer $\text{et}[[C]]$, cf. Figure 3; (iii) induction hypothesis in line four; and finally (iv) in the pre-ultimate line, the main constraint on the soundness of the inference mechanisms (1). This concludes the case.

- **CASE return(E).** By unfolding of definitions, we obtain for $\sigma \in \text{MemV}$

$$
\begin{align*}
\text{et}[[\text{return}(E)]_s^{\eta_{i+1}} t \sigma &= [S \theta] ([E] \sigma) \sigma_g \\
&= [S[\mathbf{r} \mapsto E] \theta] \sigma \\
&= [(\text{infer}[\text{return}(E)]_s T) \theta] \sigma.
\end{align*}
$$

This concludes the case.

- **CASE var x \leftarrow E \text{ in } \{C\}.** By unfolding of definitions in conjunction with application of the induction hypothesis, we obtain for $\sigma \in \text{MemV}$

$$
\begin{align*}
\text{et}[[\text{var } x \leftarrow E \text{ in } \{C\}]_s^{\eta_{i+1}} t \sigma &= \text{et}[[C]]_s^{\eta_{i+1}} t \sigma[x \mapsto [E] \sigma] \\
&\le [(\text{infer}[C]_s T) \theta] \sigma[x \mapsto [E] \sigma] \\
&= [(\text{infer}[C]_s T)[x \mapsto E] \theta] \sigma \\
&= [(\text{infer}[\text{var } x \leftarrow E \text{ in } \{C\}]_s T) \theta] \sigma
\end{align*}
$$

In the third line, we employ that due to the variable condition the local variable $x$ is distinct from all (global) variables in the domain of $\sigma$. This concludes the case.
<figure><img src="image_4.png" /><figcaption>Fig. 6. Additional Benchmark Examples</figcaption></figure>

- **CASE C;D.** By unfolding of definitions in conjunction with two applications of the induction hypothesis, we obtain for $\sigma \in \text{Mem V}$

$$ 
\begin{align*}
\text{et}[[C;D]]_S^{\eta_{i+1}} t \sigma &= \text{et}[[C]]_S^{\eta_{i+1}} (\text{et}[[D]]_S^{\eta_{i+1}} t) \sigma \\
&\leq [[(\text{infer}[C]_S T) \theta]] [[(\text{infer}[D]_S T) \theta]] \sigma \\
&= [[\text{infer}[C]_S (\text{infer}[D]_S T) \theta]] \sigma \\
&= [[(\text{infer}[C;D]_S T) \theta]] \sigma
\end{align*}
$$

Apart from applications of the induction hypothesis, in line two, we have employed monotonicity, cf. Figure 3. This concludes the case.

- **CASE if (B) {C} else {D}.** By unfolding of definitions in conjunction with two applications of the induction hypothesis, we obtain for $\sigma \in \text{Mem V}$

$$ 
\begin{align*}
\text{et}[[\text{if} (B) \{C\} \text{else } \{D\}]]_S^{\eta_{i+1}} t \sigma &= [[[B]] \sigma] \cdot \text{et}[[C]]_S^{\eta_{i+1}} t \sigma + [[[\neg B]] \sigma] \cdot \text{et}[[D]]_S^{\eta_{i+1}} t \sigma \\
&\leq [[[B]] \sigma] \cdot [[(\text{infer}[C]_S T) \theta]] + [[[\neg B]] \sigma] \cdot [[(\text{infer}[D]_S T) \theta]] \\
&= [[([B] \cdot \text{infer}[C]_S T) + [\neg B] \cdot \text{infer}[D]_S T) \theta]] \sigma \\
&= [[(\text{infer}[\text{if} (B) \{C\} \text{else } \{D\}])_S T) \theta]] \sigma
\end{align*}
$$

This concludes the case.

- **CASE while (B) {D}.** In this case, $(\text{infer}[C]_S T) \theta = U \theta$ for some term U, satisfying (i) $B \vdash \text{infer}[C]_S T \leq U$ and (ii) $\neg B \vdash T \leq U$. By induction hypothesis and monotonicity of et$[[D]]_S^\eta$, this says nothing more than that $[[U \theta]]$ is a loop-invariant for the while loop, wrt. $[[T \theta]]$ (see Figure 3).

- **CASE C <> D.** In this case, $(\text{infer}[C]_S T) \theta = U \theta$ for some term U, satisfying (i) $\vdash \text{infer}[C]_S N \leq U$ and (ii) $\vdash \text{infer}[D]_S N \leq U$. The case follows by of two applications of induction hypothesis. $\square$

## ADDITIONAL BENCHMARKS

We detail in Figures 6 and 7 the benchmarks from this paper, whose evaluation results are given in Table 2 (h).

Received 2022-11-10; accepted 2023-03-31
<table><caption>Fig. 7. Additional benchmark examples (cont'd)</caption><tr><td>def f(n):<br/>var m<br/>if (n > 0) {<br/>m := f(n - 1)<br/>};<br/>m := m + Bernoulli(1/2);<br/>return m</td><td>def every(bins):<br/>var d, number, k<br/>k := 1<br/>while (k ≤ 5) {<br/>if (Bernoulli(5-k+1/5)) {k := k + 1}<br/>number := number + 1<br/>};<br/>return number</td></tr><tr><td>(a) rec1 benchmark example.</td><td>(b) every-while benchmark example.</td></tr></table>