# Taming reluctant random walks in the positive quadrant

Jeremie Lumbroso, Marni Mishna, Yann Ponty

► To cite this version:

Jeremie Lumbroso, Marni Mishna, Yann Ponty. Taming reluctant random walks in the positive quadrant. GASCOM - 10th conference on random generation of combinatorial structures - 2016, Jun 2016, Bastia, France. pp.99 – 114, 10.1016/j.endm.2017.05.008 . hal-01291164

HAL Id: hal-01291164

https://inria.hal.science/hal-01291164

Submitted on 21 Mar 2016

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

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
TAMING RELUCTANT RANDOM WALKS IN THE POSITIVE
QUADRANT

JEREMIE LUMBROSO, MARNI MISHNA AND YANN PONTY

ABSTRACT. A lattice walk model is said to be reluctant if the defining step set has a strong drift towards the boundaries. We describe efficient random generation strategies for these walks.

# 1. INTRODUCTION

Walks on lattices are fundamental combinatorial classes. They appear in many guises particularly in formal language theory, queuing theory, and combinatorics as they naturally encode common relations. A typical lattice path model is a set of walks defined by a fixed, finite set of allowable moves (called the *step set*), and a region to which the walks are confined (typically a convex cone). The exact and asymptotic enumeration of lattice paths restricted to the first quadrant (known as quarter plane models) have been a particularly active area of study of late because of some new and interesting techniques coming from different areas of computer algebra and complex analysis [5, 7, 18, 20].

Efficient uniform random generation is useful to study the typical large scale behavior of walks under different conditions. Models which restrict walks to the upper half plane can be specified by an algebraic combinatorial grammar [11, 3]. Consequently, efficient random generation schemes can be obtained using several systematic strategies, such as recursive generation [14] and Boltzmann sampling [12].

Intriguingly, walks restricted to the first quadrant are more complex. Rare is the quarter-plane model with an algebraic generating function that cannot be trivially reformulated as a half-plane model. Overwhelmingly, the cyclic lemma, combinatorial identities and other grammar-based techniques that are so fruitful in the half-plane case, do not easily apply. Furthermore, there is only a small proportion of models whose generating function satisfies a differential equation with polynomial coefficients<fnref n="1" />, again excluding a potential source of direct, generic generation techniques [1].

Rejection sampling is the term for a general technique where one generates from a simpler superclass, and then rejects elements until an element from the desired class is obtained. In the case of lattice paths, a naive rejection strategy could use unrestricted walks as a superset. This is only practical for those quarter-plane models whose counting sequences grow essentially like those of the unrestricted walks. Such is the case when the *drift*, or vector sum of the stepset, is positive coordinate-wise. Anticipated rejection can also be used when the drift is 0, and provides a provably efficient algorithm [2]. However, any such strategy is demonstrably doomed to failure when the drift of the step set is negative in any coordinate, as the probability of generating long unconstrained walks

<footnote n="1">¹For example, amongst the 20 804 small step models with less than 6 steps in 3 dimensions, only around 150 appear to be D-finite [4].</footnote>
<figure><img src="image_2.png" /><figcaption>FIGURE 1. A RANDOM WALK WITH 18 000 STEPS IN THE QUARTERPLANE USING THE STEPSET $S = \{(1, 0), (0, 1), (-1, 0), (1, -1), (-1, -1), (-2, -1)\}$.</figcaption></figure>

which remain in the first quadrant becomes exponentially small. One strategy in the literature is to change the probability on the allowable steps, and consequently forgo the uniformity of the generation [6]. It appears then that the problem of efficient, uniform random generation algorithms for generic quarter plane lattice path models is a relatively undeveloped topic.

**Our contribution.** The main result of this paper is an efficient rejection algorithm for the uniform random generation of walks in the quarter plane. It is an application of recent results due Johnson, Mishna and Yeats [17], Garbit and Raschel [15] amongst others. It is provably efficient and straightforward to implement. It is most impressive on walks whose drift is negative in both coordinates, a property we call *reluctant*, but it also offers notable gains for any model which tends to either boundary.

More precisely, we describe a strategy in which every walk of length $n$ is generated with equal likelihood. The efficacy result holds for quarter plane models with any step set, and is easily generalized to higher dimensions. Figure 1 illustrates a walk of over 18000 steps that was generated uniformly at random for the quarter plane model with reluctant step set

$$S = \{(1, 0), (0, 1), (-1, 0), (1, -1), (-1, -1), (-2, -1)\}.$$

The probability of generating a walk of this length by rejection from the set of unrestricted sequences of steps $S^*$ is less than $\frac{5.3299^{18000}}{6} \sim 1.75 \cdot 10^{-926}$. However, with our strategy, it was generated (relatively quickly).

Rejection from an unrestricted walk is not the only competition. For the purposes of comparison, we describe a recursive strategy which requires exact enumeration results to be tabulated in advance. This has potential to be efficient, and is insensitive to the drift of the model, but does require a lot of storage. We discuss this algorithm in Section 3.2.

Our alternative sampler is based on a straightforward combinatorial interpretation of an enumerative result. Roughly, we use that for any quarter plane model, there is a
corresponding half plane model such that asymptotically, both models have the same exponential growth factor. This implies that a rejection strategy from this half plane has sub-exponential rate of rejection. The sub-exponential factors are conjectured to also match in many cases, suggesting that it is in fact a particularly efficient strategy.

Asymptotic enumerative results are recalled in the next section. A *baseline* algorithm, based on a trivial recurrence, is presented in Section 3. Section 4 describes our main rejection algorithm. Its practical implementation depends on the rationality of the slope of the half-plane model, and is discussed in Subsection 4.4. We conclude with some remarks regarding implementation aspects, along with possible extensions.

## 2. 2D LATTICE PATH BASICS

A 2D lattice path model is a combinatorial class consisting of walks on the 2D integer lattice, starting at the origin, and taking steps from some finite multi-set $S \subset \mathbb{Z}^2$ of allowable steps. In this work, we consider the restriction of such walks to the positive quadrant $Q = \mathbb{Z}^2_{\ge 0}$, although the general strategy works for a wider set of cones. We use the half-plane $H_\theta$ defined by a line through the origin:

$$H_{\theta} = \{x \sin \theta + y \cos \theta \ge 0\}.$$

For a fixed, finite step set $S \subset \mathbb{Z}^2$, a given cone $C$, and positive integer $n$ we define walks$(C, S, n)$ to be the class of walks of length $n$ starting at the origin, taking steps in $S$, and staying in $C$. Formally,

$$\text{walks}(C, S, n) = \{x_0, x_1, \dots, x_n \mid x_0 = (0,0) \wedge x_{j+1} - x_j \in S \wedge x_i \in C\}.$$

The complete class is given by walks$(C, S) = \bigcup_{n \ge 0} \text{walks}(C, S, n)$.

We use the following enumerative quantities in our analysis:

$$(1) \qquad q_n^S = |\text{walks}(Q, S, n)| \qquad h_n^S(\theta) = |\text{walks}(H_\theta, S, n)|.$$

The asymptotic regime for $q_n^S$ is always of the form

$$(2) \qquad q_n^S \sim \gamma \rho^{-n} n^{-r},$$

for real numbers $\rho$ and $r$. We refer to $\rho^{-1}$ as the exponential growth factor of the model. The asymptotic regime critically (although not exclusively) depends on the *drift* of the step set $S$, defined as $\text{drift}(S) = \sum_{s \in S} s$. A walk model walks$(Q, S)$ is said to be *reluctant* when $\text{drift}(S) = (\delta_1, \delta_2)$ with $\delta_1 < 0$ and $\delta_2 < 0$.

Reluctant models for the positive quadrant have exponential growth factors that are lower than the number of steps. It follows that the naive algorithm that performs rejection from unconstrained walks, has exponential time complexity, motivating our algorithmic contribution. Indeed, even when one of $\delta_1 < 0$ or $\delta_2 < 0$, the exponential growth factor can be less than the number of steps.

## 3. BASIC RECURSIVE RANDOM GENERATOR

The exact value of $q_n^S$ can be expressed using a recurrence. This motivates a straightforward instance of the recursive method [21, 14], where steps are simply drawn sequentially, using probabilities that depend both on the current position reached, and the number of remaining steps.
3.1. **Exact enumeration of walks.** Define $q_n^S(x, y)$ to be the number of positive suffixes of walks in walks$(Q, S, n)$ which start from the point $(x, y)$ and remain in the positive quadrant. Such suffix walks of length $n$ can be factored as a first step $(i, j) \in S$, keeping the walk in the positive quadrant, followed by another positive suffix of length $n-1$ starting at $(x + i, y + j)$. This leads to the recurrence:

$$ (3) \qquad q_n^S(x, y) = \begin{cases} \sum_{(i,j) \in S \text{ s.t. } x+i \ge 0, y+j \ge 0} q_{n-1}^S(x+i, y+j) & \text{if } n > 0, \\ 1 & \text{if } n = 0 \end{cases} $$

A quadrant walk is also the positive suffix of a walk starting at (0, 0), thus $q_n^S := q_n^S(0, 0)$. This recurrence can also be trivially adapted to handle general cones, higher dimensions, or for further constraining the end-point, e.g. to count/generate meanders, or walks ending on the diagonal.

3.2. **Algorithm and complexity analysis.** Once the cardinalities $q_n^S(x, y)$ are available, a uniform random walk is generated, by choosing one of the steps with probabilities proportional to the number of possible suffixes.

(1) **Preprocessing.** Precompute $q_{n'}^S(x, y)$ for each $n' \in [0, n]$ and $(x, y) \in [0, n \cdot a] \times [0, n \cdot b]$, where $a := \max_{(i,j) \in S} i$ and $b := \max_{(i,j) \in S} j$;

(2) **Generation.** Initially starting from (0, 0) and $n' := n$, iterate until $n' = 0$:

(a) Choose a step $(i, j) \in S$ with probability $q_{n'-1}^S(x + i, y + j)/q_{n'}^S(x, y)$;

(b) Add $(i, j)$ to the walk, update the current point $((x, y) := (x + i, y + j))$, and decrease the remaining length ($n' := n' - 1$);

**Theorem 1 (Complexity/correctness).** The random uniform generation of k 2-dimensional walks confined to the positive quadrant can be performed in $\Theta(k \cdot n + n^3)$ arithmetic operations, using storage for $\Theta(n^3)$ numbers.

*Proof.* The preprocessing stage should only be computed once in the generation of k sequences. It involves $\Theta(|S| \cdot n^{d+1})$ arithmetic operations, and requires storage for $\Theta(n^{d+1})$ large numbers. The generation of a single walk requires the generation of $\Theta(n)$ random numbers and, for each of them, their comparisons to $\Theta(|S|)$ other numbers.

An induction argument establishes the correctness of the algorithm. Assume that, for all $n' < N$ and $(x, y) \in [0, n' \cdot a] \times [0, n' \cdot b]$, the positive suffixes are uniformly generated, a fact that can be verified when $n' = 0$. Then for $n' = N$, the algorithm chooses a suitable step $(i, j) \in S$, and then recursively generates a – uniform from the induction hypothesis – suffix from the updated position. The probability of generating any such walk is therefore

$$ \mathbb{P}(w) = \frac{q_{N-1}^S(x + i, y + j)}{q_N(x, y)} \times \frac{1}{q_{N-1}^S(x + i, y + j)} = \frac{1}{q_N^S(x, y)} $$

and we conclude with the uniformity of the generation. $\square$

In practice however, the memory consumption of the algorithm grows in $\Theta(n^4)$ bits, which limits the utility of this strategy to $n < 500$. Thus the above algorithm only serves as a baseline for our alternative based on rejection.
## 4. EFFICIENT REJECTION SAMPLER FROM 1D MODELS

We recall some basics of rejection sampling for our analysis. Let $\mathcal{A}$ be a combinatorial class which contains the sub-class $\mathcal{C}$. Given a random sampler for $\mathcal{A}$, we can use a rejection strategy to make a random sampler for $\mathcal{C}$. Let $a_n$ and $c_n$ respectively count the number of elements of size $n$ in $\mathcal{A}$ and $\mathcal{C}$. Following Devroye [10, Chapter II.3], we say that class $\mathcal{A}$ efficiently covers $\mathcal{C}$ if

$$ \left( \frac{a_n}{c_n} \right) \in \mathcal{O}(n^p), $$

Here $p \ge 0$ is some constant independent of $n$. In other words, asymptotically, the expected number of elements drawn from $\mathcal{A}$ before generating an element in $\mathcal{C}$ is polynomial in $n$. Ideally $p$ is as small as possible.

### 4.1. Candidate superclass: Half-plane model.
Our algorithm arises from the surprising observation made by Johnson, Mishna, and Yeats [17], later proven by Garbit and Raschel [15]:

**Theorem 2 (Garbit and Raschel [15]).** Consider a step set $\mathcal{S}$, let $\rho(\theta)^{-1} := \lim_{n \to \infty} h_n^S(\theta)^{1/n}$ be the exponential growth factor of the half-plane model walks($H_\theta, \mathcal{S}$), and define

$$ \theta^* := \underset{0 \le \theta \le \pi/2}{\operatorname{argmax}} \rho(\theta), $$

Then the growth factor $\rho^{-1} := \lim_{n \to \infty} (q_n^S)^{1/n}$ of walks in the positive quadrant $\mathcal{Q}$ satisfies:

$$ (4) \qquad \rho = \rho(\theta^*). $$

This says that the exponential growth of the quarter-plane model is equal to the exponential growth of a superclass half-plane model. Furthermore the value of $\theta^*$ is explicitly computable.

**Corollary 3.** The combinatorial class walks($H_{\theta^*}, \mathcal{S}$) efficiently covers walks($\mathcal{Q}, \mathcal{S}$).

Next we consider the sub-exponential factors, as this gives the polynomial complexity of the rejection. On the side of the half-plane walks, the asymptotic formulas for $h_n^S(\theta)$ can be deduced from the complete generating function study of Banderier and Flajolet [3]. The sub-exponential factors are either $n^0, n^{-1/2}$, or $n^{-3/2}$, depending on the drift of the model (positive, zero and negative respectively).

For quarter-plane walks, the picture is less complete. The case of excursions for models with zero drift was described by Denisov and Wachtel [9], and from this work Duraj [13, Theorem II] was able to conclude explicit formulas for reluctant walks: Let $S(x, y) = \sum_{(i,j) \in \mathcal{S}} x^i y^j$, and let $(\alpha, \beta)$ be the unique positive critical point of $S(x, y)$. Such a point always exists, provided that $\mathcal{S}$ satisfies some non-triviality conditions. Then,

$$ (5) \qquad q_n^S \sim \gamma \rho^{-n} n^{-r}, $$

where $\rho$ and $r$ satisfy

$$ \rho = \frac{1}{S(\alpha, \beta)} \quad \text{and} \quad r = 1 + \pi \arccos \frac{S_{xy}(\alpha, \beta)}{\sqrt{S_{xx}(\alpha, \beta)S_{yy}(\alpha, \beta)}}. $$
**Data:** Reluctant step set $\mathcal{S} \subset \mathbb{Z}^2$, length $n$

**Result:** $w \in \text{walks}(\mathcal{Q}, \mathcal{S}, n)$ drawn uniformly at random

// Determine optimal slope $m = \arctan(\theta^*)$ following [17]

if $\mathcal{S}$ is singular then $m \leftarrow 0$;

else

Set $S(x, y) = \sum_{(i,j) \in \mathcal{S}} x^i y^j$;

Determine $(x, y) = (\alpha, \beta)$, the unique positive solution of

$$ \frac{d}{dx}S(x, y) = \frac{d}{dy}S(x, y) = 0; $$

if $\beta = 1$ then $m = \infty$;

else $m = \ln \alpha / \ln \beta$;

end

// Create suitable grammar $\mathcal{G}$

if $m = \infty$ then $p \leftarrow 1$ and $q \leftarrow 0$;

else if $m$ is rational then find $p, q \in \mathbb{N}$ so that $m = p/q$;

else find $p/q$, a $1/\sqrt{n}$-rational approximation to $m$;

$\mathcal{A} \rightarrow \{ip + jq : (i, j) \in \mathcal{S}\}$;

$\mathcal{G} \rightarrow \text{grammar}(\mathcal{A})$;

// Main rejection loop

repeat

| $w \rightarrow \text{UniformDraw}(\mathcal{G}, n)$

until $2\text{DMap}(w) \in \mathcal{Q}$;

**Algorithm 1:** Outline of our rejection algorithm. `UniformDraw(\mathcal{G}, n)` denotes a uniform sampler of walks of length $n$ for the grammar $\mathcal{G}$, and $2\text{DMap}(w)$ indicates the reinterpretation of $w$ as a sequence of 2D steps.

4.2. **The algorithm.** Algorithm 1 implements a classic rejection from a carefully-chosen half-plane model.

**Theorem 4 (Complexity of Algorithm 1).** Let $\mathcal{S}$ be a reluctant walk model and let $M(n)$ denote the time complexity of generating a walk for the half-plane model $H_{\theta^*}$. The expected time taken by Algorithm 1 to generate a walk in the positive quadrant is in $\Theta(M(n) \times n^{r-3/2})$.

This immediately follows from formula (5), from which we deduce that the expected number of trials is $h_n^S(\theta)/q_n^S \in \Theta(n^{r-3/2})$. For reluctant small step models, one has $3.3 < r < 7.5$. More recently, Garbit and Raschel have conjectured formulas for the sub-exponential factor in the general case, and remarkably suggest that for many (non-reluctant) models, $h_n^S(\theta)/q_n^S \in \mathcal{O}(1)$.

Next we address the efficient uniform random generation of walks in $\text{walks}(H_{\theta^*}, \mathcal{S})$.

4.3. **Half plane models as unidimensional walks.** We now describe efficient samplers for half-plane models $\text{walks}(H_{\theta}, \mathcal{S})$. Remark that walks in any half-plane can be generated as positive 1D walks, by taking 1D steps that are the orthogonal projections of those in $\mathcal{S}$ onto the half-plane boundary.

A unidimensional model is defined by a set $\mathcal{A} \subset \mathbb{R}$. The nontriviality conditions imply that $\mathcal{A}$ contains both a positive and a negative element. The associated class of walks
begins at 0, takes steps which are elements from $\mathcal{A}$ such that the sum over any prefix of the walk is nonnegative. If $\mathcal{A}$ is a multiple of a set of integers, then the class is modelled by a context free grammar, which we describe in the next section. Otherwise, the class cannot be trivially modelled by a context-free grammar, as is proven in Section 4.5.1.

Given $\mathcal{S}$ and $\theta$, we define the associated unidimensional model

$$ \mathcal{A}(\theta) = \{i \sin \theta + j \cos \theta : (i, j) \in \mathcal{S}\}. $$

The classes walks($H_\theta$, $\mathcal{S}$) and walks($\mathcal{A}(\theta)$, $\mathbb{R}_{\ge 0}$) are in a straightforward bijection.

**Remark 1** ([17]). If $\mathcal{S}$ defines a non-trivial 2D quarterplane model, then $\mathcal{A}(\theta)$ defines a non-trivial unidimensional model. Moreover if $\mathcal{S}$ is relunctant, then the drift of $\mathcal{A}(\theta)$ is negative. Finally multiplying steps by a positive constant does not affect the language of positive walks.

Two cases arise, depending on whether or not $\mathcal{A}(\theta^*)$ consists of rational-valued steps (up to rescaling). This is equivalent to asking if the slope of the boundary of the half-plane, $m = \tan(\theta^*)$ is rational.

### 4.4. Case 1: Rational projected steps.
When $\tan(\theta^*)$ is rational, the steps in $\mathcal{A}(\theta^*)$ can be scaled to be integers, as mentioned in Remark 1, therefore we consider unidimensional models $\mathcal{A}$ which consist of integer-valued steps.

Combinatorial specifications and, specifically, context-free grammars can then be used for random generation. Context-free grammars are indeed suitable to describe objects following rules which depend on a single, integer counter—and place certain, finite constraints on this counter. For the purpose of random walks, this counter may typically keep track of the height of the walk, and be constrained to always remain positive (i.e., the walk remains above the x-axis). From a grammar, random objects can be sampled using a variety of generic methods. More generally, this is equivalent to saying that grammars can describe walks that are confined within a half-plane.

To build the grammar grammar($\mathcal{A}$) for a unidimensional model defined by step set $\mathcal{A}$, we first distinguish the positive, negative and neutral steps

$$ \mathcal{A}^+ := \{a \mid a \in \mathcal{A} \text{ and } w(a) > 0\}, $$

$$ \mathcal{A}^- := \{a \mid a \in \mathcal{A} \text{ and } w(a) < 0\}, $$

$$ \mathcal{A}^0 := \{a \mid a \in \mathcal{A} \text{ and } w(a) = 0\}, $$

and define the largest upward and downward step lengths

$$ \bar{a} := \max \mathcal{A}^+ $$

$$ \bar{b} := -\min \mathcal{A}^- $$

Note that both of these lengths are positive, and are well-defined when the step set satisfies the conditions of non-triviality.
<figure><img src="image_3.png" /><figcaption>FIGURE 2. Typical decomposition of a walk in $\mathcal{D}$ with first step of height of height $a \ge k$.</figcaption></figure>

Using these three sets, and these two values we define the associated grammar $\mathbf{grammar}(\mathcal{A})$, whose terminals are given by $\mathcal{A}$, and non-terminals are defined as follows:

$$ \mathcal{P} = \mathcal{D} \times \mathcal{P}_{\text{aux}} \qquad \mathcal{L}_i = \sum_{\substack{a \in \mathcal{A} \\ w(a)=i}} a + \sum_{k=i+1}^{\min(\bar{a}, i+\bar{b})} \mathcal{L}_k \mathcal{R}_{k-i} $$

$$ \mathcal{P}_{\text{aux}} = \varepsilon + \sum_{k=1}^{\bar{a}} \mathcal{L}_k \times \mathcal{P}_{\text{aux}} \qquad \mathcal{R}_j = \sum_{\substack{b \in \mathcal{A} \\ w(b)=-j}} b + \sum_{k=j+1}^{\min(j+\bar{a}, \bar{b})} \mathcal{L}_{k-j} \mathcal{R}_k $$

$$ \mathcal{D} = \sum_{\substack{c \in S \\ w(c)=0}} c \times \mathcal{D} + \sum_{k=1}^{\max(\bar{a}, \bar{b})} \mathcal{L}_k \times \mathcal{D} \times \mathcal{R}_k \times \mathcal{D} $$

This follows from Duchon [11], Bousquet-Mélou and Ponty [8], with minor corrections to the indices that prevent the grammar from referencing undefined rules. The decomposition of a walk is unique and a schematic of a typical decomposition is presented in Figure 2.

Given the step set $\mathcal{A}$, $\mathcal{G} = \mathbf{grammar}(\mathcal{A})$ can be built in constant time (proportional to $\max(\bar{a}, \bar{b})^2$). To generate an element from a context free grammar, one either uses recursive methods [21] or Boltzmann generation [12]. The grammar here is straightforward, so most common optimizations apply [16].

**Theorem 5** (Complexity of rational half-plane sampling). Let walks$(\mathbb{R}_{\ge 0}, \mathcal{A})$ be a non-trivial unidimensional model defined by a rational multiset $\mathcal{A} \subset \mathbb{Z}$. The uniform random generation of $k$ walks of length $n$ in walks$(\mathbb{R}_{\ge 0}, \mathcal{A})$ can be performed in $O(k \cdot n \log n)$ arithmetic operations using storage for $O(1)$ numbers.

**Corollary 6.** When the step set $\mathcal{S}$ yields a rational $\mathcal{A}$ unidimensional projection, Algorithm 1 generates $k$ walks in the positive quadrant using $O(k \cdot n^{r-1/2} \log n)$ arithmetic operations, where $r$ is the exponent of the subexponential term in the asymptotics of walks$(\mathcal{Q}, \mathcal{S})$.

4.5. **Case 2: Non-rational projected steps.** When the projected step set $\mathcal{A}(\theta^*)$ contains non-rational steps, then the associated language is not context-free, and grammars can no longer be used directly. However, it is still possible to use a rational approximation of the perfect half-plane model, at the expense of the algorithmic efficiency.

4.5.1. *Contextuality of associated languages.*
**Lemma 7.** Let $S \subset \mathbb{Z}^2$ be a finite set which defines a non-trivial quarterplane model. Let $\theta^*$ be angle determined by Theorem 2 and assume furthermore that $m = \tan(\theta^*)$ is irrational. Then, the language $\mathcal{L}_m$ whose alphabet is made from the pairs $(i, j) \in S$, and the words are restricted to walks in walks($H_{\theta^*}$, $S$) is not context-free.

*Proof.* Consider two steps $a, b \in \mathcal{A}$, encoded by symbols $s_a$ and $s_b$, such that $a > 0$ and $b < 0$ and $a/b$ is irrational. The existence of such steps follows from the non-triviality of $S$, and the irrationality of $\tan(\theta^*)$. First, recall that the intersection of a context-free language and a rational language is a context-free language. If the intersection language

$$ \mathcal{L}_m^{\cap} = \{s_a^* s_b^*\} \cap \mathcal{L}_m = \{s_a^i s_b^j \mid a \cdot i - b \cdot j \ge 0\} $$

is not context-free, then neither is $\mathcal{L}_m$.

The fact that $\mathcal{L}_m^{\cap}$ is not context free can be proven using the context-free version of the pumping lemma, which states that, if $\mathcal{L}_m^{\cap}$ is context-free, then there exists a word length $p$ above which each word $w \in \mathcal{L}_m^{\cap}$ can be decomposed as $w = x.u.y.v.z$ such that $|u.y.v| \le p$, $|u.v| \ge 1$, and $\{x.u^i.y^i.v^i.z \mid i \in \mathbb{N}\} \subset \mathcal{L}_m^{\cap}$.

Let $\Delta(w) = |w|_{s_a} \cdot a - |w|_{s_b} \cdot b$ denote the (signed) final distance to the half plane, we establish the following technical lemma.

**Lemma 8.** For any $p \ge 0$, there exists a word $w^* \in \mathcal{L}_m^{\cap}$, $|w^*| > p$, such that $\Delta(w^*) < \Delta(w)$ for all $w \in \mathcal{L}_m^{\cap}$, $|w| < |w^*|$.

*Proof.* Assume that $p$ is given, and let $\Delta^{\le p}$ denote the smallest distance to the half plane of a word of length $\le p$, reached by some word $s_a^{x\bullet} s_b^{y\bullet} \in \mathcal{L}_m^{\cap}$ of length $x\bullet + y\bullet \le p$.

First we constructively show the existence of a word of length greater than $p$, whose final distance to the half-plane is smaller than $\Delta^{\le p}$. Consider the word

$$ w^{\circ} := s_a^{K \cdot x \bullet} s_b^{K \cdot y \bullet} s_b \quad \text{where} \quad K := \left\lfloor \frac{b}{\Delta^{\le p}} \right\rfloor. $$

Since both the slope and ratio $a/b$ are irrational, then $\Delta^{\le p} \ne 0$ and such a word exists. The final distance to the half plane of $w^{\circ}$ is given by:

$$ \Delta(w^{\circ}) = K \cdot \Delta^{\le p} - b = \left( \left\lfloor \frac{b}{\Delta^{\le p}} \right\rfloor - \frac{b}{\Delta^{\le p}} \right) \cdot \Delta^{\le p} < \Delta^{\le p}. $$

Consider now the smallest word $w^* \in \mathcal{L}_m^{\cap}$ such that $\Delta(w^*) < \Delta^{\le p}$. Such a word exists since $\Delta(w^{\circ}) < \Delta^{\le p}$ and clearly obeys $|w^*| > p$. Since $w^*$ is the smallest word such that $|w^*| > p$ and $\Delta(w^*) < \Delta^{\le p}$, then one has $\Delta(w^*) < \Delta^{\le |w^*|}$ which proves our claim. $\square$

Let us now investigate the possible factorizations as $x.u.y.v.z$ of the word $w^*$, whose existence is established by Lemma 8, and show that neither of them satisfies the pumping lemma. Focusing on $u$ and $v$, remark that neither of them should simultaneously feature both kinds of steps, otherwise any word $w^{[1]} = x.u^2.y.v^2.z \notin \mathcal{L}_m^{\cap}$, as it would feature at least two peaks. It follows that any satisfactory decomposition must be of the form $u = s_a^i$ and $v = s_b^j$, and the non-rationality of $a/b$ implies that $\Delta(u.v) \ne 0$. If $\Delta(u.v) < 0$, then the word $w^{[2]} = x.u^r.y.v^r.z$, $r > \lceil \Delta(w^*) / \Delta(u.v) \rceil$ is such that $\Delta(w^{[2]}) < 0$, and therefore $w^{[2]} \notin \mathcal{L}_m^{\cap}$. If $\Delta(u.v) > 0$, then let us observe that $u.v \in \mathcal{L}_m^{\cap}$ and has total length $i + j < p$, therefore Lemma 8 implies that $\Delta(u.v) > \Delta^{|w^*|}$. It follows that the word $w^{[3]} = x.u^0.y.v^0.z$ has final distance to the slope $\Delta^{|w^*|} - \Delta(u.v) < 0$ and therefore $w^{[3]} \notin \mathcal{L}_m^{\cap}$. Having found all
possible decompositions lacking in some respect, we conclude that $\mathcal{L}_m^\cap$ is not a context-free language, and neither is $\mathcal{L}_m$. $\square$

**4.6. Rational approximations.** All is not lost in the case of an irrational slope model, however, as we can define an approximation to the slope that is sufficiently close to the optimal slope to ensure polynomial-time rejection.

**Definition (δ-rational approximation).** A half-plane model $H_{\theta_r}(S)$ is a $\delta$-rational approximation of a half-plane model $H_{\theta}(S)$ if and only if $m_r := \tan \theta_r \in \mathbb{Q}$ and $|\tan \theta - \tan \theta_r| \le \delta$.

**Proposition 9.** For any model walks$(H_\theta, S)$ and $\delta > 0$ a desired precision, there exists a grammar with $O(1/\delta)$ non-terminals and $O(1/\delta^2)$ rules, which generates a $\delta$-rational approximation walks$(H_{\theta_r}, S)$ of walks$(H_\theta, S)$.

Remark that, as soon as $|\tan \theta - \tan \theta_r| > 0$, the exponential growth factor of the half-plane model becomes greater than that of the quarter-plane model, and Algorithm 1 becomes exponential on $n$. On the other hand, for any length $n \in \mathbb{N}$, setting $\delta := 1/(n+1)$ will define a model $H_{\theta_r}(S)$ which coincides with $H_\theta(S)$ on positive walks of length $n$. Indeed, the accumulated error due to the approximation of the step set remains too small to lead to the acceptance of some walk in $H_{\theta_r}(S)$ and not in $H_\theta(S)$ (and vice-versa).

**Theorem 10 (Complexity of $1/(n+1)$-rational approximation).** Let walks$(\mathbb{R}_{\ge 0}, A, n)$ be a non-trivial unidimensional model defined by a non-rational multiset $A$. The uniform random generation of $k$ walks of length $n$ in walks$(\mathbb{R}_{\ge 0}, A, n)$ can be performed in:

* $O(k \cdot n \log n + n^3)$ arithmetic operations, using storage for $O(n^3)$ large integers [14];

* $O(k \cdot n^3 \log n + n^2)$ arithmetic operations, using storage for $O(n^2)$ large integers [16];

* $O(k \cdot n^2 + n^2)$ arithmetic operations, using storage for $O(n^2)$ real values (Oracle) [12, 19].

Finally, we conjecture that a polynomial rejection is actually reached using a $(1/\sqrt{n})$-rational approximation.

**Conjecture 1.** Let $h_n(\theta)$ be the number of walks of length $n$ in an half-plane model $H_\theta(S)$, there exists an infinite sequence of angles $\{\theta_n\}_{n \le 0}$ such that:

* For all $n \ge 0$, $H_{\theta_n}(S)$ is a $(1/\sqrt{n})$-rational approximation of $H_{\theta^*}(S)$;

* The number of rejections in Algorithm 1 remains polynomial in $n$:

$$ \exists p \in \mathbb{R}, \lim_{n \to +\infty} \frac{h_n(\theta_n)}{h_n(\theta)} \in O(n^p) $$

## 5. REMARKS AND FUTURE EXTENSIONS

We have implemented this algorithm in Python, with external calls to Sage to compute $\theta^*$, and to Maple for Boltzmann Generation. Experimentally, in the case of irrational projected steps, even crude approximations led to much improved empirical complexities than both the default half-plane generation and the naive recursive generator. On the other hand, increasingly precise approximations led to an overwhelming growth in the size of the grammar, as could be expected from the asymptotic complexity. This raises interesting questions about the precise interplay between the size of the grammar and the
complexity, starting with our conjecture which we hope to address in a future version of
this work.

There are many possible optimizations, notably, anticipated rejection. We expect this
should have a positive effect on the complexity, particularly in the null-drift cases, possibly
after projection onto the targeted half-plane.

Natural extensions and generalizations. Finally, many natural extensions come to mind. Generating excursions in the quarter plane is difficult, but using our grammar-based approach it is completely straightforward. Finally, there are analogous “best hyperplane” theorems in higher dimensions, and for more general cones, and our general approach could in principle generalize to these cases.

**Acknowledgements.** We are very grateful for discussions with Kilian Raschel and Julien Courtiel.

## REFERENCES

[1] Axel Bacher, Olivier Bodini, and Alice Jacquot. Exact-size sampling for motzkin trees in linear time via boltzmann samplers and holonomic specification. In *ANALCO*, pages 52–61. SIAM, 2013.

[2] Axel Bacher and Andrea Sportiello. Anticipated rejection algorithms and the darling?mandelbrot distribution. In *AofA*, 2014.

[3] Cyril Banderier and Philippe Flajolet. Basic analytic combinatorics of directed lattice paths. *Theoret. Comput. Sci.*, 281(1-2):37–80, 2002. Selected papers in honour of Maurice Nivat.

[4] Alin Bostan, Mireille Bousquet-Mélou, Manuel Kauers, and Stephen Melczer. Three dimensional lattice walks restricted to the positive octant. *Annals of Combinatorics*, to appear.

[5] Alin Bostan and Manuel Kauers. Automatic classification of restricted lattice walks. In *DMTCS Proceedings of the 21st International Conference on Formal Power Series and Algebraic Combinatorics (FPSAC'09)*, Hagenberg, Austria, pages 203–217, 2009.

[6] Mireille Bousquet-Mélou. On the importance sampling of self-avoiding walks. *Combinatorics, Probability and Computing*, pages 1–24, 2011.

[7] Mireille Bousquet-Mélou and Marni Mishna. Walks with small steps in the quarter plane. In *Algorithmic probability and combinatorics*, volume 520 of *Contemp. Math.*, pages 1–39. Amer. Math. Soc., Providence, RI, 2010.

[8] Mireille Bousquet-Mélou and Yann Ponty. Culminating paths. *Discrete Math. Theor. Comput. Sci.*, 10(2):125–152, 2008.

[9] Denis Denisov and Vitali Wachtel. Random walks in cones. *Ann. Probab.*, 43(3):992–1044, 2015.

[10] Luc Devroye. *Nonuniform random variate generation*. Springer-Verlag, New York, 1986.

[11] Philippe Duchon. On the enumeration and generation of generalized Dyck words. *Discrete Math.*, 225(1-3):121–135, 2000. Formal power series and algebraic combinatorics (Toronto, ON, 1998).

[12] Philippe Duchon, Philippe Flajolet, Guy Louchard, and Gilles Schaeffer. Random sampling from Boltzmann principles. In Peter Widmayer et al., editor, *Automata, Languages, and Programming*, number 2380 in Lecture Notes in Computer Science, pages 501–513. Springer Verlag, 2002.

[13] Jetlir Duraj. Random walks in cones: the case of nonzero drift. *Stochastic Process. Appl.*, 124(4):1503–1518, 2014.

[14] Philippe Flajolet, Paul Zimmerman, and Bernard Van Cutsem. A calculus for the random generation of labelled combinatorial structures. *Theoretical Computer Science*, 132(1-2):1–35, 1994.

[15] Rodolphe Garbit and Kilian Raschel. On the exit time from a cone for Brownian motion with drift. *Electron. J. Probab.*, 19:no. 63, 27, 2014.

[16] Massimiliano Goldwurm. Random generation of words in an algebraic language in linear binary space. *Information Processing Letters*, 54(4):229–233, 1995.

[17] Samuel Johnson, Marni Mishna, and Karen Yeats. Towards a combinatorial understanding of lattice path asymptotics. arXiv:1305.7418, 2013.
[18] Irina Kurkova and Kilian Raschel. Random walks in $(\mathbb{Z}_+)^2$ with non-zero drift absorbed at the axes. *Bull. Soc. Math. France*, 139(3):341–387, 2011.

[19] Carine Pivoteau, Bruno Salvy, and MichÃlle Soria. Boltzmann oracle for combinatorial systems. In *Algorithms, Trees, Combinatorics and Probabilities*, pages 475–488. Discrete Mathematics and Theoretical Computer Science, 2008. Proceedings of the Fifth Colloquium on Mathematics and Computer Science. Blaubeuren, Germany. September 22-26, 2008.

[20] Kilian Raschel. Counting walks in a quadrant: a unified approach via boundary value problems. *J. Eur. Math. Soc. (JEMS)*, 14(3):749–777, 2012.

[21] Herbert S. Wilf. A unified setting for sequencing, ranking, and selection algorithms for combinatorial objects. *Advances in Mathematics*, 24:281–291, 1977.

JL: DEPARTMENT OF COMPUTER SCIENCE, PRINCETON UNIVERSITY, PRINCETON, NJ, USA; MM: DEPARTMENT OF MATHEMATICS, SIMON FRASER UNIVERSITY, CANADA; YP: CNRS – LIX AND AMIB PROJECT, ECOLE POLYTECHNIQUE AND INRIA SACLAY, FRANCE