diff --git a/src/doc/en/reference/combinat/module_list.rst b/src/doc/en/reference/combinat/module_list.rst
index 859be86ceae..ebfa7b6317f 100644
--- a/src/doc/en/reference/combinat/module_list.rst
+++ b/src/doc/en/reference/combinat/module_list.rst
@@ -25,6 +25,10 @@ Comprehensive Module list
     sage/combinat/baxter_permutations
     sage/combinat/binary_recurrence_sequences
     sage/combinat/binary_tree
+    sage/combinat/boltzmann_sampling/__init__
+    sage/combinat/boltzmann_sampling/generator
+    sage/combinat/boltzmann_sampling/grammar
+    sage/combinat/boltzmann_sampling/oracle
     sage/combinat/cartesian_product
     sage/combinat/catalog_partitions
     sage/combinat/chas/__init__
diff --git a/src/sage/combinat/all.py b/src/sage/combinat/all.py
index 8e298c898c0..ef17f199356 100644
--- a/src/sage/combinat/all.py
+++ b/src/sage/combinat/all.py
@@ -177,6 +177,9 @@
 
 from .species.all import *
 
+# Boltzmann sampling
+from .boltzmann_sampling.all import *
+
 lazy_import('sage.combinat.kazhdan_lusztig', 'KazhdanLusztigPolynomial')
 
 lazy_import('sage.combinat.degree_sequences', 'DegreeSequences')
diff --git a/src/sage/combinat/boltzmann_sampling/__init__.py b/src/sage/combinat/boltzmann_sampling/__init__.py
index e69de29bb2d..cea13eb6f7e 100644
--- a/src/sage/combinat/boltzmann_sampling/__init__.py
+++ b/src/sage/combinat/boltzmann_sampling/__init__.py
@@ -0,0 +1,8 @@
+r"""
+Boltzmann Sampling
+==================
+
+- :ref:`sage.combinat.boltzmann_sampling.grammar`
+- :ref:`sage.combinat.boltzmann_sampling.oracle`
+- :ref:`sage.combinat.boltzmann_sampling.generator`
+"""
diff --git a/src/sage/combinat/boltzmann_sampling/all.py b/src/sage/combinat/boltzmann_sampling/all.py
new file mode 100644
index 00000000000..682e6397b61
--- /dev/null
+++ b/src/sage/combinat/boltzmann_sampling/all.py
@@ -0,0 +1,8 @@
+r"""
+Boltzmann sampling features that are imported by default in the interpreter
+namespace
+"""
+
+from .generator import Generator, UnionBuilder
+from .grammar import Atom, Grammar, Product, Ref, Rule, Union
+from .oracle import OracleFromFunctions, SimpleOracle
diff --git a/src/sage/combinat/boltzmann_sampling/generator.pyx b/src/sage/combinat/boltzmann_sampling/generator.pyx
index 38bbb394b3f..827fabfc0c2 100644
--- a/src/sage/combinat/boltzmann_sampling/generator.pyx
+++ b/src/sage/combinat/boltzmann_sampling/generator.pyx
@@ -3,14 +3,15 @@ r"""Boltzmann generator for Context-free grammars.
 
 This module provides functions for generating combinatorial objects (i.e.
 objects described by a combinatorial specification, see
-:ref:`sage.combinat.bolzmann_sampling.grammar`) according to the Boltzmann
+:mod:`sage.combinat.boltzmann_sampling.grammar`) according to the Boltzmann
 distribution.
 
 Given an unlabelled combinatorial class A, the Boltzmann distribution of
-parameter x is such that an object of size n is drawn with the probability ``x^n
-/ A(x)`` where ``A(x)`` denotes the ordinary generating function of A. For
-labelled classes, this probability is set to ``x^n / (n! * A(x))`` where A(x)
-denotes the exponential generating function of A. See [DuFlLoSc04] for details.
+parameter `x` is such that an object of size n is drawn with the probability
+`\frac{x^n}{A(x)}` where `A(x)` denotes the ordinary generating function of A.
+For labelled classes, this probability is set to `\frac{x^n}{n!A(x)}` where
+`A(x)` denotes the exponential generating function of A. See [DuFlLoSc04] for
+details.
 
 By default, the objects produced by the generator are nested tuples of strings
 (the atoms). For instance ``('z', ('z', 'e', 'e'), ('z', 'e', 'e'))`` is a
@@ -22,22 +23,23 @@ generation. For instance, in order to generate Dyck words using the grammar for
 binary trees, one can use a builder that return ``""`` for each leaf ``"(" +
 left child + ")" + right child`` for each node. The builders will receive a
 tuple for each product, a string for each atom and builders for unions should be
-computed using the ``UnionBuilder`` helper. See the example below for the case
-of Dyck words.
+computed using the :func:`UnionBuilder` helper. See the example below for the
+case of Dyck words.
 
 EXAMPLES::
 
-    sage: from sage.combinat.boltzmann_sampling.grammar import *
-    sage: from sage.combinat.boltzmann_sampling.generator import *
     sage: leaf = Atom("e", size=0)
     sage: z = Atom("z")
     sage: grammar = Grammar(rules={"B": Union(leaf, Product(z, "B", "B"))})
     sage: generator = Generator(grammar)
+
     sage: def leaf_builder(_):
     ....:     return ""
+
     sage: def node_builder(tuple):
     ....:     _, left, right = tuple
     ....:     return "(" + left + ")" + right
+
     sage: generator.set_builder("B", UnionBuilder(leaf_builder, node_builder))
     sage: dyck_word, _ = generator.gen("B", (10, 20))
     sage: dyck_word  # random
@@ -49,23 +51,24 @@ the tree.
 
 EXAMPLES::
 
-    sage: from sage.combinat.boltzmann_sampling.grammar import *
-    sage: from sage.combinat.boltzmann_sampling.generator import *
     sage: leaf = Atom("e", size=0)
     sage: z = Atom("z")
     sage: grammar = Grammar(rules={"B": Union(leaf, Product(z, "B", "B"))})
     sage: generator = Generator(grammar)
+
     sage: def leaf_height(_):
     ....:     return 0
+
     sage: def node_height(tuple):
     ....:     _, left, right = tuple
     ....:     return 1 + max(left, right)
+
     sage: generator.set_builder("B", UnionBuilder(leaf_height, node_height))
     sage: height, _ = generator.gen("B", (10, 20))
     sage: height  # random
     6
 
-REFERENCES::
+REFERENCES:
 
 .. [DuFlLoSc04] Philippe Duchon, Philippe Flajolet, Guy Louchard, and Gilles
    Schaeffer. 2004. Boltzmann Samplers for the Random Generation of
@@ -177,6 +180,7 @@ cdef c_generate(int id, float weight, rules, builders, randstate rstate):
     generated = []
     cdef list todo = [(REF, weight, id)]
     cdef float r = 0.
+
     while todo:
         type, weight, args = todo.pop()
         if type == REF:
@@ -217,8 +221,9 @@ cdef c_generate(int id, float weight, rules, builders, randstate rstate):
     obj, = generated
     return obj
 
-
 cdef c_gen(int id, float weight, rules, int size_min, int size_max, int max_retry, builders):
+    """Search for a tree in a given size window. Wrapper around c_simulate and
+    c_generate."""
     cdef int nb_rejections = 0
     cdef int cumulative_rejected_size = 0
     cdef int size = -1
@@ -259,23 +264,26 @@ cdef inline identity(x):
     return x
 
 def UnionBuilder(*builders):
-    """"Helper for computing the builder of a union out of the builders of
-    its components. For instance, in order to compute the height of a binary
-    tree on the fly:
+    """Helper for computing the builder of a union out of the builders of its
+    components in Boltzmann samplers. For instance, in order to compute the
+    height of a binary tree on the fly:
 
     EXAMPLES::
 
         sage: # Grammar: B = Union(leaf, Product(z, B, B))
-        sage: from sage.combinat.boltzmann_sampling.generator import UnionBuilder
+
         sage: def leaf_builder(_):
         ....:     return 0
+
         sage: def node_builder(tuple):
         ....:     _, left, right = tuple
         ....:     return 1 + max(left, right)
+
         sage: builder = UnionBuilder(leaf_builder, node_builder)
         sage: choice_number = 0
         sage: builder((choice_number, "leaf"))
         0
+
         sage: choice_number = 1
         sage: builder((choice_number, ('z', 37, 23)))
         38
@@ -293,7 +301,9 @@ cdef inline ProductBuilder(builders):
     return build
 
 cdef make_default_builder(rule):
-    """Generate the default builders for a rule"""
+    """Generate the default builders for a rule.
+
+    For use with Boltzmann samplers :mod:`sage.combinat.boltzmann_sampling`"""
     type, __, args = rule
     if type == REF:
         return identity
@@ -312,9 +322,10 @@ class Generator:
     def __init__(self, grammar, oracle=None):
         """Make a Generator out of a grammar.
 
-        INPUT::
+        INPUT:
 
         - ``grammar`` -- a combinatorial grammar
+
         - ``oracle`` (default: None) -- an oracle for the grammar. If not
           supplied, a default generic oracle is automatically generated.
         """
@@ -338,9 +349,10 @@ class Generator:
     def set_builder(self, non_terminal, func):
         """Set the builder for a non-terminal symbol.
 
-        INPUT::
+        INPUT:
 
         - ``non_terminal`` -- string, the name of the non-terminal symbol
+
         - ``func`` -- function, the builder
         """
         symbol_id = self.name_to_id[non_terminal]
@@ -349,7 +361,7 @@ class Generator:
     def get_builder(self, non_terminal):
         """Retrieve the current builder for a non-terminal symbol.
 
-        INPUT::
+        INPUT:
 
         - ``non_terminal`` -- string, the name of the non-terminal symbol
         """
@@ -359,13 +371,15 @@ class Generator:
     def gen(self, name, window, max_retry=2000):
         """Generate a term of the grammar in a given size window.
 
-        INPUT::
+        INPUT:
 
         - ``name`` -- string, the name of the symbol of the grammar you want to
           generate
+
         - ``window`` -- pair of integers, the size of the generated object will
           be greater than the first component of the window and lower than the
           second component
+
         - ``max_retry`` (default: 2000) -- integer, maximum number of attempts.
           If no object in the size window is found in less that ``max_retry``
           attempts, the generator returns None
diff --git a/src/sage/combinat/boltzmann_sampling/grammar.py b/src/sage/combinat/boltzmann_sampling/grammar.py
index 135a18c8613..ca97d48e748 100644
--- a/src/sage/combinat/boltzmann_sampling/grammar.py
+++ b/src/sage/combinat/boltzmann_sampling/grammar.py
@@ -2,35 +2,35 @@
 Context-free grammars for Boltzmann generation.
 
 Grammars use the basic operators of the symbolic method of analytic
-combinatorics (currently ``+``, ``*`` and atoms) to specify unlabelled
-combinatorial classes. For instance, binary tree can be specified by
-``B = leaf + Z * B * B`` which, using the syntax implemented in this module,
-looks like:
+combinatorics to specify labelled or unlabelled combinatorial classes. For
+instance, binary tree can be specified by ``B = leaf + Z * B * B`` which, using
+the syntax implemented in this module, looks like:
 
-EXAMPLES::
+EXAMPLE::
 
-    sage: from sage.combinat.boltzmann_sampling.grammar import *
     sage: z = Atom("z")
     sage: leaf = Atom("leaf", size=0)
-    sage: bintree = Grammar(rules={"B": Union(leaf, Product(z, "B", "B"))})
+    sage: Grammar(rules={"B": Union(leaf, Product(z, "B", "B"))})
+    B -> Union(leaf, Product(z, B, B))
 
 Grammars are not limited to a single rule:
 
-EXAMPLES::
+EXAMPLE::
 
-    sage: from sage.combinat.boltzmann_sampling.grammar import *
     sage: z = Atom("z")
     sage: leaf = Atom("leaf", size=0)
-    sage: planetree = Grammar(rules={
+    sage: Grammar(rules={
     ....:     "T": Product(z, "S"),
     ....:     "S": Union(leaf, Product("T", "S")),
     ....: })
+    S -> Union(leaf, Product(T, S))
+    T -> Product(z, S)
 
 
-Note that at the moment, we only support unlabelled classes.
-
 AUTHORS:
+
 - Matthieu Dien (2019): initial version
+
 - Martin Pépin (2019): initial version
 """
 
@@ -47,7 +47,22 @@ class Rule(SageObject):
 
 
 class Atom(Rule):
-    """Grammar atoms (terminal symbols)."""
+    """Terminal symbol of a grammar.
+
+    EXAMPLES::
+
+        sage: z = Atom("z")
+        sage: z
+        z
+
+        sage: z = Atom("z", size=4)
+        sage: z
+        z^4
+
+        sage: x = Atom("x", size=0)
+        sage: x
+        x^0
+    """
 
     def __init__(self, name, size=1):
         """Create an Atom.
@@ -65,7 +80,6 @@ def _latex_(self):
         r"""Return the LaTeX representation of an atom.
 
         EXAMPLES::
-            sage: from sage.combinat.boltzmann_sampling.grammar import Atom
 
             sage: z = Atom("z")
             sage: latex(z)
@@ -95,8 +109,10 @@ def _latex_(self):
             return "{}^{}".format(nice_name, self.size)
 
     def _repr_(self):
-        # XXX. shall it be copy-pastable?
-        return self.name
+        if self.size != 1:
+            return "{}^{}".format(self.name, self.size)
+        else:
+            return self.name
 
     def _to_combsys(self):
         if self.size > 0:
@@ -106,10 +122,11 @@ def _to_combsys(self):
 
 
 class Ref(Rule):
-    """References to non terminal symbols.
+    """Non terminal symbols of a grammar.
 
-    Instances of this class reprensent recursive references to non-terminal
-    symbols inside grammar rules.
+    Instances of this class represent recursive references to non-terminal
+    symbols inside grammar rules. In general you should not use this class
+    directly.
     """
 
     def __init__(self, name):
@@ -125,7 +142,6 @@ def _latex_(self):
         r"""Return the LaTeX representation of a non-terminal symbol.
 
         EXAMPLES::
-            sage: from sage.combinat.boltzmann_sampling.grammar import Ref
 
             sage: latex(Ref("X"))
             X
@@ -156,10 +172,17 @@ def _to_rule(r):
 class Union(Rule):
     """Union of two or more rules.
 
-    Union(A, B, C) corresponds to the following grammar in BNF syntax:
-    _ ::= A
-        | B
-        | C
+    D = Union(A, B, C) corresponds to the following grammar in BNF syntax:
+    ``D ::= A | B | C``
+
+    EXAMPLES::
+
+        sage: Union("A", "B", "C")
+        Union(A, B, C)
+
+        sage: z = Atom("z")
+        sage: Union(z, "A")
+        Union(z, A)
     """
 
     def __init__(self, *args):
@@ -168,17 +191,6 @@ def __init__(self, *args):
         INPUT:
 
         - ``args`` -- list of strings or Rules; strings are interpreted as Refs
-
-        EXAMPLES::
-            sage: from sage.combinat.boltzmann_sampling.grammar import *
-
-            sage: A, B, C = Ref("A"), Ref("B"), Ref("C")
-            sage: Union(A, B, C)
-            Union(A, B, C)
-
-            sage: z = Atom("z")
-            sage: Union(z, "A")
-            Union(z, A)
         """
         if len(args) < 2:
             if len(args) == 0:
@@ -196,7 +208,6 @@ def _latex_(self):
         """Return the LaTeX representation of a union.
 
         EXAMPLES::
-            sage: from sage.combinat.boltzmann_sampling.grammar import *
 
             sage: A, B = Ref("A"), Ref("B")
             sage: latex(Union(A, B))
@@ -224,6 +235,15 @@ class Product(Rule):
 
     Product(A, B, C) corresponds to the following grammar in BNF syntax:
     x ::= A × B × C
+
+    EXAMPLES::
+
+        sage: Product("A", "B", "C")
+        Product(A, B, C)
+
+        sage: z = Atom("z")
+        sage: Product(z, "A")
+        Product(z, A)
     """
 
     def __init__(self, *args):
@@ -232,17 +252,6 @@ def __init__(self, *args):
         INPUT:
 
         - ``args`` -- list of strings or Rules; strings are interpreted as Refs
-
-        EXAMPLES::
-            sage: from sage.combinat.boltzmann_sampling.grammar import *
-
-            sage: A, B, C = Ref("A"), Ref("B"), Ref("C")
-            sage: Product(A, B, C)
-            Product(A, B, C)
-
-            sage: z = Atom("z")
-            sage: Product(z, "A")
-            Product(z, A)
         """
         if len(args) < 2:
             if len(args) == 0:
@@ -260,7 +269,6 @@ def _latex_(self):
         r"""Return the LaTeX representation of a product.
 
         EXAMPLES::
-            sage: from sage.combinat.boltzmann_sampling.grammar import *
 
             sage: A, B = Ref("A"), Ref("B")
             sage: latex(Product(A, B))
@@ -292,7 +300,7 @@ class Grammar(SageObject):
     """Context free grammars."""
 
     def __init__(self, rules=None, labelled=False):
-        """Create a grammar.
+        r"""Create a grammar.
 
         INPUT:
 
@@ -301,25 +309,25 @@ def __init__(self, rules=None, labelled=False):
 
         - ``labelled`` (default: False) -- whether the atoms of the grammar are
           labelled or not. In the labelled case, the generator module will draw
-          objects according to the labelled distribution (i.e. P_x[a] =
-          x^|a|/(|a|!A(x))). At the moment, we do not support mixing labelled
-          and unlabelled atoms in a grammar.
+          objects according to the labelled distribution (i.e.
+          `P_x[a] = \frac{x^{|a|}{|a|!A(x)}`).
 
         EXAMPLES::
-            sage: from sage.combinat.boltzmann_sampling.grammar import *
 
             sage: z = Atom("z")
-            sage: leaf = Atom("leaf", size=0)
-            sage: Grammar(rules={"B": Union(leaf, Product(z, "B", "B"))})
-            B -> Union(leaf, Product(z, B, B))
+            sage: eps = Atom("eps", size=0)
+
+            sage: Grammar(rules={"S": Union(eps, Product(z, "S"))})
+            S -> Union(eps, Product(z, S))
+
+            sage: Grammar(rules={"B": Union(eps, Product(z, "B", "B"))})
+            B -> Union(eps, Product(z, B, B))
 
-            sage: z = Atom("z")
-            sage: nil = Atom("nil", size=0)
             sage: g = Grammar()
             sage: g.set_rule("D", Union(z, Product(z, "S", z)))
-            sage: g.set_rule("S", Union(nil, Product("D", "S")))
+            sage: g.set_rule("S", Union(eps, Product("D", "S")))
             sage: g
-            S -> Union(nil, Product(D, S))
+            S -> Union(eps, Product(D, S))
             D -> Union(z, Product(z, S, z))
         """
         self.rules = {}
@@ -338,8 +346,7 @@ def set_rule(self, name, rule):
 
         - ``rule`` -- a Rule
 
-        EXAMPLES::
-            sage: from sage.combinat.boltzmann_sampling.grammar import *
+        EXAMPLES:
 
             sage: g = Grammar()
             sage: g.set_rule("A", Union("B", "C"))
@@ -395,7 +402,6 @@ def _latex_(self):
         r"""Return a LaTeX representation of the grammar.
 
         EXAMPLES::
-            sage: from sage.combinat.boltzmann_sampling.grammar import *
 
             sage: y, z = Atom("y"), Atom("z")
             sage: leaf = Atom("leaf", size=0)
diff --git a/src/sage/combinat/boltzmann_sampling/oracle.py b/src/sage/combinat/boltzmann_sampling/oracle.py
index f3da2a7a54c..e4e60785423 100644
--- a/src/sage/combinat/boltzmann_sampling/oracle.py
+++ b/src/sage/combinat/boltzmann_sampling/oracle.py
@@ -1,5 +1,5 @@
 """
-Various implementations of oracles for Boltzmann sampling.
+Various oracle implementations for Boltzmann sampling.
 
 Oracles are used to get (often approximate) values of generating functions.
 Thanks to the symbolic method, functionnal equations can be derived from
@@ -7,10 +7,12 @@
 genrating functions based on these equations.
 
 Currently two oracles are implemented:
-- ``SimpleOracle`` implements approximation by simple iteration of the
+
+- :class:`SimpleOracle` implements approximation by simple iteration of the
   equations.
-- ``OracleFromFunctions`` wraps an generating function given in the form of
-  a python ore sage function as an oracle.
+
+- :class:`OracleFromFunctions` wraps an generating function given in the form
+  of a python ore sage function as an oracle.
 
 AUTHORS:
 - Matthieu Dien (2019): initial version
@@ -32,8 +34,6 @@ class SimpleOracle:
     """Simple oracle for critical Boltzmann sampling based on iteration.
 
     EXAMPLES::
-        sage: from sage.combinat.boltzmann_sampling.grammar import *
-        sage: from sage.combinat.boltzmann_sampling.oracle import SimpleOracle
 
         sage: leaf = Atom("leaf", size=0)
         sage: z = Atom("z")
@@ -41,6 +41,7 @@ class SimpleOracle:
         sage: oracle = SimpleOracle(g)
         sage: oracle("z")  # abs tol 0.0001
         0.25
+
         sage: oracle("B")  # abs tol 0.0001
         2
     """
@@ -122,8 +123,6 @@ def _find_singularity(self):
     def _compute_weights(self):
         self.weights = self._find_singularity()
 
-    # XXX. No _latex_ method, is this a problem?
-
     def _repr_(self):
         return "SimpleOracle for {}".format(latex(self.grammar))
 
@@ -162,12 +161,12 @@ def __init__(self, variables, gen_funs):
           the variables from the first argument as named argument.
 
         EXAMPLES::
-            sage: from sage.combinat.boltzmann_sampling.oracle import *
 
             sage: B(z) = (1 - sqrt(1 - 4 * z)) / (2 * z)
             sage: oracle = OracleFromFunctions({"z": 1/4}, {"B": B})
             sage: oracle("z")  # abs tol 0.0000001
             0.25
+
             sage: oracle("B")  # abs tol 0.0000001
             2
         """