arXiv:1404.3875v2 [cs.DS] 28 Apr 2014

Boltzmann samplers
for random generation of lambda terms
Pierre Lescanne
University of Lyon,
École normale supérieure de Lyon,
LIP (UMR 5668 CNRS ENS Lyon UCBL INRIA)
46 allée d’Italie, 69364 Lyon, France
pierre.lescanne@ens-lyon.fr
April 29, 2014
Abstract
Randomly generating structured objects is important in testing and
optimizing functional programs, whereas generating random λ-terms is
more specifically needed for testing and optimizing compilers. For that a
tool called QuickCheck has been proposed, but in this tool the control of
the random generation is left to the programmer. Ten years ago, a method
called Boltzmann samplers has been proposed to generate combinatorial
structures. In this paper, we show how Boltzmann samplers can be developed to generate random λ-terms, but also other random data structures
like trees. These samplers rely on a critical value which parameters the
main random selector and which is exhibited here with explanations on
how it is computed. Haskell programs are proposed to show how samplers
are actually implemented.
Keywords: lambda calculus, combinatorics, functional programming,
test, random generator, Boltzmann sampler

1

Introduction

Claessen and Hughes [5] ask a fundamental question:
[H]ow would one choose a random closed λ-term with a uniform
distribution?
Actually generating random lambda terms and more specifically generating random typable λ-terms are fundamental for debugging functional programming
compilers [22], but also for answering many questions concerning λ-calculus and
functional programming. We will see that this is possible using an apparatus called Boltzmann samplers (Duchon et al. [9]) and based on the ideas of
1

statistical physics initiated by the Austrian physicist Ludwig Boltzmann (18441906). The key issue is the ability to generate random combinatorial structures
with some flexibility on the size. In this paper the combinatorial structures are
λ-terms, but also since we are interested in functional programming they are
other data structures, like binary trees (the data type Tree in [5]) or 1-2-trees.
These latter structures will allow us to introduce the concept of Boltzmann
samplers. Assume a program generating random λ-terms is built. Somewhere
the program chooses whether it generates a variable, an abstraction or an application. This is done according to a probability distribution. Why the chosen
probability distribution is respectively
variable
0.3703026

abstraction
0.25939476

application
0.3703026

will be explained further. For now let us say that there is a critical value
ρ ≈ 0.509308127, which plays a main role in counting λ-terms and that the
above values are respectively (1 − ρ2 )/2 and ρ2 . Once we know how to generate
large plain λ-terms, we are able to address the random generation of large
typable λ-terms.
This paper is based on four milestone papers.
i [8] Nicolaas Govert de Bruijn (1972). Lambda calculus notation with nameless dummies, a tool for automatic formula manipulation, with application
to the Church-Rosser theorem. Indagationes mathematicae, 34(5), 381–
392.
ii [5] Koen Claessen and John Hughes. QuickCheck: a lightweight tool for
random testing of Haskell programs. In Martin Odersky and Philip Wadler,
editors, ICFP, pages 268–279. ACM, 2000. This paper received the Most
Influential ICFP Paper Award for the year 2000.
iii [9] Philippe Duchon, Philippe Flajolet, Guy Louchard, and Gilles Schaeffer. Boltzmann samplers for the random generation of combinatorial
structures. Combinatorics, Probability & Computing, 13(4-5):577–625,
2004.
iv [25] John Tromp. Binary lambda calculus and combinatory logic. In Marcus Hutter, Wolfgang Merkle, and Paul M. B. Vitányi, editors, Kolmogorov
Complexity and Applications, volume 06051 of Dagstuhl Seminar Proceedings. Internationales Begegnungs- und Forschungszentrum fuer Informatik
(IBFI), Schloss Dagstuhl, Germany, 2006.
The first paper established what is now called de Bruijn indices which play
an important role here in counting λ-terms up to α-conversion. But Nicolaas
de Bruijn is better known outside the Functional Programming community as
one of the pioneer of modern combinatorics [7] (see [14] fourth cover pages). We
will use indirectly his combinatorics results here.

2

The above quotation of Claessen and Hughes is extracted from the introduction of the second paper [5] and is completed by1
We have chosen to put distribution under the human tester’s
control.
However unlike Claessen and Hughes, the philosophy of our paper is that we
have chosen to put the distribution under a strict mathematical control, based
on Boltzmann samplers and probability theory as introduced in the third paper [9]. For that we use a counting method introduced by John Tromp in [25].
This amazing paper deals with algorithmic complexity theory [19] and proposes
to replace Turing machines by λ-calculus in this theory. Therefore it provides
a very simple and elegant concrete definition of descriptional complexity (Kolmogorov complexity). For this Tromp proposes to represent λ-terms by binary
chains and in this framework he defines a very small self interpreter of size 210.
His coding of λ-terms in {0, 1}∗ is as follows2 :
n
b
d
λM

=
=

1n 0
c
00M

[
M
N

=

cN
b
01M

Representing λ-terms by binary chains confers to them a natural notion of size.
This counting is realistic3 as it gives variables a weight depending on the distance
from their binder. In other words a variable which is deep in the stack is heavier
(larger) than a shallow one .
Moreover our paper relies on the book of Flajolet and Sedgewick [11] which
is the reference on generating functions. Our paper is also the application of a
recent article containing results on counting binary λ-terms [13].
In this paper we speak about random generation of λ-terms. This applies
in testing and optimizing functional programming compilers [22, 21], but it is
clear that this applies to other combinatorial objects having bound variables,
like imperative programs with a block structure [26].
A survey of Boltzmann approach
Recall that a Boltzmann sampler generates uniformly random objects with a
tolerance in the size of the generated objects. In other words, the sampler
generates the objects in a cloud around a given size.
1 See

appendix for the whole quotation.
exactly since Tromp starts indices at 0 like [16] and we starts indices at 1 like de
Bruijn [8].
3 In a previous submission [12] where we chose to give de Bruijn indices size 0, a referee
said: “If the authors want to use the de Bruijn representation, another interesting experiment
could be done: rather than to count variables as size 0, they should be counted using their
unary representation. This would penalize deep lexical scoping, which is not a bad idea since
’local’ terms are much easier to understand and analyze than deep terms”.
2 Not

3

An orographic cloud4

Before a Boltzmann sampler for combinatorial structures is defined, the generating function associated with these structures is first expressed. Its coefficients count the structures by size. In our case, i.e., trees and λ-terms, the
expressions of the generating functions contain a square root of a polynomial,
whose smallest root in module is a real which we call the critical value. This
is the radius of convergence of the associated series. A Boltzmann sampler is
derived from this series and based on a value x which allows generating random
structures with a mean and a standard deviation which depend on x. Actually
if the user has a specific size in mind, she provides the desired mean and a
calculation returns the x that yields this mean. But in our case we assume that
we attempt to generate large structures of size, say at least 100, perhaps 1000
or 1000000 or more. In our domain, trees and λ-terms, the choice of x is not
any choice, but x shall be set to the critical value, for mathematical reasons
that are explained in [9]. This way we get a uniform random sampler for all big
trees or λ-terms.
Now let us focus on trees and λ-terms and assume that the structures are
distributed in kinds. For instance, for λ-terms, these kinds are variables, abstractions and applications. A random selector among these kinds is built using
probabilities computed using the critical value. Then the selector is embedded
in a sampler which calls itself recursively. But if any freedom is left to the sampler, it may generate huge structures as well as small ones. To avoid the pitfall
of generating huge objects, the sampler has an upper limit of structures it can
generate. Furthermore after the generation a sieve retains only the structures
that are large enough, yielding the tolerance threshold (the cloud) we mentioned
previously.
Structure of the paper The paper is structured as follows. Since the problem of counting λ-terms and presenting their Boltzmann samplers is not straightforward, we present first in Section 2 the generating functions for trees: binary
trees and 1-2 trees. Then we present the generating function for counting λterms in Section 3. In Section 4 we present a less efficient method for generating
λ-terms not based on Boltzmann samplers. In the following sections we address
Boltzmann samplers. First in Section 5 we present the general perspective of
Boltzmann samplers. Then in Section 6 we specialize on trees: binary trees and
1-2 trees. Section 7 is devoted to Boltzmann samplers for λ-terms. Until that
point the paper addresses only the problem of counting and generating plain
4 Public Domain Picture, origin: Wikimedia, source: U.S. National Oceanic and Atmospheric Administration

4

λ-terms. In Section 8 we survey the harder problem of counting and generating closed terms for which we cannot propose a Boltzmann sampler, since the
generating function is not defined by a unique equation.

2

Counting trees

Before presenting how to count λ-terms, we will show in this section how to
count trees. We will first show how to count binary trees, then how to count
1-2 trees aka Motzkin trees.

2.1

Binary trees

A binary tree (the data type Tree in [5]) is either a leaf or a compound tree
made of a node and two subtrees. It can be described as follows:
B=

B × • × B.

+

We assume that leaves and nodes have size 1. General techniques of generating
functions ([11] page 738) leads to the following equation for the generating
function
X
B(z) =
Bn z n
n≥0

whose coefficients count binary trees of size n:
B(z) = z + zB(z)2
or
zB(z)2 − B(z) + z = 0
which yields:
B(z) =

1−

√

1 − 4z 2
.
2z

B(z) is not defined for z ∈ R larger that 12 . We say that B(z) has a singularity
ρB which is equal to 12 . We are going to use this fact in our sampler of random
binary trees.

2.2

Motzkin Trees

We are going to consider another family of trees which have analogies with
lambda terms, namely 1-2 trees ([11] pages 81). They are also called Motzkin
trees. They have two kinds of internal nodes, unary nodes and binary nodes.
The class of Motzkin trees is described as
M =

+

•×M

5

+

M × • × M.

unrankT oo :: Int → Integer → T erm
unrankToo n k
| k == (trompoo n) = Index $ fromIntegral (n − 1)
| k <= (trompoo (n − 2)) = Abs (unrankToo (n − 2) k )
| otherwise = unrankApp (n − 2) 0 (k − trompoo (n − 2))
where unrankApp n j r | r <= tjnj = let (dv , rm) = (r − 1) 0 divMod 0 tnj
in App (unrankToo j (dv + 1))
(unrankToo (n − j ) (rm + 1))
| otherwise = unrankApp n (j + 1) (r − tjnj )
where tnj = trompoo (n − j )
tjnj = (trompoo j ) ∗ tnj
Figure 1: The unranking function of λ-terms in Haskell
which means that a Motzkin tree is either a leaf of size 1 or a tree of size n + 1
rooted on another tree of size n or a tree of size n1 + n2 + 1 rooted on two trees
of size n1 and n2 . The generating function
M (z) =

∞
X

Mn z n

n=0

where Mn counts the number of Motzkin trees of size n, is a solution of the
equation:
M (z) = z + zM (z) + zM (z)2
or
zM (z)2 − (1 − z)M (z) + z = 0.
Hence
M (z) =

1−z−

√

1 − 2z − 3z 2
.
2z

It is given by the smallest root in module of the polynomial 1 − 2z − 3z 2 , namely
ρM = 1/3 (the other root −1 is larger in module).

3

Counting lambda terms

We now consider λ-terms up to α-conversion. To count equivalence classes
modulo α-conversion it is convenient to count canonical representatives of classes
modulo α. For that the best way is to count terms with de Bruijn indices, since
those terms are unique representatives of classes modulo α.
We work with de Bruijn indices starting at 1. Sometimes we call those
indices improperly “variables”. For instance λ1 is the representative of the αconversion equivalence class of the term λx.x and λλ1 2 is the representative of
the α-conversion equivalence class of the term λx.λy.y x. We assume that the

6

size of the index i is i + 1 and the size of an abstraction is 2 and the size of
an application is 2 as well. Therefore we get the equation for the numbers5 of
terms of size n:
S∞,0

=

S∞,1 = 0,

S∞,n+2

=

1 + S∞,n +

n
X

S∞,k S∞,n−k .

k=0

This means that there is no term of size 0 or 1 and for terms of size n+2, there is
one term
Pn which is an index namely n + 1, plus S∞,n terms that are abstractions
plus k=0 S∞,k S∞,n−k terms that are applications of terms of size k on terms
of size n − k, for k from 0 to n. Sequence (S∞,n )n∈N can be found in the On-line
Encyclopedia of Integer Sequences with the entry number A114851. Its first
20 values are:
0, 0, 1, 1, 2, 2, 4, 5, 10, 14, 27, 41,
78, 126, 237, 399, 745, 1292, 2404, 4259.
Recall some results of Grygiel-Lescanne [13]. Let S∞ (z) denote the generating
function for the sequence (S∞,n )n∈N , that is
S∞ (z)

=
=

∞
X
n=0
∞
X

S∞,n z n
z n+2 +

n=0

+

∞
X

S∞,n z n+2

n=0
∞
n
X
X
n=0

!
S∞,k S∞,n−k

z n+2

k=0

S∞ (z) fulfills the equation
S∞ (z) =

z2
+ z 2 S∞ (z) + z 2 S∞ (z)2 .
1−z

which yields
S∞ (z) =

z3 − z2 − z + 1 −

√

z6 + 2 z5 − 5 z4 + 4 z3 − z2 − 2 z + 1
.
2z 2 (1 − z)

The radius of convergence of this series i.e., the smallest singularity of the associated analytic function, given as the smallest real root of the polynomial
z6 + 2 z5 − 5 z4 + 4 z3 − z2 − 2 z + 1
5

4

3

=

2

(z − 1)(z + 3z − 2z + 2z + z − 1).
5 The index ∞ is justified by the fact that S
∞,n is the limit of Sm,n when m goes to ∞.
See Section 8.

7

A computation using the mathematical software Sage [24] tells us that this
smallest root is approximately
ρS∞ = 0.509308127.
Knowledge of the singularity ρS∞ (which we will write just ρ in what follows) is
crucial for computing quantities related to the sequence (S∞,n )n∈N . Especially
it plays a key role in computing parameters of the Boltzmann sampler. The
following theorem is in [13].
Theorem 1 (Asymptotic evaluation of S∞,n )
The number of all binary λ-terms of size n satisfies
S∞,n ∼ (1/ρ)n ·

C
,
n3/2

.
.
where 1/ρ = 1.963447954 and C = 1.021874073.
It says that whereas binary words are of size 2n , binary λ-terms are approximately of size 1.963n for large values of n. For the skeptic reader, we have
checked numerically this result with good accuracy [17].

4

Unranking lambda terms

We know how to count the λ-terms, which means basically that we know how
to order the terms, assigning a number to each term (ranking). We may use
this fact the other way around and build a λ term from its rank. This function
is called unranking and is implemented in Haskell as unrankToo (see Figure 1).
The function trompoo n 6 implements S∞,n . It works as follows. Given a number
k (1 ≤ k ≤ S∞,n ) we consider which interval it belongs to. If k = S∞,n ,
then k corresponds to a de Bruijn index, namely the de Bruijn index n − 1. If
k ≤ S∞,n−2 , then k corresponds to an abstraction, namely the abstraction of the
term obtained by unranking k in the interval [1..S∞,n−2 ]. If S∞,n−2 < k < S∞,n ,
then k corresponds to an application of two terms which are then computed.
Unranking has at least three applications:
• counting typable λ-terms,
• generating random λ-terms and
• generating random typable λ-terms.

4.1

Counting typable λ-terms

We know very little about the combinatorial properties of simply typable λterms. For instance, we do not know any formula for counting them, we have
no idea of their distribution among plain λ-terms and we have no information
6 The

notation oo mimics ∞.

8

on their ratio among plain λ-terms. We propose to count how many terms of
size n are typable. We proceed as follows, we generate all the plain terms of
size n and we filter those that are typable and we count them. This process
is relatively tedious and is limited by the number of plain terms. Here are the
numbers T∞,n of typable terms of size n we computed:
n
0
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

n
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42

T∞,n
0
0
1
1
2
2
3
5
8
13
22
36
58
103
177
307
535
949
1645
2936
5207

T∞,n
9330
16613
29921
53588
96808
174443
316267
572092
1040596
1888505
3441755
6268500
11449522
20902152
38256759
70004696
128336318
235302612
432050796
793513690
1459062947
2683714350

We did not go further than 42. Notice that
S∞,42 = 7 395 529 009

and

S∞,43 = 14 023 075 765.

In words this means that there are more than seven billions of terms of size 42
and about 14 billions of terms of size 43 which we have to traverse to know how
many are simply typable. We conjecture [17] that like the S∞,n ’s, the T∞,n ’s
increase like 1.963447954n and the ratio T∞,n /S∞,n is polynomial.

4.2

Generating random λ-terms (plain and typable)

Thanks to unranking we have a very natural method for generating random
λ-terms [20] (called by Duchon et al. [9] the recursive method ). It works as
follows: to generate a random λ-terms of size n one generates first a random
number say k in the interval [1..S∞,n ] and then by unranking k one creates the
associated λ-term. The size of the term we can generate this way is limited by
the size of the number S∞,n we can manipulate.
9

To generate a random typable λ-term of size n, one repeats the process of
generating plain terms until one gets a typable λ-term. The process is inefficient
for two reasons: one has to repeat the generation of plain terms with a lot of
waste and one has to handle huge numbers. But unlike Boltzmann samplers,
unranking generates random λ-terms of a precise size. However in applications
like compiler testing, we seldom need to generate terms of a very precise size.
An approximation on the size is enough.

5

Boltzmann samplers

In this section we present the notion of Boltzmann samplers. For more detail
the reader is invited to look at [9].
A Boltzmann sampler returns a random object in a given class C according
to a given index x. By convention a Boltzmann sampler for C depending on the
index x is written ΓC(x). In Haskell we will write it gammaC gen x , because
we need a generator of random numbers to be passed as parameter.
Before specializing the notion of Boltzmann sampler to count specifically
lambda terms of size n assume we consider a generic class C of combinatorial
structures, where objects of size n are counted by numbers Cn . Assume in
addition that the class C admits the generating function C, which means that
C(x) =

∞
X

Cn xn .

n=0

Uniform probability distribution assigns to each γ ∈ Cn the probability:
PCn (γ) =

1
Cn

We want to generate random objects with some flexibility on the size. In other
words we want the objects to be generated in some cloud around a given size n,
that is so that the size N of the objects lies in some interval (1−ε)n ≤ N ≤ (1+ε)
for some factor ε > 0 which we call “tolerance”. Such a method is called the
approximate-size uniform random generation.
The Boltzmann models assign to any object γ ∈ C the following probability:
PCn ,x (γ) =

1
· x|γ| .
C(x)

The size of an object in a Boltzmann model is a random variable N . The
probability of drawing an object of size n under the model of index x is
PCn (N = n) =
This is a probability since
X
PCn (N = n) =
n≥0

Cn xn
.
C(x)

1 X
Cn xn = 1
C(x)
n≥0

10

The random variable N has a first moment and a second moment [9]:
Ex (N ) = x

C 0 (x)
C(x)

Ex (N 2 ) =

x2 C 00 (x) + xC 0 (x)
.
C(x)

and a standard deviation:
σCn (x)

=
=

p
Ex (N 2 ) − Ex (N )2
s
x2 C 00 (x) + xC 0 (x)
C 0 (x)2
− x2
C(x)
C(x)2

One builds a Boltzmann generator for a class C according to the recursive
construction of the class C .
Disjoint union
If the class is built as C = A + B, then Cn = An + Bn and C(z) = A(z) + B(z).
The probability of the occurrence of some object is
PC ,x (γ ∈ A ) =

A(x)
,
C(x)

PC ,x (γ ∈ B) =

B(x)
C(x)

Assume we have defined a Haskell data type
data KindC = IsA | IsB
A generator for a Bernoulli variable is given by the Haskell function:
bern :: Double → StdGen → (KindC , StdGen)
bern x gen = let (p, g) = randomR (0, 1) gen
in if p < A(x )/C (x ) then (IsA, g)
else (IsB , g)
The class C is represented by the type TypeC :
data TypeC = A TypeA | B TypeB
where TypeA and TypeB are two types.7 Assume two functions gammaA and
gammaB . Then
gammaC :: StdGen → Double → (C , StdGen)
gammaC gen x = let (class, g) = bern x gen
in case class of
IsA → let (a, g1) = gammaA g x
in (A a, g1)
IsB → let (b, g1) = gammaB g x
in (B b, g1)
In what follows we will also consider not only a disjoint union of two classes
C = A + B, but a disjoint union of three classes like C = A + B + D. We call
randSel the function which replaces bern in this case.
7 It can be the case that TypeA and TypeB are TypeC itself or are defined recursively using
TypeC .

11

Cartesian product
If the class is built as C = A × B, then the generating function satisfies C(z) =
A(z) · B(z), since
X
C(z) =
z |α|+|β| .
hα,βi∈A ×B

γ = hα, βi ∈ C has the probability:
PC ,x (γ) =

x|γ|
x|α|
x|β|
=
·
.
C(x)
A(x) B(x)

Assume we have defined a Haskell data type
data TypeC = Times TypeA TypeB
In this case the Boltzmann sampler is
gammaC :: Double → StdGen → (TypeC , StdGen)
gammaC x gen = let (a, g1) = gammaA gen x
(b, g2) = gammaB g1 x
in (g2, Times a b)
Sequences
C can be the class of all finite sequences of elements of a class A . Its associated
type is:
data TypeC = Nil + Times TypeA TypeC
This means that C is the solution of the equation
C =1+A ×C
where 1 denotes the empty sequence. The associated generating function is
solution of C(z) = 1 + A(z) · C(z) which yields:
C(z) =

1
.
1 − A(z)

This means that Nil is selected with probability 1/C(x) = 1 − A(x) and a non
empty sequence is selected with probability A(x). A Boltzmann sampler for
sequences is
gammaC :: Double → StdGen → (TypeC 0 , StdGen)
gammaC x gen = let (p, g) = randomR (0, 1) gen
in if p < A(x )
then let (a, g1) = gammaA x g
(c, g2) = gammaC x g1
in (Times a c, g2)
else (Nil , g)

12

The case x = ρC
Assume that the generating function we consider is of the form:
p
PC (x) − QC (x)
C(x) =
RC (x)
where PC (x), QC (x) and RC (x) are three polynomials and where ρC is such
that QC (ρC ) = 0 and where QC (x) 6= 0 and RC (x) 6= 0 for 0 ≤ x ≤ ρC . Those
properties are fulfilled for the three generating functions B(x), M (x) and S∞ (x)
we have seen. Notice that
PC (ρC )
C(ρC ) =
RC (ρC )
is finite. On another hand
C 0 (x)

=

PC0 (x)
Q0 (x)
− p C
RC (x) 2 QC (x)RC (x)
p
0
(PC (x) − QC (x))RC
(x)
−
2
RC (x)

shows that
lim C 0 (x) = ∞.

x→ρC

Hence
lim Ex (N ) = lim

x→ρC

x→ρC

xC 0 (x)
= ∞.
C(x)

Therefore if we choose x to be ρC , the size of the generated structures will be
distributed all over the natural numbers.

6
6.1

Boltzmann samplers for trees
Boltzmann samplers for Motzkin trees

Given a number x, the mean value formula (or first moment) Ex,M (N ) =
xM 0 (x)/M (x) allows us by solving the equation xM 0 (x)/M (x) = n to tell which
value of x provides a sampler with n as mean value for the size. For instance,
xM 0 (x)/M (x) = 100 returns x100 = 0.33330833286456574 it says that the sampler gammaM gen x100 (see below) will generate a random term with a mean
value 100. Similarly since x600 = 0.3333326388880542, gammaM gen x600 will
generate a random term with mean value 600. For a small value, say 10 we get
x1 = 0.3308286281723805. The standard deviation σMn (x) can be instantiated
to x10 or x100 or to x600 . Actually
σMn (x10 ) ≈ 25 σMn (x100 ) ≈ 816

σMn (x600 ) ≈ 1200.

Added to the fact that x100 and x600 are very close to the ρM = 1/3 (which is
called the critical value) it appears that a Boltzmann sampler for large values
13

of n is obtained by parameterizing gammaM with ρM , namely 1/3. See [9]
Section 7.2 for a mathematical justification of this choice.
Since there are three components in the addition, we generalize the Bernoulli
generator and then we specialize it to the case of Motzkin trees, producing a
function we call randSelM . Let us call
M0 (z)
z
=
M (z)
M (z)
the first component of (z + zM (z) + zM (z)2 )/M (z),
M1 (z)
=z
M (z)
the second component of (z + zM (z) + zM (z)2 )/M (z) and
M2 (z)
= zM (z)
M (z)
the third component of (z + zM (z) + zM (z)2 )/M (z). We define the following
Haskell functions:
m z = (1 − z − sqrt(1 − 2 ∗ z − 3 ∗ z ∗ z ))/(2 ∗ z )
p1 x = x / m x
p2 x = p1 x + x
One may notice that p2 x + x ∗ m x == 1 for every value of x, which means
that when one adds all the probabilities of the events one gets 1.
data KindMotzkin = IsLeaf | IsUnary | IsBinary
randSelM :: StdGen → Double → (KindMotzkin, StdGen)
randSelM gen x = let (p, g) = randomR (0, 1) gen
in if p < p1 x then (IsLeaf , g)
else if p < p2 x then (IsUnary, g)
else (IsBinary, g)
Then we define the data type
data MTree = MNil | MU MTree | MB MTree MTree
for Motzkin trees and we get the Boltzmann sampler for Motzkin trees:
gammaM :: StdGen → Double → (MTree, StdGen)
gammaM gen x = let (kind , g) = randSelM gen x
in case kind of
IsLeaf → (MNil , g)
IsUnary → let (t1, g1) = gammaM g x
in (MU t1, g1)
IsBinary → let (t1, g1) = gammaM g x
(t2, g2) = gammaM g1 x
in (MB t1 t2, g2)
14

ceiledGammaM :: StdGen → (EMTree, Int, StdGen)
ceiledGammaM gen =
let (kind , g0) = selOneThird gen
in case kind of
IsLeaf → (OK MNil , 1, g0)
IsUnary → let (et1, n1, g1) = ceiledGammaM g0
in case et1 of
OK t1 → if n1 <= upLimit
then (OK (MU t1), n1 + 1, g1)
else (Error , upLimit, g1)
Error → (Error , upLimit, g1)
IsBinary → let (et1, n1, g1) = ceiledGammaM g0
(et2, n2, g2) = ceiledGammaM g1
in case et1 of
OK t1 → case et2 of
OK t2 → if n1 + n2 + 1 <= upLimit
then (OK (MB t1 t2), n1 + n2 + 2, g2)
else (Error , upLimit, g2)
Error → (Error , upLimit, g2)
Error → (Error , upLimit, g2)
Figure 2: The function ceiledGammaM
An interesting feature should be noticed in this program. Indeed the function
gammaM gen x calls gammaM g x and should not terminate since we can prove
in no way that the first arguments namely gen versus g decrease. Thus the
termination is only probabilistic and depends on the choice of x. The implementation of randSelM makes gammaM gen x to terminate with probability 1,
thanks to the whole theory of Boltzmann samplers. For instance if the generator is buggy and takes its values only in the interval ( 31 , 1] the program will not
terminate.
Samplers for large Motzkin trees
Note that M (ρM ) = 1. For large Motzkin trees the right choice is to take
x = ρM = 1/3. This gives the probability
• for drawing a leaf

ρM
M (ρM )

= 31 ,

• for drawing a tree rooted by a unary node ρM = 13 ,
• for drawing a tree rooted by a binary node ρM M (ρM ) = 13 .
Actually since choosing ρM generates large trees, it is advisable to limit the
size of the tree to be generated during the process. Therefore, we consider what
Duchon et al. call ceiled random generation. We define a constant upLimit which
15

sets a size upper limit that the sampler should respect during the generation.
We then define a data type:
data EMTree = OK MTree | Error
and a selection function selOneThird which selects a kind among IsLeaf , IsUnary
and IsBinary with probability 1/3, 1/3, 1/3. The function ceiledGammaM is
defined in Figure 2.
If we want to generate a term between say a and b we generate terms with
an upLimit set to b until we get a term larger than a. Duchon et al. have shown
that the complexity of this method is linear. In other words, the generation of
a tree with a size between a and b takes a time 0(b).

6.2

Boltzmann samplers for binary trees

The binary trees are implemented in Haskell by
data BTree = BNil | BNode BTree BTree
If we compute for B(z) the means associated to some values of n, a we did above
for M (z), we get:
x10

=

0.4960783708241233

x100

=

0.4999739685017913

x600

=

0.49999930090217926

for which we obtain the standard deviations:
σBn (x10 ) ≈ 22

σBn (x100 ) ≈ 970

σBn (x600 ) ≈ 14623.

There are two kinds of binary trees:
• leafs (represented by BNil ) which correspond to
z
B0 (z)
=
B(z)
B(z)
which is the first component of (z + zB(z)2 )/B(z), and
• nodes (represented by BNode t1 t2) which correspond to
B1 (z)
= zB(z)
B(z)
which is the second component of (z + zB(z)2 )/B(z).
Let us define
b z = (1 − sqrt(1 − 4 ∗ z ∗ z ))/ (2 ∗ z )
16

q x = x /bx
We introduce
data KindBTree = IsBNil | IsBNode
and
randSelB :: StdGen → Double → (KindBTree, StdGen)
randSelB gen x = let (p, g) = randomR (0, 1) gen
in if p < q x
then (IsBNil , g)
else (IsBNode, g)
We can then introduce the sampler:
gammaB :: StdGen → Double → (BTree, StdGen)
gammaB gen x = let (kind , g) = randSelB gen x
in case kind of
IsBNil → (BNil , g)
IsBNode → let (t1, g1) = gammaB g x
(t2, g2) = gammaB g1 x
in (BNode t1 t2, g2)
Once again we see that the best value for x to generate random large binary
trees is ρB = 21 . Thus gammaB gen 0.5 generates large binary trees. For a
similar approach of the same problem with in particular a discussion of ceiled
random generation see [27]. For other methods for generating random binary
trees, see [15], pages 18-19 and [23].

7

Boltzmann samplers for lambda terms

Like for Motzkin trees we consider the equation of the generating function:
S∞ (z) =

z2
+ z 2 S∞ (z) + z 2 S∞ (z)2 .
1−z

with three components: the first corresponds to de Bruijn indices, the second to
abstractions, the third to applications. Like for Motzkin trees we have to build
a random selector among three probabilities. First we describe a data type:
data KindTerm = IsVariable | IsAbstraction | IsApplication
a function soo
soo z = (z 3 − z 2 − z + 1 − sq)/(2 ∗ z ∗ z ∗ (1 − z ))
where sq = sqrt(z 6 + 2 ∗ (z 5 ) − 5 ∗ (z 4 ) + 4 ∗ (z 3 ) − z 2 − 2 ∗ z + 1)
and two functions:
17

p1 x = x ∗ x / (1 − x ) / soo x
p2 x = p1 x + x 2
Using Sage we computed the values:
x100

=

0.5092252666102192

x500

=

0.5093048407797965

x600

=

0.5093058457062517

x1000

=

0.5093073063214039

which are the choices of x yielding mean values 100, 500, 600 and 1000 respectively.
General samplers of λ-terms
The values of the probabilities for a given x are
• pv (x) =

x2
(1−x)S∞ (x)

for variables,

• pabs (x) = x2 for abstractions,
• papp (x) = x2 S∞ (x) for applications.
For selecting among IsVariable, IsAbstraction and IsApplication we use the
following selector written in Haskell
randSel :: StdGen → Double → (KindTerm, StdGen)
randSel gen x = let (p, g) = randomR (0, 1) gen
in if p < pv x then (IsVariable, g)
else if p < pv x + pabs x
then (IsAbstraction, g)
else (IsApplication, g)
Samplers for large λ-terms
Like for trees, the best choice of x for generating large λ-terms is ρ which we
call rho in Haskell:
rho :: Double
rho = 0.5093081270242373.
whose square is
rhosquare = rho ∗ rho
which yields ρ2 = 0.25939476825293667. Notice that since ρ is a root of the
2
polynomial below the square root, S∞ (ρ) = 1−ρ
2ρ2 . The values of the probabilities
for selecting among variables, abstractions and applications are:
18

• pv (ρ) =

2ρ4
(1−ρ)(1−ρ2 )

for variables,

• pabs (ρ) = ρ2 for abstractions,
• papp (ρ) =
Let us simplify

1−ρ2
2

for applications.

2ρ4
(1−ρ)(1−ρ2 )

into

1−ρ2
2

2ρ4
1 − ρ2
−
2
(1 − ρ)(1 − ρ )
2

by computing the difference:
4ρ4 − (1 − ρ2 )2 (1 − ρ)
2(1 − ρ)(1 − ρ2 )
ρ5 + 3ρ4 − 2ρ3 + 2ρ2 + ρ − 1
2(1 − ρ)(1 − ρ2 )
0.

=
=
=

Therefore we get the result announced in the introduction namely that
• pv (ρ) =

1−ρ2
2

≈ 0.3703026 for variables,

2

• pabs (ρ) = ρ ≈ 0.25939476 for abstractions,
• papp (ρ) =

1−ρ2
2

≈ 0.3703026 for applications.

We build the random selector:
pvrho = (1 − rhosquare)
p2rho = pvrho + rhosquare
randSelRho :: StdGen → (KindTerm, StdGen)
randSelRho gen = let (p, g) = randomR (0, 1) gen
in if p < pvrho then (IsVariable, g)
else if p < p2rho
then (IsAbstraction, g)
else (IsApplication, g)
We generate a Boltzmann sampler for indices (or variables). A de Bruijn index
is a natural number i.e., a sequence of 1’s. Hence we take the Boltzmann sampler
for sequences with A(z) = z and A(ρ) = ρ.
gammaV :: StdGen → (Integer , StdGen)
gammaV gen = let (a, g) = randomR (0, 1) gen
in if a < rho
then let (n, g1) = gammaV g
in (n + 1, g1)
else (1, g)
We want to generate terms that are below a certain upLimit. Thus when the
limit is passed we generate an error. For that we create a type with error.
data ETerm = OK Term | Error
19

This method recalls the sized generation of [5]. A Boltzmann sampler for large
λ-terms which is given:
ceiledGammaSoo :: StdGen → (ETerm, Int, StdGen)
ceiledGammaSoo gen =
let (kind , g0) = randSelRho gen
in case kind of
IsVariable → let (n, g1) = gammaV g0
in if fromIntegral n <= upLimit
then (OK (Index n), fromIntegral (n + 1), g1)
else (Error , upLimit, g1)
IsAbstraction → let (et1, n1, g1) = ceiledGammaSoo g0
in case et1 of
OK t1 → if n1 + 2 <= upLimit
then (OK (Abs t1), n1 + 2, g1)
else (Error , upLimit, g1)
Error → (Error , upLimit, g1)
IsApplication → let (et1, n1, g1) = ceiledGammaSoo g0
(et2, n2, g2) = ceiledGammaSoo g1
in case et1 of
OK t1 → case et2 of
OK t2 → if n1 + n2 + 2 <= upLimit
then (OK (App t1 t2), n1 + n2 + 2, g2)
else (Error , upLimit, g2)
Error → (Error , upLimit, g2)
Error → (Error , upLimit, g2)
To generate a large plain λ-term, terms are filtered until a term that is large
enough is generated. Recall that the method is linear in time complexity. Thus
the generation of a term of size 100, 000 takes a few seconds, the generation of
a term of size one million takes three minutes and the generation of a term of
size five millions takes five minutes on a laptop.
To generate large typable λ-terms we generate λ-terms and check their typability. Currently we are able to generate random typable λ-terms of size 500.

8

Ranking and unranking closed lambda terms

Boltzmann samplers work only when one knows the generating function explicitly. Thus if one wants to generate typable closed terms, one has two solutions.
First one generates a random typable term using the Boltzmann sampler described above and if the produced term is not closed, one closes it by adding
the adequate number of abstractions. Second one uses an unranking method
similar to this described in Section 4.2. For that let us recall how one counts
terms with at most m free indices [13]. The specific case of the closed terms is
the case when the number of free indices is at most 0.
Sm,0

= Sm,1 = 0,
20

Sm,n+2

=

[m ≥ n + 1] + Sm+1,n +

n
X

Sm,k Sm,n−k .

k=0

Given a predicate P , [P (~x)] denotes the Iverson symbol, i.e., [P (~x)] = 1 if P (~x)
and [P (~x)] = 0 if ¬P (~x). Figure 3 gives a program for unranking terms with
at most m indices, where tromp m n is the Haskell function representing Sm,n .
unrankT :: Int → Int → Integer → Term
unrankT m n k
| m >= n − 1 && k == (tromp m n) = Index $ fromIntegral (n − 1) — terms 1n−1 0
| k <= (tromp (m + 1) (n − 2)) = Abs (unrankT (m + 1) (n − 2) k ) — terms 00M
| otherwise = unrankApp (n − 2) 0 (k − tromp (m + 1) (n − 2)) — terms 01M N
where unrankApp n j h
| h <= tmjtmnj = let(dv , rm) = (h − 1) 0 divMod 0 tmnj
in App (unrankT m j (dv + 1)) (unrankT m (n − j ) (rm + 1))
| otherwise = unrankApp n (j + 1) (h − tmjtmnj )
where tmnj = tromp m (n − j )
tmjtmnj = (tromp m j ) ∗ tmnj

Figure 3: unrankT function in Haskell.
S0,n counts closed terms and is computed by tromp 0 n.

9

Related works

In the introduction we cited papers that are clearly connected to this work. In a
recent work, Bacher et al. [1] propose an improved random generation of binary
trees and Motzkin trees, based on Rémy algorithm [23] (or algorithm R as Knuth
calls it [15]). Instead of growing the trees from the root, they propose like Rémy
to grow the trees from inside by an operation called grafting. It is not clear how
this can be generalized to λ-terms as this requires “to find a combinatorial
interpretation for the holonomic equations [which] is not [...] always possible,
and even for simple combinatorial objects this is not elementary” (Conclusion
of [1] page 16).
We would like also to mention papers on counting λ-terms [18, 12] and
evaluating their combinatorial proprieties namely [2, 6, 3, 4]. Another related
paper is [10] which proposes Haskell programs for enumerating structures.

10

Acknowledgments

We would like to thank Katarzyna Grygiel, Bruno Salvy and John Tromp, for
fruitful discussions. We are also indebted to Neil Sloane and the developers of

21

the On-line Encyclopedia of Integer Sequences which allowed us to know about
sequence A114851 which was the starting point of this research.

11

Conclusion

Boltzmann samplers are central tools for the uniform generation of random
structures, like trees or λ-terms. Two directions are now open for applications:
first to integrate the programs proposed here in actual testers and optimizers and
second to extend Boltzmann samplers to other kinds of programs, for instance
programs with block structures. From the theoretical point of view, more should
be known about generating functions for closed λ-terms or λ-terms with fixed
bound on the number of variables. Boltzmann samplers should be designed
for such terms which requires to extend the theory. Concerning combinatorial
properties of simply typable λ-terms many question are left open and seem to
be hard. Besides since we are interested in generating typable terms, it could
be worth to build random uniform samplers delivering directly typable terms,
for instance based on a hashing table built offline.

References
[1] A. Bacher, O. Bodini, and A. Jacquot. Efficient random sampling of binary
and unary-binary trees via holonomic equations. CoRR, abs/1401.1140,
2014.
[2] O. Bodini, D. Gardy, and B. Gittenberger. Lambda-terms of bounded unary
height. 2011 Proceedings of the Eighth Workshop on Analytic Algorithmics
and Combinatorics (ANALCO), 2011.
[3] O. Bodini, D. Gardy, B. Gittenberger, and A. Jacquot. Enumeration of
generalized BCI lambda-terms. Electr. J. Comb., 20(4):P30, 2013.
[4] O. Bodini, D. Gardy, and A. Jacquot. Asymptotics and random sampling
for BCI and BCK lambda terms. Theor. Comput. Sci., 502:227–238, 2013.
[5] K. Claessen and J. Hughes. QuickCheck: a lightweight tool for random
testing of Haskell programs. In M. Odersky and P. Wadler, editors, ICFP,
pages 268–279. ACM, 2000.
[6] R. David, K. Grygiel, J. Kozik, C. Raffalli, G. Theyssier, and M. Zaionc.
Asymptotically almost all λ-terms are strongly normalizing. Logical Methods in Computer Science, 9(1:02):1–30, 2013.
[7] N. G. de Bruijn. Asymptotic Methods in Analysis. North-Holland, Amsterdam, 1958.
[8] N. G. de Bruijn. Lambda calculus with nameless dummies, a tool for
automatic formula manipulation, with application to the Church-Rosser
22

theorem. Proc. Koninkl. Nederl. Akademie van Wetenschappen, 75(5):381–
392, 1972.
[9] P. Duchon, P. Flajolet, G. Louchard, and G. Schaeffer. Boltzmann samplers
for the random generation of combinatorial structures. Combinatorics,
Probability & Computing, 13(4-5):577–625, 2004.
[10] J. Duregård, P. Jansson, and M. Wang. Feat: functional enumeration of
algebraic types. In J. Voigtländer, editor, Proceedings of the 5th ACM
SIGPLAN Symposium on Haskell, Haskell 2012, Copenhagen, Denmark,
13 September 2012, pages 61–72. ACM, 2012.
[11] P. Flajolet and R. Sedgewick. Analytic Combinatorics. Cambridge University Press, 2008.
[12] K. Grygiel and P. Lescanne. Counting and generating lambda terms. J.
Funct. Program., 23(5):594–628, 2013.
[13] K. Grygiel and P. Lescanne. Counting terms in the binary lambda calculus.
CoRR, abs/1401.0379, 2014.
[14] D. E. Knuth. Selected Papers on Analysis of Algorithms, volume 102 of
CSLI Lecture Notes. Stanford, California: Center for the Study of Language
and Information, 2000.
[15] D. E. Knuth. The Art of Computer Programming, Volume 4, Fascicle 4:
Generating All Trees, History of Combinatorial Generation (Art of Computer Programming). Addison-Wesley Professional, 2006.
[16] P. Lescanne. From λσ to λυ, a journey through calculi of explicit substitutions. In H. Boehm, editor, Proceedings of the 21st Annual ACM Symposium on Principles Of Programming Languages, Portland (Or., USA),
pages 60–69. ACM, 1994.
[17] P. Lescanne. An exercise on streams: convergence acceleration. CoRR,
abs/1312.4917, 2013.
[18] P. Lescanne. On counting untyped lambda terms. Theor. Comput. Sci.,
474:80–97, 2013.
[19] M. Li and P. Vitányi. An introduction to Kolmogorov complexity and its
applications (2nd ed.). Springer-Verlag New York, Inc., Secaucus, NJ, USA,
1997.
[20] A. Nijenhuis and H. S. Wilf. Combinatorial algorithms for computers and
calculators. Computer science and applied mathematics. Academic Press,
New York, 1978. First ed. published in 1975 under title: Combinatorial
algorithms.

23

[21] M. Palka. Testing an optimising compiler by generating random lambda
terms. Licentiatavhandling, Department of Computer Science and Engineering, Chalmers University of Technology and Göteborg University, may
2012.
[22] M. H. Palka, K. Claessen, A. Russo, and J. Hughes. Testing an optimising
compiler by generating random lambda terms. In Proceedings of the 6th
International Workshop on Automation of Software Test, AST ’11, pages
91–97, New York, NY, USA, 2011. ACM.
[23] J.-L. Rémy. Un procédé itératif de dénombrement d’arbres binaires et son
application à leur génération aléatoire. ITA, 19(2):179–195, 1985.
[24] W. Stein et al. Sage Mathematics Software (Version 5.8). The Sage Development Team, 2012. http://www.sagemath.org.
[25] J. Tromp. Binary lambda calculus and combinatory logic. In M. Hutter, W. Merkle, and P. M. B. Vitányi, editors, Kolmogorov Complexity
and Applications, volume 06051 of Dagstuhl Seminar Proceedings. Internationales Begegnungs- und Forschungszentrum fuer Informatik (IBFI),
Schloss Dagstuhl, Germany, 2006.
[26] X. Yang, Y. Chen, E. Eide, and J. Regehr. Finding and understanding bugs
in C compilers. In M. W. Hall and D. A. Padua, editors, Proceedings of the
32nd ACM SIGPLAN Conference on Programming Language Design and
Implementation, PLDI 2011, San Jose, CA, USA, June 4-8, 2011, pages
283–294. ACM, 2011.
[27] B. Yorgey. Random binary trees with a size-limited critical Boltzmann
sampler. http://byorgey.wordpress.com/2013/04/25/, April 2013.
Here is the whole quotation of Claessen and Hughes [5]:
We have chosen to put distribution under the human tester’s
control, by defining a test data generation language (also embedded
in Haskell), and a way to observe the distribution of test cases. By
programming a suitable generator, the tester can not only control the
distribution of test cases, but also ensure that they satisfy arbitrarily
complex invariants.

24

