The article provides an extension of Boltzmann sampling to the realm of labeled combinatorial structures, enumerated by exponential generating functions and constructed using the box product. ---- [Boltzmann sampling](https://en.wikipedia.org/wiki/Boltzmann_sampler) is a general, pretty fast, technique for generate uniformely random combinatorial structures specified by structural equations. The key idea of Boltzmann sampling can be illustrated by the following example. Suppose you have a generating function $$f(x) = x + x^2 + x^2 + x^3 + x^3 + x^3 + x^3 + x^4 ...$$ corresponding to the one object of size one, two objects of size 2, four objects of size 3, etc. If you set, for example, $x = \frac{1}{2}$ then $$f\left(\frac{1}{2}\right) = \frac{1}{2} + \left(\frac{1}{2}\right)^2 + \left(\frac{1}{2}\right)^2 + \left(\frac{1}{2}\right)^3 + \left(\frac{1}{2}\right)^3 + \left(\frac{1}{2}\right)^3 + \left(\frac{1}{2}\right)^3 + \left(\frac{1}{2}\right)^4 ...$$ Now let's take $x$ such that $f(x) < \infty$, and chose an integer $k$. As you think, it quickly becomes clear that you can use this to chose a random object of size $k$. The probability equals precisely $\frac{x^k}{f(x)}$. For a given $k$ it is uniform. The elaboration of this idea in the context of structural equations (a.k.a. recursive decompositions, [symbolic method](https://en.wikipedia.org/wiki/Symbolic_method_%28combinatorics%29)) leads you to Boltzmann sampling technique. And knowing the values of generating functions together with corresponding recursive decompositions will allow you to generate a random object according to a uniform probability measure for any size. Shifting $x$ you will shift the mean size of a generating object. See [Boltzmann Samplers For The Random Generation Of Combinatorial Structures](https://www.researchgate.net/publication/2562585_Boltzmann_Samplers_For_The_Random_Generation_Of_Combinatorial_Structures) paper by Duchon, Flajolet, Louchard and Schaeffer for such elaboration.

Nice paper! Algorithm 4, line 14: remove $new$ from $x^{new}_v=x_v+r_v$. I don't understand the parallelization/scheduling part. May come back to that latter. I have an implementation of the push method here: https://github.com/maxdan94/push I will compare the speed.