$\Biggl{|}\mbox{$\>\\!$}\mathrm{e}^{\mathrm{i}t}-\sum_{j=0}^{k-1}\frac{(\mathrm{i}\mbox{$\>\\!$}t)^{j}}{j!}\Biggr{|}\leq\frac{\,|t|^{k}}{k!\,}.$
(5.11)

Now, combining the uniform expansions (5.10) and (5.11) (the latter with $k=2$
or $k=1$ as appropriate) and returning to (5.7), we obtain

$\displaystyle\sum_{\ell\in\mathbb{N}^{q}}\log\left(1+w_{\ell}\right)$
$\displaystyle=\sum_{\ell\in\mathbb{N}^{q}}\frac{z_{1}^{\ell}z_{2}}{1+z_{1}^{\ell}z_{2}}\left(\mathrm{i}\mbox{$\;\\!$}(\tilde{t}_{1}\ell+\tilde{t}_{2})-\frac{1}{2}\left(\tilde{t}_{1}\ell+\tilde{t}_{2}\right)^{2}+O(1)\left(\tilde{t}_{1}\ell+\tilde{t}_{2}\right)^{3}\right)$
$\displaystyle\quad+O(1)\sum_{\ell\in\mathbb{N}^{q}}\frac{z_{1}^{2\ell}z_{2}^{2}}{(1+z_{1}^{\ell}z_{2})^{2}}\left(\tilde{t}_{1}\ell+\tilde{t}_{2}\right)^{2}$
$\displaystyle=\mathrm{i}\mbox{$\;\\!$}\tilde{t}_{1}\sum_{\ell\in\mathbb{N}^{q}}\frac{\ell\mbox{$\;\\!$}z_{1}^{\ell}z_{2}}{1+z_{1}^{\ell}z_{2}}+\mathrm{i}\mbox{$\;\\!$}\tilde{t}_{2}\sum_{\ell\in\mathbb{N}^{q}}\frac{z_{1}^{\ell}z_{2}}{1+z_{1}^{\ell}z_{2}}-\frac{1}{2}\sum_{\ell\in\mathbb{N}^{q}}\frac{z_{1}^{\ell}z_{2}}{1+z_{1}^{\ell}z_{2}}\left(\tilde{t}_{1}\ell+\tilde{t}_{2}\right)^{2}$
$\displaystyle\quad+O(1)\sum_{\ell\in\mathbb{N}^{q}}z_{1}^{\ell}z_{2}\left(\tilde{t}_{1}\ell+\tilde{t}_{2}\right)^{3}+O(1)\sum_{\ell\in\mathbb{N}^{q}}z_{1}^{2\ell}z_{2}^{2}\left(\tilde{t}_{1}\ell+\tilde{t}_{2}\right)^{2}$
$\displaystyle=:\mathrm{i}\mbox{$\>\\!$}\tilde{t}_{1}\mbox{$\>\\!$}\varSigma_{1}+\mathrm{i}\mbox{$\>\\!$}\tilde{t}_{2}\mbox{$\;\\!$}\varSigma_{2}-\tfrac{1}{2}\mbox{$\>\\!$}\varSigma_{3}+O(1)\mbox{$\>\\!$}\varSigma_{4}+O(1)\mbox{$\>\\!$}\varSigma_{5}.$
(5.12)

According to the calibration equations (see (3.27) and (3.30)), the first two
sums in (5.12) are known exactly,

$\varSigma_{1}=\langle N\rangle,\qquad\varSigma_{2}=\langle M\rangle.$ (5.13)

Next, the error sums $\varSigma_{4}$ and $\varSigma_{5}$ can be shown to
asymptotically vanish. Indeed, using the elementary inequality $(a+b)^{r}\leq
2^{r-1}\left(a^{r}+b^{r}\right)$ ($r\geq 1$) and combining Lemma 3.2 with
formulas (3.25), (3.26) and (3.35) gives upon simple calculations the estimate

$\displaystyle 0\leq\varSigma_{4}$ $\displaystyle\leq
4\mbox{$\;\\!$}\tilde{t}_{1}^{\mbox{$\;\\!$}3}\sum_{\ell\in\mathbb{N}^{q}}\ell^{3}z_{1}^{\ell}z_{2}+4\mbox{$\;\\!$}\tilde{t}_{2}^{\mbox{$\;\\!$}3}\sum_{\ell\in\mathbb{N}^{q}}z_{1}^{\ell}z_{2}$
$\displaystyle=\frac{O(1)\mbox{$\;\\!$}\langle
M\rangle^{3/2}\mbox{$\>\\!$}z_{2}}{\langle
N\rangle^{3}\mbox{$\;\\!$}\gamma^{3+1/q}}+\frac{O(1)\mbox{$\;\\!$}z_{2}}{\langle
M\rangle^{3/2}\mbox{$\;\\!$}\gamma^{1/q}}=\frac{O(1)}{\langle
M\rangle^{1/2}}=o(1).$ (5.14)

Similarly,

$\displaystyle 0\leq\varSigma_{5}$ $\displaystyle\leq
2\mbox{$\;\\!$}\tilde{t}_{1}^{\mbox{$\;\\!$}2}\sum_{\ell\in\mathbb{N}^{q}}\ell^{2}z_{1}^{2\ell}z_{2}^{2}+2\mbox{$\;\\!$}\tilde{t}_{2}^{\mbox{$\;\\!$}2}\sum_{\ell\in\mathbb{N}^{q}}z_{1}^{2\ell}z_{2}^{2}$
$\displaystyle=\frac{O(1)\mbox{$\;\\!$}\langle
M\rangle\mbox{$\>\\!$}z_{2}^{2}}{\langle
N\rangle^{2}\mbox{$\;\\!$}\gamma^{2+1/q}}+\frac{O(1)\mbox{$\;\\!$}z_{2}^{2}}{\langle
M\rangle\mbox{$\;\\!$}\gamma^{1/q}}=O(\kappa^{1/q})=o(1).$ (5.15)

Finally, consider the sum $\varSigma_{3}$ in (5.12) which, as we will see,
provides the main contribution to (5.7). To this end, observe (cf. (3.27))
that

$\displaystyle
0\leq\sum_{\ell\in\mathbb{N}^{q}}z_{1}^{\ell}z_{2}\left(\tilde{t}_{1}\ell+\tilde{t}_{2}\right)^{2}-\varSigma_{3}$
$\displaystyle=\sum_{\ell\in\mathbb{N}^{q}}\frac{z_{1}^{2\ell}z_{2}^{2}}{1+z_{1}^{\ell}z_{2}}\left(\tilde{t}_{1}\ell+\tilde{t}_{2}\right)^{2}$
$\displaystyle\leq\sum_{\ell\in\mathbb{N}^{q}}z_{1}^{2\ell}z_{2}^{2}\left(\tilde{t}_{1}\ell+\tilde{t}_{2}\right)^{2}=\varSigma_{5}=o(1),$
(5.16)

as shown in (5.15). In turn, using the asymptotic results (3.43), (3.44) and
(3.45), and recalling the rescaling expressions (5.5), we obtain the limit

$\displaystyle\sum_{\ell\in\mathbb{N}^{q}}z_{1}^{\ell}z_{2}\left(\tilde{t}_{1}\ell+\tilde{t}_{2}\right)^{2}$
$\displaystyle=\tilde{t}_{1}^{\mbox{$\;\\!$}2}\sum_{\ell\in\mathbb{N}^{q}}\ell^{2}z_{1}^{\ell}z_{2}+2\mbox{$\;\\!$}\tilde{t}_{1}\tilde{t}_{2}\sum_{\ell\in\mathbb{N}^{q}}\ell\mbox{$\>\\!$}z_{1}^{\ell}z_{2}+\tilde{t}_{2}^{\mbox{$\;\\!$}2}\sum_{\ell\in\mathbb{N}^{q}}z_{1}^{\ell}z_{2}$
$\displaystyle\to
t_{1}^{2}+\frac{2\mbox{$\;\\!$}t_{1}t_{2}}{\sqrt{q+1}}+t_{2}^{2}\mbox{$\;\\!$},$
(5.17)

which is a quadratic form with matrix (5.2).

Thus, substituting the estimates (5.13), (5.14), (5.15), (5.16) and (5.17)
into (5.12) and returning to (5.7), we obtain

$\varphi(t_{1},t_{2})\to\exp\left\\{-\frac{1}{2}\left(t_{1}^{2}+\frac{2\mbox{$\;\\!$}t_{1}t_{2}}{\sqrt{q+1}}+t_{2}^{2}\right)\right\\}\\!,$

and the proof of Theorem 5.1 is complete. ∎

######  Corollary 5.2.

Under the hypotheses of Theorem 5.1, the following laws of large numbers hold
under the Boltzmann measure ${\mathsf{P}\mbox{$\>\\!\\!$}}_{\bm{z}}$,

$\frac{M_{\lambda}}{\langle
M\rangle}\stackrel{{\scriptstyle\mathrm{p}}}{{\longrightarrow}}1,\qquad\frac{N_{\lambda}}{\langle
N\rangle}\stackrel{{\scriptstyle\mathrm{p}}}{{\longrightarrow}}1.$ (5.18)

######  Remark 5.1.

As a curiosity, we observe that the limiting distribution in Theorem 4.1
formally conforms to Theorem 5.1 under the additional limit as $\langle
M\rangle\to\infty$. Indeed, start with the intermediate limit (4.11) (with
$\langle M\rangle$ fixed) and switch from Laplace transform to characteristic
function by formally changing $(s_{1},s_{2})$ ($s_{i}\geq 0$) to
$-\mathrm{i}\mbox{$\;\\!$}(t_{1},t_{2})$ ($t_{i}\in\mathbb{R}$). Then, bearing
in mind the normalization (5.1), we obtain

$\displaystyle-\langle
M\rangle\\!\left(1-\mathrm{e}^{\mbox{$\>\\!$}\mathrm{i}\mbox{$\>\\!$}t_{2}/\sqrt{\langle
M\rangle}}\left(1-\frac{\mathrm{i}\mbox{$\;\\!$}t_{1}q}{\sqrt{\langle
M\rangle}}\right)^{-1/q}\right)+\frac{\mathrm{i}\mbox{$\;\\!$}t_{1}\sqrt{\langle
M\rangle}}{\sqrt{q+1}}$
$\displaystyle\to-\frac{1}{2}\left(t_{1}^{2}+\frac{2\mbox{$\;\\!$}t_{1}t_{2}}{\sqrt{q+1}}+t_{2}^{2}\right)\\!,$
(5.19)

by Taylor expanding the left-hand side of (5.19) up to the second order in
parameter $\langle M\rangle^{-1/2}=o(1)$.

We finish this section by stating a theorem concerning the marginal local-type
asymptotics for the length $M_{\lambda}$ and the corresponding conditional
asymptotics for the weight $N_{\lambda}$.

######  Theorem 5.3.

Under Assumptions 3.2 and 5.1, the following distributional asymptotics hold
under the Boltzmann distribution ${\mathsf{P}\mbox{$\>\\!\\!$}}_{\bm{z}}$ on
the space $\check{\varLambda}^{q}$.

  * (a)

(Local limit theorem for $M_{\lambda}$) For any $m$ such that $m-\langle
M\rangle=O\bigl{(}\sqrt{\langle M\rangle}\mbox{$\>\\!$}\bigr{)}$,

$Q_{\bm{z}}(M_{\lambda}=m)\sim f_{\langle M\rangle}(m),$ (5.20)

where $f_{\langle M\rangle}(x)$ is the normal density with mean $\langle
M\rangle$ and standard deviation $\sqrt{\langle M\rangle}$,

$f_{\langle M\rangle}(x)=\frac{1}{\sqrt{2\pi\langle
M\rangle}}\,\exp\left\\{-\frac{(x-\langle
M\rangle)^{2}}{2\mbox{$\>\\!$}\langle M\rangle}\right\\}\\!,\qquad
x\in\mathbb{R}.$

  * (b)

(Conditional limit theorem for $N_{\lambda}$) Conditionally on $M_{\lambda}=m$
with $m-\langle M\rangle=O\bigl{(}\sqrt{\langle
M\rangle}\mbox{$\>\\!$}\bigr{)}$, the partition weight $N_{\lambda}$ is
asymptotically normal with mean $m\mbox{$\>\\!$}\langle N\rangle/\langle
M\rangle$ and variance $q\mbox{$\;\\!$}\langle N\rangle^{2}/\langle M\rangle$,
that is,

$\frac{\sqrt{\langle M\rangle}}{\sqrt{q}}\left(\frac{N_{\lambda}}{\langle
N\rangle}-\frac{m}{\langle
M\rangle}\right)\stackrel{{\scriptstyle\mathrm{d}}}{{\longrightarrow}}_{M_{\lambda}=m}\mathcal{N}(0,1).$
(5.21)

Part (a) of this theorem can be anticipated from the marginal asymptotic
normality of $M_{\lambda}$ proven in Theorem 5.1. The result in part (b) can
be formally derived from the joint asymptotic normality of the pair
$(N_{\lambda},M_{\lambda})$ by calculating the conditional density of
$N_{\lambda}$ given $M_{\lambda}$. Alternatively, one can refer to the “normal
correlation” result (see, e.g., [69, Sec. II.13, Theorem 2, pp. 303–304]). A
rigorous proof of Theorem 5.3 can be obtained by adapting the standard proof
of a “one-dimensional” local limit theorem for $N_{\lambda}$ (see, e.g., [13],
[80]). We will return to this issue in a separate publication.

###  5.2 Limit shape of Young diagrams

In this section we show that, under the slow growth condition on $\langle
M\rangle$, properly scaled Young diagrams of random partitions
$\lambda\in\check{\varLambda}^{q}$ have a limit shape given by the tail of the
gamma integral,

$\omega_{q}^{*}(x):=\frac{1}{\Gamma(1/q)}\int_{x}^{\infty}\\!u^{1/q-1}\mbox{$\;\\!$}\mathrm{e}^{-u}\,\mathrm{d}{u}=1-G_{1/q}(x),\qquad
x\geq 0,$ (5.22)

where $G_{1/q}(x)$ is the distribution function of
$\mathrm{Gamma}\mbox{$\;\\!$}(1/q)$ (see (4.1)). In particular, for $q=1$ the
definition (5.22) is reduced to

$\omega_{1}^{*}(x)=\mathrm{e}^{-x},\qquad x\geq 0.$

Specifically, set

$A:=q\mbox{$\;\\!$}\langle N\rangle/\langle M\rangle,\qquad B=\langle
M\rangle,$ (5.23)

and consider a scaled Young diagram with upper boundary

$\widetilde{Y}_{\lambda}(x)=B^{-1}\mbox{$\;\\!$}Y_{\lambda}(Ax),\qquad x\geq
0,$ (5.24)

where (see (2.2))

$Y_{\lambda}(x)=\sum_{\ell\geq x}\nu_{\ell},\qquad x\geq 0.$ (5.25)

######  Remark 5.2.

The area under the scaled Young diagram is given by

$\int_{0}^{\infty}\mbox{$\>\\!\\!$}\widetilde{Y}_{\lambda}(x)\,\mathrm{d}{x}=B^{-1}\int_{0}^{\infty}\mbox{$\>\\!\\!$}Y_{\lambda}(Ax)\,\mathrm{d}{x}=\frac{\langle
N\rangle}{AB}=\frac{1}{q}.$

Naturally, this condition is preserved by the limit shape; indeed, integrating
by parts we get

$\int_{0}^{\infty}\\!\omega_{q}^{*}(x)\,\mathrm{d}{x}=\frac{1}{\Gamma(1/q)}\int_{0}^{\infty}x^{1/q}\mbox{$\;\\!$}\mathrm{e}^{-x}\,\mathrm{d}{x}=\frac{\Gamma(1+1/q)}{\Gamma(1/q)}=\frac{1}{q}.$

First, we obtain the _expected limit shape_ result.

######  Theorem 5.4.

Under Assumptions 3.2 and 5.1, uniformly in $x\geq 0$

${\mathsf{E}}_{\bm{z}}\mbox{$\;\\!\\!$}\bigl{(}\widetilde{Y}_{\lambda}(x)\bigr{)}\to\omega_{q}^{*}(x),$
(5.26)

where the limit shape $x\mapsto\omega_{q}^{*}(x)$ is defined in (5.22).

###### Proof.

We first show that for each $x\geq 0$, the convergence (5.26) holds. By Lemma
3.1,

$\displaystyle{\mathsf{E}}_{\bm{z}}\bigl{(}Y_{\lambda}(Ax)\bigr{)}=\sum_{\ell\geq
Ax}{\mathsf{E}}_{\bm{z}}(\nu_{\ell})$ $\displaystyle=\sum_{\ell\geq
Ax}\frac{z_{1}^{\ell}z_{2}}{1+z_{1}^{\ell}z_{2}}=\sum_{\ell\geq
Ax}z_{1}^{\ell}z_{2}-\widetilde{R}_{1}(\bm{z}),$ (5.27)

where (cf. (3.27) and (3.28))

$0\leq\widetilde{R}_{1}(\bm{z})=\sum_{\ell\geq
Ax}\frac{z_{1}^{2\ell}\mbox{$\;\\!\\!$}z_{2}^{2}}{1+z_{1}^{\ell}z_{2}}\leq\sum_{\ell\in\mathbb{N}^{q}}z_{1}^{2\ell}\mbox{$\;\\!\\!$}z_{2}^{2}=\braket{M}O(\kappa^{1/q}),$
(5.28)

by virtue of (3.41) (with $r=2$). Furthermore, applying Lemma 3.4 (with
$\gamma=-\log z_{1}$ and $s=0$) and noting that $\gamma A\to 1$, we obtain

$\sum_{\ell\geq
Ax}z_{1}^{\ell}z_{2}\sim\frac{z_{2}}{q\,\gamma^{1/q}}\int_{\gamma
Ax}^{\infty}\mbox{$\>\\!\\!$}u^{1/q-1}\mbox{$\;\\!$}\mathrm{e}^{-u}\,\mathrm{d}{u}\sim\frac{\langle
M\rangle}{\Gamma(1/q)}\int_{x}^{\infty}\\!u^{1/q-1}\mbox{$\;\\!$}\mathrm{e}^{-u}\,\mathrm{d}{u},$
(5.29)

using the asymptotic formula (3.35), Thus, substituting (5.28) and (5.29) into
(5.27) gives

$\displaystyle{\mathsf{E}}_{\bm{z}}\bigl{(}\widetilde{Y}_{\lambda}(x)\bigr{)}=\frac{1}{\langle
M\rangle}\,{\mathsf{E}}_{\bm{z}}\bigl{(}Y_{\lambda}(Ax)\bigr{)}\sim\frac{1}{\Gamma(1/q)}\int_{x}^{\infty}\\!u^{1/q-1}\mbox{$\;\\!$}\mathrm{e}^{-u}\,\mathrm{d}{u}=\omega_{q}^{*}(x),$

as claimed. Finally, the uniform convergence in formula (5.26) follows by
Lemma 3.6, noting that the function $x\mapsto\omega_{q}^{*}(x)$ is continuous,
bounded, and decreasing on $[0,\infty)$. ∎

Now, we are ready to state and prove the main result of this section.

######  Theorem 5.5.

Under Assumptions 3.2 and 5.1, the rescaled Young diagrams converge to the
limit shape $y=\omega_{q}^{*}(x)$ in
${\mathsf{P}\mbox{$\>\\!\\!$}}_{\bm{z}}$-probability uniformly for $x\geq 0$,
that is,

${\mathsf{P}\mbox{$\>\\!\\!$}}_{\bm{z}}\\!\left(\lambda\in\check{\varLambda}^{q}\colon\sup_{x\geq
0}\,\bigl{|}\mbox{$\>\\!$}\widetilde{Y}_{\lambda}(x)-\omega_{q}^{*}(x)\bigr{|}>\varepsilon\right)\to
0.$

###### Proof.

By virtue of Theorem 5.4, letting
$Y^{0}_{\lambda}(x):=Y_{\lambda}(x)-{\mathsf{E}}_{\bm{z}}\bigl{(}Y_{\lambda}(x)\bigr{)}$
it suffices to check that

${\mathsf{P}\mbox{$\>\\!\\!$}}_{\bm{z}}\\!\left(\sup_{x\geq
0}\,\bigl{|}Y^{0}_{\lambda}(Ax)\bigr{|}>B\mbox{$\>\\!$}\varepsilon\right)\to
0.$ (5.30)

Put $Z_{\lambda}(x):=Y_{\lambda}(x^{-1})$ for $0\leq x\leq\infty$; in
particular, $Z_{\lambda}(0)=Y_{\lambda}(\infty)=0$,
$Z_{\lambda}(\infty)=Y_{\lambda}(0)=M_{\lambda}$. By the definition (2.2), for
any $0<s<t\leq\infty$ we have

$Z_{\lambda}(t)-Z_{\lambda}(s)=Y_{\lambda}(t^{-1})-Y_{\lambda}(s^{-1})=\sum_{t^{-1}\leq\mbox{$\>\\!$}\ell\mbox{$\>\\!$}<s^{-1}}\\!\nu_{\ell}\mbox{$\;\\!$},$

which implies that the random process $Z_{\lambda}(x)$ ($x\geq 0$) has
independent increments. Hence,
$Z^{0}_{\lambda}(x):=Z_{\lambda}(x)-{\mathsf{E}}_{\bm{z}}\bigl{(}Z_{\lambda}(x)\bigr{)}$
is a martingale with respect to the filtration
${\mathcal{F}}_{x}=\sigma\\{\nu_{\ell}\mbox{$\>\\!$},\,\ell\geq x^{-1}\\}$.
From (2.2) it is also evident that $Z^{0}_{\lambda}(x)$ is _càdlàg_ (i.e., its
paths are everywhere right-continuous and have left limits, cf. Figure 3).
Therefore, by the Doob–Kolmogorov submartingale inequality (see, e.g., [89,
Theorem 6.16, p. 101]) we obtain

$\displaystyle{\mathsf{P}\mbox{$\>\\!\\!$}}_{\bm{z}}\\!\left(\sup_{x\geq
0}\,\bigl{|}Y^{0}_{\lambda}(Ax)\bigr{|}>B\mbox{$\>\\!$}\varepsilon\right)$
$\displaystyle\equiv{\mathsf{P}\mbox{$\>\\!\\!$}}_{\bm{z}}\\!\left(\sup_{y\leq\infty}\mbox{$\>\\!\\!$}|Z^{0}_{\lambda}(y\mbox{$\>\\!$}A^{-1})|>B\mbox{$\>\\!$}\varepsilon\right)$
$\displaystyle\leq\frac{{\mathsf{Var}}_{\bm{z}}\bigl{(}Z_{\lambda}(\infty)\bigr{)}}{B^{2}\mbox{$\>\\!$}\varepsilon^{2}}=\frac{{\mathsf{Var}}_{\bm{z}}\bigl{(}Y_{\lambda}(0)\bigr{)}}{B^{2}\mbox{$\>\\!$}\varepsilon^{2}}.$
(5.31)

Recalling that $Y_{\lambda}(0)=M_{\lambda}$ and using Theorem 3.9, the right-
hand side of (5.31) is estimated by $O(\langle M\rangle^{-1})$. Thus, the
claim (5.30) follows and the proof of Theorem 5.5 is complete. ∎

Convergence of normalized Young diagrams to their limits shape is illustrated
in Figure 3 for $q=1$ and $q=2$. Random partitions were simulated using a
suitable Boltzmann sampler implemented as Algorithm 1 (see Section 6.1).

$q=2$$q=1$ $\displaystyle\omega_{1}^{*}(x)=\mathrm{e}^{-x}$
---
$\omega^{*}_{2}(x)=1-G_{1/2}(x)$

Figure 3: Illustration of convergence to the limit shape for $q=1$ and $q=2$
(in the online version shown in blue and red, respectively). The step plots
depict the upper boundary of the scaled Young diagrams (see (5.24)), while the
smooth lines represent the limit shape $\omega_{q}^{*}(x)=1-G_{1/q}(x)$ (see
(5.22)). The corresponding partitions $\lambda\in\check{\varLambda}^{q}$ were
sampled using Algorithm 1 with hyper-parameters $\braket{M}=50$ and
$\braket{N}=2.5\cdot 10^{5}$ ($q=1$) or $\braket{N}=1.25\cdot 10^{7}$ ($q=2)$;
in both cases, $\kappa=0.01$ (cf. Assumption 3.2). The respective sample
weight and length are $N_{\lambda}=236{,}\mbox{$\>\\!$}369$, $M_{\lambda}=52$
($q=1$) and $N_{\lambda}=12{,}\mbox{$\>\\!$}733{,}\mbox{$\>\\!$}323$,
$M_{\lambda}=45$ ($q=2$).

Finally, we can analyse asymptotic fluctuations of scaled Young diagrams.

######  Theorem 5.6.

Under Assumptions 3.2 and 5.1, for any $x>0$ the random value
$\widetilde{Y}_{\lambda}(x)$ is asymptotically normal with variance
$\omega_{q}^{*}(x)/\langle M\rangle$, that is,

$\widetilde{Y}^{*}_{\lambda}(x):=\frac{\sqrt{\braket{M}}\left(\widetilde{Y}_{\lambda}(x)-{\mathsf{E}}_{\bm{z}}\bigl{(}\widetilde{Y}_{\lambda}(x)\bigr{)}\right)}{\sqrt{\mbox{$\>\\!$}\omega_{\smash{q}}^{*}(x)}}\stackrel{{\scriptstyle\mathrm{d}}}{{\longrightarrow}}\mathcal{N}(0,1).$
(5.32)

###### Proof.

Consider the characteristic function of $\widetilde{Y}^{*}_{\lambda}(x)$,

$\varphi_{\bm{z}}(t;x):={\mathsf{E}}_{\bm{z}}\bigl{[}\exp\mbox{$\;\\!\\!$}\bigl{(}\mathrm{i}\mbox{$\;\\!$}t\mbox{$\;\\!$}\widetilde{Y}^{*}_{\lambda}(x)\bigr{)}\bigr{]}\qquad(t\in\mathbb{R}).$
(5.33)

Substituting the definition (5.32) and using (5.23), (5.24) and (5.25), this
is transformed as

$\varphi_{\bm{z}}(t;x)=\exp\mbox{$\;\\!\\!$}\bigl{\\{}-\mathrm{i}\mbox{$\;\\!$}\tilde{t}\,{\mathsf{E}}_{\bm{z}}\bigl{(}Y_{\lambda}(Ax)\bigr{)}\mbox{$\;\\!$}\bigr{\\}}\,{\mathsf{E}}_{\bm{z}}\bigl{[}\exp\mbox{$\;\\!\\!$}\bigl{(}\mathrm{i}\mbox{$\;\\!$}\tilde{t}\,Y_{\lambda}(Ax)\bigr{)}\bigr{]},\qquad\tilde{t}:=\frac{t}{\sqrt{\braket{M}\omega_{\smash{q}}^{*}(x)}}\mbox{$\>\\!$},$
(5.34)

and furthermore (see (5.27))

${\mathsf{E}}_{\bm{z}}\bigl{(}Y_{\lambda}(Ax)\bigr{)}=\sum_{\ell\geq
Ax}\frac{z_{1}^{\ell}z_{2}}{1+z_{1}^{\ell}z_{2}}.$ (5.35)

Next, similarly to (5.6) and (5.7) the last expectation in (5.34) is expressed
as

${\mathsf{E}}_{\bm{z}}\\!\left[\exp\left(\mathrm{i}\mbox{$\;\\!$}\tilde{t}\sum_{\ell\geq
Ax}\nu_{\ell}\right)\right]=\prod_{\ell\geq
Ax}\frac{1+z_{1}^{\ell}z_{2}\,\mathrm{e}^{\mbox{$\;\\!$}\mathrm{i}\mbox{$\>\\!$}\tilde{t}}}{1+z_{1}^{\ell}z_{2}}=\prod_{\ell\geq
Ax}\bigl{(}1+w_{\ell}(\tilde{t}\mbox{$\>\\!$})\bigr{)},$ (5.36)

where

$w_{\ell}(t):=\frac{z_{1}^{\ell}z_{2}}{1+z_{1}^{\ell}z_{2}}\,(\mathrm{e}^{\mbox{$\;\\!$}\mathrm{i}\mbox{$\>\\!$}t}-1).$
(5.37)

Choosing the principal branch of the logarithm and using (5.35) and (5.36),
from (5.34) we get

$\log\varphi_{\bm{z}}(t;x)=-\mathrm{i}\mbox{$\;\\!$}\tilde{t}\sum_{\ell\geq
Ax}\frac{z_{1}^{\ell}z_{2}}{1+z_{1}^{\ell}z_{2}}+\sum_{\ell\geq
Ax}\log\bigl{(}1+w_{\ell}(\tilde{t}\mbox{$\>\\!$})\bigr{)}.$ (5.38)

In turn, similarly to (5.12) we obtain

$\sum_{\ell\geq
Ax}\log\bigl{(}1+w_{\ell}(\tilde{t}\mbox{$\>\\!$})\bigr{)}=\Bigl{(}\mathrm{i}\mbox{$\;\\!$}\tilde{t}-\tfrac{1}{2}\mbox{$\>\\!$}\tilde{t}^{2}+O(\tilde{t}^{3})\Bigr{)}\sum_{\ell\geq
Ax}\frac{z_{1}^{\ell}z_{2}}{1+z_{1}^{\ell}z_{2}}+O(\tilde{t}^{2})\sum_{\ell\geq
Ax}\frac{z_{1}^{2\ell}z_{2}^{2}}{(1+z_{1}^{\ell}z_{2})^{2}}.$ (5.39)

As was shown in the proof of Theorem 5.4 (see (5.27), (5.28)) and (5.29)),

$\sum_{\ell\geq Ax}\frac{z_{1}^{\ell}z_{2}}{1+z_{1}^{\ell}z_{2}}\sim\langle
M\rangle\,\omega_{q}^{*}(x),\qquad\sum_{\ell\geq
Ax}\frac{z_{1}^{2\ell}z_{2}^{2}}{(1+z_{1}^{\ell}z_{2})^{2}}=\braket{M}O(\kappa^{1/q})=\braket{M}o(1).$
(5.40)

Using (5.40) and recalling the notation of $\tilde{t}$ in (5.34), from (5.39)
we get

$\sum_{\ell\geq
Ax}\log\bigl{(}1+w_{\ell}(\tilde{t}\mbox{$\>\\!$})\bigr{)}=\mathrm{i}\mbox{$\;\\!$}\tilde{t}\sum_{\ell\geq
Ax}\frac{z_{1}^{\ell}z_{2}}{1+z_{1}^{\ell}z_{2}}-\tfrac{1}{2}\mbox{$\>\\!$}t^{2}+o(1).$

Finally, returning to (5.38), we see that
$\varphi_{\bm{z}}(t;x)\to-\tfrac{1}{2}\mbox{$\>\\!$}t^{2}$, which proves the
theorem. ∎

######  Remark 5.3.

The reason for using in (5.32) the intrinsic centering
${\mathsf{E}}_{\bm{z}}\bigl{(}\widetilde{Y}_{\lambda}(x)\bigr{)}$ rather than
the limit shape value $\omega_{q}^{*}(x)$ is that the error terms in the
asymptotic estimates (5.40) are of order $\braket{M}\kappa^{1/q}$, where
$\kappa=\langle M\rangle^{q+1}/\langle N\rangle=o(1)$ (see Assumption 5.1).
Combined with the factor $\tilde{t}=O\bigl{(}\braket{M}^{-1/2}\bigr{)}$, this
produces the error bound of order $\langle
M\rangle^{1/2}\mbox{$\;\\!$}\kappa^{1/q}$, which is not guaranteed to be
small. Thus, a stronger assumption to this end is $\langle
M\rangle^{1/2}\mbox{$\>\\!$}\kappa^{1/q}=o(1)$, that is, $\langle
M\rangle^{1+3\mbox{$\>\\!$}q/2}/\langle N\rangle=o(1)$. On the other hand,
lifting any control over the length may restore the limit-shape centering; for
example, for $q=1$ (ordinary strict partitions), a central theorem of that
kind was proved in [80].

###  5.3 A joint limit theorem for the extreme parts (growing expected length)

We use the notation $\lambda_{\rm min}$ and $\lambda_{\rm max}$ (cf. Section
4.3).

######  Theorem 5.7.

Under Assumptions 3.2 and 5.1, set

$b_{q}:=\left(\frac{q\mbox{$\;\\!\\!$}\braket{M}}{\displaystyle\Gamma(1/q)}\right)^{\mbox{$\;\\!\\!$}q}\\!,\qquad
B_{q}:=\log\mbox{$\;\\!\\!$}\langle
M\rangle-\bigl{(}1-\tfrac{1}{q}\bigr{)}\log\log\mbox{$\;\\!\\!$}\langle
M\rangle-\log\Gamma(1/q),$ (5.41)

and consider the normalized versions of $\lambda_{\rm min}$ and $\lambda_{\rm
max}$ defined as follows,

$\lambda_{\rm min}^{*}:=\gamma\mbox{$\;\\!$}b_{q}\mbox{$\;\\!$}\lambda_{\rm
min},\qquad\lambda_{\rm max}^{*}:=\gamma\mbox{$\>\\!$}\lambda_{\rm
max}-B_{q}\mbox{$\;\\!$}.$ (5.42)

Then $\lambda^{*}_{\rm min}$ and $\lambda^{*}_{\rm max}$ are asymptotically
independent under the measure ${\mathsf{P}\mbox{$\>\\!\\!$}}_{\bm{z}}$ as
$\braket{N}\to\infty$ and their marginal limiting laws are given,
respectively, by a Weibull distribution with shape parameter $1/q$ and the
standard double-exponential (Gumbel) distribution,

$\displaystyle{\mathsf{P}\mbox{$\>\\!\\!$}}_{\bm{z}}(\lambda_{\rm
min}^{*}>x_{1})$
$\displaystyle\to\exp\mbox{$\;\\!\\!$}\bigl{(}-x_{1}^{1/q}\bigr{)},\qquad$
$\displaystyle x_{1}\geq 0,$ (5.43)
$\displaystyle{\mathsf{P}\mbox{$\>\\!\\!$}}_{\bm{z}}(\lambda_{\rm max}^{*}\leq
x_{2})$
$\displaystyle\to\exp\mbox{$\;\\!\\!$}\bigl{(}-\mathrm{e}^{-x_{2}}\bigr{)},$
$\displaystyle x_{2}\in\mathbb{R}.$ (5.44)

###### Proof.

For $x_{1}\geq 0$ and $x_{2}\in\mathbb{R}$, set

$\ell^{*}_{1}(x_{1}):=\min\\{\ell\in\mathbb{N}^{q}\colon\ell>x_{1}/(\gamma\mbox{$\;\\!$}b_{q})\\},\qquad\ell^{*}_{2}(x_{2}):=\min\\{\ell\in\mathbb{N}^{q}\colon\ell>(B_{q}+x_{2})\mbox{$\;\\!$}\gamma^{-1}\\}.$

Recalling the asymptotic relations (3.25) and (3.22), observe that

$\displaystyle\gamma\mbox{$\;\\!$}\ell^{*}_{1}(x_{1})$
$\displaystyle=\frac{x_{1}}{b_{q}}+O(\gamma)\sim
x_{1}\left(\frac{\Gamma(1/q)}{\displaystyle
q\braket{M}}\right)^{\mbox{$\;\\!\\!$}q}\\!,$ (5.45)
$\displaystyle\gamma\mbox{$\;\\!$}\ell^{*}_{2}(x_{2})$
$\displaystyle=B_{q}+x_{2}+O(\gamma)\sim\log\braket{M}\mbox{$\>\\!\\!$}.$
(5.46)

Like in the proof of Theorem 4.4, we have

$\displaystyle{\mathsf{P}\mbox{$\>\\!\\!$}}_{\bm{z}}(\lambda_{\rm
min}^{*}>x_{1},\,\lambda_{\rm max}^{*}\leq x_{2})$
$\displaystyle=\exp\Biggl{\\{}-\Biggl{(}\mbox{$\>\\!$}\sum_{\ell\in\mathbb{N}^{q}}^{\infty}-\sum_{\ell\geq\ell^{*}_{1}(x_{1})}^{\infty}+\sum_{\ell\geq\ell^{*}_{2}(x_{2})}^{\infty}\Biggr{)}\log\mbox{$\>\\!\\!$}\\!\left(1+z_{1}^{\ell}z_{2}\right)\Biggr{\\}}.$

Applying Lemmas 3.3 and 3.5 (with $\gamma=-\log z_{1}$ and $\eta=z_{2}$) and
using the asymptotic relation (3.35), we obtain

$-\log{\mathsf{P}\mbox{$\>\\!\\!$}}_{\bm{z}}(\lambda^{*}_{\rm
min}>x_{1},\,\lambda^{*}_{\rm max}\leq
x_{2})\sim\frac{\braket{M}}{\Gamma(1/q)}\left(\int_{0}^{\gamma\mbox{$\>\\!$}\ell^{*}_{1}(x_{1})}\\!+\int_{\gamma\mbox{$\>\\!$}\ell^{*}_{2}(x_{2})}^{\infty}\right)u^{1/q-1}\mbox{$\;\\!$}\mathrm{e}^{-u}\,\mathrm{d}{u}.$
(5.47)

Integrating by parts and using the asymptotic relation (5.45), we obtain

$\int_{0}^{\gamma\mbox{$\>\\!$}\ell^{*}_{1}(x_{1})}\mbox{$\>\\!\\!$}u^{1/q-1}\mbox{$\;\\!$}\mathrm{e}^{-u}\,\mathrm{d}{u}=q\int_{0}^{\gamma\mbox{$\>\\!$}\ell^{*}_{1}(x_{1})}\\!\mathrm{e}^{-u}\,\mathrm{d}{(u^{1/q})}\sim
q\mbox{$\;\\!$}\bigl{(}\gamma\mbox{$\;\\!$}\ell^{*}_{1}(x_{1})\bigr{)}^{1/q}\sim\frac{\Gamma(1/q)\,x_{1}^{1/q}}{\displaystyle\braket{M}}.$
(5.48)

Next, using (5.46) we get

$\displaystyle\int_{\gamma\mbox{$\>\\!$}\ell^{*}_{2}(x_{2})}^{\infty}\mbox{$\>\\!\\!$}u^{1/q-1}\mbox{$\;\\!$}\mathrm{e}^{u}\,\mathrm{d}{u}$
$\displaystyle\sim\bigl{(}\gamma\mbox{$\;\\!$}\ell^{*}_{2}(x_{2})\bigr{)}^{1/q-1}\mbox{$\;\\!$}\mathrm{e}^{-\gamma\mbox{$\>\\!$}\ell^{*}_{2}(x_{2})}$
$\displaystyle\sim
B_{q}^{1/q-1}\mbox{$\;\\!$}\mathrm{e}^{-(B_{q}+x_{2})}\sim\frac{\Gamma(1/q)\,\mathrm{e}^{-x_{2}}}{\braket{M}}.$
(5.49)

Hence, substituting (5.48) and (5.49) into (5.47) yields

$\displaystyle-\log{\mathsf{P}\mbox{$\>\\!\\!$}}_{\bm{z}}(\lambda^{*}_{\rm
min}>x_{1},\,\lambda^{*}_{\rm max}\leq x_{2})\sim
x_{1}^{1/q}+\mathrm{e}^{-x_{2}},$

which completes the proof of the theorem. ∎

######  Remark 5.4.

The necessary use of the intrinsic calibration parameter $\gamma=-\log z_{1}$
in Theorem 5.7 may be a little disappointing. This can be easily improved
under a slightly stronger condition on slow growth of $\langle M\rangle$ than
in Assumption 5.1; namely, $\kappa^{1/q}\log\langle M\rangle=o(1)$, that is,
$\langle M\rangle^{q+1}\mbox{$\;\\!$}(\log\mbox{$\>\\!\\!$}\langle
M\rangle)^{q}/\langle N\rangle=o(1)$. In this case, the normalization (5.42)
can be written more explicitly by replacing $\gamma$ with $\gamma_{0}=\langle
M\rangle/(q\mbox{$\;\\!$}\langle N\rangle)$ (see (3.34)).

######  Corollary 5.8.

Under the hypotheses of Theorem 5.7, the following law of large numbers holds,

$\frac{\langle M\rangle\mbox{$\>\\!$}\lambda_{\rm max}}{q\mbox{$\;\\!$}\langle
N\rangle\log\mbox{$\>\\!\\!$}\langle
M\rangle}\stackrel{{\scriptstyle\mathrm{p}}}{{\longrightarrow}}1,$ (5.50)

where the symbol $\stackrel{{\scriptstyle\mathrm{p}}}{{\to}}$ indicates
convergence in ${\mathsf{P}\mbox{$\>\\!\\!$}}_{\bm{z}}$-probability.

###### Proof.

Theorem 5.7 implies that $\gamma\mbox{$\>\\!$}\lambda_{\rm
max}/B_{q}\stackrel{{\scriptstyle\mathrm{p}}}{{\to}}1$, and the claim (5.50)
follows by noting that $\gamma\sim\gamma_{0}=\langle
M\rangle/(q\mbox{$\;\\!$}\langle N\rangle)$ and
$B_{q}\sim\log\mbox{$\>\\!\\!$}\langle M\rangle$. ∎

######  Remark 5.5.

Theorem 5.7 indicates that, under the condition of slow growth of
$\braket{M}$, the smallest part $\lambda_{\mathrm{min}}$ of a
${\mathsf{P}\mbox{$\>\\!\\!$}}_{\bm{z}}$-typical partition
$\lambda\in\check{\varLambda}^{q}$ “lives” on the scale
$A_{*}=(\gamma\mbox{$\;\\!$}b_{q})^{-1}\propto\langle N\rangle/\langle
M\rangle^{q+1}=\kappa^{-1}$. On the other hand, Corollary 5.8 shows that the
scale of variation of the largest part $\lambda_{\mathrm{max}}$ is given by
$A^{*}\mbox{$\>\\!\\!$}=B_{q}\mbox{$\;\\!$}\gamma^{-1}\\!\propto\\!\langle
N\rangle\log\mbox{$\>\\!\\!$}\langle M\rangle/\langle M\rangle$. This is to be
compared with the typical behaviour in the bulk of the partition “spectrum”,
where the scale of variation is given by $A\sim\gamma^{-1}\propto\langle
N\rangle/\langle M\rangle$.

######  Remark 5.6.

Continuing an asymptotic linkage between the cases of fixed or slowly growing
parameter $\langle M\rangle$, observed above in Remark 5.1, the limiting
distributions of Theorem 4.4 formally conform to Theorem 5.7 in the limit as
$\langle M\rangle\to\infty$. Indeed, using (4.43) we have

$\displaystyle-\log G_{\rm max}(x+B_{q})$
$\displaystyle=\frac{\braket{M}\mbox{$\>\\!$}\Gamma(1/q,x+B_{q})}{\Gamma(1/q)}$
$\displaystyle\sim\frac{\braket{M}}{\Gamma(1/q)}\,(x+B_{q})^{1/q-1}\mbox{$\;\\!$}\mathrm{e}^{-x-B_{q}}$
$\displaystyle\sim\frac{\braket{M}}{\Gamma(1/q)}\,(\log\mbox{$\>\\!\\!$}\langle
M\rangle)^{1-1/q}\,\mathrm{e}^{-x}\cdot\mathrm{e}^{-B_{q}}=\mathrm{e}^{-x},$

according to the definition of $B_{q}$ in (5.41). Similarly, using (4.44) and
(5.48) we have

$\displaystyle-\log G^{c}_{\rm min}(x/b_{q})$
$\displaystyle=\braket{M}\\!\left(1-\frac{\Gamma(1/q,x/b_{q})}{\Gamma(1/q)}\right)$
$\displaystyle\sim\frac{\braket{M}\mbox{$\;\\!\\!$}q}{\Gamma(1/q)}\left(\frac{x}{b_{q}}\right)^{1/q}=x^{1/q},$

by the definition of $b_{q}$ in (5.41).

######  Remark 5.7 (Possible generalization to non-integral powers).

The reader may have observed that our main results (including Theorems 4.1 and
5.1 about joint limiting laws of weight $N_{\lambda}$ and length
$M_{\lambda}$; Theorems 4.4 and 5.7 about the joint asymptotics of extreme
parts; the limit shape Theorems 5.5 and 5.6; and the cardinality Theorem 4.3)
continue to make sense _for any real power $q>0$_, although our proofs proceed
from the assumption that $q$ is integer, which makes the partition model
combinatorically well defined. The interest in non-integral powers is not new
but was mostly motivated by applications in statistical physics [3, 19, 68]. A
mathematically meaningful interpretation of partitions with non-integral power
parts may be based on the idea recently proposed by Lipnik et al. [55],
whereby the multiplicative Boltzmann structure is introduced on the hidden
“substrate” $\mathbb{N}=\\{k\\}$, from which the $q$-power parts of the random
partition $\lambda=(\ell_{i})$ are formed as $\ell=\lfloor k^{q}\rfloor$, with
any $q\in\mathbb{R}_{+}$. We will review and extend our results under this
approach in a separate paper.

##  6 Application to random sampling

_Boltzmann sampling_ is a powerful technique conceptualized, streamlined and
popularized by Duchon et al. [23] in the context of single-parameter
combinatorial structures (for multi-parametric extensions, see Bendkowski et
al. [10] and the references therein). Random integer partitions with
controlled expected weight and length provide an “exactly soluble” instance of
a two-parametric combinatorial structure, where the issues of Boltzmann
sampling implementation and efficiency can be analysed in some depth.

Specifically, in this section we discuss sampling from the Boltzmann
distribution on partition spaces $\check{\varLambda}^{q}$ (i.e., into distinct
$q$-power parts), calibrated under the predefined hyper-parameters
$\braket{N}$ and $\braket{M}$, which have the meaning of the expected weight
and length, respectively. The two controlling parameters in question are
$z_{1}$ and $z_{2}$, which are amenable to asymptotic analysis as was shown in
Section 3.3. Once these parameters are fixed, due to the mutual independence
of the multiplicities $(\nu_{\ell})$ (see Proposition 2.1 and Lemma 3.1), the
Boltzmann sampling is essentially reduced to an iterated independent testing
of potential parts $\ell=j^{q}$ via dichotomous (Bernoulli) random trials with
success probabilities
${\mathsf{P}\mbox{$\>\\!\\!$}}_{\bm{z}}(\nu_{\ell}=1)=z_{1}^{\ell}z_{2}\mbox{$\;\\!$}(1+z_{1}^{\ell}z_{2})^{-1}$.
The practical implementation of such sampling algorithms thus relies on a
random number generator $\mathtt{Ber}(p)$, in each call producing an
independent pseudo-random value $1$ or $0$ with probabilities $p$ and $1-p$,
respectively.

It is convenient to distinguish between the _free samplers_ and the _rejection
samplers_ , with the former just producing independent random realizations of
partitions under the said Boltzmann distribution, and the latter comprising
one or more rejection loops that iterate a free Boltzmann sampler until the
desired targets are met. We discuss these two versions in Sections 6.1 and
6.2, respectively.

Computer codes were implemented using the programming language C and Intel®
oneAPI DPC++ compiler, and run on a desktop CPU Intel® CoreTM i5-10600
(processor base frequency 3.30 GHz, turbo boost frequency 4.80 GHz). Numerical
calculations were carried out using MapleTM (Release 2022.1, licensed to the
University of Leeds).

###  6.1 Free sampler

In this subsection, we delineate a free Boltzmann sampler (see Algorithm 1
below) under the calibration through the hyper-parameters $\braket{N}$ and
$\braket{M}$. It should be noted that, despite an intuitive appeal of iterated
Bernoulli-type tests, there are some implementation concern that have to be
addressed. We discuss them below before presenting the algorithm.

####  6.1.1 Correcting the bias

The first issue to consider is that of choosing the control parameters $z_{1}$
and $z_{2}$ to ensure that the sampler is unbiased, that is,
${\mathsf{E}}_{\bm{z}}(N_{\lambda})=\braket{N}$ and
${\mathsf{E}}_{\bm{z}}(M_{\lambda})=\braket{M}$. Unfortunately, we can solve
this set of equations only asymptotically (see Lemma 3.1). In a “crude”
version of Algorithm 1, we use the leading terms in the asymptotics by setting
(cf. (3.25) and (3.26))

$z_{1}=\mathrm{e}^{-\gamma_{0}},\qquad
z_{2}=\frac{\braket{M}\gamma_{0}^{1/q}}{\Gamma(1+1/q)},$ (6.1)

where $\gamma_{0}=\langle M\rangle/(q\mbox{$\;\\!\\!$}\braket{N})$.
Inevitably, this causes a bias in the resulting expectations. More precisely,
the first source of this bias clearly comes from dropping the (positive)
remainder terms $R_{1}(\bm{z})$ and $R_{2}(\bm{z})$ in the approximate series
representations of the aforementioned expected values (see equations (3.27)
and (3.30)). A further error occurs when replacing the resulting series with
the corresponding integrals, using Lemma 3.2.

In fact, one can show that the overall bias due to (6.1) is always negative.
Indeed, recalling that $\Delta_{0}(\gamma)<0$ (see (3.11)101010The inequality
$\Delta_{0}(\gamma)<0$ can also be seen directly by monotonicity of the
function $x\mapsto\mathrm{e}^{-\gamma\mbox{$\>\\!$}x^{q}}$.), we have

$\displaystyle{\mathsf{E}}_{\bm{z}}(M_{\lambda})=\sum_{\ell\in\mathbb{N}^{q}}\frac{z_{1}^{\ell}z_{2}}{1+z_{1}^{\ell}z_{2}}<z_{2}\sum_{\ell\in\mathbb{N}^{q}}z_{1}^{\ell}$
$\displaystyle=z_{2}\sum_{j=1}^{\infty}\mathrm{e}^{-\gamma_{0}\mbox{$\>\\!$}j^{q}}$
$\displaystyle<z_{2}\int_{0}^{\infty}\\!\mathrm{e}^{-\gamma_{0}\mbox{$\>\\!$}x^{q}}\mathrm{d}{x}=\frac{z_{2}\,\Gamma(1/q)}{q\,\gamma_{0}^{1/q}}=\langle
M\rangle,$

according to the parameter choice (6.1). Turning to
${\mathsf{E}}_{\bm{z}}(N_{\lambda})$, recall from (3.31) that the error term
$-R_{2}(\bm{z})$ is negative, and furthermore,

$R_{2}(\bm{z})\sim
z_{2}^{2}\sum_{\ell\in\mathbb{N}^{q}}\ell\mbox{$\>\\!$}z_{1}^{2\ell}\sim\frac{z_{2}\,q\mbox{$\;\\!\\!$}\braket{M}\mbox{$\;\\!\\!$}\gamma_{0}^{1/q}}{\Gamma(1/q)}\cdot\frac{\Gamma(1+1/q)}{q\,\gamma_{0}^{1+{1/q}}}=\frac{z_{2}\braket{M}}{q\,\gamma_{0}}.$
(6.2)

On the other hand, the error due to replacing the sum
$\sum_{\ell}\ell\mbox{$\>\\!$}z_{1}^{\ell}$ in (3.30) by the corresponding
integral is bounded, according to (3.16), by
$z_{2}\,O\bigl{(}\gamma_{0}^{-1+1/q}\bigr{)}$. Since $\braket{M}$ is bounded
away from zero, it follows that the $R_{2}$-term (6.2) is dominant and,
therefore, the overall bias in targeting $\braket{N}$ is negative.

A practical recipe towards correcting the bias may be to move the error terms
$R_{1}(\bm{z})$ and $R_{2}(\bm{z})$ to the left-hand side of equations in
(3.27) and (3.30), for simplicity using their integral approximations.
Effectively, this amounts to redefining the hyper-parameters,

$\displaystyle\braket{\tilde{N}}$
$\displaystyle:=\braket{N}+z_{2}^{2}\sum_{\ell\in\mathbb{N}^{q}}\ell\mbox{$\>\\!$}z_{1}^{2\ell}\approx\braket{N}+z_{2}^{2}\int_{0}^{\infty}\\!x^{q}\,\mathrm{e}^{-2\mbox{$\>\\!$}\gamma_{0}\mbox{$\>\\!$}x^{q}}\mbox{$\>\\!$}\mathrm{d}{x}=\braket{N}+\frac{\braket{M}^{2}\gamma_{0}^{1/q}}{2^{1/q}\,\Gamma(1/q)},$
$\displaystyle\braket{\tilde{M}}$
$\displaystyle:=\braket{M}+z_{2}^{2}\sum_{\ell\in\mathbb{N}^{q}}z_{1}^{2\ell}\approx\braket{M}+z_{2}^{2}\int_{0}^{\infty}\\!\mathrm{e}^{-2\mbox{$\>\\!$}\gamma_{0}\mbox{$\>\\!$}x^{q}}\mbox{$\>\\!$}\mathrm{d}{x}=\braket{N}+\frac{q\mbox{$\;\\!\\!$}\braket{M}^{2}\gamma_{0}^{1/q}}{2^{1/q}\,\Gamma(1/q)}.$

Accordingly, we redefine
$\tilde{\gamma_{0}}=\langle\tilde{M}\rangle/(q\mbox{$\;\\!\\!$}\braket{\tilde{N}})$
and (cf. (6.1))

$\tilde{z}_{1}=\mathrm{e}^{-\tilde{\gamma_{0}}},\qquad\tilde{z}_{2}=\frac{q\mbox{$\;\\!\\!$}\braket{\tilde{M}}\tilde{\gamma_{0}}^{1/q}}{\Gamma(1/q)}.$
(6.3)

A numerical illustration of the proposed modification is presented in Table 1,
showing a significant reduction of bias. The expected values were computed
from the exact series expansions (3.27) and (3.30) using Maple.

Table 1: Expected values of $N_{\lambda}$ and $M_{\lambda}$ for $q=1$ and $q=2$ under two choices of the calibrating parameters: using the leading asymptotic terms (6.1) and after a heuristic correction (6.3). $q=1$ | $q=2$
---|---
$\vphantom{\int^{)}}\braket{N}=10^{6}$ | $\braket{\tilde{N}}\doteq 1{,}\mbox{$\>\\!$}002{,}\mbox{$\>\\!$}499.50$ | $\braket{N}=10^{7}$ | $\braket{\tilde{N}}\doteq 10{,}\mbox{$\>\\!$}315{,}\mbox{$\>\\!$}391.57$
$\vphantom{\int^{)}}\braket{M}=100$ | $\braket{\tilde{M}}\doteq 100.4999$ | $\braket{M}=50$ | $\braket{\tilde{M}}\doteq 53.1539$
$z_{1}\doteq 0.9999000$ | $\tilde{z}_{1}=0.9998998$ | $z_{1}\doteq 0.9999975$ | $\tilde{z}_{1}\doteq 0.9999974$
$z_{2}\doteq 0.010000$ | $\tilde{z}_{2}\doteq 0.0100750$ | $z_{2}\doteq 0.0892062$ | $\tilde{z}_{2}\doteq 0.0962720$
${\mathsf{E}}_{\bm{z}}(N_{\lambda})\doteq 997{,}\mbox{$\>\\!$}510.70$ | ${\mathsf{E}}_{\tilde{z}}(N_{\lambda})\doteq 999{,}\mbox{$\>\\!$}985.73$ | ${\mathsf{E}}_{\bm{z}}(N_{\lambda})\doteq 9{,}\mbox{$\>\\!$}699{,}\mbox{$\>\\!$}070.63$ | ${\mathsf{E}}_{\tilde{z}}(N_{\lambda})\doteq 9{,}\mbox{$\>\\!$}981{,}\mbox{$\>\\!$}802.71$
${\mathsf{E}}_{\bm{z}}(M_{\lambda})\doteq 99.498326$ | ${\mathsf{E}}_{\tilde{z}}(M_{\lambda})\doteq 99.992019$ | ${\mathsf{E}}_{\bm{z}}(M_{\lambda})\doteq 47.018388$ | ${\mathsf{E}}_{\tilde{z}}(M_{\lambda})\doteq 49.754495$

If the remaining (small) bias is still an issue, a further recalibration can
be carried out by a suitable refinement of the solution $\bm{z}=(z_{1},z_{2})$
to the equations (3.27) and (3.30), for instance, by using a two-dimensional
Newton–Raphson method. For a general approach to the multidimensional tuning
of parameters based on convex optimization, see Bendkowski et al. [10].

####  6.1.2 Truncation of the parts pipeline

We deal with a finitary computation, so should rule out the risk of indefinite
processing. That is to say, the algorithm must have a well-defined stopping
rule that would guarantee a finite-time termination. In a free sampler, the
sequence of productive outcomes in successive Bernoulli trials (i.e., with
sample multiplicities $\nu_{\ell}=1$ corresponding to non-zero parts) is
${\mathsf{P}\mbox{$\>\\!\\!$}}_{\bm{z}}$-a.s. finite (see Lemma 2.2). More
precisely, the last successful trial selects the largest part $\lambda_{\rm
max}$, after which the testing settles down to pure idling. The
${\mathsf{P}\mbox{$\>\\!\\!$}}_{\bm{z}}$-distribution of $\lambda_{\rm max}$
is given by (cf. Sections 4.3 and 5.3)

$\displaystyle{\mathsf{P}\mbox{$\>\\!\\!$}}_{\bm{z}}(\lambda_{\rm max}\leq L)$
$\displaystyle={\mathsf{P}\mbox{$\>\\!\\!$}}_{\bm{z}}(\nu_{\ell}\equiv 0\
\text{for all}\ \ell>L)$
$\displaystyle=\prod_{\ell>L}\frac{1}{1+z_{1}^{\ell}z_{2}}=\frac{1}{F(\bm{z})}\prod_{\ell\leq
L}(1+z_{1}^{\ell}z_{2}),$ (6.4)

where $F(\bm{z})=\prod_{\ell\in\mathbb{N}^{q}}(1+z_{1}^{\ell}z_{2})$ is the
generating function of the partition space $\check{\varLambda}^{p}$ (see
(3.2)). It is also easy to see that conditioning on $\lambda_{\rm
max}=j_{0}^{q}$ does not change the distribution of the preceding
multiplicities $\\{\nu_{j^{q}},1\leq j\leq j_{0}-1\\}$, that is, they remain
mutually independent and with Bernoulli distributions (3.1). Thus, if the
numerical value of $F(\bm{z})$ can be calculated in advance, which is a common
convention in computing known as an _oracle_ (see, e.g., [23, 30, 10]), then
we can sample the random value $\lambda_{\rm max}$ using formula (6.4) and
then sample independently the preceding candidate parts via the respective
Bernoulli trials. Unfortunately, this approach embeds a computational error
through the numerical calculation of $F(\bm{z})$, so it is not quite “exact”;
besides, convergence of the infinite product may not be fast, given that the
parameter $z_{1}$ is close to $1$ (see (3.25)). Specifically, using Lemma 3.5
one can check that the truncation error arising from a partial product up to
$\ell_{*}$ is of order
$\braket{M}(\gamma\mbox{$\>\\!$}\ell_{*})^{1/q-1}\mbox{$\>\\!$}\mathrm{e}^{-\gamma\mbox{$\>\\!$}\ell_{*}}$,
which dictates that $\ell_{*}$ be chosen much bigger than $\gamma^{-1}\sim
q\mbox{$\;\\!$}\langle N\rangle/\langle M\rangle$.

An alternative idea is to truncate the pipeline of potential parts
$\ell\in\mathbb{N}^{q}$ subject to testing at an appropriate threshold $L$
(see Section 2.5), so that the Bernoulli testing only runs over $\ell\leq L$.
A simple pragmatic solution is to choose the threshold $L$ so that the
probability of exceeding it in an indefinite free sampler is small enough,
that is,

${\mathsf{P}\mbox{$\>\\!\\!$}}_{\bm{z}}(\lambda_{\mathrm{max}}>L)\leq\delta,$
(6.5)

where the confidence tolerance $\delta>0$ can be chosen in advance to be as
small as desired. Then the corresponding threshold $L=L(\delta)$ can be
determined from a suitable limit theorem for the largest part, namely, Theorem
4.4 if $\braket{M}>0$ is fixed, or Theorem 5.7 for slow growth of
$\braket{M}$. In the former case, threshold $L$ is determined by the
asymptotic equation (see (4.43))

$\Gamma(1/q,\gamma_{0}L)=\frac{\Gamma(1/q)}{\braket{M}}\cdot\log\frac{1}{1-\delta},$
(6.6)

where, as before,
$\gamma_{0}=\braket{M}\mbox{$\>\\!\\!$}/(q\mbox{$\;\\!\\!$}\braket{N})$. In
the latter case, we obtain from (5.44)

$L=\frac{1}{\gamma_{0}}\left(B_{q}-\log\log\frac{1}{1-\delta}\right)\\!,$
(6.7)

where (see (5.41))

$B_{q}=\log\mbox{$\;\\!\\!$}\langle
M\rangle-\bigl{(}1-\tfrac{1}{q}\bigr{)}\log\log\mbox{$\;\\!\\!$}\langle
M\rangle-\log\Gamma(1/q).$

Note that for $q=1$ the bounds (6.6) and (6.7) coincide, reducing to

$L=\frac{1}{\gamma_{0}}\left(\log\mbox{$\;\\!\\!$}\langle
M\rangle-\log\log\frac{1}{1-\delta}\right)\\!.$ (6.8)

An illustration of evaluation of the threshold $L$ is presented in Table 2 for
$q=1$ and $q=2$. The equation (6.6) was solved numerically using Maple. One
can observe from the table that while the value $\delta=10^{-k}$
($k=1,2,\dots$) is decreasing geometrically, the growth of threshold $L$ is
only about linear. Intuitively, this is explained by the fact that the
confidence probability $1-\delta$ enters expressions (6.6) and (6.7) under the
logarithm, that is, as $\log\left(1-\delta\right)$. More precisely, it readily
follows from (6.7) that

$L-\frac{B_{q}}{\gamma_{0}}\sim\frac{1}{\gamma_{0}}\log\frac{1}{\delta}\qquad(\delta\to
0{+}).$

Likewise, equation (6.6) asymptotically solves to yield

$L-\frac{\log\braket{M}}{\gamma_{0}\,\Gamma(1/q)}\sim\frac{1}{\gamma_{0}}\log\frac{1}{\delta}\qquad(\delta\to
0{+}).$

Table 2: Threshold $L$ for the largest part $\lambda_{\rm{max}}$ with confidence probability $1-\delta$, calculated from expressions (6.8) ($q=1$) and (6.6) or (6.7) ($q=2$) and rounded down to the nearest $q$-th power. $\vphantom{\int^{t}}q=1$, $\braket{N}=10^{6}$, $\braket{M}=100$ |  $q=2$, $\braket{N}=10^{7}$, $\braket{M}=50$
---|---
$\vphantom{\int^{t}}\delta$ |  $L$ (6.8) | $\delta$ |  $L$ (6.6) |  $L$ (6.7)
$0.1$ | $68{,}\mbox{$\>\\!$}555$ | $0.1$ | $1{,}\mbox{$\>\\!$}890{,}\mbox{$\>\\!$}625=1375^{2}$ | $1{,}\mbox{$\>\\!$}962{,}\mbox{$\>\\!$}801=1401^{2}\vphantom{\int^{t}}$
$0.01$ | $92{,}\mbox{$\>\\!$}053$ | $0.01$ | $2{,}\mbox{$\>\\!$}762{,}\mbox{$\>\\!$}244=1662^{2}$ | $2{,}\mbox{$\>\\!$}900{,}\mbox{$\>\\!$}209=1703^{2}\vphantom{\int^{t}}$
$0.001$ | $115{,}\mbox{$\>\\!$}124$ | $0.001$ | $3{,}\mbox{$\>\\!$}636{,}\mbox{$\>\\!$}649=1907^{2}$ | $3{,}\mbox{$\>\\!$}825{,}\mbox{$\>\\!$}936=1956^{2}\vphantom{\int^{t}}$
$0.0001$ | $138{,}\mbox{$\>\\!$}154$ | $0.0001$ | $4{,}\mbox{$\>\\!$}515{,}\mbox{$\>\\!$}625=2125^{2}$ | $4{,}\mbox{$\>\\!$}743{,}\mbox{$\>\\!$}684=2178^{2}\vphantom{\int^{t}}$

According to Lemma 2.3 (with
$\tilde{\varLambda}^{\dagger}=\check{\varLambda}_{L}$), the output of a
truncated sampling algorithm follows the Boltzmann distribution with a smaller
source set $\mathbb{A}_{L}=\\{\ell\in\mathbb{N}^{q}\colon\ell\leq L\\}$, which
nonetheless approximates well the target Boltzmann distribution
${\mathsf{P}\mbox{$\>\\!\\!$}}_{\bm{z}}$ (see Lemma 2.6). One should be wary
though that truncation contributes to the negative bias (see (2.31) and
(2.32)), which may require a refined calibration through the parameters
$z_{1}$ and $z_{2}$.

Note that the confidence guarantee $1-\delta$ as discussed above is valid only
in the case of a single output instance. If the purpose of the free algorithm
is to produce an independent sample of, say, $k$ random Boltzmann partitions,
then the overall confidence probability is approximately given by
$(1-\delta)^{k}$, which may be exponentially small if $k$ is large while
$\delta$ stays fixed. A simple upper bound for the error probability is based
on the Bernoulli inequality, yielding $1-(1-\delta)^{k}\leq k\delta$. This
motivates the well-known Bonferroni correction, which amounts to choosing the
individual error probability $\delta_{0}=\delta/k$ in order to ensure the
overall error probability not exceeding $\delta$. As an example, if
$k=1{}\mbox{$\>\\!$}000$ and we would like to guarantee the overall error
probability bound $\delta=0.1$, then the individual error probability should
be taken as $\delta_{0}=0.0001$. The approximation is quite accurate here, as
the exact solution is $\delta_{0}\doteq 0.000105355$. Clearly, switching from
$\delta$ to $\delta_{0}$ leads to a higher threshold $L$. For instance, Table
2 shows that in the case $q=1$ the suitable threshold $L$ needs to double. In
general, the increase of $L$ due to multiple errors is not really dramatic
because of the logarithmic dependence on $\delta$ mentioned above.

Another unwelcome outcome of the Bernoulli testing is that it may not return
any parts at all, which has a positive
${\mathsf{P}\mbox{$\>\\!\\!$}}_{\bm{z}}$-probability even in the infinite
sequence of tests (see Remark 3.3). This is not critical, as the sampling
cycle can be repeated if necessary. However, it may be wasteful and can be
easily rectified by adopting a similar approach based on confidence.
Specifically, one can set a lower cutoff $L_{0}$ such that the run of the
sampler is terminated, and the cycle is repeated, if the Bernoulli tests fail
to select at least one (non-zero) part $\ell\leq L_{0}$. To this end, we
choose $L_{0}$ in such a way that
${\mathsf{P}\mbox{$\>\\!\\!$}}_{\bm{z}}(\lambda_{\rm min}>L_{0})\leq\delta$,
where $\delta>0$ is small enough (cf. (6.5)). Again referring to the limit
theorems regarding the smallest part, we obtain from Theorem 4.4 (cf. (6.6))

$\Gamma(1/q)-\Gamma(1/q,\gamma_{0}L_{0})=\frac{\Gamma(1/q)}{\braket{M}}\mbox{$\;\\!$}\log\frac{1}{\delta},$
(6.9)

and from Theorem 5.7 (cf. (6.7))

$L_{0}=\frac{1}{\gamma_{0}}\left(\frac{\Gamma(1/q)}{q\mbox{$\;\\!\\!$}\braket{M}}\log\frac{1}{\delta}\right)^{q}\\!.$
(6.10)

A numerical illustration of the confident lower threshold $L_{0}$ is presented
in Table 3 for $q=1$ and $q=2$. The equation (6.9) was solved numerically
using Maple. The match between the results produced via equations (6.9) and
(6.10) is quite close, especially for $q=2$. One should also observe a
significant difference between the thresholds $L_{0}$ and $L$, which underpins
a considerable computational saving due to the lower cutoff, activated
whenever the sampler fails to produce at least one positive part up to
$L_{0}$.

Table 3: Asymptotic threshold $L_{0}$ for the smallest part $\lambda_{\rm{min}}$ with confidence probability $1-\delta$, calculated from expressions (6.9) or (6.10), and rounded down to the nearest $q$-th power. $\vphantom{\int^{t}}q=1$, $\braket{N}=10^{6}$, $\braket{M}=100$ |  $q=2$, $\braket{N}=10^{7}$, $\braket{M}=50$
---|---
$\vphantom{\int^{t}}\delta$ |  $L_{0}$ (6.9) |  $L_{0}$ (6.10) | $\delta$ |  $L_{0}$ (6.9) |  $L_{0}$ (6.10)
$0.1$ | $232$ | $230$ | $0.1$ | $625=25^{2}$ | $625=25^{2}\vphantom{\int^{t}}$
$0.01$ | $471$ | $460$ | $0.01$ | $2{,}\mbox{$\>\\!$}601=51^{2}$ | $2{,}\mbox{$\>\\!$}601=51^{2}\vphantom{\int^{t}}$
$0.001$ | $715$ | $690$ | $0.001$ | $5{,}\mbox{$\>\\!$}929=77^{2}$ | $5{,}\mbox{$\>\\!$}929=77^{2}\vphantom{\int^{t}}$
$0.0001$ | $966$ | $921$ | $0.0001$ | $10{,}\mbox{$\>\\!$}816=104^{2}$ | $10{,}\mbox{$\>\\!$}609=103^{2}\vphantom{\int^{t}}$

Finally, if both cutoffs $L$ and $L_{0}$ are exercised as described above
then, due to the asymptotic independence of $\lambda_{\rm max}$ and
$\lambda_{\rm min}$, the overall confidence probability is (asymptotically)
given by $(1-\delta)^{2}=1-2\mbox{$\>\\!$}\delta+\delta^{2}$, hence the
resulting error probability is bounded by $2\mbox{$\>\\!$}\delta$.

####  6.1.3 Free sampling algorithm

A free Boltzmann sampler is presented below in pseudocode as Algorithm 1. For
simplicity, the algorithm incorporates only the upper threshold $L$ selected
in advance for a given confidence probability $1-\delta$. As discussed in
Section 6.1.2, this is essential to ensure termination of the code, but for
the sake of optimization a lower cutoff $L_{0}$ can also be included without
difficulty. As explained above, the confidence probability should be chosen
carefully to match a possibly multiple output.

Input: integer $q$, real $\langle N\rangle,\langle M\rangle,L$

Output: partition $\lambda\in\check{\varLambda}^{q}_{L}$, weight
$N_{\lambda}$, length $M_{\lambda}$

1 integer array $\lambda_{[\ ]}$;

2 real $z_{1},z_{2},\gamma_{0},q\mbox{$\>\\!$}$;

3 $\gamma_{0}\leftarrow\langle M\rangle/(q\mbox{$\;\\!\\!$}\braket{N})$;

4 $z_{1}\leftarrow\mathrm{e}^{-\gamma_{0}}$, $z_{2}\leftarrow
q\mbox{$\;\\!\\!$}\braket{M}\mbox{$\;\\!\\!$}\smash{\gamma_{0}^{1/q}}\mbox{$\>\\!\\!$}/\mbox{$\>\\!$}\Gamma(1/q)$;

5 integer $j^{*}$, $j$, $N$, $M$;

6 $j^{*}\leftarrow\lfloor L^{1/q}\rfloor$;

7 $N\leftarrow 0$, $M\leftarrow 0$;

8 for  _$j$ from $j^{*}$ to $1$ by $-1$_ do

9 $p\leftarrow
z_{1}^{j^{q}}\\!z_{2}\mbox{$\;\\!$}(1+z_{1}^{j^{q}}\\!z_{2}){{}^{-1}}$;

10 if  _$\mathtt{Ber}(p)=1$_ then

11 $N\leftarrow N+j^{q}$;

12 $M\leftarrow M+1$;

13 $\lambda_{M}\leftarrow j^{q}$;

14

15 end if

16

17 end for

18$N_{\lambda}\leftarrow N$, $M_{\lambda}\leftarrow M$;

return _$(\lambda,N_{\lambda},M_{\lambda})$_

Algorithm 1 FreeSampler ($q,\langle N\rangle,\langle M\rangle,L$)

The code structure is fairly straightforward and consists in a single cycle of
sequential Bernoulli tests over potential parts $\ell\in\mathbb{N}^{q}$. It is
convenient to do this via downward scoping in view of our convention to
enumerate the partition parts in decreasing order. Because the resulting
length of the output partition $\lambda=(\lambda_{i})$ is unknown in advance,
the space for the corresponding integer array is defined in the code as
$\lambda_{[\ ]}$, that is, through a dynamically allocated memory. Finally,
the calibration parameters $z_{1}$ and $z_{2}$ are specified using the
leading-term formulas (6.1); if desired, these can be replaced by the bias-
correcting values (6.3) or by any other, more refined choices.

By design, the output of Algorithm 1 is a random partition
$\lambda\in\check{\varLambda}^{q}_{L}=\check{\varLambda}^{q}\cap\\{(\lambda_{i})\colon\lambda_{\rm{max}}\leq
L\\}$. It has a Boltzmann distribution on the space
$\check{\varLambda}^{q}_{L}$, with expected values of weight $N_{\lambda}$ and
length $M_{\lambda}$ close to the predefined hyper-parameters $\braket{N}$ and
$\braket{M}$, respectively. As discussed in Section 2.5, this distribution
approximates (in total variation) the Boltzmann distribution on the infinite
partition space $\check{\varLambda}^{q}$, which may suffice for the sampling
purposes at hand.

####  6.1.4 Validation of Algorithm 1

The output performance of the code in Algorithm 1 was visually monitored via
the marginal histograms for the sample weight $N_{\lambda}$ and length
$M_{\lambda}$ (Figure 4), as well as by the bivariate histograms and frequency
level plots of the sample pairs $(N_{\lambda},M_{\lambda})$ (Figure 5). The
numerical illustration was carried out in the case of square parts, $q=2$
(selected for computational convenience in order to reduce the completion
time), and in two different regimes with regard to the hyper-parameter
$\braket{M}$, that is, “fixed” and “slow growth”, illustrated by
$\braket{M}=5$ ($\braket{N}=12{,}\mbox{$\>\\!$}500$) and $\braket{M}=50$
($\braket{N}=10^{7}$), yielding for the parameter $\kappa=\langle
M\rangle^{3}\mbox{$\;\\!\\!$}/\langle N\rangle$ values $\kappa=0.01$ and
$\kappa=0.0125$, respectively (cf. Assumption 5.1). The algorithm was run at a
very low confidence tolerance (error probability) $\delta=10^{-8}$ and with
the corresponding truncation value $L$ calculated using formulas (6.6) or
(6.7) according to the regime at hand, yielding $L=299^{2}=89,401$ and
$L=2{,}\mbox{$\>\\!$}903^{2}=8{,}\mbox{$\>\\!$}427{,}\mbox{$\>\\!$}409$,
respectively (cf. Table 2).

(a) $\braket{M}=5$, $\braket{N}=12{,}\mbox{$\>\\!$}500$ ($\kappa=0.01$). A dot
at the origin on the left panel shows a discrete atom
$\pi_{0}=\mathrm{e}^{-\braket{M}}\doteq 0.00674$ (see (4.15)).

(b) $\braket{M}=50$, $\braket{N}=10^{7}$ ($\kappa=0.0125$). The atom
$\pi_{0}=\mathrm{e}^{-\braket{M}}\doteq 2\cdot 10^{-22}$ is too small to be
visible.

Figure 4: Marginal histograms for the weight $N_{\lambda}$ (left) and length
$M_{\lambda}$ (right) for random samples (of size $10^{5}$ each) from the
partition space $\check{\varLambda}^{q}$ ($q=2$), simulated using a free
Boltzmann sampler as set out in Algorithm 1. Colour coding (online version):
_blue_ designates the limiting distributions under the “fixed” regime, that
is, compound Poisson-Gamma (left) and Poisson (right); _red_ indicates a
normal approximation; _magenta_ depicts a mean-gamma approximation (see
Section 4.1). In the black-and-white version, the normal curves on the left
are identifiable by a noticeable positive shift.

The empirical results (with $10^{5}$ output partitions in both cases) were
compared with the theoretical predictions from Theorems 4.1 and 5.1. The
marginal histogram plots for $N_{\lambda}$ and $M_{\lambda}$ shown in Figure 4
depict bell-shaped unimodal empirical distributions with the sample modes
noticeably shifted to the left of the calibration hyper-parameters
$\braket{N}$ and $\braket{M}$, respectively. The discrepancy between the modes
and the means is observed especially well in the weight plots (the more so for
smaller $\braket{M}$), which is essentially due to the fact that the
underlying gamma and Poisson distributions are right-skewed; for example, the
mode of $\mathrm{Gamma}\mbox{$\;\\!$}(\alpha)$ is given by
$\max\\{\alpha-1,0\\}$ whereas the mean is $\alpha$. Another (minor) reason
for a negative bias is due to a certain miscalibration, as was pointed out in
Section 6.1.1. In line with theoretical predictions, this mismatch is
vanishing with the growth of the hyper-parameters $\braket{N}$ and
$\braket{M}$, together with the improving accuracy of the normal
approximations, both for $N_{\lambda}$ and $M_{\lambda}$. It is also
interesting to note that the mean-gamma approximation
$\mathrm{Gamma}\mbox{$\;\\!$}(\braket{M})$ (see Section 4.1) nearly perfectly
matches the exact compound Poisson-Gamma distribution for $\braket{M}=50$ (see
Figure 4(b), left); for smaller values of $\braket{M}$, this approximates is
rather crude, however it still captures well the mode of the Poisson-Gamma
distribution and also its right shoulder (see Figure 4(a), left).

A remarkable exception to the unimodality of the plots in Figure 4 is the
compound Poisson-Gamma plot for the weight $N_{\lambda}$, with a relatively
small value of $\braket{M}=5$ (see Figure 4(a), left), where one can clearly
see a singularity at zero (cf. Section 4.1). This theoretical prediction is
supported by the empirical results, with an obvious excess of smaller weights.
With $\braket{M}=5$ and $\braket{N}=12{,}\mbox{$\>\\!$}500$, the local minimum
of the theoretical density $g(x)$ defined in (4.17) occurs at $x_{0}\doteq
0.10340$ with value111111Note that the asymptotic formula (4.18) gives a
pretty accurate approximation $g(x_{0})\approx 0.14334$. $g(x_{0})\doteq
0.19632$, which corresponds to weight $n_{0}=\lceil
x_{0}/\gamma_{0}\rceil=517$. If the density $g(x)$ continued to decay to the
left of $x_{0}$, this would predict the (asymptotic) probability of getting
weights smaller than $n_{0}$ (together with an empty partition) loosely
bounded by $x_{0}\mbox{$\;\\!$}g(x_{0})+\pi_{0}\doteq 0.02704$. But the actual
compound Poisson-Gamma probability is higher, $G(x_{0})\doteq 0.03120$ (see
(4.15)). The excess of “small” partitions is reminiscent of a partition
interpretation of the Bose–Einstein condensation (see [79]). As already
mentioned in Section 4.1, this is a truly finite-length phenomenon, which
vanishes as $\braket{M}\to\infty$ (cf. Figure 4(b), where $\braket{M}=50$).

(a) Bivariate frequency plot. (b) Frequency level sets.

Figure 5: Joint sampling distribution of weight $N_{\lambda}$ and length
$M_{\lambda}$ for $q=2$, $\braket{M}=50$ and $\braket{N}=10^{7}$ (cf. marginal
plots in Figure 4(b)). A random Boltzmann sample of partitions
$\lambda\in\check{\varLambda}^{q}$ (of size $10^{5}$) was simulated using
Algorithm 1.

The bivariate plots in Figure 5 (for $\braket{M}=50$) appear to be
approximately consistent with the asymptotically predicted (standardized)
confidence ellipses of the form

$\mathcal{L}_{\alpha}=\\{\bm{x}\in\mathbb{R}^{2}\colon\bm{x}\mbox{$\>\\!$}\bm{K}_{q}^{-1}\bm{x}^{\top}\leq\chi^{2}_{2}(1-\alpha)\\},$
(6.11)

where $\bm{K}_{q}^{-1}$ is the inverse covariance matrix (5.2), and
$\chi^{2}_{2}(1-\alpha)$ is the quantile of the chi-squared distribution with
two degrees of freedom, corresponding to confidence probability $1-\alpha$.
The latter distribution simplifies to an exponential distribution with mean
$2$, hence $\chi^{2}_{2}(1-\alpha)=2\log\left(1/\alpha\right)$. According to
Theorem 5.1, a sample point $(N_{\lambda}^{*},M_{\lambda}^{*})$ belongs to the
ellipse (6.11) approximately with probability $1-\alpha$, where the
standardized values $N_{\lambda}^{*},M_{\lambda}^{*}$ are defined in (5.1).
The inverse of $\bm{K}_{q}$ is easily computed,

$\bm{K}_{q}^{-1}=\frac{q+1}{q}\begin{pmatrix}1&\displaystyle\frac{-1}{\sqrt{q+1}}\\\
\displaystyle\frac{-1}{\sqrt{q+1}}&1\end{pmatrix}\\!,$

and the confidence ellipse (6.11) specializes as follows,

$N_{\lambda}^{*2}-\frac{2\mbox{$\>\\!$}N_{\lambda}^{*}M_{\lambda}^{*}}{\sqrt{q+1}}+M_{\lambda}^{*2}\leq\frac{2\mbox{$\>\\!$}q}{q+1}\log\frac{1}{\alpha}\mbox{$\>\\!$}.$

A closer inspection of the level plots in Figure 5(b) reveals some elongation
of the frequency level sets towards bigger values of the weight $N_{\lambda}$,
thus indicating a bit of discrepancy with the predicted elliptical shape. This
observation is confirmed by comparison of the marginal histograms of
$N_{\lambda}$ and $M_{\lambda}$ in Figure 4, where the latter is reasonably
symmetric while the former is noticeably skewed to the right. A heuristic
explanation of such an effect may be based on noticing from formulas (2.7)
that, while the length $M_{\lambda}$ is built by summation of multiplicities
$\nu_{\ell}$, the weight $N_{\lambda}$ involves size-biased terms
$\ell\mbox{$\>\\!$}\nu_{\ell}$, which pinpoint skewing the distribution to the
right.

A well-localized unimodal nature of the distributions behind the outputs
$N_{\lambda}$ and $M_{\lambda}$ is sometimes referred to as a _bumpy type_
[23], characterized by an asymptotically large signal-to-noise ratio (SNR) in
response to a large signal,121212The “bumpy” property is especially helpful
for sampling with rejection designed to achieve certain targets for weight and
length, due to a guaranteed asymptotically fast delivery of the output (see
Section 6.2).

$\mathrm{SNR}(X):=\frac{[{\mathsf{E}}(X)]^{2}}{{\mathsf{Var}}(X)}\to\infty,\qquad{\mathsf{E}}(X)\to\infty.$
(6.12)

Here, the notation $X$ designates a random output in question (such as the
size), which has a large expected value. Following definition (6.12) and
applying Theorem 3.9, we readily get

$\mathrm{SNR}(N_{\lambda})\sim\frac{\braket{N}^{2}}{(q+1)\braket{N}^{2}\\!/\\!\braket{M}}=\frac{\braket{M}}{q+1},\qquad\mathrm{SNR}(M_{\lambda})\sim\frac{\braket{M}^{2}}{\braket{M}}=\braket{M},$

so that under Assumption 5.1 (slow growth of $\braket{M}$) each of the
marginal SNRs tends to infinity.

In the multivariate case, the SNR is usually defined in the literature as a
scalar value,

$\mathrm{SNR}(\bm{X}):=\bm{\mu}\mbox{$\;\\!$}\bm{K}^{-1}\mbox{$\;\\!\\!$}\bm{\mu}^{\top},\qquad\bm{\mu}:={\mathsf{E}}(\bm{X}),\
\ \bm{K}:={\mathsf{Cov}}(\bm{X},\bm{X})$ (6.13)

(see, e.g., [65, Eq. (1), p. 511]). Again using Theorem 3.9, we find the
asymptotic inverse of the covariance matrix,

$\bm{K}^{-1}(\bm{z})\sim\frac{1}{q\mbox{$\;\\!\\!$}\braket{N}^{2}}\begin{pmatrix}\braket{M}&-\braket{N}\\\
-\braket{N}&\displaystyle\frac{(q+1)\mbox{$\>\\!\\!$}\braket{N}^{2}}{\braket{M}}\end{pmatrix}\\!,$

and hence

$\mathrm{SNR}(N_{\lambda},M_{\lambda})\sim\frac{(\langle N\rangle,\langle
M\rangle)}{q\mbox{$\;\\!\\!$}\braket{N}^{2}}\begin{pmatrix}\braket{M}&-\braket{N}\\\
-\braket{N}&\displaystyle\frac{(q+1)\mbox{$\>\\!\\!$}\braket{N}^{2}}{\braket{M}}\end{pmatrix}\begin{pmatrix}\braket{N}\\\\[4.79993pt]
\braket{M}\end{pmatrix}=\braket{M}\to\infty.$

However, a scalar definition (6.13) is not entirely satisfactory — for
instance, it cannot detect whether the individual components of $\bm{X}$ are
of bumpy type. As an alternative, we propose the following matrix definition,

$\mathbf{SNR}(\bm{X}):=(\bm{\mu}\mbox{$\>\\!$}\bm{K}^{-1/2})^{\mbox{$\>\\!\\!$}\top}(\bm{\mu}\mbox{$\>\\!$}\bm{K}^{-1/2})=\bm{K}^{-1/2}\mbox{$\>\\!$}(\bm{\mu}^{\mbox{$\;\\!\\!$}\top}\mbox{$\>\\!\\!$}\bm{\mu})\mbox{$\;\\!$}\bm{K}^{-1/2},$
(6.14)

where $\bm{K}^{-1/2}$ is the (unique) positive definite square root of the
inverse covariance matrix $\bm{K}^{-1}$, that is,
$\bm{K}^{-1/2}\bm{K}^{-1/2}=\bm{K}^{-1}$ [43, Theorem 7.2.6, p. 439]. In our
case, the exact expression for $\bm{K}^{-1/2}(\bm{z})$ is cumbersome (although
available), but its asymptotic version under Assumption 5.1 simplifies to

$\bm{K}^{-1/2}(\bm{z})\sim\frac{\sqrt{\braket{M}}}{\sqrt{q\mbox{$\;\\!$}(q+1)}\braket{N}}\begin{pmatrix}\sqrt{q}&-1\\\\[4.79993pt]
-1&\displaystyle\frac{(q+1)\braket{N}}{\braket{M}}\end{pmatrix}\\!.$

Hence, after straightforward calculations we obtain from (6.14)

$\mathbf{SNR}(N_{\lambda},M_{\lambda})\sim\frac{\braket{M}}{q+1}\begin{pmatrix}1&\sqrt{q}\\\\[4.79993pt]
\sqrt{q}&q\end{pmatrix}\\!,$

which tends to infinity in matrix sense.

###  6.2 Rejection sampler

The idea of a rejection sampler discussed in this section is to run a free
sampler in a loop until a prescribed target is met. For example, if the target
is set in terms of the required partition length as $M_{\lambda}=m$, with a
fixed $m\in\mathbb{N}$, then the free sampler is iterated until a partition of
exact length $m$ is obtained. Likewise, if the target is set for the partition
weight, $N_{\lambda}=n$, with a fixed $n\in\mathbb{N}$, then the free sampling
loop runs until a partition of exact weight $n$ is found. These two targets
can be imposed simultaneously, $M_{\lambda}=m$ and $N_{\lambda}=n$; here, it
is natural to design the rejection algorithm as a juxtaposition of two loops
of the free sampler, such that the internal loop runs until the length target
is met and, every time this happens, the resulting partition is checked with
regards to the weight target and is either rejected, whereby the internal loop
starts afresh, or accepted, in which case the algorithm stops.

A more general approach, leading to the so-called approximate algorithms, is
to relax the exact targets to suitable intervals (brackets). In Algorithm 2
presented below in Section 6.2.3, we give an example of a Boltzmann rejection
sampler aiming to sample a partition $\lambda\in\check{\varLambda}^{q}$
satisfying two conditions, $M_{\lambda}=m$ and $n\leq
N_{\lambda}\leq\theta\kern 0.35999ptn$, for some predefined tolerance factor
$\theta\geq 1$. Of course, if $\theta=1$, the approximate algorithm is reduced
to an exact one. An approximate target can also be considered for length,
$m\leq M_{\lambda}\leq\theta^{\prime}m$, and furthermore, such approximate
targets can be combined if desired.

Before delineating Algorithm 2, we discuss a few implementation issues arising
therein.

####  6.2.1 Parameter calibration and truncation of parts

To start with, the choice of the calibrating parameters $z_{1}$ and $z_{2}$
now follows a slightly different logic as compared to the case of a free
sampler. If only one target is set, such as $M_{\lambda}=m$, then according to
formula (2.25) the conditional Boltzmann measure
${\mathsf{P}\mbox{$\>\\!\\!$}}_{\bm{z}}(\cdot\,|\mbox{$\;\\!$}\check{\varLambda}^{q}(\mathbin{\vbox{\hbox{\scalebox{0.5}{$\bullet$}}}}\mbox{$\>\\!$},m))$
does not depend on the parameter $z_{2}$, whereas the other parameter,
$z_{1}$, can be chosen with a view on a desired expected value $\braket{N}$ of
the output weight $N_{\lambda}$, as was the case with the free sampler.
Nonetheless, in order to maximize efficiency of the sampling algorithm, the
“free” parameter $z_{2}$ should still be chosen in line with the mean
conditions (3.21) subject to the specification $\braket{M}=m$, thus aiming to
benefit from the bumpy nature of the distribution
${\mathsf{P}\mbox{$\>\\!\\!$}}_{\bm{z}}(\cdot\,|\mbox{$\;\\!$}\check{\varLambda}^{q}(\mathbin{\vbox{\hbox{\scalebox{0.5}{$\bullet$}}}}\mbox{$\>\\!$},m))$
(cf. Section 6.1.3). Similarly, if the condition $N_{\lambda}=n$ is targeted
then the conditional Boltzmann measure
${\mathsf{P}\mbox{$\>\\!\\!$}}_{\bm{z}}(\cdot\,|\mbox{$\;\\!$}\check{\varLambda}^{q}(n,\mbox{$\;\\!\\!$}\mathbin{\vbox{\hbox{\scalebox{0.5}{$\bullet$}}}}))$
does not depend on $z_{1}$, however both parameters $z_{1}$ and $z_{2}$ are
chosen to match the mean conditions (3.21) with $\braket{N}=n$. Furthermore,
if both targets are imposed, $N_{\lambda}=n$ and $M_{\lambda}=m$, then by
Lemma 2.4 the measure
${\mathsf{P}\mbox{$\>\\!\\!$}}_{\bm{z}}(\cdot\,|\mbox{$\>\\!$}\check{\varLambda}^{q}(n,m))$
is reduced to the uniform distribution on $\check{\varLambda}^{q}(n,m)$
regardless of the parameters $z_{1}$ and $z_{2}$. But, as explained above, it
is worthwhile to calibrate them in line with the mean conditions (3.21). With
these targets in mind, in what follows (including Algorithm 2) we use the
specific hyper-parameters $\braket{N}=n$ and $\braket{M}=m$. Moreover, since
avoiding bias is no longer a concern (unlike in the case of the free sampler),
using the leading-term expressions (6.1) is a perfectly satisfactory option.
Of course, the same recipe applies to the approximate sampling.

Another issue to be addressed is whether any truncation of the source of parts
is needed in the algorithm (cf. Section 6.1.2). As long as the weight target
is involved, $N_{\lambda}=n$, one can use a natural majorant $L^{*}=n$ as a
call parameter $L$ in Algorithm 1 (see Section 6.1.3), which clearly causes no
loss in confidence (i.e., $\delta=0$, see (6.5)). The same is true in the case
of an approximate target, $N_{\lambda}\in[n,\theta\kern 0.47992ptn]$, by
choosing the majorant $L^{*}=\theta\kern 0.71997ptn$. However, if only the
length target is in place then no such majorant is available and, therefore,
confidence considerations must be deployed, as discussed in Section 6.1.2.

####  6.2.2 Censoring of iterations

As was pointed out in the Introduction, in contrast to the special value
$q=1$, in the general case with $q\geq 2$ and arbitrary $m\geq 1$ there is no
guarantee for a given natural number $n\in\mathbb{N}$ to be partitionable into
a required number $m$ of $q$-power parts (unless it is covered by a solution
of the Waring problem [76]). The requirement that the parts be distinct adds
to the complexity of the question. Therefore, the space
$\check{\varLambda}^{q}(n,m)$ may well be empty and, not knowing this in
advance, the task of sampling from such a space may be “mission impossible”.

To be specific, consider sampling subject to the joint targets $M_{\lambda}=m$
and $N_{\lambda}=n$. As already indicated at the beginning of Section 6, a
general design of the corresponding sampling algorithm is based on the two
nested loops according to the separated targets, internal for $M_{\lambda}$
and external for $N_{\lambda}$. While the internal loop is certain to produce
a random partition in $\check{\varLambda}^{q}$ with exactly $m$ parts (see
more about this below), the external loop contains an inherent loose end due
to its potential failure to satisfy the weight requirement, simply because
there may be no such partitions. That is to say, although the successful
completion of the external loop will take some time by repeatedly querying the
internal loop, it may be pointless to wait for too long as there is no
certainty if that is not wasteful.

We propose to resolve this difficulty by an appropriate “censoring” of
processing time, that is, by setting a limit $t^{*}$ on the waiting time,
chosen to ensure sufficiently high confidence in the algorithm’s ability to
deliver a successful completion within the allocated time limit, of course
provided that the task is feasible (i.e., required partitions exist). More
precisely, given a confidence tolerance (significance level) $\delta\in(0,1)$,
the threshold $t^{*}$ should be such that, if the target space is non-empty,
the probability that the algorithm would not succeed by time $t^{*}$ does not
exceed $\delta$. Such an approach is akin to statistical testing of the null
hypothesis “the target is non-empty”. Under this hypothesis, the test
(implemented as a sampling algorithm) fails to produce a required partition
(Type I error) with rate bounded by $\delta$.

The choice of threshold $t^{*}$ is determined by the sampling task at hand. A
few examples of interest are as follows (assuming the number of parts $m$ to
be fixed):

  * (T1)

_Exact sampling:_ For a given $n$, attempt to sample
$\lambda\in\check{\varLambda}^{q}(n,m)$ (in other words, check
partitionability of $n$).

  * (T2)

_Multiple exact sampling:_ For a given $n$ and some $\theta>1$, attempt to
sample $\lambda\in\check{\varLambda}^{q}(k,m)$ for each integer $k$ in the
range $k\in[n,\theta\kern 0.24005ptn]$ (that is, test partitionability of each
of these numbers).

  * (T3)

_Approximate sampling:_ Same as in task (T2) but attempting to sample
$\lambda\in\bigcup_{n\leq k\leq\theta\kern
0.24005ptn}\check{\varLambda}^{q}(k,m)$ (that is, to find at least one
partitionable number in the said range).

######  Remark 6.1.

If required, tasks (T2) and (T3) could be modified to a two-sided version,
such as $\theta^{-1}n\leq k\leq\theta\kern 0.35999ptn$ or, more generally,
$\theta_{1}n\leq k\leq\theta_{2}\mbox{$\>\\!$}n$, with some
$0<\theta_{1}<1<\theta_{2}$.

First, let us look at how the internal loop performs towards its task of
sampling a partition
$\lambda\in\check{\varLambda}^{q}_{L}(\mathbin{\vbox{\hbox{\scalebox{0.5}{$\bullet$}}}}\mbox{$\>\\!$},m)=\check{\varLambda}^{q}(\mathbin{\vbox{\hbox{\scalebox{0.5}{$\bullet$}}}}\mbox{$\>\\!$},m)\cap\\{\lambda\colon\lambda_{\rm
max}\leq L\\}$ (i.e., the source of parts is truncated to $\\{\ell\leq L\\}$,
see Section 2.5). According to Theorem 4.2(a), if $\langle M\rangle>0$ is
fixed and $L\sim\theta\mbox{$\>\\!$}\langle N\rangle$ then the distribution of
$M_{\lambda}$ conditional on $\lambda_{\rm max}\leq L$ converges to a Poisson
law with mean $\mu_{\theta}=\braket{M}G_{1/q}(a_{\theta})$, where
$a_{\theta}=\theta\mbox{$\>\\!$}\langle M\rangle\mbox{$\;\\!\\!$}/q$. A
stronger result concerns a Poisson approximation (cf. Remark 4.1) under a
suitable metric, such as the total variation distance between distributions.
Namely, the said conditional distribution of $M_{\lambda}$ can be replaced by
a Poisson distribution with mean

${\mathsf{E}}_{\bm{z}}(M_{\lambda}\,|\,\lambda_{\rm{max}}\leq
L)=\sum_{\ell\leq
L}\frac{z_{1}^{\ell}z_{2}}{1+z_{1}^{\ell}z_{2}}\sim\mu_{\theta},$ (6.15)

with the error in total variation bounded by (see [8, Theorem 1, p. 474])

$\frac{1}{\mu_{\theta}}\sum_{\ell\leq
L}\left(\frac{z_{1}^{\ell}z_{2}}{1+z_{1}^{\ell}z_{2}}\right)^{\mbox{$\>\\!\\!$}2}=O(z_{2})=O(\kappa^{1/q})\to
0,$ (6.16)

according to (4.6). The advantage of such an approximation is that it holds
true even if $\braket{M}$ is slowly growing, whereby the error estimate (6.16)
is still valid.

Returning to the analysis of the internal loop, with $\braket{N}=n$ and
$\braket{M}=m$ we have

${\mathsf{P}\mbox{$\>\\!\\!$}}_{\bm{z}}(M_{\lambda}=m\,|\,\lambda_{\rm
max}\leq\theta\kern
0.35999ptn)\sim\frac{\mu_{\theta}^{m}\,\mathrm{e}^{-\mu_{\theta}}}{m!},\qquad\mu_{\theta}=m\,G_{1/q}(\theta\mbox{$\>\\!$}m/q).$
(6.17)

Hence, the probability (6.17) is bounded away from zero, and since the
attempts within the internal loop are independent, the number of internal runs
until success has geometric distribution, with the expected time to success
being bounded by a constant (depending on $m$). If $m\to\infty$ (with
$\kappa=m^{q+1}\mbox{$\>\\!\\!$}/n=o(1)$), then $\mu_{\theta}\sim m$ and,
according to (6.15) and (6.16), we have

${\mathsf{P}\mbox{$\>\\!\\!$}}_{\bm{z}}(M_{\lambda}=m\,|\,\lambda_{\rm
max}\leq\theta\kern
0.35999ptn)=\frac{\mu_{\theta}^{m}\,\mathrm{e}^{-\mu_{\theta}}}{m!}+O(\kappa^{1/q})\sim\frac{m^{m}\,\mathrm{e}^{-m}}{m!}\sim\frac{1}{\sqrt{2\pi\mbox{$\>\\!$}m}},$
(6.18)

with the help of the Stirling formula. In turn, formula (6.18) implies that
the expected number of runs of the internal loop is of order
$O\bigl{(}\sqrt{m}\mbox{$\>\\!$}\bigr{)}$, which is not particularly large for
practical implementation.

Let us now turn to tasks (T1) – (T3) and focus on probabilistic analysis of
runs of the external loop, taking for granted that $M_{\lambda}=m$ and,
automatically, $\lambda_{\rm max}\leq L^{*}$, where $L^{*}=\theta\kern
0.35999ptn$ is a majorant in $\bigcup_{n\leq k\leq\theta
n}\check{\varLambda}^{q}(k,m)$. As stipulated above, we use the hyper-
parameters $\braket{N}=n$ and $\braket{M}=m$, and specify the calibrating
parameters $z_{1}$ and $z_{2}$ according to the leading-term expressions
(6.1). Denote for short

$p_{k}^{*}:={\mathsf{P}\mbox{$\>\\!\\!$}}_{\bm{z}}(N_{\lambda}=k\,|\mbox{$\;\\!$}M_{\lambda}=m,\mbox{$\>\\!$}\lambda_{\rm
max}\leq\theta\kern 0.35999ptn),\qquad n\leq k\leq\theta\kern 0.35999ptn.$
(6.19)

for simplicity omitting reference to $n$ and $m$. Of course, if
$\check{\varLambda}^{q}(k,m)=\varnothing$ then $p_{k}^{*}=0$. The probability
(6.19) can be interpreted as the probability of successfully sampling a
partition of weight $k$ in a single run of the external loop. Due to
independence of successive runs, the number $T_{k}$ of attempts until success
for a targeted weight value $k$ follows geometric distribution,

${\mathsf{P}\mbox{$\>\\!\\!$}}_{\bm{z}}(T_{k}>t\,|\mbox{$\;\\!$}M_{\lambda}=m,\lambda_{\rm
max}\leq n)=\left(1-p_{k}^{*}\right)^{t},\qquad t\in\mathbb{N}_{0}.$ (6.20)

This includes the case $p_{k}^{*}=0$, whereby $T_{k}=\infty$
(${\mathsf{P}\mbox{$\>\\!\\!$}}_{\bm{z}}$-a.s.). Also note that $(T_{k})$ are
mutually independent for different $k$.

  * (T1)

Here, $\theta=1$, so $L^{*}=n$. Suppose that
$\check{\varLambda}^{q}(n,m)\neq\varnothing$ and let
$\lambda_{*}\in\check{\varLambda}^{q}(n,m)$, so that $N_{\lambda_{*}}=n$ and
$M_{\lambda_{*}}=m$. Then we can write

${\mathsf{P}\mbox{$\>\\!\\!$}}_{\bm{z}}(\lambda_{*}\mbox{$\;\\!$}|\mbox{$\;\\!$}M_{\lambda}=m,\mbox{$\>\\!$}\lambda_{\rm
max}\leq
n)=\frac{{\mathsf{P}\mbox{$\>\\!\\!$}}_{z}(\lambda_{*}\mbox{$\>\\!$}|\mbox{$\;\\!$}\lambda_{\rm
max}\leq n)}{{\mathsf{P}\mbox{$\>\\!\\!$}}_{z}(M_{\lambda}=m\,|\,\lambda_{\rm
max}\leq n)}.$ (6.21)

Starting with the numerator, we have

${\mathsf{P}\mbox{$\>\\!\\!$}}_{z}(\lambda_{*}\mbox{$\>\\!$}|\mbox{$\;\\!$}\lambda_{\rm
max}\leq n)=z_{1}^{n}z_{2}^{m}\prod_{\ell\leq
n}\frac{1}{1+z_{1}^{\ell}\mbox{$\;\\!\\!$}z_{2}}.$

Substituting formulas (6.1), we obtain

$z_{1}^{n}z_{2}^{m}\sim\mathrm{e}^{-m/q}\,\frac{q^{2m}(m/q)^{m+m/q}}{\bigl{(}\Gamma(1/q)\bigr{)}^{m}\mbox{$\;\\!$}n^{m/q}},$

while Lemma 3.5, with the help of the asymptotic relation (3.35), yields

$\prod_{\ell\leq
n}\frac{1}{1+z_{1}^{\ell}\mbox{$\;\\!\\!$}z_{2}}\sim\exp\bigl{(}-m\,G_{1/q}(m/q)\bigr{)}.$

Furthermore, by Theorem 4.2(a) the denominator in (6.21) is asymptotically
given by

${\mathsf{P}\mbox{$\>\\!\\!$}}_{z}(M_{\lambda}=m\,|\,\lambda_{\rm max}\leq
n)\sim\frac{m^{m}\mbox{$\>\\!$}\bigl{(}G_{1/q}(m/q)\bigr{)}^{m}\exp\bigl{(}-m\,G_{1/q}(m/q)\bigr{)}}{m!}.$
(6.22)

Hence, returning to (6.21) we obtain

$p_{n}^{*}\geq{\mathsf{P}\mbox{$\>\\!\\!$}}_{\bm{z}}(\lambda_{*}\mbox{$\;\\!$}|\mbox{$\;\\!$}M_{\lambda}=m,\lambda_{\rm
max}\leq n)\sim\frac{1}{C_{1}(m,q)\,n^{m/q}},$ (6.23)

where

$C_{1}(m,q):=\frac{1}{m!}\left(\frac{\mathrm{e}^{1/q}\,\Gamma(1/q)\,G_{1/q}(m/q)}{q^{1-1/q}\mbox{$\>\\!$}m^{1/q}}\right)^{\mbox{$\>\\!\\!$}m}\\!.$
(6.24)

Combining (6.20) and (6.23), we have, asymptotically,

${\mathsf{P}\mbox{$\>\\!\\!$}}_{\bm{z}}(T_{n}>t\,|\mbox{$\;\\!$}M_{\lambda}=m,\lambda_{\rm
max}\leq
n)=\left(1-p_{n}^{*}\right)^{t}\leq\left(1-\frac{1+o(1)}{C_{1}(m,q)\,n^{m/q}}\right)^{\mbox{$\;\\!\\!$}t}\\!.$
(6.25)

Thus, for the probability (6.25) not to exceed a predefined (small) confidence
tolerance $\delta>0$, it suffices to choose the threshold $t=t_{n}^{*}$ as
follows,

$t_{n}^{*}\simeq\frac{\displaystyle\log\delta}{\displaystyle\log\left(1-\frac{1+o(1)}{C_{1}(m,q)\,n^{m/q}}\right)}\sim
C_{1}(m,q)\,n^{m/q}\log\frac{1}{\delta}\mbox{$\>\\!$}.$ (6.26)

######  Remark 6.2.

The bound (6.26) is very conservative due to a crude estimate (6.23)
leveraging just one instance $\lambda_{*}\in\check{\varLambda}^{q}(n,m)$. If
more information was available about the size of the space
$\check{\varLambda}^{q}(n,m)$, the bound (6.26) could be reduced accordingly.
For example, if $q=1$ then it is known that $\\#\check{\varLambda}(n,m)\sim
n^{m-1}(m!\,(m-1)!)^{-1}$ (see (4.40)). Hence, formula (6.26) for the time
threshold in the case $q=1$ is replaced by a much better and more realistic
estimate,

$\tilde{t}_{n}^{*}\simeq\frac{C_{1}(m,1)}{\\#\check{\varLambda}(n,m)}\,n^{m}\log\frac{1}{\delta}\sim\frac{\mathrm{e}^{m}\left(1-\mathrm{e}^{-m}\right)^{m}(m-1)!}{m^{m}}\,n\mbox{$\>\\!$}\log\frac{1}{\delta}\mbox{$\>\\!$}.$
(6.27)

Likewise, for $q=2$ we get, using (4.42) (with $C_{m}=1$),

$\tilde{t}_{n}^{*}\simeq\frac{C_{1}(m,2)}{\\#\check{\varLambda}^{2}(n,m)}\,n^{m/2}\log\frac{1}{\delta}\sim\left(\frac{2\mbox{$\;\\!$}\mathrm{e}}{m}\right)^{\\!m/2}\Gamma(m/2)\,\bigl{(}G_{1/2}(m/2)\bigr{)}^{m}\,n\mbox{$\>\\!$}\log\frac{1}{\delta}\mbox{$\>\\!$},$
(6.28)

which again grows only linearly in $n$.

  * (T2)

For the multiple exact sampling in the range $k\in[n,\theta\kern 0.24005ptn]$,
we can just repeat the procedure in task (T1) for each $k$ in that range.
According to (6.20), the probability that the number of attempts until
success, $T_{k}$, exceeds a threshold $t_{k}^{*}$ is given by (cf. (6.20))

${\mathsf{P}\mbox{$\>\\!\\!$}}_{\bm{z}}\bigl{(}T_{k}>t_{k}^{*}\,|\mbox{$\;\\!$}M_{\lambda}=m,\lambda_{\rm
max}\leq k\bigr{)}=(1-p_{k}^{*})^{t_{k}^{*}},\qquad n\leq k\leq\theta\kern
0.35999ptn.$

Taking into account only partitionable numbers $k\in[n,\theta\kern
0.24005ptn]$ (i.e., such that $\check{\varLambda}^{q}(k,m)\neq\varnothing$
and, therefore, $p_{k}^{*}>0$) and using a Bonferroni-type inequality, the
probability of Type I error for task (T2) (i.e., that the external loop fails
for at least one such $k$) is bounded as follows,

${\mathsf{P}\mbox{$\>\\!\\!$}}_{\bm{z}}\Biggl{(}\,\bigcup_{k=n}^{\lfloor\theta
n\rfloor}\\{t_{k}^{*}<T_{k}<\infty\\}\,\Bigl{|}\mbox{$\;\\!$}M_{\lambda}=m,\lambda_{\rm
max}\leq\theta\kern
0.35999ptn\Biggr{)}\leq\sum_{k\colon\mbox{$\>\\!\\!$}p_{k}^{*}>0}(1-p_{k}^{*})^{t_{k}^{*}}.$
(6.29)

Motivated by formula (6.26), we can look for the time limits $t_{k}^{*}$ in
the form

$t_{k}^{*}\sim c\mbox{$\;\\!$}k^{m/q}.$ (6.30)

Then from (6.29) using (6.23) we get, asymptotically,

$\displaystyle\sum_{k\colon\mbox{$\>\\!\\!$}p_{k}^{*}>0}(1-p_{k}^{*})^{t_{k}^{*}}$
$\displaystyle\lesssim\sum_{k=n}^{\lfloor\theta
n\rfloor}\left(1-\frac{1}{C_{1}(m,q)\,k^{m/q}}\right)^{c\mbox{$\>\\!$}k^{m/q}.}$
$\displaystyle\leq(\theta-1)\,n\mbox{$\;\\!$}\exp\left(-\frac{c}{C_{1}(m,q)}\right)\leq\delta.$

Solving this inequality for $c$ and returning to (6.30) ultimately yields

$\displaystyle t_{k}^{*}\simeq
C_{1}(m,q)\,k^{m/q}\log\frac{(\theta-1)\mbox{$\>\\!$}n}{\delta},\qquad n\leq
k\leq\theta\kern 0.35999ptn.$ (6.31)

Thus, the time bound (6.31) follows the same formula as in a single test (cf.
(6.26)) but with a Bonferroni-type adjustment of the significance level in
order to offset the multiple testing.

######  Remark 6.3.

The same comment as in Remark 6.2 applies to task (T2). Specifically, for
$q=1$ and $q=2$ the improved formulas for the thresholds $t^{*}_{k}$ are
given, respectively, by

$\displaystyle\tilde{t}_{k}^{*}$
$\displaystyle\simeq\frac{\mathrm{e}^{m}\left(1-\mathrm{e}^{-m}\right)^{m}(m-1)!}{m^{m}}\,k\mbox{$\>\\!$}\log\frac{(\theta-1)\mbox{$\>\\!$}n}{\delta},$
(6.32) $\displaystyle\tilde{t}_{k}^{*}$
$\displaystyle\simeq\left(\frac{2\mbox{$\;\\!$}\mathrm{e}}{\pi
m}\right)^{\\!m/2}\Gamma(m/2)\,\biggl{(}\int_{0}^{m/2}\\!u^{-1/2}\mbox{$\;\\!$}\mathrm{e}^{-u}\,\mathrm{d}{u}\biggr{)}^{m}\,k\mbox{$\>\\!$}\log\frac{(\theta-1)\mbox{$\>\\!$}n}{\delta}.$
(6.33)

  * (T3)

In a single attempt, the external loop gets a partition
$\lambda\in\bigcup_{n\leq k\leq\theta\kern
0.24005ptn}\check{\varLambda}^{q}(k,m)$ with probability

$\displaystyle{\mathsf{P}\mbox{$\>\\!\\!$}}_{\bm{z}}\bigl{(}n\leq
N_{\lambda}\leq\theta\kern
0.24005ptn\>|\mbox{$\;\\!$}M_{\lambda}=m,\lambda_{\rm max}\leq\theta\kern
0.24005ptn)$ $\displaystyle=\sum_{k=n}^{\lfloor\theta n\rfloor}p_{k}^{*}\to
G_{1/q}^{\star
m}(a_{\theta}\mbox{$\>\\!$}|\mbox{$\;\\!$}a_{\theta})-G_{1/q}^{\star
m}(a_{1}|\mbox{$\;\\!$}a_{\theta}).$ (6.34)

The limit (6.34) is due to Theorem 4.2(b), where
$a_{\theta}=\theta\mbox{$\>\\!$}m/q$ and $G_{1/q}^{\star
m}(x\mbox{$\;\\!$}|\mbox{$\;\\!$}a_{\theta})$ stands for the $m$-convolution
of the $a_{\theta}$-truncated gamma distribution
$G_{1/q}(x\mbox{$\;\\!$}|\mbox{$\;\\!$}a_{\theta})$.

To circumvent the trouble of computing such a convolution, observe that in the
range $0\leq x\leq a$ the distribution function $G^{\star
m}_{\alpha}(x\mbox{$\;\\!$}|\mbox{$\;\\!$}a)$ coincides with
$G_{m\mbox{$\>\\!$}\alpha}(x)$ up to the normalization factor
$G_{\alpha}(a)^{m}$. This is obvious for $m=1$, and the general case can be
seen by induction over $m$ using the convolution formula. Indeed, denoting the
corresponding densities by $g^{\star
m}_{\alpha}(x\mbox{$\;\\!$}|\mbox{$\;\\!$}a)$ and $g_{m\alpha}(x)$,
respectively, we have by definition
$g_{\alpha}(x\mbox{$\;\\!$}|\mbox{$\;\\!$}a)=g_{\alpha}(x)/G_{\alpha}(a)$
($0\leq x\leq a$), and the induction step is carried out as follows,

$\displaystyle g^{\star m}_{\alpha}(x\mbox{$\;\\!$}|\mbox{$\;\\!$}a)$
$\displaystyle=\int_{0}^{x}\\!g^{\star(m-1)}_{\alpha}(u\mbox{$\;\\!$}|\mbox{$\;\\!$}a)\,g_{\alpha}(x-u\mbox{$\;\\!$}|\mbox{$\;\\!$}a)\,\mathrm{d}{u}$
$\displaystyle=\frac{1}{G_{\alpha}(a)^{m}}\int_{0}^{x}\\!g_{(m-1)\mbox{$\>\\!$}\alpha}(u)\,g_{\alpha}(x-u)\,\mathrm{d}{u}=\frac{g_{m\mbox{$\>\\!$}\alpha}(x)}{\bigl{(}G_{\alpha}(a)\bigr{)}^{m}},\qquad
0\leq x\leq a,$

due to the convolution property of the gamma distribution. Thus, formula
(6.34) simplifies to

$\sum_{k=n}^{\lfloor\theta
n\rfloor}p_{k}^{*}\to\frac{G_{m/q}(\theta\mbox{$\>\\!$}m/q)-G_{m/q}(m/q)}{{\bigl{(}G_{1/q}(\theta\mbox{$\>\\!$}m/q)\bigr{)}^{m}}}=:C_{3}(m,q,\theta).$
(6.35)

Now, since individual runs of the external loop are independent, the
probability that the number of attempts until success, $T$, exceeds a
threshold $t$ is given by (cf. (6.20))

${\mathsf{P}\mbox{$\>\\!\\!$}}_{\bm{z}}(T>t\,|\mbox{$\;\\!$}M_{\lambda}=m,\lambda_{\rm
max}\leq\theta\kern 0.24005ptn)=\Biggl{(}1-\sum_{k=n}^{\lfloor\theta
n\rfloor}p_{k}^{*}\Biggr{)}^{\\!t}\simeq\bigl{(}1-C_{3}(m,q,\theta)\bigr{)}^{t},$
(6.36)

on account of (6.35). Hence, in order that this probability be bounded by a
confidence tolerance $\delta>0$, we may choose the threshold $t=t^{*}$ as
follows,

$t^{*}\simeq\frac{\log\delta}{\log\left(1-C_{3}(m,q,\theta)\right)}.$ (6.37)

######  Remark 6.4.

According to formula (6.37), the threshold $t^{*}$ does not depend on $n$. It
is of interest to look at how it depends on the growth of the power $q$. To
this end, by a direct analysis of the gamma distribution (see (4.1)) one can
verify that

$G_{\alpha}(\theta\mbox{$\>\\!$}\alpha)=(\theta\mbox{$\>\\!$}\alpha)^{\alpha}\bigl{(}1+O(\alpha)\bigr{)}\qquad(\alpha\to
0{+}).$

Hence, from formula (6.35) we get

$C_{3}(m,q,\theta)\sim
1-\theta^{-m/q}\sim\frac{m\log\theta}{q}\qquad(q\to\infty),$

and then (6.37) gives

$t^{*}\simeq\frac{q\mbox{$\;\\!$}\log\mbox{$\>\\!\\!$}(1/\delta)}{m\mbox{$\>\\!$}\log\theta}.$

####  6.2.3 Rejection sampling algorithm

A stylized example of rejection sampler is presented below in pseudocode as
Algorithm 2. It is set out in a flexible way so as to be usable in exact and
approximate sampling alike, as determined by the tasks (T1)–(T3) described in
Section 6.2.2. In particular, the range parameter $\theta$ is allowed to take
the value $\theta=1$, in which case the algorithm would work towards the exact
sampling task (T1) (i.e., with a specific weight target $N_{\lambda}=n$). As
explained in Section 6.2.1, the hyper-parameters are adapted to the desired
targets, $\braket{N}=n$ and $\braket{M}=m$, and the calibrating parameters
$z_{1}$ and $z_{2}$ are set according to the simplified expressions (6.1). A
predefined time bound $t^{*}$ for the external loop is selected according to
the task at hand, as discussed in Section 6.2.2, and on account of the
required confidence probability $1-\delta$.

As briefly indicated at the start of Section 6.2, Algorithm 2 comprises an
external loop that iterates the free sampler in an internal loop (i.e.,
Algorithm 1, with the majorant $L=\theta\kern 0.24005ptn$), which delivers, in
each productive cycle, a partition $\lambda$ that meets the length target
$M_{\lambda}=m$. This continues until the trial partition $\lambda$ meets the
weight target (e.g., $N_{\lambda}=n$ in task (T1) or $n\leq
N_{\lambda}\leq\theta\kern 0.35999ptn$ in task (T3)). However, if the limit of
attempts $t^{*}$ is reached with no success then the algorithm terminates,
returning a message ‘VOID’. It remains to add that for task (T2) involving
multiple exact sampling, the algorithm should be run in an additional loop to
scan all weight values in the range $k\in[n,\theta\kern 0.24005ptn]$.

Input: integer $q,n,m$, real $\theta\geq 1$, $t^{*}$

Output: partition
$\lambda\in\check{\varLambda}^{q}(\mathbin{\vbox{\hbox{\scalebox{0.5}{$\bullet$}}}}\mbox{$\>\\!$},m)$
with $N_{\lambda}\in[n,\theta\kern 0.24005ptn]$, otherwise ‘VOID’

1 integer array $\lambda_{[\ ]}$;

2 integer $N,M,t$;

3 real $L$;

4 $L\leftarrow\theta\kern 0.24005ptn$;

5 $N\leftarrow 0$, $M\leftarrow 0$, $t\leftarrow 0$;

6 while  _$N\notin[n,\theta\kern 0.24005ptn]$ and $t\leq t^{*}$_ do

7 while  _$M\neq m$ _ do

8 $(\lambda,N,M)\leftarrow{\tt FreeSampler}(q,n,m,L)$

9 end while

10 $t\leftarrow t+1$;

11

12 end while

13if  _$t\leq t^{*}$_ then

14 $N_{\lambda}\leftarrow N$;

15 return $(\lambda,N_{\lambda})$

16 else

17 return ‘VOID’

18 end if

Algorithm 2 ReSampler ($q,n,m,\theta,t^{*}$)

Algorithm 2 can be optimized in a number of ways. Since the weight of a valid
output $\lambda$ should not exceed $\theta\kern 0.35999ptn$, it is clear that
the run of the internal loop can be terminated prior to collecting the
required number of parts $m$ if the next candidate part is too large, so that
the incremented weight will certainly exceed the majorant. Furthermore, if the
number of collected parts has already reached the target value $m$ then there
is no need to keep scanning the remaining values in the range $\ell\leq L$ and
the current run of the internal loop may be stopped without any loss. However,
to avoid bias and maintain the Boltzmann distribution of the output, the
corresponding proposal
$\lambda\in\check{\varLambda}_{L}^{q}(\mathbin{\vbox{\hbox{\scalebox{0.5}{$\bullet$}}}}\mbox{$\>\\!$},m)$
must be accepted only if the remaining candidate parts in the range $\ell\leq
L$ were to be rejected by the respective Bernoulli checks. Since individual
such checks are mutually independent, their multitude can be replaced by a
single Bernoulli trial with the corresponding product probability of failure.
An additional benefit of such an aggregated Bernoulli check is that this will
reduce the number of calls of the random number generator and hence improve
the efficiency of the sampler.

Another improvement of the code implementation in the multiple testing task
(T2) proceeds from the observation that the sequential procedure based on
separate testing of each target in the range $[n,\theta\kern 0.24005ptn]$ (see
Section 6.2.2) is apparently wasteful, because a partition of some
$k^{\prime}\in[n,\theta\kern 0.24005ptn]$ obtained whilst looking for
partitions of a different number $k$ would be discarded in that cycle of the
external loop, whereas keeping it would have helped to achieve success if an
earlier search with target $k^{\prime}$ failed, or to save time on a duplicate
job when the algorithm moves to the new target $k^{\prime}$. In practice, all
partitions (at least, the new ones) obtained in every run of the external loop
should be stored as long as they fit into the range $[n,\theta\kern
0.24005ptn]$, thus leaving dynamically fewer targets to address.

For the sake of presentational clarity, Algorithm 2 embeds iterated calls of
the free sampler (Algorithm 1), but this means that the calibrating parameters
$z_{1}$ and $z_{2}$ are recalculated at every such call, which is of course
wasteful. This drawback can be easily amended by writing out the code
explicitly. Note, however, that such an improvement would have no significant
bearing on the asymptotic estimation of the code complexity.

####  6.2.4 Complexity and performance

Building on the probabilistic analysis of the internal and external loops
carried out in Section 6.2.2, it is straightforward to estimate the _time
complexity_ of Algorithm 2, understood as the _expected number of elementary
runs to completion_.

Starting with the internal loop, in its crude (non-optimized) version each
internal run comprises $\lfloor(\theta\kern 0.35999ptn)^{1/q}\rfloor$ checks
of available parts $\ell\in\mathbb{N}^{q}$ not exceeding
$L^{*}\mbox{$\>\\!\\!$}=\theta\kern 0.35999ptn$. Combined with the estimate
(6.17) of the probability to collect $m$ parts in a single run and the
corresponding geometric distribution of the number of attempts, the complexity
of the internal loop is bounded by

$\mu_{\theta}^{-m}\mbox{$\>\\!$}m!\,\mathrm{e}^{\mbox{$\>\\!$}\mu_{\theta}}\mbox{$\>\\!$}(\theta\kern
0.35999ptn)^{1/q}.$ (6.38)

As for the external loop, its complexity depends on the task at hand. If
$\theta=1$ (which corresponds to task (T1) of exact sampling with the weight
target $N_{\lambda}=n$), then the time to completion, $T_{n}$, has geometric
distribution with parameter $p_{n}^{*}$ (see (6.20)). For simplicity dropping
a time bound $t^{*}$ (but still assuming that the space
$\check{\varLambda}^{q}(n,m)$ is non-empty, so that $p_{n}^{*}>0$), the
expected time to completion is given by
${\mathsf{E}}_{\bm{z}}(T_{n})=1/p_{n}^{*}$. With a time bound $t^{*}$, the
expectation is modified as follows,

$\displaystyle{\mathsf{E}}_{\bm{z}}(T_{n};T_{n}\leq t^{*})$
$\displaystyle=\sum_{t=1}^{\,t^{*}\\!}t\mbox{$\;\\!$}(1-p_{n}^{*})^{t-1}p_{n}^{*}+t^{*}(1-p_{n}^{*})^{t^{*}}=\frac{1-(1-p_{n}^{*})^{t^{*}}}{p_{n}^{*}}<\frac{1}{p_{n}^{*}}.$
(6.39)

However, the reduction in (6.39) is not significant, because under our
confidence-based choice of the time limit (see (6.23)), we always have
$(1-p_{n}^{*})^{t^{*}}\mbox{$\>\\!\\!$}\\!\leq\delta$. Thus, combining
formulas (6.38) and (6.39), the total complexity guarantee for task (T1) is
estimated by

$\frac{m!\,\mathrm{e}^{\mbox{$\>\\!$}\mu_{1}}}{\mu_{1}^{m}}\mbox{$\;\\!$}O\bigl{(}n^{1/q}\mbox{$\;\\!\\!$}/p_{n}^{*}\bigr{)},\qquad\mu_{1}=m\,G_{1/q}(m/q).$
(6.40)

Further specification depends on the informative lower bound for the
probability $p_{n}^{*}$. For example, a crude estimate (6.23) gives a more
explicit estimate for the complexity,

$\left(\frac{\exp\bigl{(}G_{1/q}(m/q)+1/q\bigr{)}\,\Gamma(1/q)}{q^{1-1/q}\mbox{$\>\\!$}m^{1+1/q}}\right)^{\mbox{$\>\\!\\!$}m}\mbox{$\;\\!$}O\bigl{(}n^{(m+1)/q}\bigr{)}.$
(6.41)

For $q=1$ and $q=2$, this estimate can be significantly improved by using
asymptotically exact cardinalities (4.40) and (4.42), respectively, yielding
the estimates

$\frac{(m!)^{2}\,\mathrm{e}^{2m}}{m^{2m+1}}\,O(n^{2})=O(n^{2})$ (6.42)

and

$\frac{2^{m/2}\mbox{$\>\\!$}m!\>\Gamma(m/2)\,\mathrm{e}^{3m/2}}{m^{3m/2}}\mbox{$\;\\!$}O\bigl{(}n^{3/2}\bigr{)}=O(n^{3/2}).$
(6.43)

Interestingly, the asymptotic bounds (6.42) and (6.43) do not depend on the
number of parts $m$.

For task (T2) (with some $\theta>1$), the above estimates just need to be
multiplied by the number of targeted weights,
$\lfloor(\theta-1)\mbox{$\;\\!$}n\rfloor+1=O(n)$. Finally, for task (T3) we
can use formula (6.40), but with $\mu_{1}$ changed to $\mu_{\theta}$ (see
(6.38)) and with the probability $p_{n}^{*}$ of success in a single attempt
replaced by the (asymptotic) probability (6.35) of at least one success in the
range $[n,\theta\kern 0.03607ptn]$, yielding

$\frac{m!\,\exp\bigl{(}m\,G_{1/q}(\theta\mbox{$\>\\!$}m/q)\bigr{)}}{m^{m}\mbox{$\;\\!$}\bigl{(}G_{m/q}(\theta\mbox{$\>\\!$}m/q)-G_{m/q}(m/q)\bigr{)}}\,O(n^{1/q})=O\bigl{(}m\mbox{$\>\\!$}n^{1/q}\bigr{)}.$

Table 4: Confident time thresholds $t^{*}$ for the external loop in tasks
(T1), (T2) and (T3), calculated for $q=1$ and $q=2$ with various values of
confidence tolerance $\delta$ using formulas (6.26), (6.31) and (6.37). In
both cases, the chosen values of $n$ and $m$ yield
$\kappa=m^{q+1}\mbox{$\>\\!\\!$}/n=0.01$. For tasks (T2) and (T3), the testing
range is set with the factor $\theta=1.1$. For comparison, corrected values
$\tilde{t}^{*}$ for tasks (T1) and (T2) are calculated from formulas (6.27),
(6.28) and (6.32), (6.33), respectively. $\vphantom{\int^{t}}q=1$,
$n=2{,}\mbox{$\>\\!$}500$, $m=5$
---
$\vphantom{\int_{y}^{t}}\delta$ |  $t_{n}^{*}$ (T1) |  $\tilde{t}_{n}^{*}\vphantom{\int^{T}}$ (T1) |  $t_{n}^{*}$ (T2) |  $\tilde{t}_{n}^{*}$ (T2) |  $t^{*}$ (T3)
$0.1$ | $8.603518\cdot 10^{13}\vphantom{\int^{t}}$ | $6{,}\mbox{$\>\\!$}343.202$ | $2.923424\cdot 10^{14}$ | $21{,}\mbox{$\>\\!$}553.82$ | $26.01955$
$0.01$ | $1.720704\cdot 10^{14}$ | $12{,}\mbox{$\>\\!$}686.40$ | $3.783776\cdot 10^{14}\vphantom{\int^{t}}$ | $27{,}\mbox{$\>\\!$}897.02$ | $52.03911$
$0.001$ | $2.581056\cdot 10^{14}$ | $19{,}\mbox{$\>\\!$}029.61$ | $4.644128\cdot 10^{14}\vphantom{\int^{t}}$ | $34{,}\mbox{$\>\\!$}240.22$ | $78.05866$
$0.0001$ | $3.441407\cdot 10^{14}$ | $25{,}\mbox{$\>\\!$}372.81$ | $5.504479\cdot 10^{14}\vphantom{\int^{t}}$ | $40{,}\mbox{$\>\\!$}583.43$ | $104.0782$
$q=2$, $n=12{,}\mbox{$\>\\!$}500$, $m=5$
$\vphantom{\int_{y}^{t}}\delta$ |  $t_{n}^{*}$ (T1) |  $\tilde{t}_{n}^{*}$ (T1) |  $t_{n}^{*}$ (T2) |  $\tilde{t}_{n}^{*}$ (T2) |  $t^{*}$ (T3)
$\vphantom{\int^{t}}0.1$ | $198{,}\mbox{$\>\\!$}687{,}\mbox{$\>\\!$}146$ | $41{,}\mbox{$\>\\!$}485.61$ | $814{,}\mbox{$\>\\!$}003{,}\mbox{$\>\\!$}357$ | $169{,}\mbox{$\>\\!$}962.8$ | $34.94282$
$\vphantom{\int^{t}}0.01$ | $397{,}\mbox{$\>\\!$}374{,}\mbox{$\>\\!$}292$ | $82{,}\mbox{$\>\\!$}971.22$ | $1{,}\mbox{$\>\\!$}012{,}\mbox{$\>\\!$}690{,}\mbox{$\>\\!$}503$ | $211{,}\mbox{$\>\\!$}448.4$ | $69.88564$
$\vphantom{\int^{t}}0.001$ | $596{,}\mbox{$\>\\!$}061{,}\mbox{$\>\\!$}438$ | $124{,}\mbox{$\>\\!$}456.8$ | $1{,}\mbox{$\>\\!$}211{,}\mbox{$\>\\!$}377{,}\mbox{$\>\\!$}649$ | $252{,}\mbox{$\>\\!$}934.0$ | $104.8285$
$\vphantom{\int^{t}}0.0001$ | $794{,}\mbox{$\>\\!$}748{,}\mbox{$\>\\!$}583$ | $165{,}\mbox{$\>\\!$}942.4$ | $1{,}\mbox{$\>\\!$}410{,}\mbox{$\>\\!$}064{,}\mbox{$\>\\!$}795$ | $294{,}\mbox{$\>\\!$}419.7$ | $139.7713$

To evaluate real time performance of the rejection sampler, we first need to
take a practical look at the censoring time limits $t^{*}$ in tasks (T1)–(T3)
proposed in Section 6.2.2. These are numerically illustrated in Table 4 for
$q=1$ and $q=2$, with various values of $n$ and $m$. Observe that the crude
bounds for tasks (T1) and (T2) calculated via formulas (6.26) and (6.31)
appear to be very high, especially for $q=1$ (of order $10^{14}$), casting
doubt on whether such limits are usable. In real terms, since each run of the
external loop is a simple check if $N_{\lambda}=n$ (see line 2 in Algorithm
2), we can assume for simplicity that it needs a single tick of the CPU clock.
If the algorithm is executed on a contemporary mid-range desktop PC (say, with
processor base frequency 3.30 GHz, which we used) then, under the estimate
(6.26) for task (T1) with $q=1$ and a fairly low confidence tolerance
$\delta=0.001$, the external loop alone may require up to $2.581056\cdot
10^{14}/(3.30\cdot 10^{9}\cdot 60\cdot 60)\doteq 21.72606\approx 22$ hours
until completion, which is unpleasantly long but not entirely unrealistic.
This estimate drops dramatically for $q=2$ to less than $1$ second. A steep
decreasing trend continues with larger powers;131313Keeping $m$ and
$\kappa=m^{q+1}\mbox{$\>\\!\\!$}/n$ fixed, from formulas (6.24) and (6.26) we
find
$\lim_{q\to\infty}t_{n}^{*}=(m^{m}\mbox{$\;\\!\\!$}/m!)\log\mbox{$\>\\!\\!$}(1/\delta)$.
For example, for $m=5$ as in Table 4 and $\delta=0.001$, this limiting value
specializes to $179.8895$ for example, for $q=3$, $n=62{,}\mbox{$\>\\!$}500$,
$m=5$ and $\delta=0.001$, formula (6.26) gives
$t_{n}^{*}=3{,}\mbox{$\>\\!$}358{,}\mbox{$\>\\!$}531$, leading to the maximum
execution time of up to $0.001$ second. Thus, the sampler becomes
progressively more efficient for larger $q$, even under a crude time bound. On
the other hand, as pointed out in Remarks 6.2 and 6.3, additional information
about the size of the corresponding partition spaces would allow a significant
reduction of the estimated bound as illustrated in Table 4 (by a factor
$10^{10}$ for $q=1$ and about $4{,}\mbox{$\>\\!$}790$ for $q=2$).

Let us now look at the real time computational cost due to the internal loop.
As mentioned before (cf. (6.38)), the expected numbers of internal runs until
collecting exactly $m$ parts is asymptotically given by
$\mu_{1}^{-m}\mbox{$\>\\!$}m!\,\mathrm{e}^{\mbox{$\>\\!$}\mu_{1}}$ (with
$\theta=1$), where $\mu_{1}=m\,G_{1/q}(m/q)$. Using for numerical illustration
the same values of $q$, $n$ and $m$ as in Table 4, this formula yields
$5.6993$ ($q=1$) and $5.7043$ ($q=2$). The average computing time for each of
such attempts is inversely proportional to the CPU base frequency (such as
3.30 GHz), but it involves many other important aspects such as the
operational efficiency of a random number generator, design of memory
allocation and data storage, numerical precision, coding implementation and
compiler used, and the overall architecture of the computer (e.g., the number
of cores and whether or not parallel processing was utilized). Thus, it is
impossible to estimate the actual computing time without real benchmarking.

To test the performance of the internal loop, we implemented the algorithm on
a desktop CPU as described at the beginning of Section 6, for simplicity using
a single core. Since internal runs are independent from each other and the
computational costs due to the multi-core design are negligible, we can simply
divide the average execution time on a single core by the number of cores at
disposal.

With the same values of $q$, $n$ and $m$ used above and in Table 4 (and with
$\theta=1$), the average number of sampling attempts (starting at line 2 of
Algorithm 2) was $5.6956$ for $q=1$ and $5.6404$ for $q=2$; note that these
sample averages match the expected values calculated above. Furthermore, the
program took on average $2.1036\cdot 10^{-3}\mbox{$\;\\!\\!$}$ seconds ($q=1$)
and $0.9780\cdot 10^{-4}\mbox{$\;\\!\\!$}$ seconds ($q=2$) per single
successful completion of the internal loop. The corresponding number of ticks
of the CPU clock per elementary check of a candidate part $\ell\leq n$ (see
formula (6.38)) is evaluated as $(2.1036\cdot
10^{-3}\mbox{$\;\\!\\!$}/(5.6956\cdot 2{}\mbox{$\;\\!$}500))\cdot 3.30\cdot
10^{9}\doteq 487.5258$ ($q=1$) and $(0.9780\cdot
10^{-4}\mbox{$\;\\!\\!$}/(5.6404\,\sqrt{12{}\mbox{$\;\\!$}500}\mbox{$\;\\!$}))\cdot
3.30\cdot 10^{9}=511.7854$, so it stays in the range about $450\div 550$.

However, there is a problem: if we combine the physical times benchmarked for
the internal loop with the time bounds $t^{*}_{n}$ for the external loop given
in Table 4 (say, with tolerance $\delta=0.001$), then for $q=1$ we obtain, by
converting seconds to minutes, hours, days and years, $2.1036\cdot
10^{-3}\mbox{$\;\\!\\!$}\cdot 2.581056\cdot 10^{14}\mbox{$\;\\!\\!$}/(60\cdot
60\cdot 24\cdot 365)\approx 17{,}\mbox{$\>\\!$}217$ years (!), which is
clearly impractical. For $q=2$, a similar calculation gives a more reasonable
estimate, $0.9780\cdot 10^{-4}\mbox{$\;\\!\\!$}\cdot
596{}\,061{}\mbox{$\;\\!$}438/(60\cdot 60)\approx 16$ hours. But with the
improved time bounds $\tilde{t}_{n}^{*}$ (see Table 4), we obtain much more
satisfactory estimates, $2.1036\cdot 10^{-3}\mbox{$\;\\!\\!$}\cdot
19\,029.61\approx 40$ seconds ($q=1$) and $0.9780\cdot 10^{-4}\cdot
124\,456.8\approx 12$ seconds ($q=2$).

As a concluding remark, Algorithm 2 could be used as an experimental tool for
searching satisfactory instances in additive problems of number theory such as
variants of the Waring problem. Here, iterated sampling attempts to find a
suitable instance subject to certain constraints — in the lack of prior
knowledge if such instances even exist — may be interpreted as statistical
testing of existence as the null hypothesis, under which a suitable confident
time limit $t^{*}$ can be determined. The sampling approach based on bounded
(although high) confidence bears similarity with primality testing procedure
such as the Miller–Rabin algorithm [66]. It can also be helpful in
verification of conjectures about the density of representable numbers via an
experimental analysis of the success rates. We will address such applications
in another paper.

## Acknowledgements

J.C.P. was supported by an EPSRC Doctoral Training Partnership scholarship at
the School of Mathematics, University of Leeds. The authors have benefited
from the useful discussions with Yuri V. Yakubovich (Saint Petersburg).

## References

  * [1]
  * [2]
  * [3] Agarwala, B.K. and Auluck, F.C. Statistical mechanics and partitions into non-integral powers of integers. _Math. Proc. Cambridge Philos. Soc._ 47 (1951), 207–216. (doi:10.1017/S0305004100026505; MR0042998)
  * [4] Andrews, G.E. The Theory of Partitions. Encyclopedia of Mathematics and its Applications, 2, G.-C. Rota, ed. Addison-Wesley, Reading, MA, 1976; reprinted by Cambridge University Press (Cambridge Mathematical Library), Cambridge, 1984. (doi:10.1017/CBO9780511608650; MR0557013)
  * [5] Arratia, R., Barbour, A.D. and Tavaré, S. _Logarithmic Combinatorial Structures: A Probabilistic Approach_. EMS Monographs in Mathematics. European Mathematical Society, Zürich, 2003. (doi:10.4171/000; MR2032426)
  * [6] Arratia, R. and Tavaré, S. Independent process approximations for random combinatorial structures. _Adv. Math._ 104 (1994), 90–154. (doi:10.1006/aima.1994.1022; MR1272071)
  * [7] Auluck, F.C. and Kothari, D.S. Statistical mechanics and the partitions of numbers. Proc. Cambridge Philos. Soc. 42 (1946), 272–277. (doi:10.1017/S030500410002303; MR0017682)
  * [8] Barbour, A. and Hall, P. On the rate of Poisson convergence. _Math. Proc. Cambridge Philos. Soc._ 95 (1984), 473–480. (doi:10.1017/S0305004100061806; MR0755837)
  * [9] Barbour, A.D., Holst, L. and Janson, S. Poisson Approximation. Oxford Studies in Probability, 2. Oxford Science Publications. Clarendon Press, Oxford, 1992. (MR1163825)
  * [10] Bendkowski, M., Bodini, O. and Dovgal, S. Tuning as convex optimisation: A polynomial tuner for multi-parametric combinatorial samplers. _Combin. Probab. Comput._ 31 (2022), 765–811. (doi:10.1017/S0963548321000547; MR4472289)
  * [11] Bernstein, M., Fahrbach, M. and Randall, D. Analyzing Boltzmann samplers for Bose–Einstein condensates with Dirichlet generating functions. In _2018 Proceedings of the Meeting on Analytic Algorithmics and Combinatorics (ANALCO)_ , M. Nebel and S. Wagner, eds. SIAM, Philadelphia, PA, 2018, pp. 107–117. (doi:10.1137/1.9781611975062.10; MR3773639)
  * [12] Bodini, O., Duchon, P., Jacquot, A. and Mutafchiev, L. Asymptotic analysis and random sampling of digitally convex polyominoes. In DGCI 2013: Discrete Geometry for Computer Imagery, Proceedings of the 17th IAPR International Conference (Seville, March 20–22, 2013), R. Gonzalez-Diaz et al., eds. Lecture Notes in Computer Science, 7749. Springer, Berlin, 2013, pp. 95–106. (doi:10.1007/978-3-642-37067-0_9; MR3155272)
  * [13] Bogachev, L.V. Unified derivation of the limit shape for multiplicative ensembles of random integer partitions with equiweighted parts. _Random Struct. Algorithms_ 47 (2015), 227–266. (doi:10.1002/rsa.20540; MR3382672)
  * [14] Bogachev, L.V. and Su, Z. Gaussian fluctuations of Young diagrams under the Plancherel measure. _Proc. Royal Soc. Ser. A_ 463 (2007), 1069–1080. (doi:10.1098/rspa.2006.1808; MR2310137)
  * [15] Bogachev, L.V. and Yakubovich, Yu.V. Limit shape of minimal difference partitions and fractional statistics. _Comm. Math. Phys._ 373 (2020), 1085–1131. (doi:10.1007/s00220-019-03513-5; MR4061406)
  * [16] Bogachev, L.V. and Zarbaliev, S.M. Universality of the limit shape of convex lattice polygonal lines. _Ann. Probab._ 39 (2011), 2271–2317. (doi:10.1214/10-AOP607; MR2932669)
  * [17] Borodin, A., Okounkov, A. and Olshanski, G. Asymptotics of Plancherel measures for symmetric groups. _J. Amer. Math. Soc._ 13 (2000), 481–515. (doi:10.1090/S0894-0347-00-00337-4, MR1758751)
  * [18] Bureaux, J. and Enriquez, N. Asymptotics of convex lattice polygonal lines with a constrained number of vertices. _Isr. J. Math._ 222 (2017), 515–549. (https://doi.org/10.1007/s11856-017-1599-3; MR3722260)
  * [19] Comtet, A., Leboeuf, P. and Majumdar, S.N. Level density of a Bose gas and extreme value statistics. Phys. Rev. Lett. 98 (2007), 070404, 1–4. (doi:10.1103/PhysRevLett.98.070404)
  * [20] Comtet, A., Majumdar, S.N., Ouvry, S. and Sabhapandit, S. Integer partitions and exclusion statistics: Limit shapes and the largest parts of Young diagrams. _J. Stat. Mech. Theory Exp._ (2007), no. 10, P10001, 1–13. (doi:10.1088/1742-5468/2007/10/P10001; MR2358050)
  * [21] Cramér, H. Mathematical Methods of Statistics. Princeton Mathematical Series, 9. Princeton University Press, Princeton, NJ, 1946\. (https://www.jstor.org/stable/j.ctt1bpm9r4; MR0016588)
  * [22] De Gregorio, P. and Rondoni, L. Microcanonical entropy, partitions of a natural number into squares and the Bose–Einstein gas in a box. _Entropy_ 20 (2018), 645\. (doi:10.3390/e20090645)
  * [23] Duchon, P., Flajolet, P., Louchard, G. and Schaeffer, G. Boltzmann samplers for the random generation of combinatorial structures. _Combin. Probab. Comput._ 13 (2004), 577–625. (doi:10.1017/S0963548304006315; MR2095975)
  * [24] Erdős, P. and Lehner, J. The distribution of the number of summands in the partitions of a positive integer. _Duke Math. J._ 8 (1941), 335–345. (doi:10.1215/S0012-7094-41-00826-8; MR0004841)
  * [25] Erdős, P. and Szalay, M. On the statistical theory of partitions. In _Topics in Classical Number Theory, Vol. I (Budapest, 1981)_ , G. Halász, ed., pp. 397–450. Colloquia Mathematica Societatis János Bolyai, 34. North-Holland, Amsterdam, 1984. (MR0781149)
  * [26] Erdős, P. and Turán, P. On some general problems in the theory of partitions, I. _Acta Arith._ 18 (1971), 53–62. (doi:10.4064/aa-18-1-53-62; MR0289446)
  * [27] Erlihson, M.M. and Granovsky, B.L. Limit shapes of Gibbs distributions on the set of integer partitions: The expansive case. _Ann. Inst. Henri Poincaré Probab. Stat._ 44 (2008), 915–945. (https://doi.org/10.1214/07-AIHP129; MR2453776)
  * [28] Ewens, W.J. The sampling theory of selectively neutral alleles. _Theoret. Population Biol._ 3 (1972), 87–112. (doi:10.1016/0040-5809(72)90035-4; MR0325177)
  * [29] Feller, W. An Introduction to Probability Theory and Its Applications, Volume II, 2nd ed. Wiley Series in Probability and Mathematical Statistics. Wiley, New York, 1971. (MR0270403)
  * [30] Flajolet, P., Fusy, É. and Pivoteau, C. Boltzmann sampling of unlabelled structures. In _2007 Proceedings of the Workshop on Analytic Algorithmics and Combinatorics (ANALCO)_ , D. Panario and R. Sedgewick, eds. SIAM, Philadelphia, PA, 2007, pp. 201–211. (doi:10.1137/1.9781611972979.5; MR2498128)
  * [31] Freiman, G.A. and Granovsky, B.L. Clustering in coagulation-fragmentation processes, random combinatorial structures and additive number systems: asymptotic formulae and limiting laws. _Trans. Amer. Math. Soc._ 357 (2005), 2483–2507. (doi:10.1090/S0002-9947-04-03617-7; MR2140447)
  * [32] Freiman, G.A. and Yudin, A.A. The interface between probability theory and additive number theory (local limit theorems and structure theory of set addition). In _Representation Theory, Dynamical Systems, and Asymptotic Combinatorics_ , V. Kaimanovich and A. Lodkin, eds. _Amer. Math. Soc. Transl. Ser. 2_ 217 (2006), 51–72. American Mathematical Society, Providence, RI. (doi:10.1090/trans2/217/05/trans2217-05.pdf; MR2276101)
  * [33] Fristedt, B. The structure of random partitions of large integers. _Trans. Amer. Math. Soc._ 337 (1993), 703–735. (doi:10.1090/S0002-9947-1993-1094553-1; MR1094553)
  * [34] Fulton, W. Young Diagrams, with Applications to Representation Theory and Geometry. London Mathematical Society Student Texts, 35. Cambridge University Press, Cambridge, 1997\. (doi:10.1017/CBO9780511626241, MR1464693)
  * [35] Garibaldi, U. and Scalas, E. _Finitary Probabilistic Methods in Econophysics._ Cambridge University Press, Cambridge, 2010. (doi:10.1017/CBO9780511777585; MR2839542)
  * [36] Goh, W.M.Y. and Hitczenko, P. Random partitions with restricted part sizes. _Random Struct. Algorithms_ 32 (2008), 440–462. (doi:10.1002/rsa.20191; MR2422389)
  * [37] Greiner, W., Neise, L. and Stöcker, H. _Thermodynamics and Statistical Mechanics._ Classical Theoretical Physics. Springer, New York, 1995. (doi:10.1007/978-1-4612-0827-3)
  * [38] Guy, R.K. _Unsolved Problems in Number Theory_ , 3rd ed. Problem Books in Mathematics. Springer, New York, 2003. (doi:10.1007/978-0-387-26677-0; MR2076335)
  * [39] Hardy, G.H. and Ramanujan, S. Asymptotic formulæ in combinatory analysis. _Proc. London Math. Soc. (2)_ 17 (1918), 75–115. (doi:10.1112/plms/s2-17.1.75; MR1575586). Also available online at http://ramanujan.sirinudi.org/Volumes/published/ram36.pdf (accessed 26.11.2022).
  * [40] Hardy, G.H. and Wright, E.M. _An Introduction to the Theory of Numbers_ , 6th ed. Revised by D. R. Heath-Brown and J. H. Silverman, with a foreword by A. Wiles. Oxford University Press, Oxford, 2008. (Publisher’s web page; MR2445243)
  * [41] Haselgrove, C.B. and Temperley, H.N.V. Asymptotic formulae in the theory of partitions. _Proc. Cambridge Philos. Soc._ 50 (1954), 225–241. (doi:10.1017/S0305004100029273; MR0062782)
  * [42] Hohloch, S. Optimal transport and integer partitions. _Discrete Appl. Math._ 190/191 (2015), 75–85. (doi:10.1016/j.dam.2015.04.002; MR3351722)
  * [43] Horn, R.A. and Johnson, C.R. Matrix Analysis, 2nd ed. Cambridge University Press, Cambridge, 2013. (doi:10.1017/9781139020411; MR2978290)
  * [44] Hua, L.-K. On the number of partitions of a number into unequal parts. _Trans. Amer. Math. Soc._ 51 (1942), 194–201. (doi:10.2307/1989985; MR0006195)
  * [45] Huang, K. _Statistical Mechanics_ , 2nd ed. Wiley, New York, 1987. (MR1042093)
  * [46] Hwang, H.-K. Limit theorems for the number of summands in integer partitions. _J. Combin. Theory Ser. A_ 96 (2001), 89–126. (doi:10.1006/jcta.2000.3170; MR1855788)
  * [47] Ingham, A.E. A Tauberian theorem for partitions. _Ann. Math. (2)_ 42 (1941), 1075–1090. (doi:10.2307/1970462; MR0005522)
  * [48] Khinchin, A.I. _Mathematical Foundations of Statistical Mechanics_. Translated from Russian by G. Gamow. Dover Publications, New York, 1949. Available online at https://archive.org/details/khinchin-mathematical-foundations-of-statistical-mechanics (accessed 05.12.2022). (MR0029808)
  * [49] Khinchin, A.Y. _Mathematical Foundations of Quantum Statistics_. Translation from the Russian edition (1951), I. Shapiro, ed. Graylock Press, Albany, NY, 1960. Available online at https://archive.org/details/khinchin-mathematical-foundations-of-quantum-statistics (accessed 05.12.2022). (MR0111217)
  * [50] Kingman, J.F.C. Random partitions in population genetics. _Proc. Royal Soc. London Ser. A_ 36 (1978), 11–20. (doi:10.1098/rspa.1978.0089; MR0526801)
  * [51] Kolchin, V.F., Sevast’yanov, B.A. and Chistyakov, V.P. _Random Allocations._ Translated from the Russian edition (1976), A. V. Balakrishnan, ed. Scripta Series in Mathematics. V. H. Winston & Sons/Scripta Technica, Inc., Washington, D.C., 1978 [distributed by Halsted Press/Wiley, New York]. (MR0471016)
  * [52] Landau, E. Über die Einteilung der positiven ganzen Zahlen in vier Klassen nach der Mindestzahl der zu ihrer additiven Zusammensetzung erforderlichen Quadrate. (German) [On the division of positive integers into four classes according to the minimum number of squares required for their additive composition]. Arch. Math. Phys. (3) 13 (1908), 305–312. Available online at https://ia600309.us.archive.org/31/items/archivdermathem37unkngoog/archivdermathem37unkngoog.pdf (accessed 28.10.2022).
  * [53] Lee, D.V. The asymptotic distribution of the number of summands in unrestricted $\varLambda$-partitions. _Acta Arith._ 65 (1993), 29–43. (doi:10.4064/aa-65-1-29-43; MR1239241)
  * [54] Lindsay, R.B. (XXII) A modification of Brillouin’s unified statistics. _Lond. Edinb. Dublin Philos. Mag. J. Sci. Ser. 7_ 17 (1934), 264–271. (doi:10.1080/14786443409462389)
  * [55] Lipnik, G.F., Madritsch, M.G. and Tichy, R.F. A central limit theorem for integer partitions into small powers. Preprint (2022), arXiv:2204.05592. Available online at https://arxiv.org/abs/2204.05592 (accessed 08.12.2022).
  * [56] Loève, M. _Probability Theory I_ , 4th ed. Graduate Texts in Mathematics, 45. Springer, New York, 1977. (doi:10.1007/978-1-4684-9464-8; MR0651017)
  * [57] Logan, B.F. and Shepp, L.A. A variational problem for random Young tableaux. _Adv. in Math._ 26 (1977), 206–222. (doi:10.1016/0001-8708(77)90030-5; MR1417317)
  * [58] Meinardus, G. Asymptotische Aussagen über Partitionen. (German) [Asymptotic statements about partitions] _Math. Z._ 59 (1954), 388–398. (doi:10.1007/BF01180268; MR0062781)
  * [59] _NIST Handbook of Mathematical Functions_ , F. W.J̇. Olver et al., eds. National Institute of Standards and Technology, U.S. Department of Commerce, Washington, DC; Cambridge University Press, Cambridge, 2010. (https://www.cambridge.org/catalogue/catalogue.asp?isbn=9780521192255; MR2723248) Online version:  _NIST Digital Library of Mathematical Functions_ , Release 1.1.7 of 2022-10-15, https://dlmf.nist.gov (accessed 28.10.2022).
  * [60] Novak, S.Y. Poisson approximation. _Probab. Surveys_ 16 (2019), 228–276. (doi:10.1214/18-PS318; MR3992498)
  * [61] Nuermaimaiti, R., Bogachev, L.V. and Voss, J. A generalized power law model of citations. In _18th International Conference on Scientometrics & Informetrics, ISSI 2021 (12–15 July 2021, KU Leuven, Belgium), Proceedings_, W. Glänzel et al., eds. International Society for Scientometrics and Informetrics (I.S.S.I.), 2021, pp. 843–848. Available online at https://issi2021.org/proceedings/ (accessed 08.12.2022).
  * [62] Pak, I. The nature of partition bijections I. Involutions. _Adv. in Appl. Math._ 33 (2004), 263–289. (doi:10.1016/j.aam.2003.08.007; MR2074399)
  * [63] Pitman, J. _Combinatorial Stochastic Processes. (Ecole d’Eté de Probabilités de Saint-Flour XXXII – 2002)_ , J. Picard, ed. Lecture Notes in Mathematics, 1875. Springer, Berlin, 2006. (doi:10.1007/b11601500; MR2245368)
  * [64] Pittel, B. On a likely shape of the random Ferrers diagram. Adv. in Appl. Math. 18 (1997), 432–488. (doi:10.1006/aama.1996.0523; MR1445358)
  * [65] Polcari, J. An informative interpretation of decision theory: the information theoretic basis for signal-to-noise ratio and log likelihood ratio. _IEEE Access_ 1 (2013), 509–522. (doi:10.1109/ACCESS.2013.2277930)
  * [66] Rabin, M.O. Probabilistic algorithm for testing primality. J. Number Theory 12 (1980), 128–138. (doi:10.1016/0022-314X(80)90084-0; MR0566880)
  * [67] Resnick, S.I. Extreme Values, Regular Variation and Point Processes. Applied Probability. A Series of the Applied Probability Trust, 4. Springer, New York, 1987. (doi:10.1007/978-0-387-75953-1; MR0900810)
  * [68] Roccia, J. and Leboeuf, P. Level density of a Fermi gas and integer partitions: A Gumbel-like finite-size correction. _Phys. Rev. C_ 81 (2010), 044301. (doi:10.1103/PhysRevC.81.044301)
  * [69] Shiryaev, A.N. _Probability_ , 2nd ed. Graduate Texts in Mathematics, 95. Springer, New York, 1996. (doi:10.1007/978-1-4757-2539-1; MR1368405)
  * [70] Sinaĭ, Ya.G. A probabilistic approach to the analysis of the statistics of convex polygonal lines. (Russian) _Funktsional. Anal. i Prilozhen._ 28 (1994), no. 2, 41–48. English translation: Probabilistic approach to the analysis of statistics for convex polygonal lines. _Funct. Anal. Appl._ 28 (1994), 108–113. (doi:10.1007/BF01076497; MR1283251)
  * [71] Su, Z. _Random Matrices and Random Partitions: Normal Convergence._ World Scientific Series on Probability Theory and Its Applications, 1. World Scientific, Hackensack, NJ, 2015. (doi:10.1142/9197, MR3381296)
  * [72] Temperley, H.N.V. Statistical mechanics and the partition of numbers, I. The transition of liquid helium. _Proc. Royal Soc. London Ser. A_ 199 (1949), 361–375. (doi:10.1098/rspa.1949.0143; MR0037247)
  * [73] Temperley, H.N.V. Statistical mechanics and the partition of numbers, II. The form of crystal surfaces. Proc. Cambridge Philos. Soc. 48 (1952), 683–697. (doi:10.1017/S0305004100076453; MR0053036)
  * [74] Uspensky, J.V. Asymptotic expressions of numerical functions occurring in problems of partition of numbers. (Russian) _Bull. Acad. Sci. Rus. Ser. VI_ 14 (1920), 199–218. Available online at https://mi.mathnet.ru/eng/izv/v14/p199 (accessed 08.12.2022). (zbMATH JFM 48.1162.01)
  * [75] Vaughan, R.C. Squares: Additive questions and partitions. _Int. J. Number Theory_ 11 (2015), 1367–1409. (doi:10.1142/S1793042115400096; MR3376217)
  * [76] Vaughan, R.C. and Wooley, T.D. Waring’s problem: A survey. _Number Theory for the Millennium, III (Urbana, IL, 2000)_ , M. A. Bennett et al., eds. A K Peters, Natick, MA, 2002, pp. 301–340. (MR1956283) Reprinted in _Surveys in Number Theory: Papers from the Millennial Conference on Number Theory_. CRC Press, Boca Raton, FL, 2003, pp. 285–324. (doi:10.1201/9780429258978)
  * [77] Vershik, A.M. Asymptotic combinatorics and algebraic analysis. In _Proceedings of the International Congress of Mathematicians (Zürich, 1994)_ , Vol. 2, S. D. Chatterji, ed. Birkhäuser, Basel, 1995, pp. 1384–1394. (doi:10.1007/978-3-0348-9078-6_133; MR1404040)
  * [78] Vershik, A.M. Statistical mechanics of combinatorial partitions, and their limit shapes. (Russian) _Funktsional. Anal. i Prilozhen._ 30 (1996), no. 2, 19–39. English translation: _Funct. Anal. Appl._ 30 (1996), 90–105. (doi:10.1007/BF02509449; MR1402079)
  * [79] Vershik, A.M. Limit distribution of the energy of a quantum ideal gas from the viewpoint of the theory of partitions of natural numbers. (Russian) Uspekhi Mat. Nauk 52(2) (1997), 139–146; English translation: Russian Math. Surveys 52 (1997), 379–386. (doi:10.1070/RM1997v052n02ABEH001782; MR1480142)
  * [80] Vershik, A.M., Freĭman, G.A. and Yakubovich, Yu.V. A local limit theorem for random partitions of natural numbers. (Russian) _Teor. Veroyatnost. i Primenen._ 44 (1999), 506–525; English translation: Freiman, G., Vershik, A.M. and Yakubovich, Yu.V. A local limit theorem for random strict partitions._Theory Probab. Appl._ 44 (2000), 453–468. (doi:10.1137/S0040585X97977719; MR1805818)
  * [81] Vershik, A.M. and Kerov, S.V. Asymptotics of the Plancherel measure of the symmetric group and the limit form of Young tableaux. (Russian) Dokl. Akad. Nauk SSSR 233 (1977), 1024–1027; available online at http://mi.mathnet.ru/eng/dan/v233/i6/p1024 (accessed 08.12.2022); English translation: Soviet Math. Doklady 18 (1977), 527–531. (MR0480398)
  * [82] Vershik, A.M. and Kerov, S.V. Asymptotic behavior of the maximum and generic dimensions of irreducible representations of the symmetric group. (Russian) Funktsional. Anal. i Prilozhen. 19 (1985), no. 1, 25–36; English translation: Funct. Anal. Appl. 19 (1985), 21–31. (doi:10.1007/BF01086021; MR0783703)
  * [83] Vershik, A. and Yakubovich, Yu. The limit shape and fluctuations of random partitions of naturals with fixed number of summands. _Mosc. Math. J._ 1 (2001), 457–468. (doi:10.17323/1609-4514-2001-1-3-457-468; MR1877604)
  * [84] Vershik, A. and Yakubovich, Yu. Fluctuations of the maximal particle energy of the quantum ideal gas and random partitions. _Comm. Math. Phys._ 261 (2006), 759–769. (doi:10.1007/s00220-005-1434-2; MR2197546)
  * [85] Whittaker, E.T. and Watson, G.N. A Course of Modern Analysis: An Introduction to the General Theory of Infinite Processes and of Analytic Functions; with an Account of the Principal Transcendental Functions, reprint of the 4th ed. (1927). Cambridge Mathematical Library. Cambridge University Press, Cambridge, 1996. (doi:10.1017/CBO9780511608759; MR1424469)
  * [86] Wright, E.M. Asymptotic partition formulae, III. Partitions into $k$-th powers. _Acta Math._ 63 (1934), 143–191. (doi:10.1007/BF02547353; MR1555393)
  * [87] Wright, E.M. The asymptotic expansion of the generalized Bessel function. _Proc. London Math. Soc. (2)_ 38 (1935), 257–270. (doi:10.1112/plms/s2-38.1.257; MR1576315)
  * [88] Yakubovich, Yu. Ergodicity of multiplicative statistics. J. Combin. Theory Ser. A 119 (2012), 1250–1279. (doi:10.1016/j.jcta.2012.03.002; MR2915644)
  * [89] Yeh, J. Martingales and Stochastic Analysis. Series on Multivariate Analysis, 1. World Scientific, Singapore, 1995. (doi:10.1142/2948; MR1412800)
  * [90] Yong, A. Critique of Hirsch’s citation index: A combinatorial Fermi problem. _Notices Amer. Math. Soc._ 61 (2014), 1040–1050. (doi:10.1090/noti1164; MR3241560)