Hongfei Fu

Kafshdar Amir

Goharshady

Automated Concentration Bound Analysis for Probabilistic Recurrence Relations Beating Karp's Cookbook a er Three Decades

Probabilistic recurrence relations (PRRs) are a standard formalism to analyze the runtime of randomized algorithms. In this work, we consider the classical problem of obtaining concentration bounds on the runtime of randomized algorithms modeled as PRRs. More speci cally, given a time limit 𝜅, we study the probability Pr[𝑇 ≥ 𝜅] that the runtime 𝑇 of the algorithm reaches or exceeds the time limit. The key research question in concentration bound analysis is to nd an upper-bound function upper(𝑛) ≥ Pr[𝑇 ≥ 𝜅] such that upper(𝑛) converges rapidly to zero as 𝑛 tends to in nity. The classical approach to this problem is the "cookbook" method by Karp (J. ACM 1994). For the past three decades, nding an approach to improve Karp's classical bounds has been a major open problem.

In this work, we rst demonstrate a method to obtain such upper-bound functions from a speci c family of supermartingale functions through the well-known Optional Stopping Theorem (OST). We then provide automated algorithmic approaches to synthesize such supermartingale functions in exponential and polynomial form. Our algorithms are based on the Taylor expansion and tight approximations of the expectation of moments that involve logarithmic terms. Compared with Karp's method, our algorithmic approaches nd signi cantly tighter concentration bounds for several classical randomized algorithms such as S , S and D C . Moreover, our algorithmic approaches can also nd concentration bounds for several randomized algorithms in computational geometry which adopt techniques such as randomized incremental construction. To the best of our knowledge, no previous method was applicable in these cases. Finally, we present experimental results demonstrating the practical versatility and e ciency of our approach. * Due to di erent academic cultural standards, Mainland Chinese co-authors are listed before other co-authors (Austria, Hong Kong). No comparison of contribution between the two groups is implied.

INTRODUCTION

Randomized algorithms. Randomized algorithms [START_REF] Karp | An introduction to randomized algorithms[END_REF][START_REF] Motwani | Randomized Algorithms[END_REF] are a class of algorithms that adopt randomization techniques, such as randomized pivoting [Hoare 1961a,b] and randomized incremental construction [START_REF] Kenneth | Four Results on Randomized Incremental Constructions[END_REF]]. As a central topic in theoretical computer science, many randomized algorithms have been devised in various application domains, such as sorting and selection [Hoare 1961a,b], computational geometry [START_REF] Kenneth | Four Results on Randomized Incremental Constructions[END_REF], distributed systems [Motwani and Raghavan 1995, Chapter 3][ Kleinberg and Tardos 2006, Chapter 13], graph theory [START_REF] Karger | Global Min-cuts in RNC, and Other Rami cations of a Simple Min-Cut Algorithm[END_REF]], and data streaming [START_REF] Robert | Counting Large Numbers of Events in Small Registers[END_REF]. Their advantage in comparison with deterministic algorithms is two-fold: rst, randomized algorithms typically have simpler design while their expected performance often matches or outperforms the deterministic counterparts (e.g. S [Hoare 1961a] vs. M S [Cormen et al. 2009, Chapter 2]); second, for various problems, especially when it comes to distributed settings, deterministic algorithms are impossible and the only solutions are through randomization [Motwani and Raghavan 1995, Chapter 3] [START_REF] Kleinberg | Algorithm design[END_REF]Tardos 2006, Chapter 13][Morris 1978].

Probabilistic recurrence relations. Probabilistic Recurrence Relations (PRRs) are the probabilistic extension of classical recurrence relations and are widely used in the runtime analysis of recursive randomized algorithms. The probabilities and random variables in these recurrences respectively model random samplings and the random size of an instance passed to the recursive procedures in a randomized algorithm. The runtime of various randomized algorithms, such as S [Hoare 1961a],

S [Hoare 1961b], and D C [Motwani and Raghavan 1995, Chapter 9], can be captured by PRRs. The runtime complexity analysis of a randomized algorithm is therefore reduced to that of its corresponding PRR. Thus, formal analysis of PRRs has become a central and classical problem in the analysis of randomized algorithms and has also received a renewed attention in recent years [START_REF] Bazzi | The Solution of Linear Probabilistic Recurrence Relations[END_REF]Chatterjee et al. 2017a;[START_REF] Chaudhuri | Probabilistic Recurrence Relations Revisited[END_REF][START_REF] Devdatt | Concentration of Measure for the Analysis of Randomized Algorithms[END_REF][START_REF] Karp | Probabilistic Recurrence Relations[END_REF][START_REF] Mcdiarmid | Large Deviations for Quicksort[END_REF][START_REF] Motwani | Randomized Algorithms[END_REF][START_REF] Tassarotti | Probabilistic Recurrence Relations for Work and Span of Parallel Algorithms[END_REF][START_REF] Tassarotti | Veri ed Tail Bounds for Randomized Programs[END_REF]. In formal analysis of PRRs, the two fundamental problems are Expected Runtime Analysis and Concentration Bound Analysis.

Expected runtime analysis. Expected runtime analysis aims to derive tight upper bounds on the expected runtime complexity of a given PRR [START_REF] Motwani | Randomized Algorithms[END_REF]. The expected runtime is a key performance criterion that examines for how long the randomized algorithm will execute in expectation. This is important since many randomized algorithms perform poorly in the worst case but are e cient on average. For example, the well-known S algorithm [Hoare 1961a] takes Ω(𝑛 2 ) time in the worst case but has an expected runtime of 𝑂 (𝑛 log 𝑛). Compared with the classical runtime analysis of non-probabilistic recurrences that can typically be solved using the Master theorem [Cormen et al. 2009, Chapter 4][ [START_REF] Akra | On the Solution of Linear Recurrence Equations[END_REF]], expected runtime analysis of PRRs is much more intricate and there is no straightforward solution similar to Master theorem, except for certain special cases [START_REF] Bazzi | The Solution of Linear Probabilistic Recurrence Relations[END_REF]. Instead, the expected runtime of PRRs is obtained either through the traditional methods of Harmonic series and generating functions via computer algebra [START_REF] Flajolet | Automatic Average-Case Analysis of Algorithm[END_REF][START_REF] Zimmermann | The automatic complexity analysis of divide-and-conquer algorithms[END_REF]] or through automated guess-and-check algorithms and approximation of integrals [Chatterjee et al. 2017a].

Concentration bound analysis. Concentration bound analysis [START_REF] Devdatt | Concentration of Measure for the Analysis of Randomized Algorithms[END_REF][START_REF] Karp | An introduction to randomized algorithms[END_REF]] considers upper bounds on the probability that a given PRR does not terminate within a prescribed time limit. The goal is to show that this probability converges rapidly to zero as the time limit tends to in nity. The most classical approach to this problem is that of [START_REF] Karp | Probabilistic Recurrence Relations[END_REF]]. Since then, nding automated approaches for obtaining tighter bounds than Karp's method has been a key open problem. Moreover, ad-hoc methods have been devised to improve these bounds over a speci c recurrence, such as S [START_REF] Mcdiarmid | Large Deviations for Quicksort[END_REF].

Comparison of the two analyses. Compared with the expected runtime analysis, concentration bound analysis is a more re ned problem. It examines the probability that the runtime of a given PRR deviates signi cantly from its expected value. Unsurprisingly, methods used for concentration bound analysis are typically much more involved than expected runtime analysis. For instance, expected runtime analysis allows a template-based approach through guess-and-check functions [Chatterjee et al. 2017a], while there is no such template-based approach for the concentration bound analysis in the literature. Furthermore, in some cases, such as S , it is possible to manually derive the exact expected runtime of probabilistic recurrences by applying the generating-function method, but no such approach is known when it comes to concentration bounds.

Karp's cookbook. The classical method to obtain concentration bounds was proposed in [START_REF] Karp | Probabilistic Recurrence Relations[END_REF]] where a "cookbook" formula similar to the Master theorem is proposed. This approach requires the solution to the expected runtime of a deterministic version from the PRR of concern, and then outputs a coarse concentration bound for the original PRR by a xed mathematical formula. Based on Karp's method, the work [START_REF] Chaudhuri | Probabilistic Recurrence Relations Revisited[END_REF] proposed an approach that weakens several conditions required in [START_REF] Karp | Probabilistic Recurrence Relations[END_REF]] and obtains bounds that are slightly worse than [START_REF] Karp | Probabilistic Recurrence Relations[END_REF]]. For the S algorithm, the work [START_REF] Mcdiarmid | Large Deviations for Quicksort[END_REF]] performed an involved ad-hoc manual analysis to derive the asymptotically optimal concentration bound. Recently, there has been a renewed interest in this problem. The work [START_REF] Tassarotti | Veri ed Tail Bounds for Randomized Programs[END_REF] mechanized the proof for concentration bound of S in the interactive theorem prover Coq [START_REF] Bertot | Interactive Theorem Proving and Program Development -Coq'Art: The Calculus of Inductive Constructions[END_REF].

Challenges and gaps. Although concentration bound analysis for PRRs is a classical problem that has been studied for around thirty years, there are still several key challenges. First, while the bounds in [START_REF] Karp | Probabilistic Recurrence Relations[END_REF]] were known to be coarse (as exempli ed by [START_REF] Mcdiarmid | Large Deviations for Quicksort[END_REF] and the results of this work), it has been a key open problem for almost three decades to nd general approaches that (a) are automated and (b) can signi cantly improve these classical results. All previous improvements lacked automation and were obtained on a xed benchmark, i.e. they consider only a speci c single recurrence, such as S . Second, in several application domains, such as computational geometry, the preprocessing time in each recurrence equation is randomized rather than deterministic, and the analysis involves an acyclic system of PRRs rather than a single equation. Unfortunately, existing approaches all assumed that the preprocessing time is deterministic and are therefore inapplicable. Third, all existing results consider concentration bounds that are exponentially decreasing and other forms of concentration bounds have not been thoroughly investigated.

Our contribution. In this work, we consider the problem of automated concentration bound analysis for PRRs. Our main contributions are as follows:

(1) We propose a sound approach based on the novel concept of Taylor Supermartingale Function Sequences (TSFS). We show that the existence of such functions for a given PRR leads to a concentration bound on its runtime. The main theoretical tools we use here are Markov's inequality (Theorem 3.2) and Optional Stopping Theorem (Theorem 3.3). (2) We develop algorithmic approaches to synthesize a TSFS in exponential and polynomial form for a given PRR. Our main tools in this synthesis process are Taylor expansion and tight approximations of expectation of moments with logarithmic terms. (3) Compared with Karp's cookbook method, our algorithmic approaches derive much tighter concentration bounds over classical randomized algorithms such as S , S , R S and D C . For example, in case of S , we nd a bound for Pr[𝑇 ≥ 20 • 𝑛] that is 10 6 times tighter. Moreover, for S and R S , our bounds become arbitrarily tighter than Karp's as 𝑛 grows. See Section 6 for details. (4) Furthermore, our algorithmic approaches can handle randomized algorithms that are beyond the scope of [START_REF] Karp | Probabilistic Recurrence Relations[END_REF]] and other existing methods such as [START_REF] Chaudhuri | Probabilistic Recurrence Relations Revisited[END_REF][START_REF] Tassarotti | Probabilistic Recurrence Relations for Work and Span of Parallel Algorithms[END_REF]. To be more precise, our algorithmic approaches can handle PRRs with randomized preprocessing phase in a randomized algorithm. In contrast, other approaches require the preprocessing phase to be deterministic. Moreover, our algorithmic approaches can handle systems of probabilistic recurrences such as those arising from the S LP [START_REF] Seidel | Small-dimensional linear programming and convex hulls made easy[END_REF]] and the S E D problem [START_REF] Welzl | Smallest enclosing disks (balls and ellipsoids)[END_REF]]. To the best of our knowledge, no existing method can handle concentration bounds for systems of probabilistic recurrences.

Key novel aspects.

The key novelties of our approach are two-fold. First, we establish the novel notion of Taylor Supermartingale Function Sequences (TSFS) for solving concentration bounds of PRRs. Second, to synthesize a TSFS in exponential or polynomial form, we consider the novel idea of Taylor expansion and over-approximation of expectation of moments with logarithmic terms in our algorithmic approaches. The novel aspect is that while these techniques are quite simple, they surprisingly lead to rst signi cant improvements in classical decades-old bounds for fundamental and well-studied randomized algorithms. Our approach can even handle classical randomized algorithms that are beyond the reach of all previous methods.

PRELIMINARIES

We rst review some necessary concepts in probability theory in Section 2.1. Then, in Section 2.2, we introduce Systems of Probabilistic Recurrence Relations (SPRRs) and formally de ne our problem.

Background Notions in Probability Theory

Throughout the paper, we denote by N, Z and R the sets of non-negative integers, integers and real numbers, respectively. Below, we present only a brief description of the concepts needed for this work. We refer to standard textbooks (e.g. [START_REF] Billingsley | Probability and Measure[END_REF][START_REF] Rey | A First Look at Rigorous Probability Theory[END_REF][START_REF] Williams | Probability with Martingales[END_REF]) for a more detailed treatment.

Discrete probability distributions.

A discrete probability distribution over a countable set 𝑈 is a function 𝜂 : 𝑈 → [0, 1] such that 𝑢 ∈𝑈 𝜂 (𝑢) = 1. The support of 𝜂 is de ned as supp(𝜂) := {𝑢 ∈ 𝑈 | 𝜂 (𝑢) > 0}. Informally, 𝜂 (𝑢) is the probability that a sampling from the distribution 𝜂 returns 𝑢. Probability spaces. A probability space is a triple (Ω, F , Pr) where Ω is a non-empty set (called the sample space), F is a 𝜎-algebra over Ω (i.e. a collection of subsets of Ω that contains the empty set ∅ and is closed under complementation and countable union) and Pr is a probability measure on F , i.e. a function Pr : F → [0, 1] such that (i) Pr(Ω) = 1, and (ii) for all pairwise-disjoint set-sequences 𝐴 1 , 𝐴 2 , • • • ∈ F (i.e. 𝐴 𝑖 ∩ 𝐴 𝑗 = ∅ whenever 𝑖 ≠ 𝑗) it holds that ∞ 𝑖=1 Pr(𝐴 𝑖 ) = Pr ∞ 𝑖=1 𝐴 𝑖 . Elements of F are conventionally called events. An event 𝐴 ∈ F holds almost-surely (a.s.) if Pr(𝐴) = 1.

Random variables. A random variable 𝑋 from a probability space (Ω, F , Pr) is an F -measurable function 𝑋 : Ω → R ∪ {-∞, +∞}, i.e. a function such that for all 𝑑 ∈ R ∪ {-∞, +∞}, the set {𝜔 ∈ Ω | 𝑋 (𝜔) < 𝑑 } lies in F .

Expectation. The expected value of a random variable 𝑋 from a probability space (Ω, F , Pr), denoted by E[𝑋 ], is de ned as the Lebesgue integral of 𝑋 w.r.t. the probability measure Pr, i.e., E[𝑋 ] := ∫ 𝑋 d Pr. The precise de nition of Lebesgue integral is rather technical and therefore omitted. See [Williams 1991, Chapter 5] for a formal de nition. In the case that ran 𝑋 = {𝑑 0 , 𝑑 1 , . . .} (ran 𝑋 is the range of the random variable 𝑋 when viewed as a function) is countable, then we have

E[𝑋 ] = ∞ 𝑘=0 𝑑 𝑘 • Pr(𝑋 = 𝑑 𝑘 ). Filtrations.
A ltration of a probability space (Ω, F , Pr) is an in nite sequence {F 𝑛 } 𝑛 ∈N of sigmaalgebras over Ω such that F 𝑛 ⊆ F 𝑛+1 ⊆ F for all 𝑛 ∈ N. Intuitively, a ltration models the information F 𝑛 available at a speci c time point 𝑛 ≥ 0.

Conditional expectation. Let 𝑋 be any random variable from a probability space (Ω, F , Pr) such that E[|𝑋 |] < ∞ (i.e. 𝑋 is integrable). Then, given any sub-sigma-algebra G ⊆ F , there exists a random variable from (Ω, F , Pr), denoted by E [𝑋 | G], such that:

(E1) E [𝑋 | G] is G-measurable, and (E2) E [|E [𝑋 | G] |] < ∞, and (E3) for all 𝐴 ∈ G, we have ∫ 𝐴 E [𝑋 | G] d Pr = ∫ 𝐴 𝑋 d Pr. The random variable E [𝑋 | G] is called the conditional expectation of 𝑋 given G, and is a.s. unique in the sense that if 𝑌 is another random variable satisfying (E1)-(E3), then Pr(𝑌 = E [𝑋 | G]) = 1.
We refer to [Williams 1991, Chapter 9] for details. The intuition is that E [𝑋 | G] represents the expectation of 𝑋 when assuming the information in G.

Discrete-time stochastic processes. A discrete-time stochastic process is an in nite sequence Γ = {𝑋 𝑛 } 𝑛 ∈N of random variables where the random variables 𝑋 𝑛 are all from the same probability space (Ω, F , Pr). The process Γ is adapted to a ltration {F 𝑛 } 𝑛 ∈N if for all 𝑛 ∈ N, 𝑋 𝑛 is F 𝑛 -measurable. Intuitively, the random variable 𝑋 𝑛 measures the value at the 𝑛-th step of the process Γ.

Stopping times. Given a ltration {F 𝑛 } 𝑛 ∈N with respect to a probability space (Ω, F , Pr), a stopping time 𝜌 is a random variable 𝜌 :

Ω → {0, 1, 2, 3, • • • , ∞}, such that for every 𝑛 ∈ {0, 1, 2, • • • , ∞}, the event {𝜔 | 𝜌 (𝜔) ≤ 𝑛} is measurable in F 𝑛 . In other words, {𝜔 | 𝜌 (𝜔) ≤ 𝑛} ∈ F 𝑛 .

Martingales and supermartingales

. A discrete-time stochastic process Γ = {𝑋 𝑛 } 𝑛 ∈N adapted to a ltration {F 𝑛 } 𝑛 ∈N is a martingale (resp. supermartingale) if for every 𝑛 ∈ N, E[|𝑋 𝑛 |] < ∞ and it holds a.s. that E [𝑋 𝑛+1 | F 𝑛 ] = 𝑋 𝑛 (resp. E [𝑋 𝑛+1 | F 𝑛 ] ≤ 𝑋 𝑛 ).
Intuitively, a martingale (resp. supermartingale) is a discrete-time stochastic process in which for an observer who has seen the values of 𝑋 0 , . . . , 𝑋 𝑛 , the expected value at the next step, i.e. E [𝑋 𝑛+1 | F 𝑛 ], is equal to (resp. no more than) the last observed value 𝑋 𝑛 . Also, note that in a martingale, the observed values for 𝑋 0 , . . . , 𝑋 𝑛-1 do not matter given that

E [𝑋 𝑛+1 | F 𝑛 ] = 𝑋 𝑛 . In contrast, in a supermartingale, the only requirement is that E [𝑋 𝑛+1 | F 𝑛 ] ≤ 𝑋 𝑛 and hence E [𝑋 𝑛+1 | F 𝑛 ] may depend on 𝑋 0 , . . . , 𝑋 𝑛-1 .
Also, note that F 𝑛 might contain more information than just the observations of 𝑋 𝑖 's. We refer to [Williams 1991, Chapter 10] for a deeper treatment.

Example 2.1. Consider an unbiased and discrete random walk, in which we start at a position 𝑋 0 , and at each second walk one step to either left or right with equal probability. Let 𝑋 𝑛 denote our position after 𝑛 seconds. It is easy to verify that E (𝑋 𝑛+1 |𝑋 0 , . . . , 𝑋 𝑛 ) = 1 2 (𝑋 𝑛 -1) + 1 2 (𝑋 𝑛 + 1) = 𝑋 𝑛 . Hence, this random walk is a martingale. Note that by de nition, every martingale is also a supermartingale. As another example, consider the classical gambler's ruin: a gambler starts with 𝑌 0 dollars of money and bets continuously until he loses all of his money. If the bets are unfair, i.e. the expected value of his money after a bet is less than its expected value before the bet, then the sequence {𝑌 𝑛 } 𝑛 ∈N 0 is a supermartingale. In this case, 𝑌 𝑛 is the gambler's total money after 𝑛 bets. On the other hand, if the bets are fair, then {𝑌 𝑛 } 𝑛 ∈N 0 is a martingale. See Appendix A for another illustrative example of martingales.

Systems of Probabilistic Recurrence Relations

In this work, we consider the formal analysis of Systems of Probabilistic Recurrence Relations (SPRRs). An SPRR is a collection of probabilistic recurrence equations that describe the running-time behavior of a randomized algorithm. Each recurrence equation corresponds to a procedure in the algorithm.

Systems of Probabilistic Recurrence Relations (SPRRs

). An SPRR Δ is a nite collection of probabilistic recurrence relations eq 1 , . . . , eq 𝑚 , which models the running time of a randomized algorithm with 𝑚 procedures, for which each recurrence equation eq 𝑟 (1 ≤ 𝑟 ≤ 𝑚) corresponds to a procedure in the algorithm (which we shall denote by proc 𝑟 ). Each eq 𝑟 (1 ≤ 𝑟 ≤ 𝑚) takes the following form:

𝑇 𝑟 (𝑛) = Cost 𝑟 (𝑛) + 𝑆 𝑟 𝑗=1 𝑇 𝑟 (ℎ 𝑟,𝑗 (𝑛)) + 𝑂 𝑟 𝑗=1 𝑇 𝑞 𝑟,𝑗 (𝑔 𝑟,𝑗 (𝑛)) if 𝑛 ≥ 𝑐 𝑟 + 1 0 if 𝑛 ≤ 𝑐 𝑟 . (eq 𝑟 )
Here, the variable 𝑛 is a non-negative integer modeling the size of the input passed to procedure proc 𝑟 and 𝑐 𝑟 is a non-negative integer constant below which the procedure proc 𝑟 halts immediately. The rest of the notation is as follows:

• 𝑇 𝑟 (𝑛) is a random variable modeling the running time of proc 𝑟 on an input of size 𝑛.

• Cost 𝑟 (𝑛) is a random variable that models the preprocessing time of proc 𝑟 on input size 𝑛.

• 𝑆 𝑟 is a constant that represents the number of procedure calls in proc 𝑟 to proc 𝑟 itself. Among these procedure calls, each ℎ 𝑟,𝑗 (𝑛) represents the random size of the instance passed to the 𝑗-th recursive procedure call (to proc 𝑟 itself). • 𝑂 𝑟 is a constant that represents the number of procedure calls in proc 𝑟 to procedures other than proc 𝑟 . We denote the callee in the 𝑗-th such call by 𝑞 𝑟,𝑗 . Note that 𝑞 𝑟,𝑗 ≠ 𝑟 . Moreover, 𝑔 𝑟,𝑗 (𝑛) is a random variable associated with input size 𝑛 that models the parameter size passed to the 𝑗-th call. Since ℎ 𝑟,𝑗 (𝑛) and 𝑔 𝑟,𝑗 (𝑛) both model input sizes to procedure calls, we assume that they observe discrete probability distributions over non-negative integers.

Intuitive description of SPRRs. Informally, given an initial input of size 𝑛 * and a prescribed main procedure proc 𝑟 * , an SPRR runs as follows. Initially, the stack contains only the main procedure proc 𝑟 * with the initial input size 𝑛 * . Then, the following steps are performed:

Step 1. If the current stack is emty, then the SPRR halts, otherwise it goes to Step 2.

Step 2. The SPRR pops the top of the current stack to get the current procedure proc 𝑟 and the current input size 𝑛.

Step 3. If 𝑛 ≤ 𝑐 𝑟 then the SPRR terminates the current procedure immediately and goes back to

Step 1, otherwise it goes to Step 4.

Step 4. The SPRR runs the preprocessing stage that takes a randomized amount Cost 𝑟 (𝑛) of execution time, and then goes to Step 5.

Step 5. The SPRR pushes to the top of the stack the 𝑆 𝑟 procedure calls to proc 𝑟 itself with their corresponding input sizes ℎ 𝑟,𝑗 (𝑛) (1 ≤ 𝑗 ≤ 𝑆 𝑟 ) and the 𝑂 𝑟 procedure calls to proc 𝑞 𝑟,𝑗 's with their input sizes 𝑔 𝑟,𝑗 (𝑛) (1 ≤ 𝑗 ≤ 𝑂 𝑟 ), and then goes back to Step 1. Note that in Step 5, we consider an arbitrarily-xed order for pushing the procedural calls to the stack, since we are only interested in the total execution time. From the intuition above, it is straightforward to see that the random variables 𝑇 𝑟 (𝑛) (1 ≤ 𝑟 ≤ 𝑚) satisfy the recurrence equations of the SPRR.

Below we present several examples that model classical randomized algorithms as SPRRs. Due to space limitations, detailed descriptions of the algorithms are relegated to Appendix B.

Example 2.2 (Q S

). Consider the problem of nding the 𝑑-th smallest element in an unordered array of 𝑛 distinct elements. A classical randomized algorithm for solving this problem is S [Hoare 1961b] that provides a simple randomized variant matching the best-known 𝑂 (𝑛) deterministic runtime. Since there is only one recursive procedure in S , we model the algorithm as the following SPRR with one recurrence equation (where we have 𝑚 = 1, 𝑇 1 = 𝑇 , Cost 1 (𝑛) = 𝑛, 𝑆 1 = 1 and 𝑂 1 = 0):

𝑇 (𝑛) = 𝑛 + 𝑇 (ℎ(𝑛)) .
Here, 𝑇 (𝑛) represents the number of comparisons performed by S over an input of size 𝑛, and ℎ(𝑛) is the random variable that captures the size of the remaining array that has to be searched recursively. It can be derived that ℎ(𝑛) = max{𝑖, 𝑛 -1 -𝑖} where 𝑖 is sampled uniformly from {0, . . . , 𝑛 -1}.

Example 2.3 (Q S ). Consider the classical problem of sorting an array of 𝑛 distinct elements. A well-known randomized algorithm for solving this problem is S [Hoare 1961a]. Just as before, there is only one recursive procedure in S . Thus, we model the algorithm as the following SPRR with one recurrence equation (where we have

𝑚 = 1, 𝑇 1 = 𝑇 , Cost 1 (𝑛) = 𝑛, 𝑆 1 = 2 and 𝑂 1 = 0): 𝑇 (𝑛) = 𝑛 + 𝑇 (ℎ 1 (𝑛)) + 𝑇 (ℎ 2 (𝑛)) .
Here, ℎ 1 (𝑛) and ℎ 2 (𝑛) are random variables that capture the sizes of the two sub-arrays. It is easy to see that ℎ 1 (𝑛) = 𝑖 and ℎ 2 (𝑛) = 𝑛 -1 -𝑖 where 𝑖 is uniformly sampled from {0, . . . , 𝑛 -1}.

Example 2.4 (S LP). Consider a linear programming problem in 𝑑 dimensions, i.e. with 𝑑 variables x = (x 1 , . . . , x 𝑑 ) and 𝑛 constraints a

𝑇 1 • x ≤ 𝑏 1 , . . . , a 𝑇 𝑛 • x ≤ 𝑏 𝑛 . Suppose the objective is to minimize c 𝑇 • x. S
LP [Motwani and Raghavan 1995, Chapter 9.10] is a well-known randomized incremental algorithm that handles this problem. To model the algorithm by an SPRR, we set up a system with 𝑑 equations 𝑒𝑞 1 , 𝑒𝑞 2 , . . . , 𝑒𝑞 𝑑 , where 𝑒𝑞 𝑟 corresponds to the 𝑟 -dimensional recursive subinstances of the algorithm. We have 𝑐 𝑟 = 𝑟 and for every 𝑛 ≥ 𝑟 + 1:

𝑇 𝑟 (𝑛) = Cost 𝑟 (𝑛) + 𝑇 𝑟 (ℎ 𝑟,1 (𝑛)) + 𝑇 𝑟 -1 (𝑔 𝑟,1 (𝑛)) (1)
Where ℎ 𝑟,1 (𝑛) = 𝑛 -1. Moreover, Cost 𝑟 and 𝑔 𝑟,1 are dependent and jointly observe the following discrete distribution: with probability 𝑟 𝑛 , Cost 𝑟 (𝑛) = 𝑟 • 𝑛 and 𝑔 𝑟,1 (𝑛) = 𝑛 -1, and with probability 1 -𝑟 𝑛 , Cost 𝑟 (𝑛) = 𝑟 and 𝑔 𝑟,1 (𝑛) = 0.

Semantics of SPRRs. We now present formal semantics that treat each random variable 𝑇 𝑟 (𝑛) as the accumulated cost of a stochastic process derived from the underlying SPRR. Below, we x an SPRR Δ = {eq 1 , • • • , eq 𝑚 } where each recurrence equation eq 𝑟 takes the form (eq 𝑟 ) as speci ed on Page 5. We represent a procedure call 𝜎 as an ordered pair (𝑟, 𝑛) of non-negative integers where 𝑟 represents the 𝑟 -th procedure and 𝑛 is the size of the parameter passed to this procedure call. We rst formally de ne the notion of stacks and stack transition functions.

Stacks.

A stack 𝛽 is a nite word 𝜎 1 . . . 𝜎 𝑘 (𝑘 ≥ 0) such that for each 1 ≤ 𝑗 ≤ 𝑘, the symbol 𝜎 𝑗 = (𝑟 𝑗 , 𝑛 𝑗 ) is a procedure call. The word 𝛽 is organized such that 𝜎 1 corresponds to the top of the stack and 𝜎 𝑘 is the bottom. We use 𝜀 to denote the empty stack. For a stack 𝛽 = (𝑟, 𝑛)𝛽 , we de ne Cost(𝛽) := Cost 𝑟 (𝑛). Intuitively, Cost(𝛽) is the preprocessing time of the call at the top of the stack.

Stack transition functions.

A stack transition function describes the stochastic e ect of popping the procedure call at the top of the stack and pushing all its sub-procedure calls. Formally, the stack transition function is a function 𝜅 that maps every procedure call 𝜎 = (𝑟, 𝑛) to the discrete probability distribution 𝜅 (𝜎) over stacks such that for each stack 𝛽, 𝜅 (𝜎) (𝛽) is the probability that exactly all the procedure calls in 𝛽 are pushed to the stack while respecting the order they appear in the word 𝛽 * . The exact form of 𝜅 is completely determined by the summations 𝑆 𝑟 𝑗=1 𝑇 𝑟 (ℎ 𝑟 𝑗 (𝑛)) and 𝑂 𝑟 𝑗=1 𝑇 𝑞 𝑟 𝑗 (𝑔 𝑟 𝑗 (𝑛)) in (eq 𝑟 ) so that it can be calculated exactly from the probability distributions of the random variables ℎ 𝑟 𝑗 (𝑛) and 𝑔 𝑟 𝑗 (𝑛). However, we do not need to provide a more detailed notation in order to present our results, which only depend on the conceptual notion of a stack transition function.

Based on stacks and the stack transition function, we now de ne the semantics of SPRRs through discrete-time Markov chains [Baier and Katoen 2008, Chapter 10].

Markov chain semantics of SPRRs. The semantics of an SPRR Δ is a discrete-time Markov chain (𝑆, P). The set 𝑆 of Markov chain states is the set of all stacks. For any two stacks 𝛽, 𝛽 , the transition probability P(𝛽, 𝛽 ) is determined by the following two cases: Case 1: 𝛽 = 𝜀. Then P(𝛽, 𝜀) := 1 and P(𝛽, 𝛽 ) := 0 for all other stacks 𝛽 ≠ 𝜀. Case 2: 𝛽 = 𝜎𝛽 with the top procedure call 𝜎. Then P(𝛽, 𝛽 ) := 0 if there is no stack 𝛾 such that 𝛽 = 𝛾 𝛽 , and P(𝛽, 𝛽 ) := 𝜅 (𝜎) (𝛾) if 𝛽 = 𝛾 𝛽 for some (unique) stack 𝛾. In other words, if the current stack is empty, i.e. the whole SPRR is terminated, which corresponds to Case 1, then the SPRR always stays at the empty stack. Otherwise (Case 2), the SPRR pops the top procedure call of the current stack and stochastically pushes sub-procedure calls back to the top of the stack. Given the Markov chain (𝑆, P) of an SPRR Δand an initial stack 𝛽 * † , the probability space (Ω Δ , F Δ , Pr Δ,𝛽 * ) for Δ is the one induced by the Markov chain of Δ, see [Baier and Katoen 2008, Chapter 10]) for details. A brief description of the construction of the probability space is also provided in Appendix C. We denote the expectation under the probability space (Ω Δ , F Δ , Pr Δ,𝛽 * ) as E Δ,𝛽 * , and omit the subscript Δ whenever it is clear from the context. Below we introduce several random variables under (Ω Δ , F Δ , Pr Δ,𝛽 * ) that will be used later.

Random stack next 𝛽 . We denote by next 𝛽 a random stack that observes the discrete probability distribution P(𝛽, -) and represents the next stack stochastically determined by the current stack 𝛽. Note that a random stack consists of a collection of random variables that describe all the information of the procedure calls in the stack.

Random traces. From the Markov chain semantics, a random execution of an SPRR with a main procedure call induces an in nite sequence of random stacks 𝛽 0 , 𝛽 1 , . . . such that the random stack 𝛽 0 corresponds to the (deterministic) main procedure call and each random stack 𝛽 𝑘 represents the stack at the 𝑘-th step of an execution, so that for each 𝑘 ≥ 0, 𝛽 𝑘+1 is sampled from next 𝛽 𝑘 . We call 𝛽 0 , 𝛽 1 , . . . a random trace.

Top-stack random variables. For each 𝑘 ≥ 0, we de ne a pair 𝜎 𝑘 = ( 𝑟 𝑘 , 𝑛 𝑘 ) of random variables such that 𝑟 𝑘 (resp. 𝑛 𝑘 ) represents the procedure (resp. the parameter size passed to the procedure) on the top of the stack at the 𝑘-th step of a run. Formally, given any in nite sequence of stacks 𝜔 = 𝛽 0 𝛽 1 . . . , if 𝛽 𝑘 = 𝜎𝛽 for some procedure call 𝜎 = (𝑟, 𝑛) and stack 𝛽 then we have 𝑟 𝑘 (𝜔) := 𝑟 and 𝑛 𝑘 (𝜔) := 𝑛. If the 𝑘-th stack is empty, we de ne 𝑟 𝑘 = 𝑛 𝑘 =⊥ .

Stopping time 𝜏. We de ne the stopping time random variable 𝜏 as the rst time point when the SPRR reaches the empty stack, i.e. 𝜏 (𝜔) := inf {𝑘 | 𝛽 𝑘 = 𝜀} for every run 𝜔 = {𝛽 𝑗 } 𝑗 ∈N of stacks (where inf ∅ := ∞). It is worth noting that the stopping time 𝜏 is di erent from the total execution time of the underlying randomized algorithm since it does not take the execution time Cost 𝑟 (𝑛) of the preprocessing stage into account. We now de ne the notion of accumulated cost for SPRRs.

Accumulated cost. We rst introduce a random variable C 𝑘 for each 𝑘 ≥ 0 that represents the accumulated cost in a run of Δ until the 𝑘-th step. We treat the time spent on the preprocessing stage as cost so that the total accumulated cost until a run of Δ reaches empty stack corresponds directly to the total execution time of the underlying randomized algorithm. Thus, for every run 𝜔, we de ne C 𝑘 (𝜔) := 𝑘-1 𝑗=0 Cost 𝑟 𝑗 (𝜔) ( 𝑛 𝑗 (𝜔)), where Cost 𝑟 (⊥) := 0. The accumulated cost is simply de ned as C 𝜏 . From the de nition, it is straightforward to see that the accumulated cost C represents the randomized execution time of the underlying randomized algorithm.

We are now formally de ne the concentration bound analysis problem over SPRRs, which is the central problem in our paper.

Concentration bound analysis of SPRRs. Consider an SPRR Δ and a main procedure call 𝜎 * = (𝑟 * , 𝑛 * ). Given a symbolic expression 𝜅 over 𝑛 * , the goal is to infer an upper-bound on the probability that the accumulated cost C 𝜏 is no smaller than 𝜅. Formally, our task is to infer an upper-bound on

Pr Δ,𝜎 * (C 𝜏 ≥ 𝜅) (2)
where 𝜅 is a polynomial over the symbolic input size 𝑛 * and ln 𝑛 * with symbolic coe cients.

Intuitively, this corresponds to bounding the probability that the runtime of the randomized algorithm exceeds the bound 𝜅. Moreover, we would like the upper-bound to decrease dramatically and rapidly tend to zero as 𝜅 increases.

THE GENERAL THEORETICAL APPROACH

In this section, we develop sound approaches for solving the concentration bound analysis problem over SPRRs. We rst recall some basic results from probability theory, then demonstrate the theoretical foundation of our approaches, and nally present our main theorems under practical considerations. Below we x an SPRR with 𝑚 recurrence equations where the 𝑟 -th recurrence equation eq 𝑟 is in the form (eq 𝑟 ).

Basic Results from Probability Theory

We rst recall the well-known Markov's inequality. Intuitively, this inequality provides a way to obtain concentration bounds for a non-negative random variable. 1991, C 6]). For any non-negative random variable 𝑋 and real number 𝑎 > 0, we have that Pr(𝑋 ≥ 𝑎) ≤ E(𝑋 ) 𝑎 . Using the following lemma, we can apply Markov's inequality to a random variable that is not necessarily non-negative: L 3.2. Consider any random variable 𝑋 and monotonically non-decreasing function 𝑓 : R → [0, ∞). For any real number 𝑑 such that 𝑓 (𝑑) > 0, we have Pr(𝑋 ≥ 𝑑) ≤ E[𝑓 (𝑋 )]/𝑓 (𝑑).

T 3.1 (M ' [W

P

. Since 𝑓 is monotonically non-decreasing, we have Pr(𝑋 ≥ 𝑑) ≤ Pr(𝑓 (𝑋 ) ≥ 𝑓 (𝑑)) as 𝑋 ≥ 𝑑 implies 𝑓 (𝑋 ) ≥ 𝑓 (𝑑). By Markov's inequality (Theorem 3.1), we obtain Pr(𝑓 (𝑋 ) ≥ 𝑓 (𝑑)) ≤ E[𝑓 (𝑋 )]/𝑓 (𝑑). The result follows.

Our approaches also depend on the classical Optional Stopping Theorem, which is well-known in martingale theory. For simplicity, we only present a minimal part of the theorem, which is required for our method: 1991, T 10.10]). For any discretetime supermartingale Γ = {𝑋 𝑛 } 𝑛 ∈N adapted to a ltration {F 𝑛 } 𝑛 ∈N where every random variable 𝑋 𝑛 is non-negative, and any stopping time 𝜌 w.r.t the ltration {F 𝑛 } 𝑛 ∈N that is almost surely nite (i.e. Pr(𝜌 < ∞) = 1), we have

T 3.3 (O [W
E[𝑋 𝜌 ] ≤ E[𝑋 0 ].

Taylor Supermartingale Function Sequences and Our Main Theorem

The core concept in our theoretical foundations is the novel notion of Taylor Supermartingale Function Sequences (TSFS). These sequences lead to concentration bounds of SPRRs through Markov's inequality and Optional Stopping Theorem. Note that, by Lemma 3.2, to solve the concentration bound problem for the accumulated cost random variable C, it su ces to choose a monotonically non-decreasing function 𝑓 and then bound the expected value of 𝑓 (C). In this work, we consider a speci c class of functions, called T-functions, as our choices for the monotonically non-decreasing function 𝑓 in Lemma 3.2.

The Class T of T-functions. We consider a special class of functions, called Taylor functions (T-functions), that have Taylor expansions and non-negative derivatives everywhere on [0, ∞). To be more precise, we de ne the T-function class T as the set of all functions 𝑓 : [0, ∞) → R that (i) have non-negative derivatives 𝑓 (𝑘) for every order 𝑘 ≥ 0, and (ii) can be expanded into Taylor series at every non-negative real number. More precisely, a function 𝑓 : [0, ∞) → R falls in the class T (i.e. is a T-function) i the following conditions hold: (C1) for every order 𝑘 ≥ 0, the 𝑘-th order derivative function 𝑓 (𝑘) exists and is non-negative (i.e. 𝑓 (𝑘) (𝑥) ≥ 0 for every real number 𝑥 ≥ 0); (C2) for any reals 𝑥, 𝑦 ≥ 0, we have

𝑓 (𝑥) = ∞ 𝑘=0 (𝑥-𝑦) 𝑘 𝑘!
• 𝑓 (𝑘) (𝑦), i.e. the function 𝑓 can be expanded into its Taylor series at any non-negative real number 𝑦. Note that the function class T includes many exponential functions (e.g. 𝑓 (𝑥) = exp(𝑡 • 𝑥) where 𝑡 > 0) and polynomial functions (e.g. 𝑓 (𝑥) = 𝑑 • 𝑥 𝑘 where 𝑑, 𝑘 > 0 and 𝑘 is an integer). Our approach combines the ideas of T-functions, Markov's inequality (Theorem 3.2) and Optional Stopping Theorem (Theorem 3.3) to bound the concentration probability Pr Δ,𝜎 * (C ≥ 𝜅) in (2) by rst choosing a suitable T-function 𝑓 and then bounding all its 𝑘-th order derivatives 𝑓 (𝑘) . T

(M T

). Let Δ be a given SPRR. Consider a T-function 𝑓 ∈ T and an in nite sequence Ψ of functions Ψ = Ψ 0 , Ψ 1 , Ψ 2 , . . . such that each Ψ 𝑘 maps a stack 𝛽 of Δ to a non-negative real value Ψ 𝑘 (𝛽). Suppose that the sequence Ψ satis es the following conditions: (A1) for every 𝑘 ≥ 0, we have Ψ 𝑘 (𝜀) ≥ 𝑓 (𝑘) (0); (A2) for every stack 𝛽 and 𝑘 ≥ 0, we have

Ψ 𝑘 (𝛽) ≥ E ∞ 𝑗=0 Ψ 𝑘+𝑗 (next𝛽) 𝑗! • Cost 𝑗 (𝛽) .
Then, for every stack 𝛽 and 𝑘 ≥ 0, it holds that Ψ 𝑘 (𝛽) ≥ E 𝛽 𝑓 (𝑘) (C) .

Intuitively, (A1) states that Ψ 𝑘 is an upper-bound for the 𝑘-th order derivative of a T-function 𝑓 at the empty stack, for which the empty stack corresponds to the zero point at 𝑓 (𝑘) ; (A2) speci es that the expectation of expanding the function Ψ 𝑘 into its Taylor series over the next probability distribution next 𝛽 and the preprocessing time Cost(𝛽) of the current stack 𝛽 is no more than its current value Ψ 𝑘 . In the sequel, we refer to conditions (A1) and (A2) as Taylor Supermartingale (TS) conditions, and to an in nite sequence Ψ of functions Ψ = Ψ 0 , Ψ 1 , Ψ 2 , . . . that ful ll the TS conditions as a Taylor Supermartingale Function Sequence (TSFS). We now provide a sketch of our proof for this theorem. For a detailed proof, see Appendix D.1. P S T 3.4. Fix any stack 𝛽 and consider the random trace 𝛽 0 , 𝛽 1 , 𝛽 2 , . . . for the SPRR Δ with the initial stack 𝛽 0 = 𝛽. For each 𝑘 ≥ 0, we de ne G 𝑘 as the smallest 𝜎-algebra that makes the random variables in 𝛽 0 , . . . , 𝛽 𝑘 and C 0 , . . . , C 𝑘-1 measurable. It is straightforward to observe that {G 𝑘 } 𝑘 ≥0 is a ltration. For each 𝑘, 𝑖 ≥ 0, we de ne the random variable 𝑌 𝑘,𝑖 :=

∞ 𝑗=0 Ψ 𝑖+𝑗 ( 𝛽 𝑘 ) 𝑗! • C 𝑗
𝑘 (where we assume 0 0 := 1). We prove that for every 𝑘, 𝑖 ≥ 0, the following properties hold:

(a) for every 𝑘, 𝑖 ≥ 0, 𝑌 𝑘,𝑖 is non-negative; (b) 𝑌 𝜏,𝑖 ≥ 𝑓 (𝑖) (C 𝜏 );

(c) For every 𝑖 ≥ 0, 𝑌 𝑘,𝑖 forms a supermartingale with respect to the ltration {G 𝑘 } 𝑘 ≥0 . First, note that Property (a) follows directly from the de nition of 𝑌 𝑘,𝑖 . Second, as 𝛽 𝜏 = 𝜀, Property (b) can be deduced from (A1) as follows:

𝑌 𝜏,𝑖 = ∞ 𝑗=0 Ψ 𝑖+𝑗 ( 𝛽 𝜏 ) 𝑗! • C 𝑗 𝜏 ≥ ∞ 𝑗=0 𝑓 (𝑖+𝑗 ) (0) 𝑗!
• C 𝑗 𝜏 = 𝑓 (𝑖) (C 𝜏 ). We now prove Property (c).

E[𝑌 𝑘+1,𝑖 | G 𝑘 ] = E ∞ 𝑗=0 Ψ 𝑖+𝑗 ( 𝛽 𝑘+1 ) 𝑗! • C 𝑗 𝑘+1 G 𝑘 By the de nition of 𝑌 𝑘+1,𝑖 = E ∞ 𝑗=0 Ψ 𝑖+𝑗 ( 𝛽 𝑘+1 ) 𝑗! • (C 𝑘 + Cost( 𝛽 𝑘 )) 𝑗 G 𝑘 (By the de nition of C 𝑘+1 ) = E ∞ 𝑗=0 Ψ 𝑖+𝑗 ( 𝛽 𝑘+1 ) 𝑗! • 𝑗 𝑙=0 𝑗 𝑙 • C 𝑙 𝑘 • Cost( 𝛽 𝑘 ) 𝑗-𝑙 G 𝑘 (By the Binomial Theorem) = E ∞ 𝑗=0 𝑗 𝑙=0 Ψ 𝑖+𝑗 ( 𝛽 𝑘+1 ) 𝑙!( 𝑗-𝑙)! • C 𝑙 𝑘 • Cost( 𝛽 𝑘 ) 𝑗-𝑙 G 𝑘 By 𝑗 𝑙 = 𝑗! 𝑙!•( 𝑗-𝑙)! = E ∞ 𝑙=0 C 𝑙 𝑘 𝑙! • ∞ 𝑗=𝑙 Ψ 𝑖+𝑗 ( 𝛽 𝑘+1 ) ( 𝑗-𝑙)! • Cost( 𝛽 𝑘 ) 𝑗-𝑙 G 𝑘 (By rearrangement of summation) = ∞ 𝑙=0 C 𝑙 𝑘 𝑙! • E ∞ 𝑗=𝑙 Ψ 𝑖+𝑗 ( 𝛽 𝑘+1 ) ( 𝑗-𝑙)! • Cost( 𝛽 𝑘 ) 𝑗-𝑙 G 𝑘 (By ( * ) below) ≤ ∞ 𝑙=0 C 𝑙 𝑘 𝑙! • Ψ 𝑖+𝑙 ( 𝛽 𝑘 ) = 𝑌 𝑘,𝑖 (By the condition (A2) and (★) in the conditions (E1)-(E3))
The ( * ) includes properties of conditional expectations such as linearity, positivity, cMON and 'taking out what is known' as in [Williams 1991, Chapter 9.7, Page 88] (which we extend to nonnegative random variables). The omitted details at (★) follow the basic de nition of conditional expectation, i.e. (E1)-(E3), to relate Ψ 𝑖+𝑗 ( 𝛽 𝑘+1 ) and Ψ 𝑖+𝑗 (next

𝛽 𝑘

). Then, by applying the Optional Stopping Theorem (Theorem 3.3), we have that E[𝑌 𝜏,𝑖 ] ≤ E[𝑌 0,𝑖 ] for every 𝑖 ≥ 0. The result follows from the fact that 𝑌 𝜏,𝑖 ≥ 𝑓 (𝑖) (C 𝜏 ) (see Property (b)) and 𝑌 0,𝑖 = Ψ 𝑖 (𝛽).

Application of Theorem 3.4. Theorem 3.4 provides a way to bound the expected value E 𝛽 [𝑓 (C 𝜏 )] by Ψ 0 (𝛽). Combined with Lemma 3.2, we obtain a general procedure for deriving concentration bounds for SPRRs as follows:

• We nd an appropriate T-function 𝑓 ∈ F.

• We nd a TSFS Ψ = Ψ 0 , Ψ 1 , Ψ 2 , . . . to obtain the bound Ψ 𝑘 (𝛽) for E 𝛽 𝑓 (𝑘) (C 𝜏 ) .

• We apply Lemma 3.2 to derive the concentration bound Pr 𝜎 * (C 𝜏 ≥ 𝜅) ≤

E 𝜎 * [𝑓 ( C 𝜏 ) ] 𝑓 (𝜅)
≤ Ψ 0 (𝜎 * ) 𝑓 (𝜅) . However, the general procedure above is infeasible to realize algorithmically, since one does not x the form of a T-function 𝑓 and the TS conditions (A1) and (A2) are de ned for all stacks and an in nite sequence of functions. In other words, one would need to consider an in nite number of stacks and functions. To overcome this di culty, in the following we consider practical versions of Theorem 3.4 for exponential and polynomial T-functions (Theorems 3.5 and 3.6). Based on these, we then develop our algorithmic approaches in Section 5.

Composition Theorems for Exponential and Polynomial T-Functions

Instead of applying Theorem 3.4 directly, we consider simpli ed versions (called composition theorems) of Theorem 3.4 over SPRRs. In many scenarios, concentration bounds are exponentially-and polynomially-decreasing. Hence, we focus on exponential and polynomial T-functions since they correspond to these two common types of concentration bounds. In detail, we devise composition theorems for exponential and polynomial T-functions respectively. Informally, the composition theorems reduce the synthesis of an in nite length TSFS into synthesizing a nite number of functions. In Section 5, we will illustrate our algorithmic approach based on these two theorems. Below, we will introduce our two composition theorems for exponential and polynomial T-functions respectively.

The Composition Theorem for exponential T-functions is stated in Theorem 3.5. Intuitively, the theorem considers 2 • 𝑚 functions 𝑓 1 , . . . , 𝑓 𝑚 and 𝑡 1 , . . . , 𝑡 𝑚 where functions 𝑓 𝑟 and 𝑡 𝑟 correspond to the 𝑟 -th recurrence equation 𝑒𝑞 𝑟 in the input SPRR, and imposes constraints on these functions. If all the constraints are satis ed, then one could compose from these 2 • 𝑚 functions a TSFS that ful lls the TS conditions (A1) and (A2) w.r.t some exponential T-function 𝑓 (𝑥) = exp(𝑡 • 𝑥) (𝑡 > 0). Thus, the theorem allows our algorithm to focus on the synthesis of only a nite number of functions, i.e. 2 • 𝑚 functions, over the recurrence equations of the input SPRR.

T 3.5 (C T F E D B
). Consider 2 • 𝑚 functions 𝑓 1 , . . . , 𝑓 𝑚 and 𝑡 1 , . . . , 𝑡 𝑚 both from N into [0, ∞). Suppose that the following condition holds for every 1 ≤ 𝑟 ≤ 𝑚 and 1 ≤ 𝑛 ≤ 𝑛 * :

exp(𝑡 𝑟 (𝑛) • 𝑓 𝑟 (𝑛)) ≥ E exp 𝑡 𝑟 (𝑛) • Cost(𝑛) + 𝑆 𝑟 𝑖=1 𝑓 𝑟 (ℎ 𝑟,𝑖 (𝑛)) + 𝑂 𝑟 𝑖=1 𝑓 𝑞 𝑟,𝑖 (𝑔 𝑟,𝑖 (𝑛)) ( ‡)
Then, the following sequence

{Ψ 𝑖 } 𝑖 ≥0 is a TSFS w.r.t the exponential T-function 𝑓 (𝑥) := exp(𝑡 ★ • 𝑥): Ψ 𝑖 (𝛽) := 𝑡 𝑖 ★ • exp(𝑡 ★ • 𝑘 𝑗=1 𝑓 𝑟 𝑖 (𝑛 𝑖 ))
for any 𝑖 ≥ 0 and stack 𝛽 = (𝑟 1 , 𝑛 1 ) . . . (𝑟 𝑘 , 𝑛 𝑘 ) where 𝑡 ★ := min 1≤𝑟 ≤𝑚 min 1≤𝑛 ≤𝑛 * 𝑡 𝑟 (𝑛) is the minimum value over all function values. P S . We will verify conditions (A1) and (A2) in our main theorem (Theorem 3.4). (A1) directly follows from the non-negativity of 𝑓 1 , 𝑓 2 , . . . , 𝑓 𝑚 's. To verify (A2), rstly, due to the convexity of 𝑡 𝑟 (𝑛), we could apply condition ( ‡) to the top element of 𝛽, and then by the de nition of next 𝛽 , we conclude that (A2) holds. Please refer to Appendix D for details.

Thus, once we synthesize appropriate functions 𝑓 1 , . . . , 𝑓 𝑚 and 𝑡 1 , . . . , 𝑡 𝑚 that satisfy the conditions speci ed in the theorem, we could construct a TSFS that leads to an exponentially-decreasing concentration bound by Theorem 3.4.

We now illustrate our Composition Theorem for polynomial T-functions. Again, the theorem reduces the synthesis of TS functions derived from some polynomial T-function into the synthesis of a nite number of functions. In detail, the theorem considers 𝑚 • (𝑘 + 1) functions 𝑓 𝑟,𝑖 (1 ≤ 𝑟 ≤ 𝑚, 0 ≤ 𝑖 ≤ 𝑘) for which each 𝑓 𝑟,𝑖 (𝑛) informally represent an over-approximation of E[𝑇 𝑟 (𝑛) 𝑖 ] (i.e. the expectation of the higher moment 𝑇 𝑟 (𝑛) 𝑖 ), and imposes two groups of constraints for these functions. One group of constraints is set up for every recurrence equation (𝑒𝑞 𝑟 ) in the given SPRR and every degree 1 ≤ 𝑖 ≤ 𝑘 , and the other group of constraints are boundary conditions for 𝑓 𝑟,0 for each recurrence equation (𝑒𝑞 𝑟 ) in the SPRR. If all the constraints are satis ed, then our theorem can construct a TSFS with respect to a polynomial T-function

𝑓 (𝑥) = 𝑥 𝑘 . T 3.6 (C T P D B
). Let Δ be an SPRR with 𝑚 recurrence equations and 𝑘 be a parameter that represents the degree of desired monomials. Consider 𝑚 • (𝑘 + 1) functions 𝑓 𝑟,𝑖 (1 ≤ 𝑟 ≤ 𝑚, 0 ≤ 𝑖 ≤ 𝑘), each of which maps natural numbers to non-negative reals. Suppose that for every 1 ≤ 𝑟 ≤ 𝑚 and 1 ≤ 𝑖 ≤ 𝑘, the following constraints are satis ed:

∀𝑛.{1 ≤ 𝑛 ≤ 𝑛 * } ⇒ 𝑓 𝑟,𝑖 (𝑛) 𝑖! ≥ E I Cost 𝑟 (𝑛) 𝑗 𝑗! • 𝑆 𝑟 𝑢=1 𝑓 𝑟,𝑝𝑢 (ℎ 𝑟,𝑢 (𝑛)) 𝑝 𝑢 ! • 𝑂 𝑟 𝑢=1 𝑓 𝑞𝑟,𝑢 ,𝑒𝑢 (𝑔 𝑟,𝑢 (𝑛)) 𝑒 𝑢 ! (3) ∀𝑛.{1 ≤ 𝑛 ≤ 𝑛 * } ⇒ 𝑓 𝑟,0 (𝑛) ≥ 1 (4)
where

I := 0 ≤ 𝑗, 𝑝 1 , . . . , 𝑝 𝑆 𝑟 , 𝑒 1 , . . . , 𝑒 𝑂 𝑟 ≤ 𝑖, 𝑆 𝑟 𝑢=1 𝑝 𝑢 + 𝑂 𝑟 𝑢=1 𝑒 𝑢 + 𝑗 = 𝑖 . Then, the following se- quence of functions Ψ 0 , Ψ 1 , . . . satis es condition (A1)-(A2) of Theorem 3.4, with respect to the T- function 𝑓 (𝑥) = 𝑥 𝑘 , given 𝛽 represented by 𝛽 = (𝑟 1 , 𝑛 1 ) . . . (𝑟 𝑡 , 𝑛 𝑡 ):
• For 0 ≤ 𝑖 ≤ 𝑘, we de ne:

Ψ 𝑖 (𝛽) := 𝑘! • 0≤𝑝 1 ,••• ,𝑝 𝑡 ≤𝑘-𝑖, 𝑡 𝑢=1 𝑝 𝑢 =𝑘-𝑖 𝑡 𝑢=1 𝑓 𝑟𝑢 ,𝑝𝑢 (𝑛 𝑢 ) 𝑝 𝑢 ! .
• For 𝑖 ≥ 𝑘 + 1, we de ne: Ψ 𝑖 (𝛽) := 0. P S . We verify conditions (A1) and (A2) in our main theorem (Theorem 3.4). ( A1) is straightforward by non-negativity of 𝑓 𝑖,𝑗 's. To verify (A2), we rst apply condition (3) to the top element of 𝛽, and then through the de nition of next 𝛽 , we conclude that (A2) holds. Please refer to Appendix D for details.

AN ILLUSTRATION OF OUR ALGORITHMIC METHOD ON QUICKSELECT

Before presenting our synthesis algorithms in detail, we rst illustrate our approach over a representative example, namely the S SPRR of Example 2.2. A formal description of our algorithms will be provided in Section 5. Consider the S SPRR in Example 2.2. We aim to derive an exponentially decreasing concentration bound for this example.

Upper-bound on expected runtime. Suppose that we are provided with a function 𝑓 (𝑛) := 6 • 𝑛 that serves as an over-approximation of E[𝑇 (𝑛)]. This function can be automatically obtained using previous approaches for expected runtime analysis, such as [Chatterjee et al. 2017a]. In this work, we assume that such upper-bounds are given as part of the input.

Synthesis goal. Since the only recurrence equation in this example is

𝑇 (𝑛) = 𝑛 + 𝑇 (ℎ(𝑛)),
in order to apply our composition theorem for exponential concentration bounds (Theorem 3.5), it su ces to synthesize a function 𝑡 (𝑛) that satis es the following condition:

∀𝑛 (1 ≤ 𝑛 ≤ 𝑛 * ) ⇒ E[exp(𝑡 (𝑛) • (𝑛 + 𝑓 (ℎ(𝑛))))] ≤ exp(𝑡 (𝑛) • 𝑓 (𝑛)).
Intuition behind the constraint. Intuitively, the LHS represents the expectation of exp(𝑡 (𝑛) • 𝑓 (𝑛)) after the one-step expansion of the recurrence equation. Therefore, the inequality above encodes a supermartingale condition for the exponential term exp(𝑡 (𝑛) • 𝑓 (𝑛)).

Rephrasing the inequality. By dividing both sides of the condition by exp(𝑡 (𝑛) • 𝑓 (𝑛)), we obtain this equivalent formulation:

∀𝑛 (1 ≤ 𝑛 ≤ 𝑛 * ) ⇒ E[exp(𝑡 (𝑛) • 𝑋 (𝑛))] ≤ 1 where 𝑋 (𝑛) := 𝑛 + 6 • ℎ(𝑛) -6 • 𝑛 is a random variable.
Se ing up a template. Following the intuition above, our approach chooses a template for 𝑡 (𝑛) and sets 𝑡 (𝑛) := 𝜆/𝑛 where 𝜆 is an unknown coe cient that should be synthesized. Now the condition becomes:

∀𝑛.(1 ≤ 𝑛 ≤ 𝑛 * ) ⇒ E[exp ((𝜆/𝑛) • (𝑛 + 6 • ℎ(𝑛) -6 • 𝑛))] ≤ 1.
Intuitively, we use 𝑡 (𝑛) as a scaling factor that allows us to bound the range of 𝑡 (𝑛) • 𝑋 (𝑛) by constants.

Binary searching. Next, our approach performs a binary search to nd the optimal value of 𝜆. The binary search is performed with the help of a veri cation procedure that, given a guess for 𝜆, veri es whether the condition above holds.

Veri cation procedure. Suppose that our current guess is 𝜆 = 1. To verify whether the condition holds for 𝜆 = 1, our approach expands the exponential term in the expectation above through Taylor's expansion and over-approximates the remainder. To over-approximate the remainder, our algorithm rstly needs to bound 𝑡 (𝑛) • 𝑋 (𝑛) by constants.

Synthesizing the Bound 𝐷. Our approach rstly bounds 𝑋 (𝑛) by 𝐵(𝑛) := 3.67 • 𝑛, i.e., we have that -𝐵(𝑛) ≤ 𝑋 (𝑛) ≤ 𝐵(𝑛). By plugging in 𝑡 (𝑛) = 1/𝑛, we have that the bound 𝐷 for 𝑡 (𝑛) • 𝑋 (𝑛), i.e. -3.67 ≤ 𝑡 (𝑛) • 𝑋 (𝑛) ≤ 3.67 This leads to the following strengthened condition that soundly approximates the original one:

1 + 𝐿 𝑖=1 𝑛 -𝑖 𝑖! • E 𝑋 (𝑛) 𝑖 + 𝑛 -(𝐿+1) •exp(3.67) (𝐿+1)! • E 𝑋 (𝑛) 𝐿+1 ≤ 1. ( 5
)
In this inequality, 𝐿 is the depth of the Taylor expansion and is always odd.

Over-approximating expectation terms. To solve this, we still need to over-approximate the higher moment E 𝑋 (𝑛) 𝑖 . We introduce a moment approximation algorithm that can derive tight upper and lower approximations up to any prescribed precision. Suppose we are to over-approximate E[𝑋 (𝑛) 2 ]. Our approach rst expands E[𝑋 (𝑛) 2 ] by applying Newton's binomial theorem to the de nition of 𝑋 (𝑛):

E[𝑋 (𝑛) 2 ] = 36 • E[ℎ(𝑛) 2 ] -60 • 𝑛 • E[ℎ(𝑛)] + 37 • 𝑛 2 . Now, it su ces to over-approximate E[ℎ(𝑛) 2 ] and under-approximate E[ℎ(𝑛)]. To over-approximate E[ℎ(𝑛) 2 ],
our approach rst expands this expectation as a summation that iterates over all of the possible elements in the support of ℎ(𝑛)'s distribution. Recall that ℎ(𝑛) = max{𝑖, 𝑛 -1 -𝑖} where 𝑖 is uniformly distributed over {0, . . . , 𝑛 -1}. Then, it soundly removes oors and ceilings from expressions such as 𝑛 2 . As expected, in the removal of oors and ceilings, our approach distinguishes between even and odd 𝑛. Below, we present the even case:

E[ℎ(𝑛) 2 ] ≤ 1 𝑛 • 𝑛-1 𝑖=𝑛/2 𝑖 2 + 𝑛-1 𝑖=𝑛/2 𝑖 2 .
Next, by applying our moment approximation algorithm that derives upper and lower bounds for any summation of the form 𝑛-1 𝑖=0 𝑖 𝑎 ln 𝑏 𝑖 up to any prescribed precision, our approach obtains that (i)

E[ℎ(𝑛) 2 ] is exactly 7 12 • 𝑛 2 -3 4 • 𝑛 + 1 6 , and (ii) E[ℎ(𝑛)] is exactly 3 8 • 𝑛 -1 4 . As a consequence, it over-approximate E[𝑋 (𝑛) 2 ] as 13 • 𝑛 2 + 3 • 𝑛 + 6.
Checking the strengthened inequality. Finally, after approximating all the higher moments that appear in Equation ( 5), our algorithm decides that the inequality in (5) holds for the guess 𝜆 = 1 and returns true ‡ .

Final results. Our approach then continues with the binary search and obtains an optimal value of 𝜆 * ≈ 1.364 for 𝜆. Thus, it successfully synthesizes 𝑡 (𝑛) := 1.364 𝑛 . By plugging 𝑡 (𝑛) back into our core theorem for exponential bounds (Theorem 3.5), our approach synthesizes the following concentration bound:

Pr[𝑇 (𝑛) ≥ 𝛼 • 𝑛] ≤ exp(-1.364 • 𝛼 + 8.184).
Comparison with Karp's cookbook. In comparison, [START_REF] Karp | Probabilistic Recurrence Relations[END_REF]] provides Pr(𝑇 (𝑛) ≥ 𝛼 • 𝑛) ≤ (3/4) 𝛼-4 . Thus, our approach produces a much tighter bound than Karp's method. For 𝛼 = 20, our method synthesizes the bound 5.089 • 10 -9 , while Karp's method outputs 0.01. In other words, our method beats [START_REF] Karp | Probabilistic Recurrence Relations[END_REF]] by a ratio of 1.984 • 10 6 .

Extension to polynomially decreasing bounds. Our approach is not limited to exponential concentration bounds. If we consider the constraint in Theorem 3.6 instead of Theorem 3.5, then every step of the algorithm above can be adapted with minimal changes in order to synthesize a polynomially decreasing concentration bound. See Section 5.4 for details. As an example, our approach is able to synthesize the following bound for S LP with dimension 𝑑 = 2:

Pr[C 𝜏 ≥ 𝛼 • 8 • 𝑛 * ] ≤ 9.33 • 𝛼 -7
In the constrast, since the problem involves stochastic cost fucntion and a system of probabilistic recurrence relations, it is beyond the theoretical framework of Karp's method.

ALGORITHMIC APPROACHES

In this section, we provide algorithmic approaches for automatic synthesis of a TSFS for exponential and polynomial T-functions. Our algorithmic approaches are based on Composition Theorems (Theorems 3.5 and 3.6) for exponential and polynomial T-functions. Additionally, we make extra assumptions on our input SPRRs. While these assumptions mean that our algorithmic approaches are not as general as the theoretical foundation presented in the previous section, they nevertheless enable us to obtain e cient algorithms that can handle real-world SPRRs and obtain tighter concentration bounds. We rst recall the syntax of SPRRs.

SPRRs. Recall that an SPRR Δ consists of 𝑚 recurrence equations for which the 𝑟 -th equation (1 ≤ 𝑟 ≤ 𝑚) in the SPRR is of the following format:

𝑇 𝑟 (𝑛) = Cost 𝑟 (𝑛) + 𝑆 𝑟 𝑗=1 𝑇 𝑟 (ℎ 𝑟,𝑗 (𝑛)) + 𝑂 𝑟 𝑗=1 𝑇 𝑞 𝑟,𝑗 (𝑔 𝑟,𝑗 (𝑛)) if 𝑛 ≥ 𝑐 𝑟 + 1 0 if 𝑛 ≤ 𝑐 𝑟 (eq 𝑟 )
where 𝑇 𝑟 represents the 𝑟 -th procedure, 𝑐 𝑟 is the non-negative integer constant for the threshold of termination of procedure 𝑇 𝑟 , Cost 𝑟 (𝑛) is the random amount of preprocessing time, 𝑇 𝑟 (ℎ 𝑟,𝑗 (𝑛))'s (1 ≤ 𝑗 ≤ 𝑆 𝑟 ) are the 𝑆 𝑟 procedure calls to 𝑇 𝑟 itself with random sizes ℎ 𝑟,𝑗 (𝑛), and 𝑇 𝑞 𝑟,𝑗 (𝑔 𝑟,𝑗 (𝑛))'s (1 ≤ 𝑗 ≤ 𝑂 𝑟 ) are the procedure calls to procedures 𝑇 𝑞 𝑟,𝑗 other than 𝑇 𝑟 with random sizes 𝑔 𝑟,𝑗 (𝑛). We assume that the 𝑚-th procedure is the main procedure and model the initial input size by 𝑛 * . Since many randomized algorithms have expected runtimes in the form of polynomials in the input size 𝑛 and its logarithm ln 𝑛, we consider concentration bounds that also involve these expressions. In this work, we refer to these polynomial-like expressions as pseudo-polynomials.

Pseudo-polynomials. A (univariate) pseudo-polynomial (in the input-size variable 𝑛) is simply a polynomial over 𝑛 and ln 𝑛. Formally, it is a function of the form 𝑝 (𝑛) = 𝑖,𝑗 ∈Z 𝑑 𝑖 𝑗 • 𝑛 𝑖 • ln 𝑗 𝑛 where 𝑑 𝑖,𝑗 's are real numbers. We also call an expression of the form 𝑑 𝑖 𝑗 • 𝑛 𝑖 • ln 𝑗 𝑛 a pseudo-monomial. Our goal is then to automatically derive exponentially and polynomially decreasing concentration ‡ Checking the inequality is not straightforward. We also provide a dedicated procedure for this step. The details are removed in order to keep this illustration high-level. bounds in the form of Pr Δ,𝜎 * (C 𝜏 ≥ 𝜅) where the expression 𝜅 is assumed to be a pseudo-monomial.

Assumptions on the SPRR structure. While our general theory (Section 3) is applicable to all SPRRs, for the sake of obtaining more e cient algorithms, we assume the input SPRRs in this section satisfy the following conditions:

• Non-mutuality: We prohibit mutual recursion and assume that the call graph among our procedures is acyclic, except for potentially the recursive calls made by a procedure to itself. Formally, for every 1 ≤ 𝑟 ≤ 𝑚 and 1 ≤ 𝑗 ≤ 𝑂 𝑟 , we assume that 𝑞 𝑟,𝑗 is in {1, . . . , 𝑟 -1}. • Ordinariness: We assume that the recursive calls are ordinary. In other words, in each procedure that calls itself recursively, we assume that either there is only one such call, or there are two such calls which correspond to a divide-and-conquer approach and the sum of sizes in the two calls are exactly 𝑛 -1. Formally, we either have 𝑆 𝑟 = 1, or 𝑆 𝑟 = 2 and ℎ 𝑟,1 (𝑛) +ℎ 𝑟,2 (𝑛) = 𝑛 -1. Note that most real-world SPRRs are ordinary. For example, S has only one recursive call to itself, whereas S ts into our divide-and-conquer category, as it makes two recursive calls to itself and the sum of sizes passed to the two calls is 𝑛 -1.

Assumptions on probability distributions. Any algorithmic approach to concentration bounds for SPRRs has to handle probability distributions in the recurrences. Aspects of a recurrence that can be probabilistic include the preprocessing cost function Cost 𝑟 (𝑛) and the sizes passed to recursive function calls, namely ℎ 𝑟,𝑗 (𝑛) and 𝑔 𝑟,𝑗 (𝑛). In the sequel, we restrict our attention to only those distributions that are used in our real-world SPRRs.

Formally, we assume that Cost 𝑟 (𝑛) is sampled from a discrete distribution with a nite support whose every element and the probability assigned to each element in this distribution is a pseudopolynomial with rational coe cients.

As mentioned before, if 𝑆 𝑟 = 2, then ℎ 𝑟,2 (𝑛) = 𝑛 -1ℎ 𝑟,1 (𝑛). Hence, we only have to specify probability distributions for 𝑔 𝑟,𝑗 (𝑛)'s and ℎ 𝑟,1 (𝑛). In the sequel, we assume that the random variable ℎ 𝑟,1 (𝑛) is indepedent with others. For ℎ 𝑟,1 (𝑛), our algorithm supports three types of distribution:

• A discrete distribution with a nite support whose every element is eithor an a ne expression with natural coe cients, or in the form (𝑛 + 𝑏)/2 for some integer 𝑏. For example, the support may contain 𝑛 -1, 2, 𝑛 2 or 𝑛-1 2 . The probability assigned to each element is a pseudo-polynomial.

• A uniform distribution over {0, 1, . . . , 𝑛 -1}, • max{𝑖, 𝑛 -𝑖 -1} where 𝑖 is distributed uniformly over {0, 1, . . . , 𝑛 -1}. For 𝑔 𝑟,𝑗 (𝑛)'s, we suppose that they observe as the same type of discrete distributions as the rst case for ℎ 𝑟,1 (𝑛).

We rst present some prerequisites of our algorithms in Sections 5.1 and 5.2 and then proceed with the main algorithms in Sections 5.3 and 5.4.

Approximations by Pseudo-Polynomials

Approximating higher-order moments. In several of our algorithmic steps in the sequel, we will need to approximate the higher-order moments of 𝑓 (𝑖), where 𝑓 is a pseudo-polynomial and either 𝑖 observes a uniform distribution or 𝑖 = max{ 𝑗, 𝑛 -1 -𝑗 } and 𝑗 observes a uniform distribution. To achieve this, we need to tightly approximate the summation

𝑛-1 𝑖=0 𝑖 𝑎 • ln 𝑏 𝑖 (6)
by a pseudo-polynomial. We now provide a short overview of how this approximation can be performed. For more details, see Appendix E. Our approximation method is M A (𝐿 𝑇 , 𝑎, 𝑏, 𝑛) to approximate (6). It approximates the summation through integrals and Taylor expansions. The parameters 𝑎 and 𝑏 represent the degrees for 𝑖 and ln 𝑖, respectively, and 𝐿 𝑇 is a parameter chosen by the user that represents the degree of the approximation. To avoid corner cases, our algorithm assumes 𝑛 ≥ exp(𝑏/𝐿 𝑇 ) + 1. The algorithm's goal is to synthesize pseudo-polynomial upper and lower bounds for (6) with constant additive error. To approximate (6), our algorithm uses its corresponding integral ∫ 𝑛 𝑑 𝑥 𝑎 • ln 𝑏 𝑥 d𝑥, since the latter could be automatically computed as a pseudo-polynomial through integration by parts. Then, our algorithm would measure the di erence between the summation and integral through Taylor's expansion with degree 𝑎 + 𝐿 𝑇 at each term. The algorithm controls the additive error for minor terms using the Riemann Zeta function.

Approximating logarithms of a ne expressions. We now present a proposition that helps us approximate ln 𝑣 (𝑎 • 𝑥 + 𝑏) by a pseudo-polynomial: P 5.1. For every 𝑎, 𝑏, 𝑥 > 0, the following inequalities hold:

ln 𝑎 + ln 𝑥 ≤ ln(𝑎 • 𝑥 + 𝑏) ≤ ln 𝑎 + ln 𝑥 + 𝑏 𝑎 •𝑥
Moreover, for every 𝑎 > 0, 𝑏 < 0, and 𝑥 ≥ -2•𝑏 𝑎 , we have:

ln 𝑎 + ln 𝑥 + 2•𝑏 𝑎 •𝑥 ≤ ln(𝑎 • 𝑥 + 𝑏) ≤ ln 𝑎 + ln 𝑥 + 𝑏 𝑎 •𝑥 .
A proof is provided in Appendix G. Using this proposition, to approximate ln 𝑣 (𝑎 • 𝑥 + 𝑏), we can rst approximate ln(𝑎 • 𝑥 + 𝑏) by the proposition above and then calculate its 𝑣-th power (and expand it by applying the binomial theorem).

Checking Non-positivity of Pseudo-polynomials

In both of our algorithms (Sections 5.4 and 5.3), we need to check whether a given pseudo-polynomial 𝑃 is non-positive for every 𝑛 ≥ 𝑛 0 . We now provide a sketch of our algorithm for handling this problem. A detailed description is provided in Appendix F.

Our algorithm rst obtains an equivalent form 𝑃 1 (𝑛) of 𝑃 (𝑛) by dividing everything by the term with the largest order. Then, it naively tries to iterate over every 𝑛 ≥ 𝑛 0 and evaluate 𝑃 1 (𝑛). This naive approach will of course take in nite time. Thus, we propose a condition for ending the iteration. The condition considers an over-approximation 𝑃 2 (𝑛) of 𝑃 1 (𝑛) by removing all nonconstant terms with negative coe cients. It also tries to nd a constant 𝑛 𝑒 such that 𝑃 2 (𝑛) decreases monotonically for every 𝑛 ≥ 𝑛 𝑒 . Intuitively, this means that if 𝑃 2 (𝑛 1 ) ≤ 0 for some 𝑛 1 ≥ 𝑛 𝑒 , then 𝑃 2 (𝑛) ≤ 0 for every 𝑛 ≥ 𝑛 1 . Combining this with the fact that 𝑃 2 is an over-approximation of 𝑃 1 , we can conclude that 𝑃 1 (𝑛) and 𝑃 (𝑛) are also non-positive for every 𝑛 ≥ 𝑛 1 . Hence, we can safely terminate the iteration after 𝑛 1 . Note that our algorithm is complete in the sense that if 𝑃 (𝑛) is always non-positive, the algorithm will eventually terminate and report this property.

Our Algorithm for Synthesizing Exponentially Decreasing Bounds

In this section, we provide a detailed explanation of the algorithm used in Section 4 to obtain concentration bounds for S . While the main ideas are the same as in Section 4, we now present them in a more general format, so that they are applicable to a wider family of SPRRs. To facilitate switching between the concrete example and the abstract explanation of the algorithm, our paragraph names in this section correspond to those of Section 4.

Sketch of the algorithm. Our algorithmic approach to exponentially decreasing bounds is based on Theorem 3.5. Recall that to apply this theorem, it su ces to synthesize 2 •𝑚 functions 𝑓 1 , . . . , 𝑓 𝑚 and 𝑡 1 , . . . , 𝑡 𝑚 that satisfy the condition in Equation ( ‡) below. Then, the bound can easily be computed by relying on the T-function 𝑓 (𝑥) = exp(𝑡 ★ • 𝑥) where 𝑡 ★ := min 1≤𝑟 ≤𝑚 min 1≤𝑛 ≤𝑛 * 𝑡 𝑟 (𝑛). We assume that the 𝑡 𝑟 's and 𝑓 𝑟 's are in a pseudo-monomial form and use a template-based method for our synthesis. As we will see, the core of the problem is to nd a value 𝜆 𝑟 , which is the coe cient in template of 𝑡 𝑟 , for each equation 1 ≤ 𝑟 ≤ 𝑚 in our SPRR. Since the SPRR is assumed to be non-mutual, we can nd these values in order from 𝜆 1 to 𝜆 𝑚 . For each 𝑟 , our goal is to nd the largest possible 𝜆 𝑟 that satis es ( ‡). Therefore, if we have a veri cation algorithm that, given 𝑡 1 , . . . , 𝑡 𝑟 and 𝑓 1 , . . . , 𝑓 𝑟 decides whether ( ‡) is satis ed, then we can nd 𝜆 𝑟 by a simple binary search. Below, we present the details of each step separately.

Intuition behind the condition ( ‡). Let us consider Theorem 3.5 again. We are looking for 2 • 𝑚 functions 𝑓 1 , . . . , 𝑓 𝑚 , 𝑡 1 , . . . , 𝑡 𝑚 : N → [0, ∞) such that:

exp(𝑡 𝑟 (𝑛) • 𝑓 𝑟 (𝑛)) ≥ E exp 𝑡 𝑟 (𝑛) • Cost 𝑟 (𝑛) + 𝑆 𝑟 𝑖=1 𝑓 𝑟 (ℎ 𝑟,𝑖 (𝑛)) + 𝑂 𝑟 𝑖=1 𝑓 𝑞 𝑟,𝑖 (𝑔 𝑟,𝑖 (𝑛)) ( ‡)
Intuitively, 𝑓 𝑟 is an over-approximation of the expected running time of the 𝑟 -th procedure and 𝑡 𝑟 is a scaling factor. Note that the 𝑓 𝑟 's can be obtained by previous automated methods for bounding expected runtimes, such as the one in [Chatterjee et al. 2017a]. Hence, our focus is on synthesizing the 𝑡 𝑟 's.

Rephrasing the condition. If we de ne

𝑋 𝑟 (𝑛) := -𝑓 𝑟 (𝑛) + Cost 𝑟 (𝑛) + 𝑆 𝑟 𝑖=1 𝑓 𝑟 (ℎ 𝑟,𝑖 (𝑛)) + 𝑂 𝑟 𝑖=1 𝑓 𝑞 𝑟,𝑖 (𝑔 𝑟,𝑖 (
𝑛)), then ( ‡) can simply be written as

E [exp(𝑡 𝑟 (𝑛) • 𝑋 𝑟 (𝑛))] ≤ 1.
As we will see in the sequel, we need to reason about the range of 𝑋 𝑟 in our veri cation algorithm.

Se ing up templates. In the rst step, the algorithm sets up a monomial template for every 𝑡 𝑟 and 𝑓 𝑟 . In other words, it assumes that each of the functions is of the following form:

• 𝑡 𝑟 (𝑛) = 𝜆 𝑟 • 𝑛 -𝑎 𝑟 • ln -𝑏 𝑟 𝑛 • 𝑓 𝑟 (𝑛) = 𝑤 𝑟 • 𝑛 𝑢 𝑟 • ln 𝑣 𝑟 𝑛
We assume that the values of 𝑎 𝑟 , 𝑏 𝑟 , 𝑢 𝑟 , 𝑣 𝑟 , and 𝑤 𝑟 are given and the goal is to synthesize 𝜆 𝑟 for every 1 ≤ 𝑟 ≤ 𝑚. This does not a ect the generality of our approach, since 𝑢 𝑟 , 𝑣 𝑟 and 𝑤 𝑟 can be automatically synthesized by previous methods that focus on bounding expected runtimes, such as [Chatterjee et al. 2017a], and 𝑎 𝑟 and 𝑏 𝑟 are typically small integers and one can simply try all possible combinations in which they are bounded by a small threshold 𝑐. Moreover, the algorithm also sets up a template for the range of 𝑋 𝑟 :

• 𝐵 𝑟 (𝑛) = 𝛾 𝑟 • 𝑛 𝑝 𝑟 • ln 𝑞 𝑟 𝑛
In other words, we assume that the range of 𝑋 𝑟 (𝑛) is bounded by a pseudo-polynomial function 𝐵 𝑟 (𝑛). Since 𝑋 𝑟 does not depend on any of the 𝑡 𝑗 's, computing the constants 𝛾 𝑟 , 𝑝 𝑟 and 𝑞 𝑟 can be automated. However, this is not the main focus of our algorithm. Hence, we have relegated the details to Appendix H.

Binary searching. Our goal is to synthesize 𝑡

𝑟 's that satisfy E [exp(𝑡 𝑟 (𝑛) • 𝑋 𝑟 (𝑛))] ≤ 1 for every 𝑟 . Each 𝑡 𝑟 is of the form 𝑡 𝑟 (𝑛) = 𝜆 𝑟 • 𝑛 -𝑎 𝑟 • ln -𝑏 𝑟 𝑛,
where the only unknown is 𝜆 𝑟 . Note that our SPRR is non-mutual. Therefore, the algorithm rst synthesizes 𝜆 1 , then uses it to synthesize 𝜆 2 and so on. When synthesizing 𝜆 𝑟 , we can assume that 𝜆 1 , . . . , 𝜆 𝑟 -1 are already known. The synthesis process itself is quite simple. The algorithm relies on a veri cation procedure that, given values for 𝜆 1 , . . . , 𝜆 𝑟 , decides whether ( ‡) holds. Using this veri cation procedure, the algorithm performs a binary search to nd the largest possible value for 𝜆 𝑟 .

Veri cation procedure. This is the most involved part of our algorithm. Given the values for 𝜆 1 , . . . , 𝜆 𝑟 , the goal is to decide whether E [exp(𝑡 𝑟 (𝑛) • 𝑋 𝑟 (𝑛))] ≤ 1 for every 𝑛 ≥ 𝑐 𝑟 + 1. It is well-known that there exists 𝜉 ∈ [0, 𝑥] such that

exp(𝑥) = 1 + 𝐿 𝑖=1 𝑥 𝑖 𝑖! + exp(𝜉) (𝐿+1)! 𝑥 𝐿+1
If -𝐵 ≤ 𝑥 ≤ 𝐵 and 𝐿 is an odd number, then 𝑥 𝐿+1 is always non-negative. Hence, we have:

exp(𝑥) ≤ 1 + 𝐿 𝑖=1 𝑥 𝑖 𝑖! + exp(𝐵) (𝐿+1)! 𝑥 𝐿+1 .
Let us apply this inequality to exp(𝑡 𝑟 (𝑛) • 𝑋 𝑟 (𝑛)). Suppose that -𝐷 𝑟 ≤ 𝑡 𝑟 (𝑛) • 𝑋 𝑟 (𝑛) ≤ 𝐷 𝑟 for some 𝐷 𝑟 > 0. Then, we have:

exp(𝑡 𝑟 (𝑛) • 𝑋 𝑟 (𝑛)) ≤ 1 + 𝐿 𝑖=1 𝑡 𝑟 (𝑛) 𝑖 𝑖! • 𝑋 𝑟 (𝑛) 𝑖 + exp(𝐷 𝑟 ) •𝑡 𝑟 (𝑛) 𝐿+1 (𝐿+1)! • 𝑋 𝑟 (𝑛) 𝐿+1 .
Now, we can strengthen our target inequality E [exp(𝑡 𝑟 (𝑛) • 𝑋 𝑟 (𝑛))] ≤ 1 to:

E 1 + 𝐿 𝑖=1 𝑡 𝑟 (𝑛) 𝑖 𝑖! • 𝑋 𝑟 (𝑛) 𝑖 + exp(𝐷 𝑟 ) •𝑡 𝑟 (𝑛) 𝐿+1 (𝐿+1)! • 𝑋 𝑟 (𝑛) 𝐿+1 ≤ 1,
which is equivalent to:

𝐿 𝑖=1 𝑡 𝑟 (𝑛) 𝑖 𝑖! • E 𝑋 𝑟 (𝑛) 𝑖 + exp(𝐷 𝑟 ) •𝑡 𝑟 (𝑛) 𝐿+1 (𝐿+1)! • E 𝑋 𝑟 (𝑛) 𝐿+1 ≤ 0.
by subtracting 1 on both side. Our veri cation procedure tries to establish the last inequality. To achieve this goal, our algorithm performs three steps:

(1) Synthesize the bound 𝐷 𝑟 .

(2) Over-approximate each expectation term E 𝑋 𝑟 (𝑛) 𝑖 for large enough 𝑛 by a pseudo-polynomial expression 𝐴 𝑟,𝑖 (𝑛). § (3) Check whether the strengthened inequality

𝐿 𝑖=1 𝑡 𝑟 (𝑛) 𝑖 𝑖! • 𝐴 𝑟,𝑖 (𝑛) + exp(𝐵 𝑟 ) •𝑡 𝑟 (𝑛) 𝐿+1 (𝐿+1)!
• 𝐴 𝑟,𝐿+1 (𝑛) ≤ 0 holds for every 𝑛 ≥ max{ℓ, 𝑐 𝑟 + 1}.

Step

(𝑛) = 𝛾 𝑟 • 𝜆 𝑟 • 𝑛 𝑝 𝑟 -𝑎 𝑟 • ln 𝑞 𝑟 -𝑏 𝑟 𝑛.
It is easy to see that we can use the maximum of 𝑉 (𝑛) as 𝐷 𝑟 . Due to the monotonicity of pseudomonomials, the algorithm considers the following cases:

• If 𝑝 𝑟 > 𝑎 𝑟 or 𝑝 𝑟 = 𝑎 𝑟 and 𝑞 𝑟 > 𝑏 𝑟 , then the algorithm outputs 𝐷 𝑟 = +∞.

• If 𝑝 𝑟 ≤ 𝑎 𝑟 and 𝑞 𝑟 ≤ 𝑏 𝑟 , then 𝑉 (𝑛) reaches its maximum at 𝑛 = 𝑐 𝑟 + 1 and the algorithm outputs 𝐷 𝑟 = 𝑉 (𝑐 𝑟 + 1).

• If 𝑝 𝑟 ≤ 𝑎 𝑟 -1 and 𝑞 𝑟 ≥ 𝑏 𝑟 + 1, then 𝑉 (𝑛) reaches its maximum at max 𝑐 𝑟 + 1,

𝑞 𝑟 -𝑏 𝑟 𝑎 𝑟 -𝑝 𝑟 . Hence, the algorithm outputs 𝐷 𝑟 = 𝑉 max 𝑐 𝑟 + 1, 𝑞 𝑟 -𝑏 𝑟 𝑎 𝑟 -𝑝 𝑟 .
Step (2). Over-approximating expectation terms. In this step, the algorithm over-approximates each expectation term E 𝑋 𝑟 (𝑛) 𝑖 by a pseudo-polynomial expression 𝐴 𝑟,𝑖 (𝑛). Recall that

𝑋 𝑟 (𝑛) = Cost 𝑟 (𝑛) + 𝑆 𝑟 𝑖=1 𝑓 𝑟 (ℎ 𝑟,𝑖 (𝑛)) + 𝑂 𝑟 𝑖=1 𝑓 𝑞 𝑟,𝑖 (𝑔 𝑟,𝑖 (𝑛)) -𝑓 𝑟 (𝑛).
The algorithm expands E 𝑋 𝑟 (𝑛) 𝑖 by simply applying the binomial theorem to the RHS of the equality above. This leads to:

E 𝑋 𝑟 (𝑛) 𝑖 = I 𝑖! • (-1) 𝑏 •𝑓 𝑟 (𝑛) 𝑏 𝑏! • E[Cost 𝑟 (𝑛) 𝑗 ] 𝑗! • E 𝑆 𝑟 𝑢=1 𝑓 𝑝𝑢 𝑟 (ℎ 𝑟,𝑢 (𝑛)) 𝑝 𝑢 ! • 𝑂 𝑟 𝑢=1 E 𝑓 𝑞𝑢 𝑞𝑟,𝑢 (𝑔 𝑟,𝑢 (𝑛))
𝑞 𝑢 ! § Since the approximation does not hold for small 𝑛, our algorithm rst nds an integer ℓ such that the approximation is sound for every 𝑛 ≥ ℓ and then simply iterates over 𝑐 𝑟 + 1 ≤ 𝑛 ≤ ℓ, evaluates the expression and trivially veri es that it satis es the condition.

where I denotes 0 ≤ 𝑏, 𝑗, 𝑝 1 , . . . , 𝑝 𝑆 𝑟 , 𝑞 1 , . . . , 𝑞 𝑂 𝑟 ≤ 𝑖, 𝑏 + 𝑗 + 𝑆 𝑟 𝑢=1 𝑝 𝑢 + 𝑂 𝑟 𝑢=1 𝑞 𝑢 = 𝑖 . Each term is treated separately. In order to obtain an over-approximation for the LHS, the algorithm over/underapproximates each term of the sum on the RHS based on the parity of 𝑏 in (-1) 𝑏 .

For each term, the algorithm has to estimate (i)

E Cost 𝑟 (𝑛) 𝑗 • 𝑂 𝑢 𝑢=1 𝑓 𝑞𝑢 𝑞𝑟,𝑢 (𝑔 𝑟,𝑢 (𝑛)) 𝑞 𝑢 !
, and (ii)

E 𝑆 𝑟 𝑢=1 𝑓 𝑝𝑢 𝑟 (ℎ 𝑟,𝑢 (𝑛)) 𝑝 𝑢 !
. To approximate (i), our algorithm performs the following steps: (Recall that all random variables in (i) observe discrete distributions):

• It iterates over every possible combination of values for Cost 𝑟 (𝑛) and the 𝑔 𝑟,𝑢 (𝑛)'s in their respective discrete distributions and approximates Cost 𝑟 (𝑛) 𝑗 and each 𝑓 𝑞 𝑟,𝑢 (𝑔 𝑟,𝑢 ) (𝑛) 𝑞 𝑢 by a pseudo-polynomial -Cost 𝑟 (𝑛) 𝑗 can represented exactly and straightforwardly by a pseudo-polynomial.

-Whenever an element in the support involves (𝑛 + 𝑏)/2 , our algorithm will distinguish between even and odd 𝑛, and transform the oor/ceiling term by an a ne expression precisely. -Now, each 𝑔 𝑟,𝑢 (𝑛) is in the form of an a ne expression 𝑝 •𝑛 +𝑞 (𝑝, 𝑞 ∈ Q). Plugging this into the template of 𝑓 , the algorithm derives a combination of terms of the form ln 𝑣 (𝑝 • 𝑛 + 𝑞).

It then applies Proposition 5.1 and approximates these terms by pseudo-polynomials. Note that if 𝑞 < 0, the proposition only applies for 𝑛 ≥ -2•𝑞 𝑝 . In this case, the algorithm sets ℓ := max{ℓ, -2•𝑞 𝑝 } and manually checks every 𝑐 𝑟 + 1 ≤ 𝑛 ≤ ℓ in order to handle the corner case. We now focus on how the algorithm estimates (ii). The treatment depends on whether 𝑆 𝑟 = 1 or 𝑆 𝑟 = 2, as well as the probability distributions. Based on the assumptions made at the beginning of this section, the algorithm considers the following six cases:

• 𝑆 𝑟 = 1 ℎ 𝑟,1 (𝑛) has a discrete distribution, we follow exactly the same procedure for approximating (i). ℎ 𝑟,1 (𝑛) is distributed uniformly: In this case, the algorithm needs to approximate

1 𝑛 • 𝑛-1 𝑖=0 𝑖 𝑢 𝑟 •𝑝 1 ln 𝑣 𝑟 •𝑝 1 𝑖 Relying on M A (𝐿 𝑇 , 𝑢 𝑟 •𝑝 1 , 𝑣 𝑟 •𝑝 1 , 𝑛)
, the algorithm computes pseudo-polynomial bounds for this term. Note that M A requires the assumption 𝑛 ≥ 1+ exp(

𝑣 𝑟 •𝑝 1
𝐿 𝑇 ) , thus the algorithm sets ℓ := max{ℓ, 1 + exp(

𝑣 𝑟 •𝑝 1 𝐿 𝑇 ) } and handles each 𝑐 𝑟 + 1 ≤ 𝑛 ≤ ℓ by evaluation. -ℎ 𝑟,1 (𝑛) = max{𝑖, 𝑛 -1 -𝑖} where 𝑖 is uniformly distributed: In this case, the algorithm should approximate 1 𝑛 • 𝑛-1 𝑖= (𝑛-1)/2 𝑖 𝑢 𝑟 •𝑝 1 • ln 𝑣 𝑟 •𝑝 1 𝑖 + 𝑛-1 𝑖= (𝑛-1)/2 𝑖 𝑢 𝑟 •𝑝 1 • ln 𝑣 𝑟 •𝑝 1 𝑖 .
It distinguishes between even and odd 𝑛 to remove oors and ceils, and then handles the two summations separately and exactly as in the previous case, applying M A . Also we need to incorporate corner cases and updates ℓ.

• 𝑆 𝑟 = 2 and ℎ 𝑟,1 (𝑛) + ℎ 𝑟,2 (𝑛) = 𝑛 -1 : In all the cases below, the algorithm over-estimates ℎ 𝑟,2 (𝑛) 𝑎 𝑟 • ln 𝑏 𝑟 ℎ 𝑟,2 (𝑛) as ℎ 𝑟,2 (𝑛) 𝑎 𝑟 • ln 𝑏 𝑟 𝑛. This is sound since it leads to an over-approximation of 𝑋 𝑟 (𝑛) and thus strengthens the condition.

ℎ 𝑟,1 (𝑛) has a discrete distribution. This case is handled similarly to when 𝑆 𝑟 = 1. The only di erence is that now the algorithm has to consider both recursive calls and should consider:

(𝑝 • 𝑛 + 𝑞) 𝑢 𝑟 •𝑝 1 • ln 𝑣 𝑟 •𝑝 1 (𝑝 • 𝑛 + 𝑞) • (𝑛 -1 -𝑝 • 𝑛 -𝑞) 𝑢 𝑟 •𝑝 2 • ln 𝑣 𝑟 •𝑝 2 𝑛.
ℎ 𝑟,1 (𝑛) is distributed uniformly: In this case, the algorithm has to approximate

1 𝑛 • 𝑛-1 𝑖=0 𝑖 𝑢 𝑟 •𝑝 1 • ln 𝑣 𝑟 •𝑝 1 𝑖 • (𝑛 -𝑖 -1) 𝑢 𝑟 •𝑝 2 • ln 𝑣 𝑟 •𝑝 1 𝑛.
To handle this expression, the algorithm rst applies binomial expansion and then uses M A as in the case where 𝑆 𝑟 = 1., and considers corner cases similar to previous cases.

ℎ 𝑟,1 (𝑛) = max{𝑖, 𝑛 -𝑖 -1} where 𝑖 is distributed uniformly: The algorithm should approximate

1 𝑛 • 𝑛-1 𝑖= (𝑛-1)/2 𝑖 𝑢 𝑟 •𝑝 1 • ln 𝑣 𝑟 •𝑝 1 𝑖 • (𝑛 -𝑖 -1) 𝑢 𝑟 •𝑝 2 • ln 𝑣 𝑟 •𝑝 1 𝑛 + 𝑛-1 𝑖= (𝑛-1)/2 𝑖 𝑢 𝑟 •𝑝 1 • ln 𝑣 𝑟 •𝑝 1 𝑖 • (𝑛 -𝑖 -1) 𝑢 𝑟 •𝑝 2 • ln 𝑣 𝑟 •𝑝 1 𝑛 .
It considers each sum separately and applies binomial expansion followed by M A as in the previous case.

Step (3). Checking the strengthened inequality. Finally, after nding the over-approximations in the previous steps, the algorithm considers the strengthened inequality

𝐿 𝑖=1 𝑡 𝑟 (𝑛) 𝑖 𝑖! • 𝐴 𝑟,𝑖 (𝑛) + exp(𝐵 𝑟 ) •𝑡 𝑟 (𝑛) 𝐿+1 (𝐿+1)! • 𝐴 𝑟,𝐿+1 (𝑛) ≤ 0
and computes it symbolically. The goal of this step is to decide whether this inequality holds for every 𝑛 ≥ max{ℓ, 𝑐 𝑟 + 1}. This step is achieved as outlined in Section 5.2 and Appendix F.

Our Algorithm for Synthesizing Polynomially Decreasing Bounds

Sketch of the algorithm. Our algorithm for synthesis of polynomial bounds is quite similar to the algorithm in Section 5.3, except that it considers a di erent set of constraints, i.e. the constraints of Theorem 3.6. To apply this theorem, we would need to synthesize 𝑚 • (𝑘 + 1) functions 𝑓 1,0 , 𝑓 where I := 0 ≤ 𝑗, 𝑝 1 , . . . , 𝑝 𝑆 𝑟 , 𝑒 1 , . . . , 𝑒 𝑂 𝑟 ≤ 𝑖, 𝑆 𝑟 𝑢=1 𝑝 𝑢 + 𝑂 𝑟 𝑢=1 𝑒 𝑢 + 𝑗 = 𝑖 . As in the previous algorithm, we start by setting up a template for each function 𝑓 𝑟,𝑖 . The core of the problem is then to nd values for coe cients 𝜆 𝑟,𝑖 that appear in the templates for 𝑓 𝑟,𝑖 . Given the non-mutuality of the input SPRR, we can solve the 𝜆 𝑟,𝑖 's in lexicographic order. As before, we design a veri cation algorithm that gets values for 𝜆 1,0 , 𝜆 1,1 , . . . , 𝜆 1,𝑘 , 𝜆 2,0 , . . . , 𝜆 𝑟,𝑖 as input and decides whether (7) is violated. We then rely on this veri cation procedure to nd each 𝜆 𝑟,𝑖 using binary search.

Intuition behind the functions and templates. In the constraints above, each 𝑓 𝑟,𝑖 (𝑛) intuitively serves as an over-approximation of E[𝑇 𝑟 (𝑛) 𝑖 ]. This, together with the fact that most randomized algorithms considered in this work have pseudo-polynomial runtimes, naturally leads to the idea of establishing pseudo-monomial templates for each 𝑓 𝑟,𝑖 .

Se ing up templates. In the rst step, the algorithm sets up a pseduo-monomial template for each 𝑓 𝑟,𝑖 . In other words, it assumes that each 𝑓 𝑟,𝑖 is in the following form:

• 𝑓 𝑟,𝑖 (𝑛) = 𝜆 𝑟,𝑖 • 𝑛 𝑢 𝑟,𝑖 • ln 𝑣 𝑟,𝑖 𝑛. As in Section 5.3, we assume that 𝑢 𝑟,𝑖 and 𝑣 𝑟,𝑖 are given and focus on synthesizing 𝜆 𝑟,𝑖 . Note that 𝜆 𝑟,𝑖 must be positive. As discussed at the end of Section 3, our algorithm will synthesize the bound 𝜆 𝑟,𝑖 • 𝜅 -𝑘 • (𝑛 * ) 𝑢 𝑟 * ,𝑘 • (ln 𝑛 * ) 𝑣 𝑟 * ,𝑘 . Given that 𝜅 is a pseudo-polynomial over 𝑛 * , this bound is polynomially decreasing.

Binary searching. Given that the SPRR is assumed to be non-mutual, we can solve the 𝜆 𝑟,𝑖 's in lexicographic order 𝜆 1,0 , 𝜆 1,1 , . . . , 𝜆 1,𝑘 , 𝜆 2,0 , . . . , 𝜆 2,𝑘 , . . . , 𝜆 𝑚,𝑘 . Assuming that we have a veri cation procedure (details below) that, given values for 𝜆 1,0 , . . . , 𝜆 𝑟,𝑖 , decides whether (7) holds, the algorithm performs a binary search on 𝜆 𝑟,𝑖 . Since no upper-bound is known apriori, it rst keeps doubling its guess until an upper-bound is found and then proceeds with classical binary search.

Veri cation procedure. As input, our veri cation procedure gets a guessed value for 𝜆 𝑟,𝑖 , as well as xed values for all 𝜆 𝑟 ,𝑖 's that come lexicographically before it. It should then check whether they satisfy

𝑓 𝑟,𝑖 (𝑛) 𝑖! ≥ E I 𝐶𝑜𝑠𝑡 𝑟 (𝑛) 𝑗 𝑗! • 𝑆 𝑟 𝑢=1 𝑓 𝑟,𝑝𝑢 (ℎ 𝑟,𝑢 (𝑛)) 𝑝 𝑢 ! • 𝑂 𝑟 𝑢=1 𝑓 𝑞𝑟,𝑢 ,𝑒𝑢 (𝑔 𝑟,𝑢 (𝑛)) 𝑒 𝑢 !
.

where

I := 0 ≤ 𝑗, 𝑝 1 , . . . , 𝑝 𝑆 𝑟 , 𝑒 1 , . . . , 𝑒 𝑂 𝑟 ≤ 𝑖, 𝑆 𝑟 𝑢=1 𝑝 𝑢 + 𝑂 𝑟 𝑢=1 𝑒 𝑢 + 𝑗 = 𝑖 .
The process is similar to that of Section 5.3. The algorithm iterates over every summation term (in I) and over-approximates the term by a pseudo-polynomial. It then checks whether the strengthened inequality obtained by replacing each term with its over-approximation holds. Since checking this inequality is done in the exact same manner as in Section 5.3, we only focus on the over-approximation step.

Over-approximations. Based on linearity of expectation, for each summation term, the algorithm has to over-approximate (i) E 𝐶𝑜𝑠𝑡 𝑟 (𝑛) , by pseudo-polynomials. (i) is approximated using the same procedure as in Section 5.3. To approximate (ii), the algorithm considers six cases, based on the number of recursive calls to the 𝑟 -th procedure and the distribution of ℎ 𝑟,1 (𝑛).

• 𝑆 𝑟 = 1 ℎ 𝑟,1 (𝑛) follows a discrete distribution, we follow the same procedure for approximating (i).

ℎ 𝑟,1 (𝑛) is distributed uniformly, The algorithm considers 1 𝑛 • 𝑛-1 𝑖=0 𝑖 𝑢 𝑟 •𝑝 1 ln 𝑣 𝑟 •𝑝 1 𝑖. We could also follow the same procedure as we do for our algorithm in exponentially decreasing bounds.

ℎ 𝑟,1 (𝑛) = max{𝑖, 𝑛 -𝑖 -1}, where 𝑖 is distributed uniformly: The algorithm considers

1 𝑛 • 𝑛-1 𝑖= (𝑛-1)/2 𝑖 𝑢 𝑟 •𝑝 1 ln 𝑣 𝑟 •𝑝 1 𝑖 + 𝑛-1 𝑖= (𝑛-1)/2 𝑖 𝑢 𝑟 •𝑝 1 ln 𝑣 𝑟 •𝑝 1 𝑖
and applies the method in the previous case to each sum separately. • 𝑆 𝑟 = 2 : In the cases below, we always have ℎ 𝑟,2 (𝑛) = 𝑛 -1ℎ 𝑟,1 (𝑛). Moreover, the algorithm over-approximates ℎ 𝑟,2 (𝑛) 𝑎 𝑟 • ln 𝑏 𝑟 ℎ 𝑟,2 (𝑛) by ℎ 𝑟,2 (𝑛) 𝑎 𝑟 • ln 𝑏 𝑟 𝑛.

ℎ 𝑟,1 (𝑛) follows a discrete distribution: this follows exactly the same procedure as we do for our algorithm in exponentially decreasing bounds. ℎ 𝑟,1 (𝑛) is distributed uniformly: In this case, the algorithm should over-approximate

1 𝑛 • 𝑛-1 𝑖=0 𝑖 𝑢 𝑟 •𝑝 1 • ln 𝑣 𝑟 •𝑝 1 𝑖 • (𝑛 -𝑖 -1) 𝑢 𝑟 •𝑝 2 • ln 𝑣 𝑟 •𝑝 1 𝑛.
The over-approximation is obtained by rst applying a binomial expansion and then applying M A to each term. ℎ 𝑟,1 (𝑛) = max{𝑖, 𝑛 -𝑖 -1}, where 𝑖 is distributed uniformly: In this case, the algorithm combines the ideas from previous cases. It has to nd an upper-bound for

1 𝑛 • 𝑛-1 𝑖= (𝑛-1)/2 𝑖 𝑢 𝑟 •𝑝 1 • ln 𝑣 𝑟 •𝑝 1 𝑖 • (𝑛 -𝑖 -1) 𝑢 𝑟 •𝑝 2 • ln 𝑣 𝑟 •𝑝 1 𝑛 + 𝑛-1 𝑖= (𝑛-1)/2 𝑖 𝑢 𝑟 •𝑝 1 • ln 𝑣 𝑟 •𝑝 1 𝑖 • (𝑛 -𝑖 -1) 𝑢 𝑟 •𝑝 2 • ln 𝑣 𝑟 •𝑝 1 𝑛 .
It considers each sum separately and rst applies a binomial expansion. Then, it applies the M A procedure to each term.

Checking the strengthened inequality. Finally, the veri cation algorithm checks the strengthened inequality obtained by replacing every term in ( 7) by its over-approximation as a pseudopolynomial. The process is the same as in Section 5.3 and relies on the algorithm of Section 5.2

Benchmark SPRR Base case threshold(s)

L1D 𝑇 (𝑛) = 𝑛 + 𝑇 (U (𝑛)) 𝑐 = 1 L2D 𝑇 (𝑛) = 𝑛 • ln 𝑛 + 𝑇 (U (𝑛)) 𝑐 = 1 S 𝑇 (𝑛) = 𝑛 + 𝑇 (MU (𝑛)) 𝑐 = 1 S S 𝑇 1 (𝑛) = 𝑛 + 𝑇 1 (MU (𝑛)) 𝑇 2 (𝑛) = 𝑇 1 (𝑛) + 𝑇 2 ( (𝑛 -1)/2 ) + 𝑇 2 ( (𝑛 -1)/2 ) 𝑐 1 = 1, 𝑐 2 = 1 S 𝑇 (𝑛) = 𝑛 + 𝑇 (ℎ 1 (𝑛)) + 𝑇 (𝑛 -1 -ℎ 1 (𝑛)) where ℎ 1 (𝑛) is sampled from U (𝑛) 𝑐 = 1 R S 𝑇 (𝑛) = 1 + 𝑇 (MU (𝑛)) 𝑐 = 1 C 𝑇 (𝑛) = 1 + 𝑇 1 𝑒 : 𝑛 -1; 1 -1 𝑒 : 𝑛 𝑐 = 1 R 𝑇 (𝑛) = {0.5 : -0.05; 0.5 : 0.05} + 𝑇 (𝑛 -1) 𝑐 = 0 S LP in dimension 𝑑 𝑇 1 (𝑛) = 𝑛 𝑇 𝑟 (𝑛) = 𝑟 𝑛 : 𝑟 • 𝑛 + 𝑇 𝑟 -1 (𝑛 -1); 1 -𝑟 𝑛 : 𝑟 + 𝑇 𝑟 (𝑛 -1) (2 ≤ 𝑟 ≤ 𝑑) 𝑐 𝑟 = 𝑟 for every 1 ≤ 𝑟 ≤ 𝑑 D in dimension 𝑑 𝑇 1 (𝑛) = 𝑛 𝑇 𝑟 (𝑛) = 1 + 𝑇 𝑟 (𝑛 -1) + 𝑇 𝑟 -1 ( { 𝑟 𝑛 : 𝑛 -1; 1 -𝑟 𝑛 : 0}) (2 ≤ 𝑟 ≤ 𝑑) 𝑐 𝑟 = 𝑟 for every 1 ≤ 𝑟 ≤ 𝑑
Table 1. Our Benchmark SPRRs. We use the notation {𝑝 1 : expr 1 , . . . , 𝑝 𝑘 : expr 𝑘 } to denote that with probability 𝑝 𝑖 the value is expr 𝑖 . MU (𝑛) denotes max{𝑖, 𝑛 -1 -𝑖} where 𝑖 is uniformly distributed.

(Appendix F).

EXPERIMENTAL RESULTS

In this section, we provide experimental results over SPRRs obtained from various classical randomized algorithms. These results demonstrate that our approach is able to obtain much tighter bounds than [START_REF] Karp | Probabilistic Recurrence Relations[END_REF]] in practice. Moreover, it can also handle a wider family of SPRRs.

Benchmarks. We evaluated the two algorithmic approaches of [START_REF] Welzl | Smallest enclosing disks (balls and ellipsoids)[END_REF]]. We also considered a manually-crafted algorithm, S S , that sorts an array using a divide-and-conquer method in which in each iteration the median element is obtained by S and used as a pivot. Due to space constraints, a detailed description of these algorithms is provided in Appendix B. Table 1 provides a summary of their SPRRs. The initial procedure 𝑟 * is always the last procedure in the SPRR and the initial input size is denoted by 𝑛 * . Implementation and environment. We implemented our algorithms in C++ used the miracl [START_REF] Miracl Ltd | [END_REF] These examples show that our approach handles a wider set of SPRRs than previous methods. • Polynomially decreasing concentration bounds: Table 4 shows our experimental results using the algorithm of Section 5.4 to obtain polynomially decreasing bounds. In all benchmarks, we let 𝑢 𝑟,𝑖 = 𝑣 𝑟,𝑖 = 𝑖 in our templates for every 1 ≤ 𝑟 ≤ 𝑚 and 1 ≤ 𝑖 ≤ 𝑘. To the best of our knowledge, we are providing the rst algorithm for synthesizing polynomial bounds. As such, there are no previous methods to compare against. Moreover, previous methods for exponential bounds do not yield any results on these benchmarks. We believe this is probably because the best possible concentration bounds are polynomially decreasing.

( 1 2 ) 𝛼 -2 3 • 𝑛 2 • 𝑛 exp(1.540 • (3.0 -𝛼)) 2.16 3.5 • 𝑛 2.5 • 𝑛 exp(1.774 • (3.5 -𝛼)) 2.08 L2D Pr[ C𝜏 ≥ 𝛼 • 𝑛 * • ln𝑛 * ] 2.5 •𝑛 •ln 𝑛 1.5 • 𝑛 •ln𝑛 𝜆 𝑛•ln 𝑛 exp(1.008 • (2.5 -𝛼)) 3.09 ( 1 2 ) 𝛼 -2 3 •𝑛 •ln𝑛 2 •𝑛 •ln𝑛 exp(1.338 • (3.0 -𝛼)) 3.10 3.5 •𝑛 •ln𝑛 2.5 •𝑛 •ln𝑛 exp(1.189 • (3.5 -𝛼)) 3.10 S Pr[ C𝜏 ≥ 𝛼 • 𝑛 * ] 5.0 • 𝑛 3 • 𝑛 𝜆 𝑛 exp(1.007 • (5.0 -𝛼)) 1.92 ( 3
D (𝑑 = 3) 10 • 𝑛 * Pr[ C 𝜏 ≥ 𝛼 • 10 • 𝑛 * ] 𝑘 = 3 1.487 • 𝛼 -3 0.15 𝑘 = 4 2.125 • 𝛼 -4 0.37 𝑘 = 5 3.318 • 𝛼 -5 1.75 D (𝑑 = 4) 41 • 𝑛 * Pr[ C 𝜏 ≥ 𝛼 • 41 • 𝑛 * ] 𝑘 = 2 1.137 • 𝛼 -2 0.21 𝑘 = 3 1.439 • 𝛼 -3 0.28 𝑘 = 4 2.000 • 𝛼 -4 3.40

E ciency.

As shown in Tables 2 and4, our algorithms are extremely e cient in practice and can synthesize much tighter bounds in a matter of a few seconds. In our experiments, the maximum runtime was only 3.40 seconds, while the average runtime was 1.38 seconds.

RELATED WORKS

In this section, we compare our approach with the most relevant results in the literature.

Expected runtime analysis. Expected runtime analysis of probabilistic recurrences has been widely considered in e.g. [START_REF] Bazzi | The Solution of Linear Probabilistic Recurrence Relations[END_REF]Chatterjee et al. 2017a;[START_REF] Flajolet | Automatic Average-Case Analysis of Algorithm[END_REF][START_REF] Zimmermann | The automatic complexity analysis of divide-and-conquer algorithms[END_REF]]. However, methods for expected runtime analysis cannot be adapted to concentration bound analysis, since it is a much more intricate problem.

Karp's cookbook. For the concentration bound analysis of probabilistic recurrences, the most well-known and related result is the classical method of [START_REF] Karp | Probabilistic Recurrence Relations[END_REF]], where it is proven that concentration bounds of a wide class of probabilistic recurrences can be directly computed by a mathematical formula. Compared with [START_REF] Karp | Probabilistic Recurrence Relations[END_REF]], our algorithmic approaches can derive much tigher concentration bounds over probabilistic recurrences for classical randomized algorithms including S , S , R S and D C . Furthermore, our algorithmic approaches can handle systems of probabilistic recurrences with randomized preprocessing time, such as S LP and S E D , which are beyond the scope of Karp's method and other existing approaches such as [START_REF] Chaudhuri | Probabilistic Recurrence Relations Revisited[END_REF][START_REF] Mcdiarmid | Large Deviations for Quicksort[END_REF][START_REF] Tassarotti | Probabilistic Recurrence Relations for Work and Span of Parallel Algorithms[END_REF][START_REF] Tassarotti | Veri ed Tail Bounds for Randomized Programs[END_REF].

Previous results on concentration bounds. The work [START_REF] Chaudhuri | Probabilistic Recurrence Relations Revisited[END_REF] weakens several conditions required by [START_REF] Karp | Probabilistic Recurrence Relations[END_REF]], but obtains bounds that are theoretically worse than [START_REF] Karp | Probabilistic Recurrence Relations[END_REF]]. Moreover, [START_REF] Chaudhuri | Probabilistic Recurrence Relations Revisited[END_REF] cannot handle probabilistic recurrences with randomized preprocessing time and systems of probabilistic recurrences. Thus, the comparison between our approach and [START_REF] Karp | Probabilistic Recurrence Relations[END_REF]] carries over to [START_REF] Chaudhuri | Probabilistic Recurrence Relations Revisited[END_REF]. The work [START_REF] Mcdiarmid | Large Deviations for Quicksort[END_REF] performs an involved ad-hoc manual analysis to derive the asymptotically optimal concentration bound for S . In contrast, our approach is automated for a wide class of probabilistic recurrences and derives a concentration bound for the S that is slightly worse than [START_REF] Mcdiarmid | Large Deviations for Quicksort[END_REF]. The work [START_REF] Tassarotti | Veri ed Tail Bounds for Randomized Programs[END_REF] mechanizes the proof of the concentration bound for S in Coq, thus it is incomparable to our approach. Finally, the recent results in [START_REF] Tassarotti | Probabilistic Recurrence Relations for Work and Span of Parallel Algorithms[END_REF]] introduce a new approach for solving work and span runtime of parallel randomized algorithms. The method of [START_REF] Tassarotti | Probabilistic Recurrence Relations for Work and Span of Parallel Algorithms[END_REF]] is manual and can neither handle randomized preprocessing time nor systems of probabilistic recurrences. Moreover, our automated approach can derive much tighter concentration bounds for S against this method: our automatic approach generate the concentration bound Pr[𝑇 (𝑛) ≥ 11 • 𝑛 * ln 𝑛 * + 12 • 𝑛 * ] ≤ exp(-5.615 ln 𝑛 * -13.476), while this method outputs Pr[𝑇 (𝑛) ≥ 11 • 𝑛 * ln 𝑛 * + 12 • 𝑛 * ] ≤ exp(-0.24 • ln 𝑛 * -1.47).

Non-probabilistic Recurrence

Relations. There many results on the formal analysis of classical non-probabilistic recurrence relations in the literature. For example, [START_REF] Wegbreit | Mechanical Program Analysis[END_REF]] considered solving recurrence relations through di erence equations and generating functions. [START_REF] Zimmermann | The automatic complexity analysis of divide-and-conquer algorithms[END_REF]] also considered transforming them into di erence equations. Grobauer [START_REF] Grobauer | Cost Recurrences for DML Programs[END_REF]] uses generating recurrence relations from DML for the worst-case analysis. The COSTA project [START_REF] Albert | Termination and Cost Analysis with COSTA and its User Interfaces[END_REF][START_REF] Albert | Automatic Inference of Upper Bounds for Recurrence Relations in Cost Analysis[END_REF][START_REF] Albert | Cost Analysis of Java Bytecode[END_REF] transforms Java bytecode into recurrence relations and solves them through ranking functions. The PURRS tool [START_REF] Bagnara | PURRS: Towards Computer Algebra Support for Fully Automatic Worst-Case Complexity Analysis[END_REF] addresses nite linear recurrences (with bounded summation), and some restricted linear in nite recurrence relations (with unbounded summation). Recently, new approaches are proposed for solving non-probabilistic recurrence relations through compositional approaches [START_REF] Farzan | Compositional Recurrence Analysis[END_REF][START_REF] Kincaid | Compositional recurrence analysis revisited[END_REF], abstract interpretation [START_REF] Kincaid | Non-linear reasoning for invariant synthesis[END_REF] and templates [START_REF] Breck | Templates and recurrences: better together[END_REF]. Since these approaches do not consider probability, they cannot be applied to the formal analysis of probabilistic recurrence relations.

Probabilistic programs. There are also many relevant results in formal analysis of probabilistic programs, such as martingale-based approaches of [Chatterjee and Fu 2017;[START_REF] Chatterjee | Termination Analysis of Probabilistic Programs Through Positivstellensatz's[END_REF]Chatterjee et al. , 2018aChatterjee et al. ,b, 2017b;;[START_REF] María | Probabilistic Termination: Soundness, Completeness, and Compositionality[END_REF][START_REF] Fu | Termination of Nondeterministic Probabilistic Programs[END_REF][START_REF] Huang | New Approaches for Almost-Sure Termination of Probabilistic Programs[END_REF][START_REF] Huang | Modular veri cation for almost-sure termination of probabilistic programs[END_REF][START_REF] Kura | Tail Probabilities for Randomized Program Runtimes via Martingales for Higher Moments[END_REF][START_REF] Mciver | A new proof rule for almost-sure termination[END_REF]Wang et al. 2020b[START_REF] Wang | Cost analysis of nondeterministic probabilistic programs[END_REF]. These approaches rst synthesize linear/polynomial (ranking) supermartingales through constraint solving, and then use concentration inequalities to derive concentration bounds of probabilistic programs. Compared with these results, our approaches synthesize Taylor supermartingale functions in exponential and polynomial form by Taylor expansion and tight approximation of expectation of moments with logarithm, thus they are completely di erent, both in terms of the problem they solve and the employed techniques.

CONCLUSION

In this work, we considered the problem of formal analysis of concentration bounds of probabilistic recurrence relations. We rst proposed a general theoretical approach through Markov's inequality and Optional Stopping Theorem. The core notion in our general approach is that of Taylor supermaringale function sequences, which establish supermartingale conditions through Taylor series. Then, we introduced two algorithmic approaches that synthesize exponential and polynomial supermartingales for deriving the concentration bounds. Our experimental results demonstrate that our algorithmic approaches can obtain much tighter concentration bounds than existing approaches and can also handle a wider variety of recurrence relations.

Since the two methods are rather dissimilar, an interesting future direction would be to explore the integration of our approaches and the classical cookbook method by Karp in order to obtain even tighter bounds. Another interesting direction is to automate the process of obtaining SPRRs from probabilistic programs and randomized algorithms, similar to how the COSTA project [START_REF] Albert | Termination and Cost Analysis with COSTA and its User Interfaces[END_REF][START_REF] Albert | Automatic Inference of Upper Bounds for Recurrence Relations in Cost Analysis[END_REF][START_REF] Albert | Cost Analysis of Java Bytecode[END_REF] translates Java bytecode into non-probabilistic recurrences. While this work provides automated approaches for concentration bound analysis of SPRRs, the translation from a program to an SPRR is still performed manually. The SPRRs for classical algorithms are well-known and widely studied, but there is still no automated method to obtain SPRRs for arbitrary programs.

A AN ILLUSTRATIVE EXAMPLE OF A MARTINGALE

Example A.1 (Pólya's Urn [START_REF] Hosam | Pólya urn models[END_REF]). As a more interesting example, consider an urn that initially contains 𝑅 0 red and 𝐵 0 blue marbles (𝑅 0 + 𝐵 0 > 0). At each step, we take one marble from the urn, chosen uniformly at random, look at its color and then add two marbles of that color to the urn. Let 𝐵 𝑛 , 𝑅 𝑛 and 𝑀 𝑛 respectively be the number of red, blue and all marbles after 𝑛 steps. Also, let 𝛽 𝑛 = 𝐵 𝑛 𝑀 𝑛 and 𝜌 𝑛 = 𝑅 𝑛 𝑀 𝑛 be the proportion of marbles that are blue (resp. red) after 𝑛 steps. Let F 𝑛 model the observations until the 𝑛-th step. The process described above leads to the following equations:

𝑀 𝑛+1 = 1 + 𝑀 𝑛 , E [𝐵 𝑛+1 | F 𝑛 ] = E [𝐵 𝑛+1 | 𝐵 1 , . . . , 𝐵 𝑛 ] = 𝐵 𝑛 𝑀 𝑛 • (𝐵 𝑛 + 1) + 𝑅 𝑛 𝑀 𝑛 • 𝐵 𝑛 , E [𝑅 𝑛+1 | F 𝑛 ] = E [𝑅 𝑛+1 | 𝐵 1 , . . . , 𝐵 𝑛 ] = 𝑅 𝑛 𝑀 𝑛 • (𝑅 𝑛 + 1) + 𝐵 𝑛 𝑀 𝑛 • 𝑅 𝑛 .
Note that we did not need to care about observing 𝑅 𝑖 's, 𝑀 𝑖 's, 𝛽 𝑖 's or 𝜌 𝑖 's, because they can be uniquely computed in terms of 𝐵 𝑖 's. More generally, an observer can observe only 𝐵 𝑖 's, or only 𝑅 𝑖 's, or only 𝛽 𝑖 's or 𝜌 𝑖 's and can then compute the rest using this information. Based on the equations above, we have:

E [𝛽 𝑛+1 | F 𝑛 ] = 𝐵 𝑛 𝑀 𝑛 • 𝐵 𝑛 + 1 𝑀 𝑛 + 1 + 𝑀 𝑛 -𝐵 𝑛 𝑀 𝑛 • 𝐵 𝑛 𝑀 𝑛 + 1 = 𝐵 𝑛 𝑀 𝑛 = 𝛽 𝑛 , E [𝜌 𝑛+1 | F 𝑛 ] = 𝑅 𝑛 𝑀 𝑛 • 𝑅 𝑛 + 1 𝑀 𝑛 + 1 + 𝑀 𝑛 -𝑅 𝑛 𝑀 𝑛 • 𝑅 𝑛 𝑀 𝑛 + 1 = 𝑅 𝑛 𝑀 𝑛 = 𝜌 𝑛 .
Hence, both {𝛽 𝑛 } 𝑛 ∈N 0 and {𝜌 𝑛 } 𝑛 ∈N 0 are martingales. Informally, this means that the expected proportion of blue marbles in the next step is exactly equal to their observed proportion in the current step. This might be counter-intuitive. For example, consider a state where 99% of the marbles are blue. Then, it is more likely that we will add a blue marble in the next state. However, this is mitigated by the fact that adding a blue marble changes the proportions much less dramatically than adding a red marble.

B EXAMPLES OF (SYSTEMS OF) PROBABILISTIC RECURRENCE RELATIONS

Example B.1 (Q S ).
Consider the problem of nding the 𝑑-th smallest element in an unordered array of 𝑛 distinct elements. A classical randomized algorithm for solving this problem is S [Hoare 1961b]. It begins by choosing a pivot element 𝑢 of the array uniformly at random. It then compares all the other 𝑛 -1 elements of the array with 𝑢 and divides them into two parts: (i) those that are smaller than 𝑢 and (ii) those that are larger. Suppose that there are 𝑑 elements in part (i). If 𝑑 < 𝑑 -1, then the algorithm recursively searches for the (𝑑 -𝑑 -1)-th smallest element of part (ii). If 𝑑 = 𝑑 -1, the algorithm terminates by returning 𝑢 as the desired answer. Finally, if 𝑑 > 𝑑 -1, the algorithm recursively nds the 𝑑-th smallest element in part (i). Note that the classical median selection algorithm is a special case of S . While there is an involved deterministic algorithm for solving the same problem, more involved linear-time non-randomized algorithms exist for the same problem, S provides a simple randomized variant matching the best-known 𝑂 (𝑛) deterministic runtime. Since there is only one recursive procedure in S , we model the algorithm as the following SPRR with one recurrence equation (where we have

𝑚 = 1, 𝑇 1 = 𝑇 , Cost 1 (𝑛) = 𝑛 -1, 𝑆 1 = 1 and 𝑂 1 = 0): 𝑇 (𝑛) = 𝑛 + 𝑇 (ℎ(𝑛)) .
Here, 𝑇 (𝑛) represents the number of comparisons performed by S over an input of size 𝑛, and ℎ(𝑛) is the random variable that captures the size of the remaining array that has to be searched recursively. It can be derived that ℎ(𝑛) = max{𝑖 -1, 𝑛 -𝑖} where 𝑖 is sampled uniformly from {1, . . . , 𝑛}.

Example B.2 (Q S

). Consider the classical problem of sorting an array of 𝑛 distinct elements. A well-known randomized algorithm for solving this problem is S . Similar to S , S begins by choosing a pivot element 𝑢 of the array uniformly at random. It then compares all the other 𝑛 -1 elements of the array with 𝑢 and divides them into two sub-arrays: (i) those that are smaller than 𝑢 and (ii) those that are larger. The algorithm then runs recursively on both sub-arrays. Compared with its deterministic counterpart, M S , S provides a simple randomized variant whose expected runtime matches the best-known deterministic runtime 𝑂 (𝑛 log 𝑛).

Just as before, there is only one recursive procedure in S . Thus, we model the algorithm as the following SPRR with only one recurrence equation (where we have

𝑚 = 1, Cost 1 (𝑛) = 𝑛 -1, 𝑆 1 = 2, 𝑂 1 = 0): 𝑇 (𝑛) = 𝑛 + 𝑇 (ℎ 1 (𝑛)) + 𝑇 (ℎ 2 (𝑛)) .
Here, ℎ 1 (𝑛) and ℎ 2 (𝑛) are random variables that capture the sizes of the two sub-arrays. It is easy to see that ℎ 1 (𝑛) = 𝑖 and ℎ 2 (𝑛) = 𝑛 -1 -𝑖 where 𝑖 is uniformly sampled from {0, . . . , 𝑛 -1}.

Example B.3 (S LP).

Consider the classic linear programming problem in 𝑑 dimensions, with 𝑑 variables x = (x 1 , . . . , x 𝑑 ) and 𝑛 constraints a 𝑇 1 • x ≤ 𝑏 1 , . . . , a 𝑇 𝑛 • x ≤ 𝑏 𝑛 , and the objective function to be minimized is c 𝑇 • x. S LP is a famous randomized incremental algorithm for linear programming. The algorithm is as follows, it picks a constraint a 𝑇 𝑖 • x ≤ 𝑏 𝑖 from all 𝑛 constraints uniformly at random and recursively solves the linear programming instance without the chosen constraint. Given the optimum 𝑥 * produced by recursion, the algorithm analyzes two cases: Case 1. If 𝑥 * still satis es the ignored constraint a 𝑖 • x ≤ 𝑏 𝑖 , then 𝑥 * is still optimal. The checking required 𝑑 steps of computation. Case 2. Otherwise, if 𝑥 * violates the ignored constraint, then the new optimum must lie in the hyperplane a 𝑖 • x = 𝑏 𝑖 , in this case, we could project all other constraints onto this hyperplane and recursively solves the subinstance. Note that since the hyperplane a 𝑖 •x = 𝑏 𝑖 has dimension 𝑑 -1, the subinstance is a linear programming in 𝑑 -1 dimensions. The recursion stops when the number of constraints 𝑛 equals to the dimension 𝑑, where in this case the optimum could be directly computed. Note that the probability of Case 2 above is no more than 𝑑 𝑛 , since there is no more than 𝑑 tight constraints in a 𝑑-dimensional LP instance. To model the algorithm above by SPRRs, we could set up a system with 𝑑 relations 𝑒𝑞 1 , 𝑒𝑞 2 , . . . , 𝑒𝑞 𝑑 , where 𝑒𝑞 𝑟 corresponds to the 𝑟 -dimensional recursion subinstaces in the algorithm. We have that 𝑐 𝑟 = 𝑟 and for every 𝑛 ≥ 𝑟 + 1:

𝑇 𝑟 (𝑛) = Cost 𝑟 (𝑛) + 𝑇 𝑟 (ℎ 𝑟,1 (𝑛)) + 𝑇 𝑟 -1 (𝑔 𝑟,1 (𝑛)) (9)
Where Cost 𝑟 (𝑛) observes as a discrete distribution: with probability 1 -𝑟 𝑛 , its value is 𝑑 and with probability 𝑟 𝑛 , its value is 𝑟 • 𝑛, which corresponds the cost for checking and projection in the two cases above, ℎ 𝑟,1 (𝑛) deterministically equals to 𝑛 -1, and the 𝑔 𝑟,1 (𝑛) also observes as discrete distribution: with probability 1 -𝑟 𝑛 its value is 0 and with probability 𝑟 𝑛 its value is 𝑛 -1, since in Case 1 discussed above, there is no need to recursively calculate a 𝑑 -1-dimensional subinstance while in the Case 2 we need to solve such a subinstance with 𝑛 -1 constraints.

Example B.4 (D C

). Consider the D C algorithm [Motwani and Raghavan 1995, Chapter 9] to compute the diameter of an input nite set 𝑆 of threedimensional points. Depending on the metric, i.e. Eucledian or 𝐿 1 , we obtain two di erent recurrence relations. For Euclidean, we have the following SPRR:

𝑇 (𝑛) = 𝑛 • ln 𝑛 + 𝑇 (ℎ(𝑛)) (10)
with 𝑐 = 0, and where ℎ(𝑛) observes as U (𝑛).For the 𝐿 1 metric, the SPRR is as follows:

𝑇 (𝑛) = 𝑛 + 𝑇 (ℎ(𝑛)) (11)
with 𝑐 = 0, and where ℎ(𝑛) observes as U (𝑛). Note that here we use deterministic versions for the subroutine H I as required in D C , see [START_REF] Franco | Finding the Intersection of n Half-Spaces in Time O(n log n)[END_REF] for the Euclidean case and [Motwani and Raghavan 1995, Problem 9.6, Page 276] for the 𝐿 1 case.

Example B.5 (R S

). Consider Sherwood's R S algorithm (cf. [Mc-Connell 2001, Chapter 9]). The algorithm checks whether an integer value 𝑑 is present within the index range [𝑖, 𝑗] (0 ≤ 𝑖 ≤ 𝑗) in an integer array 𝑎𝑟 which is sorted in increasing order and is without duplicate entries. The algorithm outputs either the index of the item or -1 if 𝑑 is not present in the index range [𝑖, 𝑗]. The SPRR for this example contains only one recurrence equation:

𝑇 (𝑛) = 1 + 𝑇 (ℎ(𝑛)) (12)
with 𝑐 = 0 and ℎ(𝑛) = max{𝑖, 𝑛 -1 -𝑖}, where 𝑖 is uniformly distributed over {0, 1, . . . , 𝑛 -1}.

Example B.6 (S E D ). Consider the S E D problem from computational geometry: given 𝑛 points in a two-dimensional plane, the goal is to nd a smallest circle that covers all the points. Note that any such circle can be identi ed by three input points that are on the circle. We call these the canonical points . To solve this problem, Welzl [START_REF] Welzl | Smallest enclosing disks (balls and ellipsoids)[END_REF]] has proposed a randomized algorithm that runs in linear expected time and searches for suitable canonical points. The algorithm S E D is split into three procedures. The main procedure, proc 3 , chooses an input point 𝑥 uniformly at random. It then guesses that the chosen point is not canonical and tries to form a smallest enclosing disk by nding three canonical points amongst the other 𝑛 -1 points. If this fails, then the algorithm knows that 𝑥 is canonical and calls proc 2 , which is a procedure that gets one canonical point, in this case 𝑥, as input and tries to nd two other canonical points. This leads to the following recurrence:

𝑇 3 (𝑛) = 1 + 𝑇 3 (𝑛 -1) + 𝑇 2 (𝑔 31 (𝑛)) (𝑒𝑞 3 )
where 𝑇 3 (𝑛) is the running time of the main procedure, proc 3 , and 𝑇 2 (𝑛) is the running time of proc 2 . Moreover, 𝑔 31 (𝑛) is the size of the random instance passed to proc 2 and is equal to zero with probability 𝑛-3 𝑛 , i.e. when the algorithm could successfully nd three canonical points other than 𝑥, and equal to 𝑛 -1 with probability 3 𝑛 , i.e. when 𝑥 must be canonical. The procedure proc 2 is in turn very similar to proc 3 . It chooses a point uniformly at random and throws it out, trying to nd the canonical points among the other 𝑛 -1. If this fails, it calls proc 1 , which is a procedure that, given a xed pair of canonical points, nds the nal canonical point (or reports that it is impossible to do so).

𝑇 2 (𝑛) = 1 + 𝑇 2 (𝑛 -1) + 𝑇 1 (𝑔 21 (𝑛)) (𝑒𝑞 2 )
As before, 𝑔 21 (𝑛) is the random size of the instance passed to proc 1 . It equals 𝑛 -1 with probability 2 𝑛 and zero with probability 𝑛-2 𝑛 . Finally, proc 1 is a simple linear search that can be modeled by this recurrence:

𝑇 1 (𝑛) = 1 + 𝑇 1 (𝑛 -1) (𝑒𝑞 1 )
The canonical points are not necessarily unique.

Putting everything together, we obtain the following SPRR with three procedures to describe the runntime behavior of S E D :

𝑇 1 (𝑛) = 1 + 𝑇 1 (𝑛 -1) (𝑒𝑞 1 ) 𝑇 2 (𝑛) = 1 + 𝑇 2 (𝑛 -1) + 𝑇 1 (𝑔 21 (𝑛)) (𝑒𝑞 2 ) 𝑇 3 (𝑛) = 1 + 𝑇 3 (𝑛 -1) + 𝑇 2 (𝑔 31 (𝑛)) (𝑒𝑞 3 ) Example B.7 (C C R ).
We consider two network scenarios in which 𝑛 clients are trying to get access to a network channel. This problem is also called Resource-Contention Resolution [Kleinberg and Tardos 2006, Chapter 13]. In this problem, if more than one client tries to access the channel, then no client can access it, and if exactly one client requests access to the channel, then the request is granted. In the concurrent setting, the clients share one variable, which is the number of clients which have not yet been granted access. Also in this scenario, once a client gets an access the client does not request for access again. For this problem, we obtain an over-approximating recurrence relation

𝑇 (𝑛) = 1 + 𝑇 (ℎ(𝑛)) (13)
where ℎ(𝑛) observes as a discrete distribution with probability 1 𝑒 to be 𝑛 -1, and probability 1 -1 𝑒 to be 𝑛.

Example B.8 (S Q S

). Consider a sorting algorithm which selects the median through the S algorithm. The SPRR for this example is:

𝑇 1 (𝑛) = 𝑛 + 𝑇 1 (ℎ 1 (𝑛)) (14)
𝑇 2 (𝑛) = 𝑇 1 (𝑛) + 𝑇 2 ( (𝑛 -1)/2 ) + 𝑇 2 ( (𝑛 -1)/2 ) (15)
where 𝑇 1 is the SPRR for S (See. Example B.1).

Example B.9 (R ). Consider the execution of a Robot dead reckoning approach, as used in probabilistic robotics. Dead reckoning is an approach for position estimation starting from a known xed position at time 𝑡 = 0. In each time unit, there is 0.5 probability that the estimation error is -0.05, and 0.5 probability that the estimation error is +0.05. We use 𝑇 (𝑛) to model the total accumulated error in 𝑛 time units. Thus, the recurrence relation is:

𝑇 (𝑛) = Cost(𝑛) + 𝑇 (𝑛 -1)
with 𝑐 = 0. Cost(𝑛) observes as a discrete distribution with 0.5 probability of being -0.05 and 0.5 probability of being +0.05.

C FORMAL CONSTRUCTION OF THE PROBABILITY SPACE IN AN SPRR'S

SEMANTICS Given an initial stack 𝛽 * ∥ , the probability space (Ω Δ , F Δ , Pr Δ,𝛽 * ) for an SPRR Δ follows standard de nitions (cf. [Baier and Katoen 2008, Chapter 10]) so that (i) the sample space Ω Δ is de ned as the set of all in nite sequences (conventionally called runs) of stacks, (ii) F Δ is the smallest sigma-algebra generated from the set of all cylinder sets for which a cylinder set is the set of all runs (of stacks) sharing some common nite pre x, and (iii) Pr Δ,𝛽 * is the unique probability measure such that for every cylinder set 𝐴 speci ed by the common nite pre x 𝛽 0 . . . 𝛽 𝑘 , we have Pr . 2018, P 3.1]). Let 𝑋 be any non-negative random variable from a probability space (Ω, F , Pr) and G be a sub-sigma-algebra of F (i.e., G ⊆ F ). Then there exists a random variable E [𝑋 | G] from (Ω, F , Pr) (called a conditional expectation of 𝑋 w.r.t G) that ful lls the conditions (E1) and (E3) above. Moreover, the conditional expectation is almost-surely unique, i.e., for any random varaibles 𝑌 , 𝑍 from (Ω, F , Pr) that both ful ll (E1) and (E3) we have that Pr(𝑌 = 𝑍 ) = 1.

To prove Theorem 3.4, we need the extension of conditional expectation to non-negative random variables that are not necessarily integrable, as is illustrated in Theorem D.1. One issue that arises in extending the traditional notion of conditional expectation is whether the basic properties for conditional expectation still hold for non-negative random variables that are not necessarily integrable. Below we show that several basic properties we utilize in the proof of Theorem 3.4 is still valid upon the extension.

Below we x a probability space (Ω, F , Pr) and a sub-sigma-algebra G ⊆ F . The following properties hold for non-negative random variables 𝑋, 𝑌 , 𝑍 and 𝑋 0 , . . . , 𝑋 𝑛 , . . . from (Ω, F , Pr).

(

E4) E[E [𝑋 | G]] = E[𝑋 ]
. This follows directly from (E3) where we take 𝐴 to be the whole sample space Ω. ). This follows exactly the same lines under 'Proof of (e)' on [Williams 1991, Page 89]

(E5) If 𝑋 is G-measurable, then E [𝑋 | G] = 𝑋 a.
. (E9) ('Taking out what is known') If 𝑍 is G-measurable, then E [𝑍 • 𝑋 | G] = 𝑍 • E [𝑋 | G] a.s.
The proof is exactly the same as under 'Proof of (j)' on [Williams 1991, Page 90] for the non-negative case. Below we present the proof for our main theorem.

Theorem 3.4. Let Δ be a given SPRR. Consider a T-function 𝑓 ∈ T and an in nite sequence Ψ of functions Ψ = Ψ 0 , Ψ 1 , Ψ 2 , . . . such that each Ψ 𝑘 maps a stack 𝛽 of Δ to its corresponding non-negative real value Ψ 𝑘 (𝛽). Suppose that the sequence Ψ satis es the following conditions: (A1) for every 𝑘 ≥ 0, we have Ψ 𝑘 (𝜀) ≥ 𝑓 (𝑘) (0); (A2) for every stack 𝛽 and 𝑘 ≥ 0, we have

Ψ 𝑘 (𝛽) ≥ E ∞ 𝑗=0 Ψ 𝑘+𝑗 (next𝛽) 𝑗! • Cost 𝑗 (𝛽) .
Then for every stack 𝛽 and 𝑘 ≥ 0, it holds that Ψ 𝑘 (𝛽) ≥ E 𝛽 𝑓 (𝑘) (C) . F T 3.4. Fix an initial stack 𝛽 * and consider the random trace 𝛽 0 , 𝛽 1 , 𝛽 2 , . . . for the SPRR Δ with the initial stack 𝛽 0 = 𝛽 * . For each 𝑘 ≥ 0, we de ne G 𝑘 as the smallest sigma-algebra that makes the random variables in 𝛽 0 , . . . , 𝛽 𝑘 and C 0 , . . . , C 𝑘-1 measurable. It is straightforward to observe that {G 𝑘 } 𝑘 ≥0 is a ltration.

For each 𝑘, 𝑖 ≥ 0, we de ne the random variable

𝑌 𝑘,𝑖 := ∞ 𝑗=0 Ψ 𝑖+𝑗 ( 𝛽 𝑘 ) 𝑗! • C 𝑗 𝑘
where we assume 0 0 := 1. We prove that for every 𝑘, 𝑖 ≥ 0, the following properties hold:

(a) for every 𝑘, 𝑖 ≥ 0, 𝑌 𝑘,𝑖 is non-negative; (b) 𝑌 𝜏,𝑖 ≥ 𝑓 (𝑖) (C 𝜏 );

(c) For every 𝑖 ≥ 0, we have 𝑌 𝑘,𝑖 forms a supermartingale with respect to the ltration {G 𝑘 } 𝑘 ≥0 . First, we have that the property (a) follows directly from the de nition of 𝑌 𝑘,𝑖 . Second, as 𝛽 𝜏 = 𝜀, the property (b) can be deduced from (A1):

𝑌 𝜏,𝑖 = ∞ 𝑗=0 Ψ 𝑖+𝑗 ( 𝛽 𝜏 ) 𝑗! • C 𝑗 𝜏 ≥ ∞ 𝑗=0 𝑓 (𝑖+𝑗 ) (0) 𝑗! • C 𝑗 𝜏 = 𝑓 (𝑖) (C 𝜏 ).
Below we prove the property (c). By Theorem D.1, we can take a conditional expectation E[𝑌 𝑘+1,𝑖 | G 𝑘 ] of 𝑌 𝑘+1,𝑖 w.r.t G 𝑘 and derive that the following equalities hold (a.s.):

E[𝑌 𝑘+1,𝑖 | G 𝑘 ] = E ∞ 𝑗=0 Ψ 𝑖+𝑗 ( 𝛽 𝑘+1 ) 𝑗! • C 𝑗 𝑘+1 G 𝑘 = E ∞ 𝑗=0 Ψ 𝑖+𝑗 ( 𝛽 𝑘+1 ) 𝑗! • (C 𝑘 + Cost( 𝛽 𝑘 )) 𝑗 G 𝑘 = E ∞ 𝑗=0 Ψ 𝑖+𝑗 ( 𝛽 𝑘+1 ) 𝑗! • 𝑗 𝑙=0 𝑗 𝑙 • C 𝑙 𝑘 • Cost( 𝛽 𝑘 ) 𝑗-𝑙 G 𝑘 = E ∞ 𝑙=0 C 𝑙 𝑘 𝑙! • ∞ 𝑗=𝑙 Ψ 𝑖+𝑗 ( 𝛽 𝑘+1 ) ( 𝑗-𝑙)! • Cost( 𝛽 𝑘 ) 𝑗-𝑙 G 𝑘 = § By applying (E6)-(Linearity) and (E8)-(cMON) § ∞ 𝑙=0 E C 𝑙 𝑘 𝑙! • ∞ 𝑗=𝑙 Ψ 𝑖+𝑗 ( 𝛽 𝑘+1 ) ( 𝑗-𝑙)! • Cost( 𝛽 𝑘 ) 𝑗-𝑙 G 𝑘 = § By applying (E9)-('Taking out what is known') § ∞ 𝑙=0 C 𝑙 𝑘 𝑙! • E ∞ 𝑗=𝑙 Ψ 𝑖+𝑗 ( 𝛽 𝑘+1 ) ( 𝑗-𝑙)! • Cost( 𝛽 𝑘 ) 𝑗-𝑙 G 𝑘 .
Furthermore, we have that the following holds (a.s.): Consider any event 𝐴 ∈ G 𝑘 such that 𝐴 is the disjoint union of 𝑑 cylinder sets de ned by 𝑑 distinct nite pre xes 𝛼 1 , . . . , 𝛼 𝑑 for which each 𝛼 𝑗 (1 ≤ 𝑗 ≤ 𝑑) is a nite sequence of stacks with length 𝑘. We use P(𝛽, (𝛽 , 𝑐)) to denote the probability that (the next stack, the cost triggered) is equal to (𝛽 , 𝑐) given the current stack 𝛽. By the de nition of integration, we have that

E ∞ 𝑗=𝑙 Ψ 𝑖+𝑗 ( 𝛽 𝑘+1 ) ( 𝑗-𝑙)! • Cost( 𝛽 𝑘 ) 𝑗-𝑙 G 𝑘 = § By applying (E6)-(Linearity) and (E8)-(cMON) § ∞ 𝑗=𝑙 E Ψ 𝑖+𝑗 ( 𝛽 𝑘+1 ) ( 𝑗-𝑙)! • Cost( 𝛽 𝑘 ) 𝑗-𝑙 G 𝑘 = §
∫ 𝐴 𝛽 ,𝑐 1 ( 𝛽 𝑘 , 𝛽 𝑘+1 ,Cost( 𝛽 𝑘 ))=(𝛽,𝛽 ,𝑐) • Ψ 𝑖+𝑗 (𝛽 ) ( 𝑗 -𝑙)! • 𝑐 𝑗-𝑙 dPr 𝛽 * = 𝑑 𝑗=1 1 𝛼 𝑗 ↓=𝛽 • Pr 𝛽 * Cyl(𝛼 𝑗 ) •       𝛽 ,𝑐 P(𝛽, (𝛽 , 𝑐)) • Ψ 𝑖+𝑗 (𝛽 ) ( 𝑗 -𝑙)! • 𝑐 𝑗-𝑙       = ∫ 𝐴 1 𝛽 𝑘 =𝛽 •       𝛽 ,𝑐 P(𝛽, (𝛽 , 𝑐)) • Ψ 𝑖+𝑗 (𝛽 ) ( 𝑗 -𝑙)! • 𝑐 𝑗-𝑙       dPr 𝛽 *
where Cyl(𝛼 𝑗 ) is the cylinder set with the common pre x 𝛼 𝑗 and 𝛼 𝑗 ↓ denotes the last stack of 𝛼 𝑗 .

Thus by the de nition of conditional expectation, we have that

E[ 𝛽 ,𝑐 1 ( 𝛽 𝑘 , 𝛽 𝑘+1 ,Cost( 𝛽 𝑘 ))=(𝛽,𝛽 ,𝑐) • Ψ 𝑖+𝑗 (𝛽 ) ( 𝑗 -𝑙)! •𝑐 𝑗-𝑙 | G 𝑘 ] = 1 𝛽 𝑘 =𝛽 •       𝛽 ,𝑐 P(𝛽, (𝛽 , 𝑐)) • Ψ 𝑖+𝑗 (𝛽 ) ( 𝑗 -𝑙)! • 𝑐 𝑗-𝑙       (a.s.).
It follows that

E       ∞ 𝑗=𝑙 Ψ 𝑖+𝑗 ( 𝛽 𝑘+1 ) ( 𝑗 -𝑙)! • Cost( 𝛽 𝑘 ) 𝑗-𝑙 G 𝑘       = ∞ 𝑗=𝑙 𝛽 E 𝛽 ,𝑐 1 ( 𝛽 𝑘 , 𝛽 𝑘+1 ,Cost( 𝛽 𝑘 ))=(𝛽,𝛽 ,𝑐) • Ψ 𝑖+𝑗 (𝛽 ) ( 𝑗-𝑙)! • 𝑐 𝑗-𝑙 G 𝑘 = ∞ 𝑗=𝑙 𝛽 1 𝛽 𝑘 =𝛽 •       𝛽 ,𝑐 P(𝛽, (𝛽 , 𝑐)) • Ψ 𝑖+𝑗 (𝛽 ) ( 𝑗 -𝑙)! • 𝑐 𝑗-𝑙       = ∞ 𝑗=𝑙 𝛽 1 𝛽 𝑘 =𝛽 • E Ψ 𝑖+𝑗 (next 𝛽 ) • Cost(𝛽) 𝑗-𝑙 ( 𝑗 -𝑙)! = ∞ 𝑗=𝑙 E       Ψ 𝑖+𝑗 (next 𝛽 𝑘 ) • Cost( 𝛽 𝑘 ) 𝑗-𝑙 ( 𝑗 -𝑙)!       = E       ∞ 𝑗=𝑙 Ψ 𝑖+𝑗 (next 𝛽 𝑘 ) • Cost( 𝛽 𝑘 ) 𝑗-𝑙 ( 𝑗 -𝑙)!       ≤ Ψ 𝑖+𝑙 ( 𝛽 𝑘 ) (a.s.).

Thus we have

E[𝑌 𝑘+1,𝑖 | G 𝑘 ] = ∞ 𝑙=0 C 𝑙 𝑘 𝑙! • E ∞ 𝑗=𝑙 Ψ 𝑖+𝑗 ( 𝛽 𝑘+1 ) ( 𝑗-𝑙)! • Cost( 𝛽 𝑘 ) 𝑗-𝑙 G 𝑘 ≤ ∞ 𝑙=0 C 𝑙 𝑘 𝑙! • Ψ 𝑖+𝑙 ( 𝛽 𝑘 ) = 𝑌 𝑘,𝑖
Then by Optional Stopping Theorem (Theorem 3.3), we have that E[𝑌 𝜏,𝑖 ] ≤ E[𝑌 0,𝑖 ] for every 𝑖 ≥ 0. Finally, the theorem follows from the fact that 𝑌 𝜏,𝑖 ≥ 𝑓 (𝑖) (C 𝜏 ) (see property (b) above) and 𝑌 0,𝑖 = Ψ 𝑖 (𝛽 0 ). D.2 Full Proof of Theorem 3.5 P . Since all 𝑓 1 , . . . , 𝑓 𝑚 are non-negative, we have that Ψ 𝑖 (𝜖) ≥ 𝑡 𝑖 ★ ≥ 𝑓 (𝑖) (0), thus (A1) is satis ed. We verify (A2) by the de nition of next 𝛽 . If 𝑛 1 ≤ 𝑐 𝑟 1 , then next 𝛽 is deterministically equal to (𝑟 2 , 𝑛 2 ) • • • (𝑟 𝑘 , 𝑛 𝑘 ), and Cost(𝛽) = 0, thus for every 𝑖 ≥ 0, we have

Ψ 𝑖 (𝛽) = 𝑡 𝑖 ★ • exp 𝑡 ★ • 𝑘 𝑗=1 𝑓 𝑟 𝑖 (𝑛 𝑖 ) ≥ 𝑡 𝑖 ★ • exp 𝑡 ★ • 𝑘 𝑗=2 𝑓 𝑟 𝑖 (𝑛 𝑖 ) = E 𝑗 ≥0 Ψ 𝑘+𝑗 (next 𝛽 ) 𝑗! • Cost(𝛽) 𝑗 .
Otherwise (i.e., 𝑛 1 > 𝑐 𝑟 1 ), the random stack next 𝛽 consists of (𝑟 1 , ℎ 𝑟 1 ,1 (𝑛 1 )), . . . , (𝑟 1 , ℎ 𝑟 1 ,𝑆 𝑟 1 (𝑛 1 )) and (𝑞 𝑟 1 ,1 , 𝑔 𝑟 1 ,1 (𝑛 1 )), . . . , (𝑞 𝑟 1 ,𝑂 𝑟 1 , 𝑔 𝑟 1 ,𝑂 𝑟 1 (𝑛 1 )), followed by (𝑟 2 , 𝑛 2 ) . . . (𝑟 𝑘 , 𝑛 𝑘 ). Note that Cost(𝛽) =

Cost 𝑟 1 (𝑛 1 ). Thus, for every 𝑖 ≥ 0, we have:

Ψ 𝑖 (𝛽) = 𝑡 𝑖 ★ • exp(𝑡 • 𝑘 𝑖=1 𝑓 𝑟 𝑖 (𝑛 𝑖 )) = 𝑡 𝑖 ★ • exp(𝑡 ★ • 𝑘 𝑖=2 𝑓 𝑟 𝑖 (𝑛 𝑖 )) • exp(𝑡 ★ • 𝑓 𝑟 1 (𝑛 1 )) ≥ 𝑡 𝑖 ★ • exp(𝑡 ★ • 𝑘 𝑖=2 𝑓 𝑟 𝑖 (𝑛 𝑖 )) • E exp 𝑡 ★ • Cost 𝑟 1 (𝑛 1 ) + 𝑆 𝑟 𝑖=1 𝑓 𝑟 (ℎ 𝑟,𝑖 (𝑛 1 )) + 𝑂 𝑟 𝑖=1 𝑓 𝑞 𝑟,𝑖 (𝑔 𝑟,𝑖 (𝑛 1 )) = E 𝑗 ≥0 𝑡 𝑖+𝑗 ★ /𝑗! • Cost 𝑟 1 (𝑛 1 ) 𝑗 • exp(𝑡 ★ • ( 𝑘 𝑖=2 𝑓 𝑟 𝑖 (𝑛 𝑖 ) + 𝑆 𝑟 𝑖=1 𝑓 𝑟 (ℎ 𝑟,𝑖 (𝑛 1 )) + 𝑂 𝑟 𝑖=1 𝑓 𝑞 𝑟,𝑖 (𝑔 𝑟,𝑖 (𝑛 1 )))) = E 𝑗 ≥0 Ψ 𝑖+𝑗 (next 𝛽 ) 𝑗! • Cost(𝛽) 𝑗
where the third inequality holds by 𝑡 ★ ≤ 𝑡 𝑟 1 (𝑛 1 ) and the convexity of exp(•). In detail, we de ne 𝑋 𝑟 (𝑛) := Cost 𝑟 (𝑛)+ 𝑆 𝑟 𝑖=1 𝑓 𝑟 (ℎ 𝑟,𝑖 (𝑛))+ 𝑂 𝑟 𝑖=1 𝑓 𝑞 𝑟,𝑖 (𝑔 𝑟,𝑖 (𝑛))-𝑓 𝑟 (𝑛) to be the exponent of ( ‡) after dividing the both sides by the left hand side. Then ( ‡) could be reformulated as E[exp(𝑡 𝑟 (𝑛) • 𝑋 𝑟 (𝑛))] ≤ 1. Note that 𝑒 (𝑡) := E[exp(𝑡 • 𝑋 𝑟 (𝑛))] is a convex function of 𝑡, and 𝑒 (0) = 1, thus, if 𝑒 (𝑡 𝑟 (𝑛)) ≤ 1, for any 𝑡 ★ ≤ 𝑡 𝑟 (𝑛), 𝑒 (𝑡 ★ ) ≤ 1. Thus the lemma follows.

D.3 Full Proof of Theorem 3.6 P . We rstly verify (A1). Since 𝜃 (𝑖) (0) = 𝑘! if and only if 𝑖 = 𝑘, otherwise 𝜃 (𝑖) (0) = 0. For 𝑖 ≠ 𝑘, the (A1) is satis ed by non-negativity of 𝑓 's, for 𝑖 = 𝑘, (A1) holds by (4). We then verify (A2), for 𝑖 ≥ 𝑘 + 1, since Ψ 𝑖 (𝛽) = 0 for any stack 𝛽, the conditions naturally holds. For 0 ≤ 𝑖 ≤ 𝑘, consider 𝛽 := (𝑟 1 , 𝑛 1 ) . . . (𝑟 𝑡 , 𝑛 𝑡 ), if 𝑛 1 ≤ 𝑐 𝑟 1 then next 𝛽 deterministically equal to (𝑟 2 , 𝑛 2 ) . . . (𝑟 𝑡 , 𝑛 𝑡 ) and Cost(𝛽) = 0, hence, (A2) holds since

Ψ 𝑖 (𝛽) ≥ Ψ 𝑖 (next 𝛽 ), otherwise, next 𝛽 = (𝑟 2 , 𝑛 2 ) (𝑟 3 , 𝑛 3 ) • • • (𝑟 𝑘 , 𝑛 𝑘 ) (𝑟 1 , ℎ 𝑟 1 ,1 (𝑛 1 )) • • • (𝑟 1 , ℎ 𝑟 1 ,𝑆 𝑟 1 (𝑛 1 )) (𝑞 𝑟 1 ,1 , 𝑔 𝑟 1 ,1 (𝑛 1 )) • • • (𝑞 𝑟 1 ,𝑂 𝑟 1 , 𝑔 𝑟 1 ,𝑂 𝑟 1 (𝑛 1 )) and Cost(𝛽) = Cost 𝑟 1 (𝑛 1 )
. Hence, we have that:

Ψ 𝑖 (𝛽) = 𝑘! • 0≤𝑝 1 ,••• ,𝑝 𝑡 ≤𝑘-𝑖, 𝑡 𝑢=1 𝑝 𝑢 =𝑘-𝑖 𝑡 𝑢=1 𝑓 (𝑝𝑢 ) 𝑟𝑢 (𝑛 𝑢 ) 𝑝 𝑢 ! = 𝑘! • 𝑘-𝑖 𝑝 1 =0 𝑓 (𝑝 1 ) 𝑟 1 (𝑛 1 ) 𝑝 1 ! • ★ 𝑡 𝑢=2 𝑓 (𝑝𝑢 ) 𝑟𝑢 (𝑛 𝑢 ) 𝑝 𝑢 ! ≥ 𝑘! • 𝑘-𝑖 𝑝 1 =0 I Cost 𝑟 (𝑛) 𝑗 𝑗! • 𝑆 𝑟 𝑢=1 𝑓 (𝑝 𝑢 ) 𝑟 (ℎ 𝑟,𝑢 (𝑛)) 𝑝 𝑢 ! • 𝑂 𝑟 𝑢=1 𝑓 (𝑒 𝑢 ) 𝑞𝑟,𝑢 (𝑔 𝑟,𝑢 (𝑛)) 𝑒 𝑢 ! • I 𝑡 𝑢=2 𝑓 (𝑝𝑢 ) 𝑟𝑢 (𝑛 𝑢 ) 𝑝 𝑢 ! = E 𝑘-𝑖 𝑗=0 Cost 𝑟 1 (𝑛 1 ) 𝑗 𝑗! • Ψ 𝑖+𝑗 (next 𝛽 ) = E 𝑘-𝑖 𝑗=0 Ψ 𝑖+𝑗 (next 𝛽 ) 𝑗! • Cost(𝛽) 𝑗
where s

I := 0 ≤ 𝑗, 𝑝 1 , . . . , 𝑝 𝑆 𝑟 , 𝑒 1 , . . . , 𝑒 𝑂 𝑟 ≤ 𝑖, 𝑆 𝑟 𝑢=1 𝑝 𝑢 + 𝑂 𝑟 𝑢=1 𝑒 𝑢 + 𝑗 = 𝑝 1 I := 0 ≤ 𝑝 2 , • • • , 𝑝 𝑡 ≤ 𝑘 -𝑖 -𝑝 1 , 𝑡 𝑢=2 𝑝 𝑢 = 𝑘 -𝑖 -𝑝 1
Thus, (A2) holds.

E DETALED DESCRIPTION OF ALGMOMENT(𝐿 𝑇 , 𝐴, 𝐵, 𝑁 ): APPROXIMATING HIGH-ORDER MOMENTS WITH PSEUDO-POLYNOMIALS In the veri cation procedure of our algorithm, to approximate the high-order moments of 𝑓 (𝑖), where 𝑓 is a pseudo-monomial and 𝑖 observes as an U () or MU () distribution, we need to tightly approximate the summation 𝑛-1 𝑖=0 𝑖 𝑎 • ln 𝑏 𝑖 with pseudo-polynomials. In this part, we devise an algorithm A M (𝐿 𝑇 , 𝑎, 𝑏, 𝑛) to solve this problem through integrals and taylor expansions, where 𝑎, 𝑏 represents the degree for 𝑖 and ln 𝑖, and 𝐿 𝑇 represents the degree of approximation. i.e. our algorithm will approximate

𝑛-1 𝑖=0 𝑖 𝑎 • ln 𝑏 𝑖 (16)
with approximation degree 𝐿 𝑇 . The goal of this algorithm is to synthesize pseudo-polynomial upper and lower bounds for (18) with additive constant error. To brie y introduce the algorithm, to approximate (18), our algorithm naturally chooses to use its corresponding integral ∫ 𝑛 𝑑 𝑥 𝑎 • ln 𝑏 𝑥 d𝑥, since the latter could be automatically computed as a pseudo-polynomial through integration by parts. Then our algorithm would measure the di erence between the summation and integral through Taylor's expansion in each term of summation with degree 𝑎 + 𝐿 𝑇 , and controls the additive error for minor terms by Riemann-Zeta functions. To avoid corner cases, our algorithm assumes that 𝑛 ≥ exp(𝑏/𝐿 𝑇 ) + 1. De ne 𝑑 := exp(𝑏/𝐿 𝑇 ) . It rstly splits the summation into two parts:

𝑛-1 𝑖=0 𝑖 𝑎 • ln 𝑏 = 𝑑-1 𝑖=0 𝑖 𝑎 • ln 𝑏 𝑖 + 𝑛-1 𝑖=𝑑 𝑖 𝑎 • ln 𝑏 𝑖 (17)
where the rst part could be calculated trivially. We now focus on how to approximate the second term.

𝑛-1 𝑖=𝑑 𝑖 𝑎 • ln 𝑏 𝑖 (18)

We will rstly introduce mathematical backgrounds of this algorithm, and then give an overview of this algorithm, and nally provide algorithmic details.

Background. Below we present mathematical backgrounds that our algorithm heavily relies on. (Integration of 𝑖 𝑎 • 𝑙𝑛 𝑏 𝑖.) There is an algorithm that given 𝑎, 𝑏, 𝑑, where 𝑎, 𝑏 represents the degree of 𝑖 and ln 𝑖 and 𝑑 is the lower bound in integration, and computes the integral ∫ 𝑛 𝑑 𝑥 𝑎 • ln 𝑏 𝑥d𝑥 with pseudo-polynoimal. This directly follows from integration by parts. (Derivatives of 𝑖 𝑎 • ln 𝑏 𝑖.) There is an algorithm that given 𝑎, 𝑏, 𝑘 and symbolically expands the 𝑘-th derivative of 𝑖 𝑎 • ln 𝑏 𝑖 , i.e. (𝑖 𝑎 • ln 𝑏 𝑖) (𝑘) with pseudo-polynoimal. Below, we will use 𝛽 𝑎,𝑏,𝑘,𝑗 to represent the coe cient of 𝑖 𝑎-𝑘 • ln 𝑗 𝑖 in the 𝑘-th derivative of 𝑖 𝑎 • ln 𝑏 𝑖. Formally, we have that:

(𝑖 𝑎 • ln 𝑏 𝑖) (𝑘) = 𝑏 𝑗=0 𝛽 𝑎,𝑏,𝑘,𝑗 • 𝑖 𝑎-𝑘 • ln 𝑗 𝑖 (19)
This follows straightforwardly by the multiplication rule for derivatives. (Monotonicity of 𝑖 𝑎 • ln 𝑏 𝑖 for 𝑎 ≤ -1 and 𝑏 ≥ 0.) . For every 𝑎 ≤ -1 and 𝑏 ≥ 0, 𝑖 𝑎 • ln 𝑏 𝑖 is monotonically increasing in (0, exp(-𝑏/𝑎)] and monotonically decreasing in [exp(-𝑏/𝑎), +∞). It follows from the derivative fo 𝑖 𝑎 • ln 𝑏 𝑖.

(Constant Approximation of 𝑛-1 𝑖=𝑑 𝑖 𝑎 • ln 𝑏 𝑖 for 𝑎 ≤ -2.) For every 𝑎 ≤ -2 and 𝑏 ≥ 0, we could approximate the summation through 𝑏-th derivative of Riemman-Zeta function at 𝑥 = -𝑎, i.e. 𝜁 (𝑏) (-𝑎). Formally, we have that 0

≤ 𝑛-1 𝑖=𝑑 𝑖 𝑎 • ln 𝑏 𝑖 ≤ 𝜁 (𝑏) (-𝑎) -𝑑-1 𝑖=0 𝑖 -𝑎 • ln 𝑏 𝑖. (20)
This directly follows from the identity below:

𝜁 (𝑏) (-𝑎) = ∞ 𝑖=0 𝑖 𝑎 • ln 𝑏 𝑖 (𝑎, 𝑏 ∈ R, 𝑎 < -1)
Details. Below we will introduce details of A M . To approximate (18), Our algorithm naturally chooses to use its corresponding integral ∫ 𝑛 𝑑 𝑥 𝑎 • ln 𝑏 𝑥 d𝑥, since by Fact 1, the latter could be automatically computed as a pseudo-polynomial through integration by parts. The rest of our algorithm will carefully measure the di erence between these two terms. To achieve this, our algorithm will do the following 3 steps.

Step 1. (Subtraction and Reorganization). Our algorithm will rstly subtract one from the other, and by splitting the integral into 𝑛 parts and reorganizing the terms appropriately, we derive the formula below:

∫ 𝑛 𝑑 𝑥 𝑎 • ln 𝑏 𝑥 d𝑥 -𝑛-1 𝑖=𝑑 𝑖 𝑎 • ln 𝑏 𝑖 (21) = 𝑛-1 𝑖=𝑑 ∫ 𝑖+1 𝑖 𝑥 𝑎 • ln 𝑏 𝑥 d𝑥 -𝑛-1 𝑖=𝑑 𝑖 𝑎 • ln 𝑏 𝑖 (22) = 𝑛-1 𝑖=𝑑 ∫ 1 0 (𝑖 + 𝑥) 𝑎 ln 𝑏 (𝑖 + 𝑥) -𝑖 𝑎 ln 𝑏 𝑖 d𝑥 (23)
Step 2. (Approximating the Integration). Then, we will focus on ( 23) and try to bound this term from below and above by pseudo-polynomials. In general, our algorithm will expand the term (𝑖 + 𝑥) 𝑎 ln 𝑏 (𝑖 + 𝑥) -𝑖 𝑎 ln 𝑏 𝑖 within the integral and and do substitutions. This is achieved through two substeps.

Step 2(a). (Expanding through Taylor's Expansion). Firstly, our algorithm expands this term through Taylor's Expansion with degree 𝑎 + 𝐿 𝑇 -1 and Lagrange's Remainder:

(𝑖 + 𝑥) 𝑎 ln 𝑏 (𝑖 + 𝑥) -𝑖 𝑎 ln 𝑏 𝑖 = 𝑎+𝐿 𝑇 -1 𝑘=1 (𝑡 𝑎 ln 𝑏 𝑡 ) (𝑘 ) | 𝑡 =𝑖 𝑘! • 𝑥 𝑘 + (𝑡 𝑎 ln 𝑏 𝑡 ) (𝑎+𝐿 𝑇 ) | 𝑡 =𝜉 𝑥,𝑖 (𝑎+𝐿 𝑇 )! • 𝑥 𝑎+𝐿 𝑇 (24)
where for each 𝜉 𝑖,𝑥 is bounded in [0, 1], i.e. 0 ≤ 𝜉 𝑖,𝑥 ≤ 1.

Step 2(b). (Substitution). In the following, our algorithm would do a chain of substitutions. Firstly, by plugging in (19), we derive an equivalent form of (24):

𝑎+𝐿 𝑇 -1 𝑘=1 𝑏 𝑗=0 𝛽 𝑎,𝑏,𝑘,𝑗 𝑘! • 𝑖 𝑎-𝑘 • ln 𝑗 𝑖 • 𝑥 𝑘 + 𝛽 𝑎,𝑏,𝑎+𝐿 𝑇 ,𝑗 (𝑎+𝐿 𝑇 )! • 𝑖 + 𝜉 𝑖,𝑥 -𝐿 𝑇 • ln 𝑏 𝑖 + 𝜉 𝑖,𝑥 • 𝑥 𝑎+𝐿 𝑇 (25)
Third, by plugging in ( 25) into ( 23) and changing the order of summation, we derive that:

∫ 1 0 𝑎+𝐿 𝑇 -1 𝑘=1 𝑏 𝑗=0 𝛽 𝑎,𝑏,𝑘,𝑗 𝑘! • 𝑛-1 𝑖=𝑑 𝑖 𝑎-𝑘 • ln 𝑗 𝑖 • 𝑥 𝑘 + 𝛽 𝑎,𝑏,𝑎+𝐿 𝑇 ,𝑗 (𝑎 + 𝐿 𝑇 )! • 𝑛-1 𝑖=𝑑 𝑖 + 𝜉 𝑖,𝑥 -𝐿 𝑇 •ln 𝑏 𝑖 +𝜉 𝑖,𝑥 •𝑥 𝑎+𝐿 𝑇 d𝑥
Finally, by the identity that ∫ 1 0 𝑥 𝑘 d𝑥 = 1/(𝑘 + 1). We derive that:

𝑎+𝐿 𝑇 -1 𝑘=1 𝑏 𝑗=0 𝛽 𝑎,𝑏,𝑘,𝑗 (
𝑘+1)! • 𝑛-1 𝑖=𝑑 𝑖 𝑎-𝑘 • ln 𝑗 𝑖 + 𝑏 𝑗=0 𝛽 𝑎,𝑏,𝑎+𝐿 𝑇 ,𝑗 (𝑎+𝐿 𝑇 +1)! • 𝑛-1 𝑖=𝑑 𝑖 + 𝜉 𝑖,𝑥 -𝐿 𝑇 •ln 𝑗 𝑖 +𝜉 𝑖,𝑥 𝑆𝐴 𝑆𝐵
Step 3. (Bounding through Subinstances and Riemaman-Zeta Function). We will now focus on the equality above, note that in this formula, 𝛽 𝑎,𝑏,𝑘,𝑗 are constants, thus, it su ces to separately bound 𝑆𝐴 and 𝑆𝐵 from below and above by pseudo-polynomials.

Bounding 𝑆𝐴. Note that bounding 𝑆𝐴 could be directly solved through the subinstace A M (𝐿 𝑇 , 𝑑, 𝑎 -𝑘, 𝑗). Bounding 𝑆𝐵. Bounding 𝑆𝐵 is more complicated, it is solved through three substeps. Below we will illustrate how to bound 𝑆𝐵 from above, the procedure of getting 𝑆𝐵's lower bound is similar.

-Firstly, we replace 𝛽 𝑎,𝑏,𝑎+𝐿 𝑇 ,𝑗 with max{𝛽 𝑎,𝑏,𝑎+𝐿 𝑇 ,𝑗 , 0}, then it su ces to overapproximate the summation 𝑛-1 𝑖=𝑑 𝑖 + 𝜉 𝑖,𝑥 -𝐿 𝑇 •ln 𝑗 𝑖 +𝜉 𝑖,𝑥 .

-Then, by Fact 3 and our assumption that 𝑑 ≥ exp(𝑏/𝐿 𝑇 ), we have that 𝑖 + 𝜉 𝑖,𝑥

-𝐿 𝑇 • ln 𝑗 𝑖 +𝜉 𝑖,𝑥 ≤ 𝑖 -𝐿 𝑇 • ln 𝑏 𝑖.
-Finally, from Fact 4, we have that

𝑛-1 𝑖=𝑑 𝑖 + 𝜉 𝑖,𝑥 -𝐿 𝑇 •ln 𝑗 𝑖 +𝜉 𝑖,𝑥 ≤ 𝑛-1 𝑖=𝑑 𝑖 -𝐿 𝑇 • ln 𝑏 𝑖 ≤ 𝜁 (𝑏) (𝐿 𝑇 ) -𝑑-1 𝑖=0 𝑖 -𝐿 𝑇 • ln 𝑏 𝑖.
In summary, we could bound 𝑆𝐵 from above by: 𝑏 𝑗=0 max{𝛽 𝑎,𝑏,𝑎+𝐿 𝑇 ,𝑗 , 0}

(𝑎 + 𝐿 𝑇 + 1)! • 𝜁 (𝑏) (𝐿 𝑇 ) - 𝑑-1 𝑖=0 𝑖 -𝐿 𝑇 • ln 𝑏 𝑖
and 𝑆𝐵 could also be similarly bounded from below by: 𝑏 𝑗=0 min{𝛽 𝑎,𝑏,𝑎+𝐿 𝑇 ,𝑗 , 0}

(𝑎 + 𝐿 𝑇 + 1)! • 𝜁 (𝑏) (𝐿 𝑇 ) - 𝑑-1 𝑖=0 𝑖 -𝐿 𝑇 • ln 𝑏 𝑖
By further plugging in the corresponding pseudo-polynomial for the integral ∫ 𝑛 𝑑 𝑥 𝑎 • ln 𝑏 𝑥d𝑥, our algorithm could successfully derive a tight approximation of (18) with additive constant error.

F DETAILED DESCRIPTION OF ALGCHECK(𝑃, 𝑁 0 ): CHECKING THE NON-POSITIVITY OF PSEUDO-POLYNOMIAL 𝑃 In the veri cation procedure of our algorithm, we need to check whether the given pseudopolynomial 𝑃 is non-positive for every 𝑛 ≥ 𝑛 0 . Below we x a univariate pseudo-polynomial 𝑃 (𝑛) := 𝑎,𝑏 ∈Z 𝑐 𝑎,𝑏 • 𝑛 𝑎 • ln 𝑏 𝑛. To brie y introduce our algorithm, it rstly obtains an equivalent form 𝑃 1 (𝑛) by extracting the term with largest order and dividing the polynomial with the largest terms. Then, our algorithm follows the naive idea that try to enumerate every 𝑛 ≥ 𝑛 0 and evaluate 𝑃 1 (𝑛). However, the enumeration involves in nite steps. Thus, we propose a condition that decides the termination of enumeration. The condition considers an over-approximation 𝑃 2 (𝑛) of 𝑃 1 (𝑛) by removing all non-constant terms with negative coe cients, and a constant 𝑛 𝑒 that serves as an underapproximation of monotonically decreasing interval of 𝑃 2 (𝑛). i.e., 𝑃 2 (𝑛) will monotonically decreasing for every 𝑛 ≥ 𝑛 𝑒 . Intuitively, the condition states that if 𝑃 2 (𝑛 1 ) ≤ 0 for some 𝑛 1 ≥ 𝑛 𝑒 , then 𝑃 2 (𝑛) ≤ 0 for every 𝑛 ≥ 𝑛 1 , combining with the fact that 𝑃 2 (𝑛) is an over-approximation of 𝑃 1 (𝑛), we could conclude that 𝑃 1 (𝑛) is also non-positive for every 𝑛 ≥ 𝑛 1 , thus we could safely terminate the enumeration. This algorithm also applies to check non-positivity only for odd/even 𝑛's, by simply enumerating only odd/even numbers 𝑛 ≥ 𝑛 0 rather than every 𝑛 ≥ 𝑛 0 . Please also note that our algorithm is complete, in the sense that if 𝑃 1 (𝑛) is always non-positive, the algorithm will eventually terminate.

Our algorithm relies on the following monotonicity property of 𝑖 𝑎 • ln 𝑏 𝑖.

Monotonicity of 𝑖 𝑎 • ln 𝑏 𝑖. For every 𝑎 ≤ -1 and 𝑏 ≥ 0, 𝑖 𝑎 • ln 𝑏 𝑖 is monotonically increasing in (0, exp(-𝑏/𝑎)] and monotonically decreasing in [exp(-𝑏/𝑎), +∞). This follows directly through computing the derivative. In detail, the algorithm is as the follows.

Step 1. (Extracting and Dividing the Major Term). Our algorithm rstly nds the term with largest order. In detail, our algorithm will scan each term in the polynomial and nds the term with highest degree of 𝑛, if there are several terms with the same highest degree of 𝑛, our algorithm chooses the term among them with the highest degree of ln 𝑛. Formally, our algorithm nds the term 𝑛 𝑎 * • ln 𝑏 * 𝑛 with 𝑐 𝑎 * ,𝑏 * ≠ 0 and highest 𝑎 * , if there are more than one terms such that their degree of 𝑛 are all equal to 𝑎 * , we choose the largest 𝑏 * among them. Then, our algorithm divides each term in 𝑃 (𝑛) by 𝑛 𝑎 * • ln 𝑏 * 𝑛, and derives the pseudo-polynomial 𝑃 1 (𝑛) below:

𝑃 1 (𝑛) = 𝑐 𝑎 * ,𝑏 * + 𝑎,𝑏≠𝑎 * ,𝑏 * 𝑐 𝑎,𝑏 • 𝑛 𝑎-𝑎 * • ln 𝑏-𝑏 * 𝑛
Note that by our choice of 𝑎 * , 𝑏 * , the degree of 𝑛 and ln 𝑛 in each non-constant term in 𝑃 1 (𝑛), i.e., 𝑎 -𝑎 * , 𝑏 -𝑏 * either satis es 𝑎 -𝑎 * ≤ -1, or satis es 𝑎 -𝑎 * = 0 and 𝑏 -𝑏 * ≤ -1, also deciding non-positivity of 𝑃 (𝑛) is equivalent to deciding non-positivity of 𝑃 1 (𝑛).

Step 2. (Bounded Checking). In general, to check whether 𝑃 1 (𝑛) is always non-positive for every 𝑛 ≥ 𝑛 0 , our algorithm follows the naive idea that try to enumerate every 𝑛 ≥ 𝑛 0 and evaluate 𝑃 1 (𝑛). However, the enumeration involves in nite steps. Thus, we propose a condition (Proposition F.2) that decides the termination of enumeration. The condition considers an over-approximation 𝑃 2 (𝑛) of 𝑃 1 (𝑛) by removing all non-constant terms with negative coe cients, and a constant 𝑛 𝑒 that serves as an underapproximation of monotonically decreasing interval of 𝑃 2 (𝑛). i.e., 𝑃 2 (𝑛) will monotonically decreasing for every 𝑛 ≥ 𝑛 𝑒 . Intuitively, the condition states that if 𝑃 2 (𝑛 1 ) ≤ 0 for some 𝑛 1 ≥ 𝑛 𝑒 , then 𝑃 2 (𝑛) ≤ 0 for every 𝑛 ≥ 𝑛 1 , combining with the fact that 𝑃 2 (𝑛) is an over-approximation of 𝑃 1 (𝑛), we could conclude that 𝑃 1 (𝑛) is also non-positive for every 𝑛 ≥ 𝑛 1 , thus we could safely terminate the enumeration. Please also note that our algorithm is complete, in the sense that if 𝑃 1 (𝑛) is always non-positive, the algorithm will eventually terminate. We will rstly introduce two propositions deciding the termination of enumeration. P F.1. Let 𝑃 2 (𝑛) to be an over-approximation of 𝑃 1 (𝑛) by removing all non-constant terms with negative coe cients.i.e. P . Note that every coe cient for non-constant terms in 𝑃 2 (𝑛) is non-negative. Thus the proposition directly follows from the monotonicity of the function 𝑖 𝑎 • ln 𝑏 𝑖. P F.2. If 𝑃 2 (𝑛 1 ) ≤ 0 for some 𝑛 1 ≥ 𝑛 𝑒 , then 𝑃 2 (𝑛 1 ) ≤ 0 for every 𝑛 1 ≥ 𝑛 𝑒 .

P

. The proposition directly follows from the monotonicity of 𝑃 2 (𝑛).

In detail, our algorithm follows the 3 substeps below.

Step 2(a). (Checking Constants). If 𝑐 𝑎 * ,𝑏 * > 0, then we could directly conclude that 𝑃 1 (𝑛) > 0 for some su ciently large 𝑛. This is due to the fact that the degree of 𝑛 and ln 𝑛 in each non-constant term in 𝑃 1 (𝑛), i.e., 𝑎 -𝑎 * , 𝑏 -𝑏 * either satis es 𝑎 -𝑎 * ≤ -1, or satis es 𝑎 -𝑎 * = 0 and 𝑏 -𝑏 * ≤ -1, thus is negligible for when 𝑛 is su ciently large.

Step 2(b). (Calculating the lower bound 𝑛 𝑒 ). Following Proposition F.1, we calculate the 𝑛 𝑒 by scanning every element in 𝑃 2 (𝑛), and trivially evaluate 𝑃 1 (𝑛) for every 𝑛 0 ≤ 𝑛 ≤ 𝑛 𝑒 -1, if 𝑃 1 (𝑛) > 0 for some 1 ≤ 𝑛 ≤ 𝑛 𝑒 , then we could conclude that 𝑃 1 (𝑛) > 0 for some 𝑛 ≥ 𝑛 0 .

Step 2(c). (Enumeration from 𝑛 𝑒 ). Our algorithm then tries every 𝑛 ≥ 𝑛 𝑒 , and do the following judgements:

• If 𝑃 1 (𝑛) > 0, then we could conclude that 𝑃 1 (𝑛) > 0 for some 𝑛 ≥ 𝑛 0 and terminates the enumeration.

• If 𝑃 2 (𝑛) ≤ 0, then through Proposition F.2, we could conclude that 𝑃 2 (𝑛) ≤ 0 for every 𝑛 ≥ 𝑛 0 , combining with the fact that 𝑃 2 (𝑛) is an over-approximation of 𝑃 1 (𝑛), we could conclude that 𝑃 1 (𝑛) is also non-positive for every 𝑛 ≥ 𝑛 1 , and terminates the enumeration. The soundness is straightforward. The completeness of our algorithm is stated below: T F.3 (S ). If 𝑃 1 (𝑛) ≤ 0 for every 𝑛 ≥ 𝑛 0 , then our algorithm will terminates.

P

. We just need to check whether step 2(c) will eventually terminate. Since we pass the check at step 2(a), we have that 𝑐 𝑎 * ,𝑏 * < 0. Then due to the fact that the degree of 𝑛 and ln 𝑛 in each non-constant term in 𝑃 1 (𝑛), i.e., 𝑎 -𝑎 * , 𝑏 -𝑏 * either satis es 𝑎 -𝑎 * ≤ -1, or satis es 𝑎 -𝑎 * = 0 and 𝑏 -𝑏 * ≤ -1, thus is negligible for when 𝑛 is su ciently large. Thus 𝑃 2 (𝑛) will be negative for some su ciently large 𝑛, thus the algorithm will terminate in nite steps. The goal of this algorithm is to synthesize a pseudo-monomial 𝐵 𝑟 (𝑛) such that -𝐵 𝑟 (𝑛) ≤ 𝑋 𝑟 (𝑛) ≤ 𝐵 𝑟 (𝑛) for every 𝑛 ≥ 𝑛 0 .

The organization of this section is as follows. We rstly brie y introduce the mathematical background for this algorithm, and then gives a simpli ed algorithm for a restricted setting. Finally we incrementally add features to this algorithm.

Mathematical Backgrounds. Here we present the mathematical backgrounds that our algorithm relies on.

The rst is a proposition that approximates ln(𝑎 • 𝑥 + 𝑏) with 𝑎 ≥ 0 for positive 𝑥, through pseudo-polynomials from above and below. To achieve this, rst note that the case where 𝑎 = 0 or 𝑏 = 0 is trivial. We will next show how to solve case for 𝑎 > 0 and 𝑏 ≠ 0, the following proposition splits the case where 𝑏 > 0 or 𝑏 < 0. First we consider corner cases, if 𝑢 = 0 and 𝑣 = 0, the proposition obviously holds. Then if 𝑢 ≥ 1 and 𝑣 = 0, we have that 𝑓 (𝑥) + 𝑓 (𝑛 -1 -𝑥) = 𝑢 • (𝑥 𝑢-1 -(𝑛 -1 -𝑥) 𝑢-1 ), thus 𝑓 (𝑥) + 𝑓 (𝑛 -1 -𝑥) < 0 for every 0 ≤ 𝑥 ≤ (𝑛 -1)/2, and 𝑓 (𝑥) + 𝑓 (𝑛 -1 -𝑥) > 0 for every (𝑛 -1)/2 ≤ 𝑥 ≤ 𝑛 -1, thus the proposition holds.

If 𝑢 ≥ 1 and 𝑣 ≥ 1. we have that 𝑓 (𝑥) + 𝑓 (𝑛 -1 -𝑥) = 𝑥 𝑢-1 • ln 𝑣-1 𝑥 • (𝑢 • ln 𝑥 +𝑣) -(𝑛 -1 -𝑥) 𝑢-1 • ln 𝑣-1 (𝑛 -1 -𝑥) • (𝑢 • ln(𝑛 -1 -𝑥) + 𝑣), thus 𝑓 (𝑥) + 𝑓 (𝑛 -1 -𝑥) < 0 for every 0 ≤ 𝑥 ≤ (𝑛 -1)/2, and 𝑓 (𝑥) + 𝑓 (𝑛 -1 -𝑥) > 0 for every (𝑛 -1)/2 ≤ 𝑥 ≤ 𝑛 -1, thus the proposition holds.

A Simplifed Version of Our Algorithm. At the beginning, we present our algorithm in a restricted setting, we assume that:

• Every random variable in 𝑋 𝑟 (𝑛) observes as a discrete distribution.

• The support of all discrete distributions does not involve oors and ceils, thus each element in the support is a linear expression 𝑎 • 𝑛 + 𝑏 where 𝑎 ∈ N and 𝑏 ∈ Z. • The number of self-recursion calls is 1, i.e. 𝑆 𝑟 = 1 Our algorithm is as the follows.

Step 1. (Enumerating the Support). Since every random variable observes as a discrete distribution, our algorithm rstly enumerate every possible combination of these random variables. Formally, we enumerate every 𝜔 in the sample space and considers: Below our algorithm will try to bound 𝑋 𝑟 (𝑛) (𝜔) for some xed 𝜔, by enumeration and simply taking the maximum for all 𝜔's, we could derive the bound 𝐵 𝑟 (𝑛).

Step 2. (Approximating 𝑋 𝑟 (𝑛) (𝜔) to Pseudo-Polynomial from Above and Below). Note that since each element in the support is a linear expression 𝑎 • 𝑛 + 𝑏 where 𝑎 ∈ N and 𝑏 ∈ Z, and recall that 𝑓 𝑟 (𝑛) is a pseudo-monomial, each term in 𝑋 𝑟 (𝑛) (𝜔) is in the form 𝜆 • (𝑎 • 𝑛 + 𝑏) 𝑢 • ln 𝑣 (𝑐 • 𝑛 + 𝑑), for real 𝜆, natural numbers 𝑎, 𝑐, 𝑢, 𝑣 and integers 𝑏, 𝑑. First note that it is straightforward to expand (𝑎 • 𝑛 + 𝑏) 𝑢 as pseudo-polynomials through Newton's binomial Theorem. Then, by combining Proposition H.1, our algorithm could approximate ln(𝑐 • 𝑛 + 𝑑) from above and below, and then through Newton's binomial theorem, we could expand the 𝑣-th power of the approximate form of ln(𝑐 •𝑛+𝑑) as a pseudo-polynomial. During this procedure, we also records a value ℓ as the maximum value of -2•𝑑 𝑐 among all terms with 𝑑 < 0, which represents the range that our proposition for approximation (Proposition H.1) does not apply.

Below we represent the pseudo-polynomial lower and upper bounds for 𝑋 𝑟 (𝑛) (𝜔) as 𝐿 𝜔 𝑟 (𝑛) and 𝑅 𝜔 𝑟 (𝑛). Please note that this bound only applies to 𝑛 ≥ ℓ.

Step 3. (Ge ing the Monomial Bound). In this step, we will further overapproximate -𝐿 𝜔 𝑟 (𝑛) and 𝑅 𝜔 𝑟 (𝑛) as pseudo-monomials, by taking the maximum of these two monomials, we derive the monomial bound 𝐵 𝜔 𝑟 (𝑛) for 𝑋 𝑟 (𝑛) (𝜔) such that -𝐵 𝜔 𝑟 (𝑛) ≤ 𝑋 𝑟 (𝑛) (𝜔) ≤ 𝐵 𝜔 𝑟 (𝑛). Below we illustrate the procedure that overapproximating 𝑅 𝜔 𝑟 (𝑛) as pseudo-monomials, the overapproximation for -𝐿 𝜔 𝑟 (𝑛) exactly follows the same procedure. We represent the pseudo-bound synthesis result for 𝑅 𝜔 𝑟 (𝑛) as 𝜃 • 𝑛 𝑢 𝑟,𝜔 • ln 𝑣 𝑟,𝜔 𝑛.

Below we represent 𝑅 𝜔 𝑟 (𝑛) as 𝑅 𝜔 𝑟 (𝑛) := 𝑎,𝑏 ∈Z 𝑐 𝑎,𝑏 • 𝑛 𝑎 • ln 𝑏 𝑛. In detail, this is achieved through 3 substeps.

Step 3(a). (Extracting the Major Term). Our algorithm rstly nds the term with largest order. In detail, our algorithm will scan each term in the polynomial and nds the term with highest degree of 𝑛, if there are several terms with the same highest degree of 𝑛, our algorithm chooses the term among them with the highest degree of ln 𝑛. Formally, our algorithm nds the term 𝑛 𝑎 * • ln 𝑏 * 𝑛 with 𝑐 𝑎 * ,𝑏 * ≠ 0 and highest 𝑎 * , if there are more than one terms such that their degree of 𝑛 are all equal to 𝑎 * , we choose the largest 𝑏 * among them. After nding the major term 𝑛 𝑎 * • ln 𝑏 * 𝑛, we could derive that 𝑢 𝑟,𝜔 = 𝑎 * and 𝑣 𝑟,𝜔 = 𝑏 * .

  s. This follows from (E1), (E3) and the almost-sure uniqueness in Theorem D.1. (E6) (Linearity) E [𝑎 • 𝑋 + 𝑏 • 𝑌 | G] = 𝑎 • E [𝑋 | G] +𝑏 • E [𝑌 | G] a.s. for all 𝑎, 𝑏 ≥ 0. This follows from the linearity of integration. (E7) (Positivity) E [𝑋 | G] ≥ 0 (a.s.). This follows from exactly the same proof under 'A Positivity Result' on [Williams 1991, Page 87]. (E8) (cMON) If we have 𝑋 𝑛 ≤ 𝑋 𝑛+1 (a.s.) for all 𝑛 ≥ 0 and lim 𝑛→∞ 𝑋 𝑛 = 𝑋 (a.s.), then it holds that E [𝑋 𝑛 | G] ≤ E [𝑋 𝑛+1 | G] (a.s.) and lim 𝑛→∞ E [𝑋 𝑛 | G] = E [𝑋 | G] (a.s.

  By enumerating all possible stacks for 𝛽 𝑘 , 𝛽 𝑘+1 and values for Cost( 𝛽 𝑘 ) § ∞ 𝑗=𝑙 E 𝛽,𝛽 ,𝑐 1 ( 𝛽 𝑘 , 𝛽 𝑘+1 ,Cost( 𝛽 𝑘 ))=(𝛽,𝛽 ,𝑐) • Ψ 𝑖+𝑗 (𝛽 ) ( 𝑗-𝑙)! • 𝑐 𝑗-𝑙 G 𝑘 = § By applying (E6)-(Linearity) and (E8)-(cMON) § ∞ 𝑗=𝑙 𝛽 E 𝛽 ,𝑐 1 ( 𝛽 𝑘 , 𝛽 𝑘+1 ,Cost( 𝛽 𝑘 ))=(𝛽,𝛽 ,𝑐) • Ψ 𝑖+𝑗 (𝛽 ) ( 𝑗-𝑙)! • 𝑐 𝑗-𝑙 G 𝑘

𝑃 2 (

 2 𝑛) := 𝑐 𝑎 * ,𝑏 * + 𝑎,𝑏≠𝑎 * ,𝑏 * max{0, 𝑐 𝑎,𝑏 } • 𝑛 𝑎-𝑎 * • ln 𝑏-𝑏 * 𝑛 Let 𝑛 𝑒 := max{exp(-(𝑏 -𝑏 * )/(𝑎 -𝑎 * )) | (𝑎, 𝑏) ≠ (𝑎 * , 𝑏 * ) and 𝑐 𝑎,𝑏 ≠ 0} be a constant. Then we have that 𝑃 2 (𝑛) monotonically decreases for every 𝑛 ≥ 𝑛 𝑒 .

G

  PROOF OF THE PROPOSITION THAT APPROXIMATES LOGARITHM Proof of Proposition 5.1: P. We rstly prove the case for 𝑏 > 0. The lower bound is trivial since 𝑏 > 0. For the upper bound, we have that:ln(𝑎 • 𝑥 + 𝑏) -ln(𝑎 • 𝑥) = ln 1 + 𝑏 𝑎 • 𝑥 ≤ 𝑏 𝑎 • 𝑥For the case 𝑏 < 0, the upper bound is also derive through the same procedure as the case 𝑏 > 0. For the lower bound, rst note that ln(𝑎 • 𝑥) = ln(𝑎 • 𝑥 + 𝑏 -𝑏), and:ln((𝑎 • 𝑥 + 𝑏) -𝑏) -ln(𝑎 • 𝑥 + 𝑏) = ln 1 -𝑏 𝑎 • 𝑥 + 𝑏 ≤ -𝑏 𝑎 • 𝑥 + 𝑏Thus we have that: ln(𝑎 • 𝑥 + 𝑏) ≥ ln 𝑎 + ln 𝑥 + 𝑏 𝑎 • 𝑥 + 𝑏 Since we assume that 𝑥 ≥ -2•𝑏 𝑎 , we have that 2•𝑏 𝑎 •𝑥 < 𝑏 𝑎 •𝑥+𝑏 , thus we have that: ln(𝑎 • 𝑥 + 𝑏) ≥ ln 𝑎 + ln 𝑥 + 2 • 𝑏 2 • 𝑥 H A METHOD THAT AUTOMATICALLY CALCULATING THE BOUND 𝐵 𝑅 (𝑁 ) In this part, we propose an algorithm that automatically calculates 𝐵 𝑟 (𝑛) (the "Setting up Template" paragraph in Section 5.3). We will rstly propose a simpli ed version of our algortihm under restricted settings, and then incrementally add features to support the whole setting for our algorithm. Firstly recall that 𝑋 𝑟 (𝑛) = Cost 𝑟 (𝑛) + 𝑆 𝑟 𝑖=1 𝑓 𝑟 (ℎ 𝑟,𝑖 (𝑛)) + 𝑂 𝑟 𝑖=1 𝑓 𝑞 𝑟,𝑗 (𝑔 𝑟,𝑗 (𝑛)) -𝑓 𝑟 (𝑛)

P H. 1 .

 1 For every 𝑎 > 0, 𝑏 > 0, and 𝑥 > 0, we have that:ln 𝑎 + ln 𝑥 ≤ ln(𝑎 • 𝑥 + 𝑏) ≤ ln 𝑎 + ln 𝑥 + 𝑏 𝑎 • 𝑥 Also, for every 𝑎 > 0, 𝑏 < 0, and 𝑥 ≥ -2•𝑏 𝑎 , we have that:ln 𝑎 + ln 𝑥 + 2 • 𝑏 𝑎 • 𝑥 ≤ ln(𝑎 • 𝑥 + 𝑏) ≤ ln 𝑎 + ln 𝑥 + 𝑏 𝑎 • 𝑥 P. We rstly prove the case for 𝑏 > 0. The lower bound is trivial since 𝑏 > 0. For the upper bound, we have that:ln(𝑎 • 𝑥 + 𝑏) -ln(𝑎 • 𝑥) = ln 1 + 𝑏 𝑎 • 𝑥 ≤ 𝑏 𝑎 • 𝑥For the case 𝑏 < 0, the upper bound is also derive through the same procedure as the case 𝑏 > 0. For the lower bound, rst note that ln(𝑎 • 𝑥) = ln(𝑎 • 𝑥 + 𝑏 -𝑏), and:ln((𝑎 • 𝑥 + 𝑏) -𝑏) -ln(𝑎 • 𝑥 + 𝑏) = ln 1 -𝑏 𝑎 • 𝑥 + 𝑏 ≤ -𝑏 𝑎 • 𝑥 + 𝑏 Thus we have that: ln(𝑎 • 𝑥 + 𝑏) ≥ ln 𝑎 + ln 𝑥 + 𝑏 𝑎 • 𝑥 + 𝑏 Since we assume that 𝑥 ≥ -2•𝑏 𝑎 , we have that 2•𝑏 𝑎 •𝑥 < 𝑏 𝑎 •𝑥+𝑏 , thus we have that: ln(𝑎 • 𝑥 + 𝑏) ≥ ln 𝑎 + ln 𝑥 + 2 • 𝑏 2 • 𝑥 P H.2.For any pseudo-monomial 𝑓 (𝑛) := 𝑛 𝑢 • ln 𝑣 𝑛 with 𝑢, 𝑣 ∈ N and 𝑢 ≥ 1, 𝑓 (𝑥) + 𝑓 (𝑛 -1 -𝑥) achieves its maximum at 𝑥 = 0 and minimum at 𝑥 = 𝑛-1 2 , for 0 < 𝑥 < 𝑛 -1. P . Computing the derivative of 𝑓 (𝑥), we have that𝑓 (𝑥) = 𝑢 • 𝑥 𝑢-1 ln 𝑣 𝑥 + 𝑣 • 𝑥 𝑢-1 • ln 𝑣-1 𝑥 = 𝑥 𝑢-1 • ln 𝑣-1 𝑥 • (𝑢 • ln 𝑥 + 𝑣)

  𝑋 𝑟 (𝑛) (𝜔) = Cost 𝑟 (𝑛) (𝜔) + 𝑆 𝑟 𝑖=1 𝑓 𝑟 (ℎ 𝑟,𝑖 (𝑛) (𝜔)) + 𝑂 𝑟 𝑖=1 𝑓 𝑞 𝑟,𝑖 (𝑔 𝑟,𝑖 (𝑛) (𝜔)) -𝑓 𝑟 (𝑛)

  𝑗 • 𝑂 𝑢

	𝑢=1	𝑓 𝑞𝑟,𝑢 ,𝑒𝑢 (𝑔 𝑟,𝑢 (𝑛)) 𝑒 𝑢 !	, and (ii) E	𝑆 𝑟 𝑢=1	𝑓 𝑟,𝑝𝑢 (ℎ 𝑟,𝑢 (𝑛)) 𝑝 𝑢 !

Table 2 .

 2 Symbolic Exponentially Decreasing Concentration Bounds. For each benchmark, we provide results based on several di erent over-approximations 𝑓 (𝑛) of the expected runtime.

	4 ) 𝛼 -4

Table 3 .

 3 Concrete Exponentially Decreasing Concentration Bounds. Table2shows our symbolic exponential bounds, which were obtained using the algorithm of Section 5.3. It also provides Karp's bounds for the same SPRRs, demonstrating that our bounds are much tighter. To further illustrate this point, we provide concrete bounds (with xed values for the parameters 𝛼 and

	Experimental results. Tables 2-4 illustrate our experimental results over the benchmarks.
	• Exponentially decreasing concentration bounds:

𝛽) in Table 3. As shown in this table, our bounds consistently beat previous results. Notably, in some cases, our bound gets in nitely tighter than Karp's bound as 𝑛 * grows. Additionally, in two cases (S S and R ), Karp's method is not applicable. This is because they either involve more than one recurrence equation or have randomized preprocessing times. * Pr[ C 𝜏 ≥ 𝛼 • 8 • 𝑛 * ] * Pr[ C 𝜏 ≥ 𝛼 • 36 • 𝑛 * ]

Table 4 .

 4 Polynomially Decreasing Concentration Bounds.

  Δ,𝛽 * (𝐴) = 0 if 𝛽 0 ≠ 𝛽 * , and Pr Δ,𝛽 * (𝐴) = 𝑘-1 𝑗=0 P(𝛽 𝑗 , 𝛽 𝑗+1 ) if 𝛽 0 = 𝛽 * . In[START_REF] Agrawal | Lexicographic ranking supermartingales: an e cient approach to termination of probabilistic programs[END_REF], conditional expectation is extended to non-negative random variables that are not necessarily integrable. A random variable 𝑋 is non-negative if 𝑋 (𝜔) ≥ 0 for all elements 𝜔 in the sample space.

	D DETAILED PROOF IN SECTION 3
	D.1 Full proof of Theorem 3.4
	T	D.1 ([A

Step 3(b). (Considering Corner Cases). We enumerate every 𝑛 0 ≤ 𝑛 ≤ ℓ, evaluate the coe cient 𝑋 𝑟 (𝑛) (𝜔) 𝑛 𝑎 * •ln 𝑏 * 𝑛 , take the maximum among them as the lower bound 𝜃 0 for the coe cient 𝜃 .

Step 3(c). (Binary Search). Finnaly, we binary search 𝜃 , the veri cation proceudre of binary search simply calls the procedure A C (𝑅 𝜔 𝑟 (𝑛) -𝜃 •𝑛 𝑎 * •ln 𝑏 * 𝑛, max{ℓ, 𝑛 0 }) to check the non-positivity, thus concluding whether the currenct choice of 𝜃 leads to a bound. By taking the maximum of 𝜃 and 𝜃 0 to tackle with corner cases, we derive the nal monomial bound for 𝑅 𝜔 𝑟 (𝑛). Below, we will then try to remove these restrictions one by one, thus giving the complete version of our algorithm.

Supporting Floors and Ceils. To support oors and ceils, same as our method for in E S and P S , we just need to case analysis whether 𝑛 is even or odd and takes the maximum.

Supporting U

(𝑛) and MU (𝑛) Distribution. To support U (𝑛) distribution, since we are only interested in the bound, it su ces to consider the case when the value is 0 or 𝑛 -1, thus reducing to previous cases. To support MU (𝑛) distribution, rst by case analysis whether 𝑛 is even or odd, we remove the oors in this distribution, and then we still only need to consider two endpoints in this distribuiton.

Supporting the Case 𝑆 𝑟 = 2. In this case, after the processing of 𝑋 𝑟 (𝑛) in the case 𝑆 𝑟 = 2 of the algorithm E S , if ℎ 𝑟,1 (𝑛) observes as discrete distribution, we could directly reduce to previous cases. Otherwise, if ℎ 𝑟,1 (𝑛) observes as U (𝑛) or MU (𝑛), since 𝑓 𝑟 (𝑛) is an over-approximation of expected running time, and Cost 𝑟 (𝑛) + 𝑂 𝑟 𝑗=1 𝑇 𝑞 𝑟,𝑗 (𝑔 𝑟,𝑗 (𝑛)) ≥ 1, we have that 𝑓 𝑟 (𝑛) = Ω(𝑛), thus the Proposition H.2 applies, in this case, we only need to consider the case when the value is 𝑓 (0) + 𝑓 (𝑛 -1) or 2 • 𝑓 ( 𝑛-1

2 ), thus reducing to previous cases.