Increasing Diamonds

Olivier Bodini, Matthieu Dien, Xavier Fontaine, Antoine Genitrini,
Hsien-Kuei Hwang

► To cite this version:

Olivier Bodini, Matthieu Dien, Xavier Fontaine, Antoine Genitrini, Hsien-Kuei Hwang. Increasing Diamonds. 12th. Latin American Theoretical INformatics Symposium, Apr 2016, Ensenada, Mexico. 10.1007/978-3-662-49529-2_16 . hal-01331227

HAL Id: hal-01331227

https://hal.sorbonne-universite.fr/hal-01331227v1

Submitted on 2 Mar 2026

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

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

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

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Increasing Diamonds *

Olivier Bodini†, Matthieu Dien‡, Xavier Fontaine†, Antoine Genitrini‡ and Hsien-Kuei Hwang§

Olivier.Bodini@lipn.univ-paris13.fr, Xavier.Fontaine@polytechnique.edu, Matthieu.Dien@lip6.fr, Antoine.Genitrini@lip6.fr, hkhwang@stat.sinica.edu.tw.

A class of diamond-shaped combinatorial structures is studied whose enumerating generating functions satisfy differential equations of the form $f'' = G(f)$, for some function $G$. In addition to their own interests and being natural extensions of increasing trees, the study of such DAG-structures was motivated by modelling executions of series-parallel concurrent processes; they may also be used in other digraph contexts having simultaneously a source and a sink, and are closely connected to a few other known combinatorial structures such as trees, cacti and permutations. We explore in this extended abstract the analytic-combinatorial aspect of these structures, as well as the algorithmic issues for efficiently generating random instances.

## 1 Introduction

Simple combinatorial structures that are both mathematically tractable and physically useful in different modeling purposes have received much attention in the literature. Typical representative examples include the simply-generated family of trees characterized by the functional equation (see [13])

$$f = zG(f),$$

and the varieties of increasing trees by the differential equation (see [3])

$$f' = G(f).$$

Due to their simplicity, these tree models also appeared naturally under various guises in many areas. Three simple prototypical cases are given in the following table.

<footnote>*This research was partially supported by the ANR MetACOnc project ANR-15-CE40-0014.</footnote>

<footnote>†Laboratoire d'Informatique de Paris-Nord, CNRS UMR 7030 - Institut Galilée - Université Paris-Nord, 99, avenue Jean-Baptiste Clément, 93430 Villetaneuse, France.</footnote>

<footnote>‡Sorbonne Universités, UPMC Univ Paris 06, CNRS, LIP6 UMR 7606, 4 place Jussieu 75005 Paris.</footnote>

<footnote>§Institute of Statistical Science, Academia Sinica, Taipei 115, Taiwan.</footnote>
<table><thead><tr><th>G</th><th>f = zG(f)</th><th>f' = G(f)</th></tr></thead><tbody><tr><td>1 + z<sup>2</sup></td><td>Binary tree</td><td>Binary increasing tree</td></tr><tr><td>(1 + z)<sup>2</sup></td><td>(Catalan tree)</td><td>(Binary search tree)</td></tr><tr><td>exp(z)</td><td>Cayley tree</td><td>Recursive tree</td></tr><tr><td>1/(1-z)</td><td>Planted (ordered) tree</td><td>Plane-oriented recursive tree</td></tr></tbody></table>

In particular, binary trees have long been studied in the computer science literature (see Knuth's book [9]) and a compilation of 214 combinatorial objects leading to the same enumerating Catalan numbers can be found in Stanley's recent book [14]. On the other hand, binary increasing trees are isomorphic to binary search trees, which represent another class of fundamental data structures with a huge number of variants; they are also closely related to Quicksort in Algorithms, to Yule-Harding models in Phylogenetics, to random permutations in Combinatorics, Rényi's car-parking problem in Applied Probability, and to Eden model in Statistical Physics, to name just a few; see [8, 6] for more information.

We explore in this paper another class of combinatorial structures, which we call *increasing diamonds*: they are labelled, directed acyclic graphs (DAGs) with a source and a sink such that the labels along any path are increasing; see Figure 1 for an illustration of two different diamonds. In standard symbolic notation (see [8]), increasing diamonds can be described as

$$ \mathcal{F} = \mathcal{Z}^{\square} + \mathcal{Z}^{\square} \star \mathcal{G}(\mathcal{F}) \star \mathcal{Z}^{\square}. \qquad (1) $$

where $\mathcal{G}$ is some functional operation specifying possible degrees and construction rules, and the two symbols $\square$ and $\star$ represent the smallest and the largest labels, respectively. This equation then translates into the differential equation satisfied by the enumerating generating function<fnref n="1" />

$$ f''(z) = G(f(z)), \text{ with } f(0) = 0 \text{ and } f'(0) = 1. \qquad (2) $$

Here $f(z) = \sum_{n \ge 1} a_n z^n / n!$, where $a_n$ enumerates the number of increasing diamonds with $n$ labels.

<figure><img src="image_2.png" /><figcaption>**Figure 1:** A binary (left) and a ternary (right) increasing diamonds of size 9 and 14, respectively.</figcaption></figure>

We study in this paper three simple representative cases, and focus on asymptotic enumeration and random generation. The following table lists the dominant term in the corresponding asymptotic approximation in each case.

<footnote n="1">¹We limit our discussion in this paper to the situation when $f'(0) = 1$ for simplicity.</footnote>
<table><thead><tr><th>G</th><th>OEIS</th><th>other description</th><th>a<sub>n</sub>/n! ~</th><th>error term of order</th></tr></thead><tbody><tr><td>exp(z)</td><td>A000111</td><td>Euler or up/down numbers</td><td>2n(2/π)<sup>n</sup></td><td>exponential</td></tr><tr><td>(1 + z)<sup>2</sup></td><td>A007558</td><td>shifts 2 places left after squared</td><td>6nρ<sub>binary</sub><sup>-n-2</sup></td><td>exponential</td></tr><tr><td>1 + z<sup>m</sup> (m ≥ 3)</td><td>-</td><td>-</td><td>C<sub>m</sub>n<sup>-m-3/(m-1)</sup>ρ<sub>m-ary</sub><sup>-n-m-2/(m-1)</sup></td><td>polynomial</td></tr><tr><td>1/(1-z)</td><td>A032035</td><td>triangular cacti with bridges</td><td>ρ<sub>plane</sub><sup>1-n</sup>n<sup>2</sup>√2logn</td><td>logarithmic</td></tr><tr><td>1/(1-z)<sup>3</sup></td><td>A001147</td><td>double factorial</td><td>2<sup>n</sup>√πn</td><td>polynomial</td></tr></tbody></table>

Here OEIS stands for Sloane's Online Encyclopedia of Integer Sequences, $C_m$ is a constant (see Theorem 5), and $\rho_{\text{binary}}$, $\rho_{\text{m-ary}}$ and $\rho_{\text{plane}}$ are three constants given in (9), (12) and (14), respectively. While most properties are expected to be similar to those of increasing trees (see [3]), the higher order derivative introduces more technical difficulties, as visible from the less common asymptotic order produced when $G = 1/(1-z)$.

Structurally, increasing diamonds are bipolar digraphs with a downward planarity; they are also special cases of series-parallel graphs and are more expressive than quasi-trees in [2]. Since DAGs with a unique source and a unique sink appear naturally in many concrete applications, our increasing diamonds may be of potential use in modelling structural parameters or problem complexity in these contexts. Typical examples include: partial orders and their linear extensions, computational processes and their executions in parallel computing, network or data flows, food-webs, register sharing, machine learning, streaming analysis, grid computing, etc.

To be useful for modelling concrete structures in applications, we need either more precise statistical properties or more efficient generation algorithms for random increasing diamonds. The former will be addressed elsewhere, and for the latter, we will focus on the by now popular Boltzmann sampling algorithm proposed in [7], which was recently extended in [5] to deal with the situation of first-order differential equations. We develop further techniques to handle the second-order differential equations.

On the other hand, the type of differential equations we study in this paper ($f'' = G(f)$) also emerges naturally in other contexts, notably in a recent paper by Kuba and Panholzer [10] on multi-labeled increasing trees and hook-length formulae; see also their earlier paper [12]. While the equations are the same, our combinatorial structures here are different and to some extent more natural, and such a difference is reflected by the initial conditions: we focus on $f(0) = 0$ and $f'(0) = 1$ whereas they deal with $f(0) = f'(0) = 0$. Also we will derive asymptotic expansions. Along the same direction, Kuba and Panholzer examined in [11] another class of tree structures whose exponential generating function satisfies $f^{(m)} = (m-1)!e^{mf}$, where $m \ge 2$, which coincides with our model when $m = 2$. They studied in detail some shape characteristics in such random trees. Further connections can be made between such bucket trees and our increasing diamonds.

The paper is organized as follows. In the next section, we analyze the three classes of increasing diamonds in detail. Then Section 3 is devoted to the development of algorithmic tools for generating efficiently random diamonds that rely on the notion of uniform Boltzmann sampling.

## 2 Exact enumeration and asymptotics

In this section, we first discuss the general solution of the differential equation $f'' = G(f)$ subject to the initial conditions $f(0) = 0$ and $f'(0) = 1$ (other initial conditions can be dealt with in a similar
manner), and then concentrate our discussion on a few special cases for which we will derive more precise asymptotic approximations.

## 2.1 General solution of $f'' = G(f)$.

Multiplying both sides of (2) by $2f'$, we obtain $2f''f' = 2G(f)f'$, which implies that

$$f'(z)^2 = f'(0)^2 + \int_0^z 2G(f(t))f'(t)dt = 1 + 2\mathcal{G}(f(z)),$$

where $\mathcal{G}(z) := \int_0^z G(t)dt$. Thus $f'(z) = \pm\sqrt{1 + 2\mathcal{G}(f(z))}$, or

$$\pm \int_0^{f(z)} \frac{1}{\sqrt{1 + 2\mathcal{G}(t)}} dt = z. \qquad (3)$$

**Lemma 1.** The solution to the differential equation $f'' = G(f)$ with $f(0) = 0$ and $f'(0) = 1$ is given by

$$\int_0^{f(z)} \frac{1}{\sqrt{1 + 2\mathcal{G}(t)}} dt = z. \qquad (4)$$

*Proof.* By expanding the left-hand side of the second equation in (3) as a Taylor series, we exclude the negative solution and conclude (4). □

The solution (4) is, although implicit, useful in our asymptotic analysis even when no further simplification is possible. First recall a useful property when $f$ blows up near the dominant singularity, which is readily modified from Lemma 1 of [3].

**Lemma 2.** Given an entire function $G$, the dominant real positive singularity of the function $f(z)$, solution to $Y'' = G(Y)$ with $Y(0) = 0$ and $Y'(0) = 1$, is given by

$$\rho = \int_{0}^{\infty} \frac{dt}{\sqrt{1 + 2 \int_{0}^{t} G(v)dv}},$$

provided that the integral converges.

From the brief discussion in Introduction, we see that the coefficient $a_n$ of $f$ is well-approximated by $a_n \sim Cn!\rho^{-n}n^{\alpha}(\log n)^{\beta}$, and from this observation we expect that the singularity analysis of Flajolet and Odlyzko (see [8]) will be useful in such an analysis, as in [3].

## 2.2 Non-plane (unordered) increasing diamonds

We discuss in detail the class of non-plane increasing diamonds, which can be decomposed as sets of increasing diamonds:

$$\mathcal{F} = \mathcal{Z}^{\square} + \mathcal{Z}^{\square} \star \text{SET}(\mathcal{F}) \star \mathcal{Z}^{\square},$$

so that the corresponding exponential generating function satisfies $f'' = e^f$. The two diamonds in Figure 1 may be regarded, neglecting the order of subtrees, as two instances of non-plane increasing diamonds.
By (4) with $G(z) = e^z$ and $\mathcal{G}(z) = e^z - 1$, we see that the exponential generating function $f$ of $a_n$ has the solution

$$f(z) = -\log(1 - \sin z), \qquad (5)$$

and the number $a_n$ of such increasing diamonds with $n$ labels starts with

$$\{a_n\}_{n \ge 1} = \{1, 1, 1, 2, 5, 16, 61, 272, 1385, 7936, 50521, 353792, \dots\},$$

which coincides with A000111 in Sloane's OEIS, where many other structures with identical enumerating sequence are also given (alternating permutations, zig-zag posets, some increasing trees, etc.). This shows the richness and usefulness of the equation $f'' = e^f$ in combinatorial objects.

Note that, by the differential equation $f''' = f'f''$ (obtained by the differentiation of $f'' = e^f$), we have the recurrence relation

$$a_n = \sum_{2 \le k < n} \binom{n-3}{k-2} a_k a_{n-k} \quad (n \ge 3),$$

which is useful for numerical purposes.

**Theorem 3.** The number $a_n$ of non-plane increasing diamonds with $n$ labels satisfies

$$a_n = \frac{2^{n+1} (n-1)!}{\pi^n} \sum_{j=-\infty}^{+\infty} \frac{1}{(1+4j)^n}. \qquad (6)$$

It is less obvious that the right-hand side represents an integer.

*Proof.* By (5), we have $f'(z) = \tan z + \sec z$, which has only simple poles at $z = (2k + \frac{1}{2})\pi$. By standard expansion for meromorphic functions ([8, Ch. IV]), we obtain the expansion (6), which is not only an asymptotic expansion (expressible as Hurwitz's zeta function)

$$a_n = \frac{2^{n+1} (n-1)!}{\pi^n} \sum_{j \ge 0} \left( \frac{1}{(1+4j)^n} + \frac{(-1)^n}{(4j+3)^n} \right),$$

but also an identity for $n \ge 1$. $\square$

Another exactly solvable case is when $G(z) = (1 - z)^{-3}$. In this case, we have the surprisingly simple solution (cf. [10])

$$f(z) = 1 - \sqrt{1 - 2z}, \qquad (7)$$

leading to the simple expression for the total number of size-$n$ diamonds

$$a_n = (2n - 3)!! = \frac{(2n - 2)!}{2^{n-1}(n-1)!} \quad (n \ge 1).$$

However, exact solutions as (7) and (5) are exceptional rather than commonplace, and different techniques are needed in most cases as we will see below.
## 2.3 m-ary increasing diamonds

Consider now increasing diamonds in which the degrees of nodes are limited to $m \ge 2$; see Figure 1 for a binary and a ternary diamond. In this case, we have the specification

$$ \mathcal{F} = \mathcal{Z}^{\square} + \mathcal{Z}^{\square} \star \mathcal{F}^{m} \star \mathcal{Z}^{\square}, \qquad (8) $$

which leads to the differential equation $f'' = 1 + f^m$ with $f(0) = 0$ and $f'(0) = 1$. Closed-form solutions are possible when $m = 2$ and $m = 3$ (in terms of elliptic integrals), but they are not simple. So we present only the solution for $m = 2$ and derive asymptotic approximation for $m \ge 3$ (in a slightly more general formulation).

**Binary increasing diamonds and Weierstrass's $\wp$-function.** From (4), we see that $f$ satisfies the equation

$$ \int_{0}^{f(z)} \frac{1}{\sqrt{1 + 2t + \frac{2}{m+1}t^{m+1}}} dt = z. $$

When $m = 2$, we can express the solution in terms of Weierstrass's elliptic function $\wp$ (see [1]), which is defined periodically over a lattice that contains one double pole in a corner of each cell. Thus, by construction,

$$ \wp(z; \omega_1, \omega_2) = \frac{1}{z^2} + \sum_{(k,l) \in \mathbb{Z}^2 \setminus \{(0,0)\}} \left( \frac{1}{(z + k\omega_1 + l\omega_2)^2} - \frac{1}{(k\omega_1 + l\omega_2)^2} \right), $$

where $\omega_1$ and $\omega_2$ are the periods of $\wp$.

**Theorem 4.** The exponential generating function of the number of binary increasing diamonds can be expressed as

$$ f(z) = 6\wp(z - \rho; -\frac{1}{3}, -\frac{1}{36}) \quad \text{where } \rho := \int_{0}^{\infty} \frac{dt}{\sqrt{1 + 2t + \frac{2}{3}t^3}}, \qquad (9) $$

and the number of size-n binary increasing diamonds is given by

$$ a_n = 6 \frac{(n+1)!}{\rho^{n+2}} \sum_{(k,l) \in \mathbb{Z}^2} \frac{1}{\left(1 + \frac{k\omega_1}{\rho} + \frac{l\omega_2}{\rho}\right)^{n+2}}, \qquad (10) $$

where $\omega_1$ and $\omega_2$ are computable constants.

Asymptotically, $a_n \sim 6(n+1)!\rho^{-n-2}$, with an exponentially small error. Note that, by starting with the initial conditions $f(0) = f'(0) = 0$, we then obtain the bi-labelled increasing trees defined in [10], which corresponds to the sequence A144849 in OEIS.

*Proof. (Sketch)* The $\wp$-function satisfies the differential equation

$$ \wp'^2(z) = 4\wp^3(z) - g_2\wp(z) - g_3, $$

and we need only to identify the corresponding parameters.

By Lemma 2, we first determine the dominant singularity $\rho$; then from the series expansion of $\wp$, we deduce (10) by a direct application of singularity analysis (see [8]). □

Although few cases lead to closed-form expressions in terms of known functions, it is not difficult to derive asymptotic approximations based on complex analysis and singularity analysis, as already highlighted in the classical paper [3].
**Polynomial varieties of increasing diamonds.** As in [3], we consider the polynomial varieties of increasing diamonds, which are characterized by $G(z)$ being a polynomial, say

$$G(z) = \sum_{0 \le j \le m} b_j z^j, \qquad (11)$$

where $m \ge 2$ and $b_m > 0$. For simplicity, we may assume that $G(z) \ne z^l H(z^k)$, for some $k \ge 2$ and $l \ge 0$, namely, $G$ is aperiodic.

Then by (2), the dominant singularity is given by

$$\rho = \int_{0}^{\infty} \frac{1}{\sqrt{1 + 2 \sum_{0 \le j \le m} \frac{b_j}{j+1} t^{j+1}}} dt, \qquad (12)$$

which is absolutely convergent since $m \ge 2$. Then we apply the same idea used in [3], and obtain

$$\rho - z = \int_{f(z)}^{\infty} \frac{1}{\sqrt{1 + 2 \sum_{0 \le j \le m} \frac{b_j}{j+1} t^{j+1}}} dt \\ = \frac{\sqrt{m+1}}{(m-1)\sqrt{2b_m}} f^{-\frac{m-1}{2}} - \frac{b_{m-1}\sqrt{m+1}}{mb_m^{3/2}} f^{-\frac{m+1}{2}} + \dots,$$

as $z \to \infty$. Then by inverting, we get

$$f(z) = \left( \frac{(m-1)\sqrt{b_m}}{\sqrt{2(m+1)}} \right)^{-\frac{2}{m-1}} (\rho - z)^{-\frac{2}{m-1}} \left( 1 + O\left( |\rho - z|^{\frac{2}{m-1}} \right) \right),$$

as $z \sim \rho$, the justification following also standard line. We then deduce by the singularity analysis the following asymptotic approximation.

**Theorem 5.** Assume that $G$ is a polynomial given in (11) and $a_n > 0$ for $n \ge n_0$ for some $n_0 > 0$. Then the number of increasing diamonds with $n$ labels in a polynomial variety satisfies

$$a_n = \left( \frac{\sqrt{2(m+1)}}{(m-1)\sqrt{b_m}} \right)^{\frac{2}{m-1}} \frac{n^{-\frac{m-3}{m-1}}}{\Gamma(\frac{2}{m-1})} \rho^{-n-\frac{2}{m-1}} \left( 1 + O\left( n^{-\frac{4}{m-1}} \right) \right),$$

for $m \ge 2$, where $\rho$ is given in (12).

Note that the asymptotic estimates here are independent of the initial conditions.

In the special case when $m = 3$, it is possible to express $f$ in terms of Jacobi elliptic functions, but the expression is messy.

## 2.4 Plane increasing diamonds

We now focus on plane (ordered) increasing diamonds, which are described by

$$\mathcal{F} = \mathcal{Z}^{\square} + \mathcal{Z}^{\square} \star \mathrm{SEQ}(\mathcal{F}) \star \mathcal{Z}^{\square}, \qquad (13)$$

leading to the differential equation $f'' = \frac{1}{1-f}$ with the initial conditions $f(0) = 0$ and $f'(0) = 1$.

The analysis of such diamonds is more involved and the asymptotic expansion we obtain has a much poorer convergence rate: instead of exponential or polynomial, the terms are now in decreasing powers of $\log n$.
**Theorem 6.** The number of plane increasing diamonds with n labels satisfies

$$a_n = \frac{n! \rho^{1-n}}{n^2 \sqrt{2 \log n}} \left( \sum_{0 \le k < K} \frac{P_k(\log \log n)}{(\log n)^k} + \mathcal{O} \left( \frac{(\log \log n)^K}{(\log n)^K} \right) \right),$$

where

$$\rho := \int_{0}^{\infty} \frac{1}{\sqrt{1 - 2 \log(1 - t)}} dt = \frac{\sqrt{e}}{2} \int_{0}^{\infty} v^{-\frac{1}{2}} e^{-v} dv \approx 0.6556795424 \dots, \quad (14)$$

and the $P_k$'s are computable polynomials (of degree k).

In particular, $P_0(x) = 1$ and $P_1(x) = \frac{1}{8}(x - 3 - 2\gamma + \log 2 + 2\log \rho)$.

The method of proof is the same as above, details being omitted here. The first few terms of $a_n$ are

$$\{1, 1, 1, 3, 13, 77, 573, 5143, 54025, 650121, 8817001, 133049339, \dots\},$$

and corresponds to sequence A032035 in OEIS, which also enumerates increasing rooted (2,3)-cacti with $n-1$ nodes. Note that $f_1 = f' - 1$ satisfies the differential equation $f_1' = e^{f_1+f_1^2/2}$.

# 3 Random generation via Boltzmann samplers

## 3.1 Boltzmann Samplers for the differential classes

The Boltzmann sampling technique was first proposed in the seminal paper [7], and has been widely developed and extended since then. It captures the features any successful algorithm must have: simple, efficient and easily extensible.

In this subsection, we briefly recall this technique for labeled structures.

**Definition 7.** A Boltzmann sampler of parameter $x > 0$ is an algorithm that draws an object $\alpha$ of size $|\alpha|$ in a given combinatorial class $\mathcal{A}$ with the probability $P_x(\alpha) = \frac{x^{|\alpha|}}{|\alpha|!A(x)}$.

<figure><img src="image_3.png" /><figcaption>**Figure 2:** A random diamond of size 591 satisfying $f'' = 1 + f^3$, with $f(0) = 0$ et $f'(0) = 1$.</figcaption></figure>

Note that the output size $N$ of a Boltzmann sampler is a random variable with the law $P_x(N = n) = a_n x^n / (n!A(x))$, and the expectation of $N$ is $E_x(N) = xA'(x)/A(x)$. Here $x$ is a free variable. To generate an object of size $n$, one can choose the parameter $x$ to be the solution of the saddle point equation $E_x(N) = n$. With this choice, it is possible to devise a linear-time algorithm to generate a random instance by repeated use of trial-and-rejection until reaching an output of size in $[(1 - \epsilon)n, (1 + \epsilon)n]$ (referred to as an approximate-size algorithm).

This universal method is not only very efficient but also fully automatizable. What we need is a complete symbolic (recursive or not) description of the class in order to construct a sampler. Indeed, Boltzmann samplers for the neutral and atomic classes $\mathcal{E}$ and $\mathcal{Z}$ are trivial, and from there general procedures exist for constructing more complex samplers through elementary operations such as addition, multiplication, cycle, set, etc. We refer the reader to the original paper [7] for more details. On the other hand, the Boltzmann sampler for the box-operator of two classes was addressed in [5].

Note that Boltzmann samplers does not return a labeled object, but only the unlabeled skeleton. To complete the process, a labeling algorithm is needed.
## 3.2 Boltzmann samplers for second-order differential classes

It is natural to divide the problem into two cases, one in which the differential equation is induced by the general shape specification $\mathcal{F}'' = \phi(\mathcal{Z}, \mathcal{F})$, where $\mathcal{F}''$ denotes the class of objects of $\mathcal{F}$ in which two nodes are pointed, and the other by $\mathcal{F}'' = \phi(\mathcal{F})$.

Before considering these two issues, we recall some basic and classical properties. First, the box product and the derivative operator are linked together by the fact that $C = A^\top * B$ entails that $C' = A' \times B$. Secondly, we know how to get a sampler of parameter $x$ for $\mathcal{F}$ by just using a sampler of $\mathcal{F}'$. This surprising result is obtained by multiplying the Boltzmann parameter $x$ by a suitable continuous random variable $u$ in $[0, 1]$. Indeed, this yields the following algorithm described in [4], which can also be derived from results in [5].

<caption>**Algorithm 1:** $\Gamma_x \mathcal{F}$ from $\Gamma \mathcal{F}'$</caption>

1: **if** Bernoulli($f(0)/f(x)$) **then**
2:     **return** an object of size 0
3: **else**
4:     Draw $U \in [0, 1]$ with the density $\delta_x(u) = f'(ux)x/(f(x) - f(0)) \cdot \mathbf{1}_{[0,1]}(u)$
5:     Draw $\gamma' = \Gamma_{U_x \mathcal{F}}$
6:     **return** $\gamma'$ where the bud is replaced by an atom.
7: **end if**

In line 5, the object contains what is called a bud in Species Theory. It can be seen as a hole, that is the reason why it is replaced by an atom (in line 6).

**General Case** $\mathcal{F}'' = \phi(\mathcal{Z}, \mathcal{F})$. We consider now the case $\mathcal{F}'' = \phi(\mathcal{Z}, \mathcal{F})$, which can be dealt with by applying twice Algorithm 1. But this requires to draw two continuous random variables $U$ and $V$, and use only their product $UV$. Clearly, this can be factored by calculating directly the random variable $S = UV$. This gives the following algorithm for which the proof is similar to that of $\mathcal{F}'' = \phi(\mathcal{Z}, \mathcal{F})$ in [5].

<caption>**Algorithm 2:** $\Gamma_x \mathcal{F}$ generates an object in $\mathcal{F}$ from a sampling in $\mathcal{F}''$</caption>

1: Draw $W \in [0, 1]$ uniformly
2: **if** $W < \frac{f(0)}{f(x)}$ **then**
3:     **return** an object of size 0
4: **else if** $W < \frac{f(0) + xf'(0)}{f(x)}$ **then**
5:     **return** an object of size 1
6: **else**
7:     Draw $S \in [0, 1]$ according to the density $\delta_x(s) = \frac{x^2(1-s)f''(sx)}{f(x) - xf'(0) - f(0)} \mathbf{1}_{[0,1]}(s)$
8:     Draw $\gamma''$ using $\Gamma_{S_x \mathcal{F}''}$
9:     **return** replace the buds in $\gamma''$ by two atoms.
10: **end if**

**Particular Case** $\mathcal{F}'' = \phi(\mathcal{F})$. We consider here the special case where $\phi$ does not explicitly depend on $\mathcal{Z}$. The Algorithm from [4] can be amended to deal with uniform continuous random variables rather than non-uniform random variables that are hard to simulate.

Classical Boltzmann samplers $\Gamma$ are parametrized by $x$, so the sampler draws an object $\alpha$ in $\mathcal{A}$ with probability $P_x(\alpha) = x^{|\alpha|}/(|\alpha|!A(x))$. But in the case of functional equations where $x$ is not explicit (such as $\mathcal{F}' = \phi(\mathcal{F})$), it has been observed in [4] that it is preferable to deal with another parameter $\tau = f(x)$. In this case, the output is distributed as $P(N = n) = a_n f^{-1}(\tau)^n/(n!\tau)$. It is nevertheless
always a Boltzmann sampler but with a different parametrization. To avoid confusion, we then indicate $\Gamma_{[\tau]}\mathcal{F}$ instead of $\Gamma_x\mathcal{F}$. Thus we can now give an algorithm similar to Algorithm 1 that uses only uniform random variables.

<caption>**Algorithm 3:** $\Gamma_{[\tau]}\mathcal{F}$ generates an object in $\mathcal{F}$ from a sampling in $\mathcal{F}'$</caption>

1: **if** Bernoulli($f(0)/\tau$) **then**
2:     **return** an object of size 0
3:     **else**
4:         Draw $U$ uniformly $\in [0, 1]$
5:         $\tau_{new} \leftarrow U\tau + (1 - U)f(0)$
6:         Draw $\gamma'$ using $\Gamma_{[\tau_{new}]\mathcal{F}}'$
7:         **return** $\gamma'$ where we replace the bud by an atom.
8:     **end if**

In order to apply twice this procedure (because $\mathcal{F}'' = (\mathcal{F}')'$), we need to obtain $\Gamma_{[\tau_{new}]\mathcal{F}}'$ by using the Boltzmann sampler of $\mathcal{F}'$. For this, let $y = f^{-1}(\tau_{new})$. We have $f'(vy) = Vf'(y) + (1-V)f'(0)$ where $V$ is a uniform random variable on $[0, 1]$. Since we are looking for an algorithm $\Gamma_{[\tau_{new}]\mathcal{F}}'$ where $\tau_{new} = f(y)$, we need an expression of $f(vy)$ in function of $\tau_{new}$. But since the differential equation $f''(z) = \phi(f(z))$ can be integrated (by multiplying both sides by $f'(z)$), we then get $f'(z) = g(f(z))$, where $\frac{1}{2}g^2$ is the primitive of $\phi$ such that $f(0) = \frac{f'(0)^2}{2}$. Then we get the expression $f(vy) = g^{-1}(Vg(f(y)) + (1-V)f'(0))$. Finally, we obtain the following algorithm.

<caption>**Algorithm 4:** $\Gamma_{[\tau_0]}\mathcal{F}$ generates an object of $\mathcal{F}$ following the Boltzmann distribution of parameter $x = f^{-1}(\tau_0)$, from a sampler of $\mathcal{F}'' = \Phi(\mathcal{F})$</caption>

1: **if** Bernoulli($f(0)/\tau_0$) **then**
2:     **return** an object of size 0
3:     **else**
4:         Draw $U$ uniformly on $[0, 1]$
5:         $\sigma \leftarrow U\tau_0 + (1 - U)f(0)$
6:         **if** Bernoulli($f'(0)/g(\sigma)$) **then**
7:             **return** an object of size 1
8:         **else**
9:             Draw $V$ uniformly on $[0, 1]$
10:         $\tau \leftarrow g^{-1}(Vg(\sigma) + (1 - V)f'(0))$
11:         Draw $\gamma''$ using $\Gamma_{[\tau]}\mathcal{F}'' = \Gamma_{[\tau]}\Phi(\mathcal{F})$
12:         **return** $\gamma''$ where the buds are replaced by two atoms.
13:     **end if**
14: **end if**

In contrast to the previous algorithm, we do not need here to draw random variables with complicated laws. This very simple sampler is easily implemented for testing purposes. It remains to analyze its complexity. As already discussed above, the dominant singularity $\rho$ of $f$ is of the form $(1 - z/\rho)^{-\alpha}$ for some $\alpha > 0$. This ensures the following theorem.

**Theorem 8.** *Algorithm 4 provides a Boltzmann sampler, and its approximate-size version gives a linear time algorithm for drawing uniformly at random a diamond of type $\mathcal{F}'' = \Phi(\mathcal{F})$, where $\Phi$ is a polynomial.*

We implemented this algorithm in Java, and obtained the following table, which synthesizes benchmarks computed on a laptop (1.5 GHz CPU and 4G RAM). The examples we tested consist of ternary diamonds $f'' = 1 + f^3$ with initial conditions $f(0) = 0$ and $f'(0) = 1$, and with size tolerance set at 10 percent. One of such diamonds is depicted in Figure 2. We observe that the timing results are consistent with our analysis.
<table><tr><td>Size n</td><td>10</td><td>100</td><td>1000</td><td>5000</td><td>10000</td><td>50000</td><td>100000</td><td>150000</td></tr><tr><td>τ<sub>0</sub></td><td>8.73</td><td>80.44</td><td>794</td><td>3972</td><td>7941</td><td>39752</td><td>79559</td><td>119086</td></tr><tr><td>Time (ms)</td><td>1</td><td>7</td><td>66</td><td>322</td><td>668</td><td>3887</td><td>7098</td><td>9812</td></tr></table>

## References

[1] M. Abramowitz and I.A. Stegun. *Handbook of Mathematical Functions: with Formulas, Graphs, and Mathematical Tables*. Dover Publications, 2012.

[2] Ei Ando, Toshio Nakata, and Masafumi Yamashita. Approximating the longest path length of a stochastic DAG by a normal distribution in linear time. *J. Discrete Algorithms*, 7(4):420–438, 2009.

[3] F. Bergeron, P. Flajolet, and B. Salvy. Varieties of Increasing Trees. In *CAAP*, volume 581 of LNCS, pages 24–48. Springer, 1992.

[4] O. Bodini. *Autour de la génération aléatoire sous modèle de Boltzmann*. Habilitation thesis, UPMC, 2010.

[5] O. Bodini, O. Roussel, and M. Soria. Boltzmann samplers for first-order differential specifications. *Discrete Applied Mathematics*, 160(18):2563–2572, 2012.

[6] Hua-Huai Chern, María-Inés Fernández-Camacho, Hsien-Kuei Hwang, and Conrado Martínez. Psi-series method for equality of random trees and quadratic convolution recurrences. *Random Structures Algorithms*, 44(1):67–108, 2014.

[7] P. Duchon, P. Flajolet, G. Louchard, and G. Schaeffer. Boltzmann samplers for the random generation of combinatorial structures. *Combinatorics, Probability & Computing*, 13(4-5):577–625, 2004.

[8] P. Flajolet and R. Sedgewick. *Analytic Combinatorics*. Camb. Univ. Press, 2009.

[9] Donald E. Knuth. *The Art of Computer Programming. Vol. 1*. Addison-Wesley, Reading, MA, 1997. Fundamental algorithms, Third edition [of MR0286317].

[10] M. Kuba and A. Panholzer. Combinatorial families of multilabelled increasing trees and hook-length formulas. *Discrete Mathematics*, 339:227–254, 2016.

[11] Markus Kuba and Alois Panholzer. A combinatorial approach to the analysis of bucket recursive trees. *Theor. Comput. Sci.*, 411(34-36):3255–3273, 2010.

[12] Markus Kuba and Alois Panholzer. Bilabelled increasing trees and hook-length formulae. *European J. Combin.*, 33(2):248–258, 2012.

[13] A. Meir and J. W. Moon. On the altitude of nodes in random trees. *Canad. J. ath.*, 30(5):997–1015, 1978.

[14] R. Stanley. *Catalan Numbers*. Cambridge University Press, 2015.