diff --git a/src/doc/en/reference/combinat/module_list.rst b/src/doc/en/reference/combinat/module_list.rst
index b469a7e30f4..a4f05a0fdc7 100644
--- a/src/doc/en/reference/combinat/module_list.rst
+++ b/src/doc/en/reference/combinat/module_list.rst
@@ -291,6 +291,7 @@ Comprehensive Module List
     sage/combinat/rsk
     sage/combinat/schubert_polynomial
     sage/combinat/set_partition
+    sage/combinat/set_partition_iterator
     sage/combinat/set_partition_ordered
     sage/combinat/sf/all
     sage/combinat/sf/character
diff --git a/src/sage/combinat/combinat_cython.pxd b/src/sage/combinat/combinat_cython.pxd
index c27be7b0bed..40cae00a781 100644
--- a/src/sage/combinat/combinat_cython.pxd
+++ b/src/sage/combinat/combinat_cython.pxd
@@ -2,7 +2,4 @@ from sage.libs.gmp.all cimport mpz_t
 
 cdef mpz_stirling_s2(mpz_t s, unsigned long n, unsigned long k)
 
-cdef list from_word(list w, list base_set)
-
 cdef list convert(Py_ssize_t* f, Py_ssize_t n)
-
diff --git a/src/sage/combinat/combinat_cython.pyx b/src/sage/combinat/combinat_cython.pyx
index 1f8c56cc16d..c0905491ac2 100644
--- a/src/sage/combinat/combinat_cython.pyx
+++ b/src/sage/combinat/combinat_cython.pyx
@@ -1,3 +1,4 @@
+# cython: binding=True
 """
 Fast computation of combinatorial functions (Cython + mpz)
 
@@ -16,12 +17,16 @@ AUTHORS:
   Lyndon words, and perfect matchings
 """
 
-cimport cython
-
 from cysignals.memory cimport check_allocarray, sig_free
 
 from sage.libs.gmp.all cimport *
 from sage.rings.integer cimport Integer
+from sage.misc.lazy_import import LazyImport
+
+set_partition_iterator = LazyImport('sage.combinat.set_partition_iterator', 'set_partition_iterator', deprecation=35741)
+set_partition_iterator_blocks = LazyImport('sage.combinat.set_partition_iterator', 'set_partition_iterator_blocks', deprecation=35741)
+linear_extension_iterator = LazyImport('sage.combinat.posets.linear_extension_iterator', 'linear_extension_iterator', deprecation=35741)
+
 
 cdef void mpz_addmul_alt(mpz_t s, mpz_t t, mpz_t u, unsigned long parity):
     """
@@ -199,133 +204,6 @@ def lyndon_word_iterator(Py_ssize_t n, Py_ssize_t k):
         while a[i] == n - 1:
             i -= 1
 
-#####################################################################
-## Set partition iterators
-
-@cython.wraparound(False)
-@cython.boundscheck(False)
-cdef list from_word(list w, list base_set):
-    cdef list sp = []
-    cdef Py_ssize_t i
-    cdef Py_ssize_t b
-    for i in range(len(w)):
-        b = <Py_ssize_t> (w[i])
-        x = base_set[i]
-        if len(sp) <= b:
-            sp.append([x])
-        else:
-            sp[b].append(x)
-    return sp
-
-@cython.wraparound(False)
-@cython.boundscheck(False)
-def set_partition_iterator(base_set):
-    """
-    A fast iterator for the set partitions of the base set, which
-    returns lists of lists instead of set partitions types.
-
-    EXAMPLES::
-
-        sage: from sage.combinat.combinat_cython import set_partition_iterator
-        sage: list(set_partition_iterator([1,-1,x]))
-        [[[1, -1, x]],
-         [[1, -1], [x]],
-         [[1, x], [-1]],
-         [[1], [-1, x]],
-         [[1], [-1], [x]]]
-    """
-    cdef list base = list(base_set)
-
-    # Knuth, TAOCP 4A 7.2.1.5, Algorithm H
-    cdef Py_ssize_t N = len(base)
-    # H1: initialize
-    cdef list a = [0] * N
-    if N <= 1:
-        yield from_word(a, base)
-        return
-
-    cdef list b = [1] * N
-    cdef Py_ssize_t j
-    cdef Py_ssize_t last = N - 1
-    while True:
-        # H2: visit
-        yield from_word(a, base)
-        if a[last] == b[last]:
-            # H4: find j
-            j = N - 2
-            while a[j] == b[j]:
-                j -= 1
-            # H5: increase a_j
-            if j == 0:
-                break
-            a[j] += 1
-            # H6: zero out a_{j+1},...,a_{n-1}
-            b[last] = b[j] + int(a[j] == b[j])
-            j += 1
-            while j < N - 1:
-                a[j] = 0
-                b[j] = b[last]
-                j += 1
-            a[last] = 0
-        else:
-            # H3: increase a_{n-1}
-            a[last] += 1
-
-@cython.wraparound(False)
-@cython.boundscheck(False)
-def _set_partition_block_gen(Py_ssize_t n, Py_ssize_t k, list a):
-    r"""
-    Recursively generate set partitions of ``n`` with fixed block
-    size ``k`` using Algorithm 4.23 from [Rus2003]_.
-    ``a`` is a list of size ``n``.
-
-    EXAMPLES::
-
-        sage: from sage.combinat.combinat_cython import _set_partition_block_gen
-        sage: a = list(range(3))
-        sage: for p in _set_partition_block_gen(3, 2, a):
-        ....:     print(p)
-        [0, 1, 0]
-        [0, 1, 1]
-        [0, 0, 1]
-    """
-    cdef Py_ssize_t i
-    if n == k:
-        yield a
-        return
-
-    for i in range(k):
-        a[n-1] = i
-        for P in _set_partition_block_gen(n-1, k, a):
-            yield P
-        a[n-1] = n-1
-    if k > 1:
-        a[n-1] = k-1
-        for P in _set_partition_block_gen(n-1, k-1, a):
-            yield P
-        a[n-1] = n-1
-
-@cython.wraparound(False)
-@cython.boundscheck(False)
-def set_partition_iterator_blocks(base_set, Py_ssize_t k):
-    """
-    A fast iterator for the set partitions of the base set into the
-    specified number of blocks, which returns lists of lists
-    instead of set partitions types.
-
-    EXAMPLES::
-
-        sage: from sage.combinat.combinat_cython import set_partition_iterator_blocks
-        sage: list(set_partition_iterator_blocks([1,-1,x], 2))
-        [[[1, x], [-1]], [[1], [-1, x]], [[1, -1], [x]]]
-    """
-    cdef list base = list(base_set)
-    cdef Py_ssize_t n = len(base)
-    cdef list a = list(range(n))
-    # TODO: implement _set_partition_block_gen as an iterative algorithm
-    for P in _set_partition_block_gen(n, k, a):
-        yield from_word(<list> P, base)
-
 ## Perfect matchings iterator
 
 def perfect_matchings_iterator(Py_ssize_t n):
@@ -410,292 +288,6 @@ cdef list convert(Py_ssize_t* f, Py_ssize_t n):
             ret.append((i, f[i]))
     return ret
 
-#####################################################################
-## Linear extension iterator
-
-from copy import copy
-def _linear_extension_prepare(D):
-    r"""
-    The preprocessing routine in Figure 7 of "Generating Linear
-    Extensions Fast" by Preusse and Ruskey.
-
-    INPUT:
-
-    - ``D``, the Hasse diagram of a poset
-
-    OUTPUT:
-
-    - a triple ``(le, a, b)``, where ``le`` is the first linear
-      extension, and ``a`` and ``b`` are lists such that ``a[i]`` and
-      ``b[i]`` are minimal elements of ``D`` after removing ``a[:i]``
-      and ``b[:i]``.
-
-    TESTS::
-
-        sage: from sage.combinat.combinat_cython import _linear_extension_prepare
-        sage: D = Poset({ 0:[1,2], 1:[3], 2:[3,4] })._hasse_diagram
-        sage: _linear_extension_prepare(D)
-        ([0, 1, 2, 3, 4], [1, 3], [2, 4])
-
-    """
-    dag_copy = copy(D) # this copy is destroyed during preparation
-    le = []
-    a = []
-    b = []
-
-    # the preprocessing routine found in Figure 7 of
-    # "Generating Linear Extensions Fast" by
-    # Pruesse and Ruskey
-    while dag_copy.num_verts() != 0:
-        # find all the minimal elements of dag_copy
-        minimal_elements = dag_copy.sources()
-        if not minimal_elements:
-            raise ValueError("the digraph must be acyclic to have linear extensions")
-        elif len(minimal_elements) == 1:
-            le.append(minimal_elements[0])
-            dag_copy.delete_vertex(minimal_elements[0])
-        else:
-            ap = minimal_elements[0]
-            bp = minimal_elements[1]
-            a.append(ap)
-            b.append(bp)
-            le.append(ap)
-            le.append(bp)
-            dag_copy.delete_vertex(ap)
-            dag_copy.delete_vertex(bp)
-
-    return (le, a, b)
-
-@cython.wraparound(False)
-@cython.boundscheck(False)
-cdef void _linear_extension_switch(list _le, list _a, list _b, list _is_plus, Py_ssize_t i):
-    """
-    This implements the ``Switch`` procedure described on page 7
-    of "Generating Linear Extensions Fast" by Pruesse and Ruskey.
-
-    If ``i == -1``, then the sign is changed.  Otherwise, then
-    ``_a[i]`` and ``_b[i]`` are transposed.
-
-    """
-    cdef Py_ssize_t a_index, b_index
-    if i == -1:
-        _is_plus[0] = not _is_plus[0]
-    else:
-        a = _a[i]
-        b = _b[i]
-        a_index = _le.index(a)
-        b_index = _le.index(b)
-        _le[a_index] = b
-        _le[b_index] = a
-        _b[i] = a
-        _a[i] = b
-
-@cython.wraparound(False)
-@cython.boundscheck(False)
-cdef bint _linear_extension_right_a(_D, list _le, list _a, list _b, Py_ssize_t i):
-    """
-    Return ``True`` if and only if ``_a[i]`` is incomparable with the
-    element to its right in ``_le`` and the element to the right is
-    not ``_b[i]``.
-
-    This is the ``Right`` function described on page 8 of
-    "Generating Linear Extensions Fast" by Pruesse and Ruskey.
-
-    ::
-
-        sage: D = Poset({ 0:[1,2], 1:[3], 2:[3,4] })._hasse_diagram             # not tested
-        sage: _linear_extension_right_a(D, [0, 1, 2, 4, 3], [1, 4], [2, 3], 0)  # not tested
-        False
-        sage: _linear_extension_right_a(D, [0, 1, 2, 4, 3], [1, 4], [2, 3], 1)  # not tested
-        False
-
-    """
-    cdef Py_ssize_t yindex
-    x = _a[i]
-    yindex = _le.index(x) + 1
-    if yindex >= len(_le):
-        return False
-    y = _le[yindex]
-    return y != _b[i] and _D.are_incomparable(x, y)
-
-@cython.wraparound(False)
-@cython.boundscheck(False)
-cdef bint _linear_extension_right_b(_D, list _le, list _a, list _b, Py_ssize_t i):
-    """
-    Return True if and only if ``_b[i]`` is incomparable with the
-    elements to its right in ``_le``.
-
-    This is the ``Right`` function described on page 8 of
-    "Generating Linear Extensions Fast" by Pruesse and Ruskey.
-
-    ::
-
-        sage: D = Poset({ 0:[1,2], 1:[3], 2:[3,4] })._hasse_diagram             # not tested
-        sage: _linear_extension_right_b(D, [0, 1, 2, 4, 3], [1, 4], [2, 3], 0)  # not tested
-        False
-        sage: _linear_extension_right_b(D, [0, 1, 2, 4, 3], [1, 4], [2, 3], 1)  # not tested
-        False
-
-    """
-    cdef Py_ssize_t yindex
-    x = _b[i]
-    yindex = _le.index(x) + 1
-    if yindex >= len(_le):
-        return False
-    y = _le[yindex]
-    return _D.are_incomparable(x, y)
-
-@cython.wraparound(False)
-@cython.boundscheck(False)
-def _linear_extension_gen(_D, list _le, list _a, list _b, list _is_plus, Py_ssize_t i):
-    """
-    This a Python version of the GenLE routine found in Figure 8
-    of "Generating Linear Extensions Fast" by Pruesse and Ruskey.
-
-    TESTS::
-
-        sage: from sage.combinat.combinat_cython import _linear_extension_prepare, _linear_extension_gen
-        sage: D = Poset({ 0:[1,2], 1:[3], 2:[3,4] })._hasse_diagram
-        sage: le, a, b = _linear_extension_prepare(D)
-        sage: [e for e in _linear_extension_gen(D, le, a, b, [True], len(a)-1)]
-        [[0, 2, 1, 3, 4]]
-
-    """
-    cdef int mra, mrb, mla
-    cdef Py_ssize_t index, index1
-    cdef bint typical
-    if i == -1:
-        return
-
-    for e in _linear_extension_gen(_D, _le, _a, _b, _is_plus, i-1):
-        yield e
-    mrb = 0
-    typical = False
-    while _linear_extension_right_b(_D, _le, _a, _b, i):
-        mrb += 1
-        # move_right
-        index = _le.index(_b[i])
-        index1 = index + 1
-        _le[index] = _le[index1]
-        _le[index1] = _b[i]
-        if _is_plus[0]:
-            yield _le[:]
-
-        for e in _linear_extension_gen(_D, _le, _a, _b, _is_plus, i-1):
-            yield e
-        mra = 0
-        while _linear_extension_right_a(_D, _le, _a, _b, i):
-            typical = True
-            mra += 1
-            # move_right
-            index = _le.index(_a[i])
-            index1 = index+1
-            _le[index] = _le[index1]
-            _le[index1] = _a[i]
-            if _is_plus[0]:
-                yield _le[:]
-
-            for e in _linear_extension_gen(_D, _le, _a, _b, _is_plus, i-1):
-                yield e
-
-        if typical:
-            _linear_extension_switch(_le, _a, _b, _is_plus, i-1)
-            if _is_plus[0]:
-                yield _le[:]
-
-            for e in _linear_extension_gen(_D, _le, _a, _b, _is_plus, i-1):
-                yield e
-            if mrb % 2 == 1:
-                mla = mra - 1
-            else:
-                mla = mra + 1
-            for _ in range(mla):
-                # move_left
-                index = _le.index(_a[i])
-                index1 = index-1
-                _le[index] = _le[index1]
-                _le[index1] = _a[i]
-                if _is_plus[0]:
-                    yield _le[:]
-
-                for e in _linear_extension_gen(_D, _le, _a, _b, _is_plus, i-1):
-                    yield e
-
-    if typical and (mrb % 2 == 1):
-        # move_left
-        index = _le.index(_a[i])
-        index1 = index-1
-        _le[index] = _le[index1]
-        _le[index1] = _a[i]
-        if _is_plus[0]:
-            yield _le[:]
-    else:
-        _linear_extension_switch(_le, _a, _b, _is_plus, i-1)
-        if _is_plus[0]:
-            yield _le[:]
-    for e in _linear_extension_gen(_D, _le, _a, _b, _is_plus, i-1):
-        yield e
-    for _ in range(mrb):
-        # move_left
-        index = _le.index(_b[i])
-        index1 = index-1
-        _le[index] = _le[index1]
-        _le[index1] = _b[i]
-        if _is_plus[0]:
-            yield _le[:]
-
-        for e in _linear_extension_gen(_D, _le, _a, _b, _is_plus, i-1):
-            yield e
-
-
-def linear_extension_iterator(D):
-    """
-    Iterate over the linear extensions of the poset.
-
-    The list ``_le`` keeps track of the current linear extensions.  The
-    boolean variable ``is_plus`` keeps track of the "sign".
-
-    INPUT:
-
-    - ``D``, the Hasse diagram of a poset.
-
-    .. WARNING::
-
-        It is assumed that ``D`` is not modified while the linear
-        extensions are generated.
-
-    EXAMPLES::
-
-        sage: from sage.combinat.combinat_cython import linear_extension_iterator
-        sage: D = Poset({ 0:[1,2], 1:[3], 2:[3,4] })._hasse_diagram
-        sage: list(linear_extension_iterator(D))
-        [[0, 1, 2, 3, 4],
-         [0, 2, 1, 3, 4],
-         [0, 2, 1, 4, 3],
-         [0, 2, 4, 1, 3],
-         [0, 1, 2, 4, 3]]
-
-        sage: D = posets.BooleanLattice(3)._hasse_diagram
-        sage: len(list(linear_extension_iterator(D)))
-        48
-
-        sage: D = posets.AntichainPoset(9)._hasse_diagram
-        sage: len(list(linear_extension_iterator(D))) == factorial(9)           # long time
-        True
-    """
-    _le, _a, _b = _linear_extension_prepare(D)
-    _max_pair = len(_a) - 1
-    _is_plus = [True] # this is modified by _linear_extension_switch
-
-    yield _le[:]
-    for e in _linear_extension_gen(D, _le, _a, _b, _is_plus, _max_pair):
-        yield e
-    _linear_extension_switch(_le, _a, _b, _is_plus, _max_pair)
-    if _is_plus[0]:
-        yield _le[:]
-    for e in _linear_extension_gen(D, _le, _a, _b, _is_plus, _max_pair):
-        yield e
-
 #####################################################################
 ## Set partition composition
 
@@ -711,7 +303,7 @@ def set_partition_composition(tuple sp1, tuple sp2):
         sage: sp1 = ((1,-2),(2,-1))
         sage: sp2 = ((1,-2),(2,-1))
         sage: p, c = set_partition_composition(sp1, sp2)
-        sage: (SetPartition(p), c) == (SetPartition([[1,-1],[2,-2]]), 0)
+        sage: (SetPartition(p), c) == (SetPartition([[1,-1],[2,-2]]), 0)                # optional - sage.combinat
         True
     """
     cdef int num_loops = 0  # The number of loops removed
diff --git a/src/sage/combinat/diagram_algebras.py b/src/sage/combinat/diagram_algebras.py
index bd281fdbb77..8bef929e9c6 100644
--- a/src/sage/combinat/diagram_algebras.py
+++ b/src/sage/combinat/diagram_algebras.py
@@ -31,9 +31,10 @@
 from sage.structure.unique_representation import UniqueRepresentation
 from sage.combinat.combinat import bell_number, catalan_number
 from sage.structure.global_options import GlobalOptions
-from sage.combinat.combinat_cython import (set_partition_iterator, perfect_matchings_iterator,
+from sage.combinat.combinat_cython import (perfect_matchings_iterator,
                                            set_partition_composition)
 from sage.combinat.set_partition import SetPartitions, AbstractSetPartition
+from sage.combinat.set_partition_iterator import set_partition_iterator
 from sage.combinat.symmetric_group_algebra import SymmetricGroupAlgebra_n
 from sage.combinat.permutation import Permutations
 from sage.graphs.graph import Graph
diff --git a/src/sage/combinat/posets/hasse_cython.pyx b/src/sage/combinat/posets/hasse_cython.pyx
index 0e4d13ff033..8f7a869fc17 100644
--- a/src/sage/combinat/posets/hasse_cython.pyx
+++ b/src/sage/combinat/posets/hasse_cython.pyx
@@ -1,4 +1,4 @@
-# -*- coding: utf-8 -*-
+# cython: binding=True
 r"""
 Some fast computations for finite posets
 """
@@ -11,146 +11,14 @@ Some fast computations for finite posets
 # (at your option) any later version.
 #                  https://www.gnu.org/licenses/
 # ****************************************************************************
-from cysignals.signals cimport sig_check
-from cysignals.memory cimport sig_malloc, sig_free
-
-from sage.libs.flint.fmpz cimport *
-from sage.libs.flint.fmpz_mat cimport *
-
-from sage.matrix.matrix_integer_dense cimport Matrix_integer_dense
-
-from sage.rings.integer_ring import ZZ
-from sage.matrix.matrix_space import MatrixSpace
 from sage.categories.finite_enumerated_sets import FiniteEnumeratedSets
+from sage.misc.lazy_import import LazyImport
 from sage.sets.recursively_enumerated_set import RecursivelyEnumeratedSet_forest
 
-
-cpdef Matrix_integer_dense moebius_matrix_fast(list positions):
-    """
-    Compute the Möbius matrix of a poset by a specific triangular inversion.
-
-    INPUT:
-
-    a list of sets of integers describing the poset, as given by the
-    lazy attribute ``_leq_storage`` of Hasse diagrams.
-
-    OUTPUT:
-
-    a dense matrix
-
-    EXAMPLES::
-
-        sage: from sage.combinat.posets.hasse_cython import moebius_matrix_fast
-        sage: D = [{0,1},{1}]
-        sage: moebius_matrix_fast(D)
-        [ 1 -1]
-        [ 0  1]
-        sage: P = posets.TamariLattice(5)
-        sage: H = P._hasse_diagram
-        sage: D = H._leq_storage
-        sage: moebius_matrix_fast(D)
-        42 x 42 dense matrix over Integer Ring (...)
-    """
-    cdef Matrix_integer_dense A
-    cdef Py_ssize_t n = len(positions)
-    cdef Py_ssize_t i
-    cdef int j, k
-    A = Matrix_integer_dense.__new__(Matrix_integer_dense,
-                                     MatrixSpace(ZZ, n, n), None, None, None)
-    fmpz_mat_one(A._matrix)
-    cdef Py_ssize_t *pos_lens = <Py_ssize_t*> sig_malloc(n*sizeof(Py_ssize_t))
-    cdef int **pos_array = <int**> sig_malloc(n*sizeof(int*))
-    cdef Py_ssize_t jind, kind
-    for i in range(n):
-        pos_lens[i] = len(positions[i])
-        pos_array[i] = <int*> sig_malloc(pos_lens[i]*sizeof(int))
-        for jind,j in enumerate(positions[i]):
-            pos_array[i][jind] = j
-
-    for i in range(n - 1, -1, -1):
-        sig_check()
-        for jind in range(pos_lens[i]):
-            j = pos_array[i][jind]
-            if j != i:
-                for kind in range(pos_lens[j]):
-                    k = pos_array[j][kind]
-                    fmpz_sub(fmpz_mat_entry(A._matrix, i, k),
-                             fmpz_mat_entry(A._matrix, i, k),
-                             fmpz_mat_entry(A._matrix, j, k))
-    for i in range(n):
-        sig_free(pos_array[i])
-    sig_free(pos_array)
-    sig_free(pos_lens)
-    return A
-
-
-cpdef Matrix_integer_dense coxeter_matrix_fast(list positions):
-    """
-    Compute the Coxeter matrix of a poset by a specific algorithm.
-
-    INPUT:
-
-    a list of sets of integers describing the poset, as given by the
-    lazy attribute ``_leq_storage`` of Hasse diagrams.
-
-    OUTPUT:
-
-    a dense matrix
-
-    EXAMPLES::
-
-        sage: from sage.combinat.posets.hasse_cython import coxeter_matrix_fast
-        sage: D = [{0,1},{1}]
-        sage: coxeter_matrix_fast(D)
-        [ 0 -1]
-        [ 1 -1]
-        sage: P = posets.TamariLattice(5)
-        sage: H = P._hasse_diagram
-        sage: D = H._leq_storage
-        sage: coxeter_matrix_fast(D)
-        42 x 42 dense matrix over Integer Ring (...)
-    """
-    cdef Matrix_integer_dense A
-    cdef Py_ssize_t n = len(positions)
-    cdef Py_ssize_t i
-    cdef int j, k
-    A = Matrix_integer_dense.__new__(Matrix_integer_dense,
-                                     MatrixSpace(ZZ, n, n), None, None, None)
-    fmpz_mat_one(A._matrix)
-    cdef Py_ssize_t *pos_lens = <Py_ssize_t*> sig_malloc(n*sizeof(Py_ssize_t))
-    cdef int **pos_array = <int**> sig_malloc(n*sizeof(int*))
-    cdef Py_ssize_t jind, kind
-    for i in range(n):
-        pos_lens[i] = len(positions[i])
-        pos_array[i] = <int*> sig_malloc(pos_lens[i]*sizeof(int))
-        for jind,j in enumerate(positions[i]):
-            pos_array[i][jind] = j
-
-    for i in range(n - 1, -1, -1):
-        sig_check()
-        for jind in range(pos_lens[i]):
-            j = pos_array[i][jind]
-            if j != i:
-                for kind in range(pos_lens[j]):
-                    k = pos_array[j][kind]
-                    fmpz_sub(fmpz_mat_entry(A._matrix, k, i),
-                             fmpz_mat_entry(A._matrix, k, i),
-                             fmpz_mat_entry(A._matrix, k, j))
-    for i in range(n):
-        sig_check()
-        for jind in range(pos_lens[i]):
-            j = pos_array[i][jind]
-            if j != i:
-                for k in range(n):
-                    fmpz_add(fmpz_mat_entry(A._matrix, i, k),
-                             fmpz_mat_entry(A._matrix, i, k),
-                             fmpz_mat_entry(A._matrix, j, k))
-    for i in range(n):
-        sig_free(pos_array[i])
-    sig_free(pos_array)
-    sig_free(pos_lens)
-    fmpz_mat_scalar_mul_si(A._matrix, A._matrix, -1)
-    return A
+coxeter_matrix_fast = LazyImport('sage.combinat.posets.hasse_cython_flint', 'coxeter_matrix_fast',
+                                 deprecation=35741)
+moebius_matrix_fast = LazyImport('sage.combinat.posets.hasse_cython_flint', 'moebius_matrix_fast',
+                                 deprecation=35741)
 
 
 class IncreasingChains(RecursivelyEnumeratedSet_forest):
diff --git a/src/sage/combinat/posets/hasse_cython_flint.pyx b/src/sage/combinat/posets/hasse_cython_flint.pyx
new file mode 100644
index 00000000000..0b33eed0b3f
--- /dev/null
+++ b/src/sage/combinat/posets/hasse_cython_flint.pyx
@@ -0,0 +1,149 @@
+# cython: binding=True
+r"""
+Some fast computations for finite posets using FLINT matrices
+"""
+# ****************************************************************************
+#       Copyright (C) 2020 Frédéric Chapoton <chapoton@unistra.fr>
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 2 of the License, or
+# (at your option) any later version.
+#                  https://www.gnu.org/licenses/
+# ****************************************************************************
+from cysignals.signals cimport sig_check
+from cysignals.memory cimport sig_malloc, sig_free
+
+from sage.libs.flint.fmpz cimport *
+from sage.libs.flint.fmpz_mat cimport *
+from sage.matrix.matrix_integer_dense cimport Matrix_integer_dense
+from sage.matrix.matrix_space import MatrixSpace
+from sage.rings.integer_ring import ZZ
+
+
+cpdef Matrix_integer_dense moebius_matrix_fast(list positions):
+    """
+    Compute the Möbius matrix of a poset by a specific triangular inversion.
+
+    INPUT:
+
+    a list of sets of integers describing the poset, as given by the
+    lazy attribute ``_leq_storage`` of Hasse diagrams.
+
+    OUTPUT:
+
+    a dense matrix
+
+    EXAMPLES::
+
+        sage: from sage.combinat.posets.hasse_cython_flint import moebius_matrix_fast
+        sage: D = [{0,1},{1}]
+        sage: moebius_matrix_fast(D)
+        [ 1 -1]
+        [ 0  1]
+        sage: P = posets.TamariLattice(5)
+        sage: H = P._hasse_diagram
+        sage: D = H._leq_storage
+        sage: moebius_matrix_fast(D)
+        42 x 42 dense matrix over Integer Ring (...)
+    """
+    cdef Matrix_integer_dense A
+    cdef Py_ssize_t n = len(positions)
+    cdef Py_ssize_t i
+    cdef int j, k
+    A = Matrix_integer_dense.__new__(Matrix_integer_dense,
+                                     MatrixSpace(ZZ, n, n), None, None, None)
+    fmpz_mat_one(A._matrix)
+    cdef Py_ssize_t *pos_lens = <Py_ssize_t*> sig_malloc(n*sizeof(Py_ssize_t))
+    cdef int **pos_array = <int**> sig_malloc(n*sizeof(int*))
+    cdef Py_ssize_t jind, kind
+    for i in range(n):
+        pos_lens[i] = len(positions[i])
+        pos_array[i] = <int*> sig_malloc(pos_lens[i]*sizeof(int))
+        for jind,j in enumerate(positions[i]):
+            pos_array[i][jind] = j
+
+    for i in range(n - 1, -1, -1):
+        sig_check()
+        for jind in range(pos_lens[i]):
+            j = pos_array[i][jind]
+            if j != i:
+                for kind in range(pos_lens[j]):
+                    k = pos_array[j][kind]
+                    fmpz_sub(fmpz_mat_entry(A._matrix, i, k),
+                             fmpz_mat_entry(A._matrix, i, k),
+                             fmpz_mat_entry(A._matrix, j, k))
+    for i in range(n):
+        sig_free(pos_array[i])
+    sig_free(pos_array)
+    sig_free(pos_lens)
+    return A
+
+
+cpdef Matrix_integer_dense coxeter_matrix_fast(list positions):
+    """
+    Compute the Coxeter matrix of a poset by a specific algorithm.
+
+    INPUT:
+
+    a list of sets of integers describing the poset, as given by the
+    lazy attribute ``_leq_storage`` of Hasse diagrams.
+
+    OUTPUT:
+
+    a dense matrix
+
+    EXAMPLES::
+
+        sage: from sage.combinat.posets.hasse_cython_flint import coxeter_matrix_fast
+        sage: D = [{0,1},{1}]
+        sage: coxeter_matrix_fast(D)
+        [ 0 -1]
+        [ 1 -1]
+        sage: P = posets.TamariLattice(5)
+        sage: H = P._hasse_diagram
+        sage: D = H._leq_storage
+        sage: coxeter_matrix_fast(D)
+        42 x 42 dense matrix over Integer Ring (...)
+    """
+    cdef Matrix_integer_dense A
+    cdef Py_ssize_t n = len(positions)
+    cdef Py_ssize_t i
+    cdef int j, k
+    A = Matrix_integer_dense.__new__(Matrix_integer_dense,
+                                     MatrixSpace(ZZ, n, n), None, None, None)
+    fmpz_mat_one(A._matrix)
+    cdef Py_ssize_t *pos_lens = <Py_ssize_t*> sig_malloc(n*sizeof(Py_ssize_t))
+    cdef int **pos_array = <int**> sig_malloc(n*sizeof(int*))
+    cdef Py_ssize_t jind, kind
+    for i in range(n):
+        pos_lens[i] = len(positions[i])
+        pos_array[i] = <int*> sig_malloc(pos_lens[i]*sizeof(int))
+        for jind,j in enumerate(positions[i]):
+            pos_array[i][jind] = j
+
+    for i in range(n - 1, -1, -1):
+        sig_check()
+        for jind in range(pos_lens[i]):
+            j = pos_array[i][jind]
+            if j != i:
+                for kind in range(pos_lens[j]):
+                    k = pos_array[j][kind]
+                    fmpz_sub(fmpz_mat_entry(A._matrix, k, i),
+                             fmpz_mat_entry(A._matrix, k, i),
+                             fmpz_mat_entry(A._matrix, k, j))
+    for i in range(n):
+        sig_check()
+        for jind in range(pos_lens[i]):
+            j = pos_array[i][jind]
+            if j != i:
+                for k in range(n):
+                    fmpz_add(fmpz_mat_entry(A._matrix, i, k),
+                             fmpz_mat_entry(A._matrix, i, k),
+                             fmpz_mat_entry(A._matrix, j, k))
+    for i in range(n):
+        sig_free(pos_array[i])
+    sig_free(pos_array)
+    sig_free(pos_lens)
+    fmpz_mat_scalar_mul_si(A._matrix, A._matrix, -1)
+    return A
diff --git a/src/sage/combinat/posets/hasse_diagram.py b/src/sage/combinat/posets/hasse_diagram.py
index 7e1f1c6e4c6..30477e09c10 100644
--- a/src/sage/combinat/posets/hasse_diagram.py
+++ b/src/sage/combinat/posets/hasse_diagram.py
@@ -1,4 +1,3 @@
-# -*- coding: utf-8 -*-
 r"""
 Hasse diagrams of posets
 
@@ -17,19 +16,21 @@
 # ****************************************************************************
 from __future__ import annotations
 
+from collections import deque
+
+from sage.arith.misc import binomial
+from sage.combinat.posets.hasse_cython import IncreasingChains
 from sage.graphs.digraph import DiGraph
-from sage.matrix.constructor import matrix
-from sage.rings.integer_ring import ZZ
-from sage.rings.finite_rings.finite_field_constructor import GF
-from sage.misc.lazy_attribute import lazy_attribute
 from sage.misc.cachefunc import cached_method
-from sage.arith.misc import binomial
+from sage.misc.lazy_attribute import lazy_attribute
+from sage.misc.lazy_import import lazy_import
 from sage.misc.rest_index_of_methods import gen_rest_table_index
-from sage.combinat.posets.hasse_cython import (moebius_matrix_fast,
-                                               coxeter_matrix_fast,
-                                               IncreasingChains)
+from sage.rings.integer_ring import ZZ
 
-from collections import deque
+lazy_import('sage.combinat.posets.hasse_cython_flint',
+            ['moebius_matrix_fast', 'coxeter_matrix_fast'])
+lazy_import('sage.matrix.constructor', 'matrix')
+lazy_import('sage.rings.finite_rings.finite_field_constructor', 'GF')
 
 
 class LatticeError(ValueError):
@@ -126,7 +127,7 @@ def linear_extensions(self):
             sage: list(H.linear_extensions())
             [[0, 1, 2, 3], [0, 2, 1, 3]]
         """
-        from sage.combinat.combinat_cython import linear_extension_iterator
+        from sage.combinat.posets.linear_extension_iterator import linear_extension_iterator
         return linear_extension_iterator(self)
 
     def greedy_linear_extensions_iterator(self):
diff --git a/src/sage/combinat/posets/linear_extension_iterator.pyx b/src/sage/combinat/posets/linear_extension_iterator.pyx
new file mode 100644
index 00000000000..5eb101b32e2
--- /dev/null
+++ b/src/sage/combinat/posets/linear_extension_iterator.pyx
@@ -0,0 +1,288 @@
+# cython: binding=True
+r"""
+Fast linear extension iterator
+"""
+cimport cython
+
+from copy import copy
+def _linear_extension_prepare(D):
+    r"""
+    The preprocessing routine in Figure 7 of "Generating Linear
+    Extensions Fast" by Preusse and Ruskey.
+
+    INPUT:
+
+    - ``D``, the Hasse diagram of a poset
+
+    OUTPUT:
+
+    - a triple ``(le, a, b)``, where ``le`` is the first linear
+      extension, and ``a`` and ``b`` are lists such that ``a[i]`` and
+      ``b[i]`` are minimal elements of ``D`` after removing ``a[:i]``
+      and ``b[:i]``.
+
+    TESTS::
+
+        sage: from sage.combinat.posets.linear_extension_iterator import _linear_extension_prepare
+        sage: D = Poset({ 0:[1,2], 1:[3], 2:[3,4] })._hasse_diagram
+        sage: _linear_extension_prepare(D)
+        ([0, 1, 2, 3, 4], [1, 3], [2, 4])
+
+    """
+    dag_copy = copy(D) # this copy is destroyed during preparation
+    le = []
+    a = []
+    b = []
+
+    # the preprocessing routine found in Figure 7 of
+    # "Generating Linear Extensions Fast" by
+    # Pruesse and Ruskey
+    while dag_copy.num_verts() != 0:
+        # find all the minimal elements of dag_copy
+        minimal_elements = dag_copy.sources()
+        if not minimal_elements:
+            raise ValueError("the digraph must be acyclic to have linear extensions")
+        elif len(minimal_elements) == 1:
+            le.append(minimal_elements[0])
+            dag_copy.delete_vertex(minimal_elements[0])
+        else:
+            ap = minimal_elements[0]
+            bp = minimal_elements[1]
+            a.append(ap)
+            b.append(bp)
+            le.append(ap)
+            le.append(bp)
+            dag_copy.delete_vertex(ap)
+            dag_copy.delete_vertex(bp)
+
+    return (le, a, b)
+
+@cython.wraparound(False)
+@cython.boundscheck(False)
+cdef void _linear_extension_switch(list _le, list _a, list _b, list _is_plus, Py_ssize_t i):
+    """
+    This implements the ``Switch`` procedure described on page 7
+    of "Generating Linear Extensions Fast" by Pruesse and Ruskey.
+
+    If ``i == -1``, then the sign is changed.  Otherwise, then
+    ``_a[i]`` and ``_b[i]`` are transposed.
+
+    """
+    cdef Py_ssize_t a_index, b_index
+    if i == -1:
+        _is_plus[0] = not _is_plus[0]
+    else:
+        a = _a[i]
+        b = _b[i]
+        a_index = _le.index(a)
+        b_index = _le.index(b)
+        _le[a_index] = b
+        _le[b_index] = a
+        _b[i] = a
+        _a[i] = b
+
+@cython.wraparound(False)
+@cython.boundscheck(False)
+cdef bint _linear_extension_right_a(_D, list _le, list _a, list _b, Py_ssize_t i):
+    """
+    Return ``True`` if and only if ``_a[i]`` is incomparable with the
+    element to its right in ``_le`` and the element to the right is
+    not ``_b[i]``.
+
+    This is the ``Right`` function described on page 8 of
+    "Generating Linear Extensions Fast" by Pruesse and Ruskey.
+
+    ::
+
+        sage: D = Poset({ 0:[1,2], 1:[3], 2:[3,4] })._hasse_diagram             # not tested
+        sage: _linear_extension_right_a(D, [0, 1, 2, 4, 3], [1, 4], [2, 3], 0)  # not tested
+        False
+        sage: _linear_extension_right_a(D, [0, 1, 2, 4, 3], [1, 4], [2, 3], 1)  # not tested
+        False
+
+    """
+    cdef Py_ssize_t yindex
+    x = _a[i]
+    yindex = _le.index(x) + 1
+    if yindex >= len(_le):
+        return False
+    y = _le[yindex]
+    return y != _b[i] and _D.are_incomparable(x, y)
+
+@cython.wraparound(False)
+@cython.boundscheck(False)
+cdef bint _linear_extension_right_b(_D, list _le, list _a, list _b, Py_ssize_t i):
+    """
+    Return True if and only if ``_b[i]`` is incomparable with the
+    elements to its right in ``_le``.
+
+    This is the ``Right`` function described on page 8 of
+    "Generating Linear Extensions Fast" by Pruesse and Ruskey.
+
+    ::
+
+        sage: D = Poset({ 0:[1,2], 1:[3], 2:[3,4] })._hasse_diagram             # not tested
+        sage: _linear_extension_right_b(D, [0, 1, 2, 4, 3], [1, 4], [2, 3], 0)  # not tested
+        False
+        sage: _linear_extension_right_b(D, [0, 1, 2, 4, 3], [1, 4], [2, 3], 1)  # not tested
+        False
+
+    """
+    cdef Py_ssize_t yindex
+    x = _b[i]
+    yindex = _le.index(x) + 1
+    if yindex >= len(_le):
+        return False
+    y = _le[yindex]
+    return _D.are_incomparable(x, y)
+
+@cython.wraparound(False)
+@cython.boundscheck(False)
+def _linear_extension_gen(_D, list _le, list _a, list _b, list _is_plus, Py_ssize_t i):
+    """
+    This a Python version of the GenLE routine found in Figure 8
+    of "Generating Linear Extensions Fast" by Pruesse and Ruskey.
+
+    TESTS::
+
+        sage: from sage.combinat.posets.linear_extension_iterator import _linear_extension_prepare, _linear_extension_gen
+        sage: D = Poset({ 0:[1,2], 1:[3], 2:[3,4] })._hasse_diagram
+        sage: le, a, b = _linear_extension_prepare(D)
+        sage: [e for e in _linear_extension_gen(D, le, a, b, [True], len(a)-1)]
+        [[0, 2, 1, 3, 4]]
+
+    """
+    cdef int mra, mrb, mla
+    cdef Py_ssize_t index, index1
+    cdef bint typical
+    if i == -1:
+        return
+
+    for e in _linear_extension_gen(_D, _le, _a, _b, _is_plus, i-1):
+        yield e
+    mrb = 0
+    typical = False
+    while _linear_extension_right_b(_D, _le, _a, _b, i):
+        mrb += 1
+        # move_right
+        index = _le.index(_b[i])
+        index1 = index + 1
+        _le[index] = _le[index1]
+        _le[index1] = _b[i]
+        if _is_plus[0]:
+            yield _le[:]
+
+        for e in _linear_extension_gen(_D, _le, _a, _b, _is_plus, i-1):
+            yield e
+        mra = 0
+        while _linear_extension_right_a(_D, _le, _a, _b, i):
+            typical = True
+            mra += 1
+            # move_right
+            index = _le.index(_a[i])
+            index1 = index+1
+            _le[index] = _le[index1]
+            _le[index1] = _a[i]
+            if _is_plus[0]:
+                yield _le[:]
+
+            for e in _linear_extension_gen(_D, _le, _a, _b, _is_plus, i-1):
+                yield e
+
+        if typical:
+            _linear_extension_switch(_le, _a, _b, _is_plus, i-1)
+            if _is_plus[0]:
+                yield _le[:]
+
+            for e in _linear_extension_gen(_D, _le, _a, _b, _is_plus, i-1):
+                yield e
+            if mrb % 2 == 1:
+                mla = mra - 1
+            else:
+                mla = mra + 1
+            for _ in range(mla):
+                # move_left
+                index = _le.index(_a[i])
+                index1 = index-1
+                _le[index] = _le[index1]
+                _le[index1] = _a[i]
+                if _is_plus[0]:
+                    yield _le[:]
+
+                for e in _linear_extension_gen(_D, _le, _a, _b, _is_plus, i-1):
+                    yield e
+
+    if typical and (mrb % 2 == 1):
+        # move_left
+        index = _le.index(_a[i])
+        index1 = index-1
+        _le[index] = _le[index1]
+        _le[index1] = _a[i]
+        if _is_plus[0]:
+            yield _le[:]
+    else:
+        _linear_extension_switch(_le, _a, _b, _is_plus, i-1)
+        if _is_plus[0]:
+            yield _le[:]
+    for e in _linear_extension_gen(_D, _le, _a, _b, _is_plus, i-1):
+        yield e
+    for _ in range(mrb):
+        # move_left
+        index = _le.index(_b[i])
+        index1 = index-1
+        _le[index] = _le[index1]
+        _le[index1] = _b[i]
+        if _is_plus[0]:
+            yield _le[:]
+
+        for e in _linear_extension_gen(_D, _le, _a, _b, _is_plus, i-1):
+            yield e
+
+
+def linear_extension_iterator(D):
+    """
+    Iterate over the linear extensions of the poset.
+
+    The list ``_le`` keeps track of the current linear extensions.  The
+    boolean variable ``is_plus`` keeps track of the "sign".
+
+    INPUT:
+
+    - ``D``, the Hasse diagram of a poset.
+
+    .. WARNING::
+
+        It is assumed that ``D`` is not modified while the linear
+        extensions are generated.
+
+    EXAMPLES::
+
+        sage: from sage.combinat.posets.linear_extension_iterator import linear_extension_iterator
+        sage: D = Poset({ 0:[1,2], 1:[3], 2:[3,4] })._hasse_diagram
+        sage: list(linear_extension_iterator(D))
+        [[0, 1, 2, 3, 4],
+         [0, 2, 1, 3, 4],
+         [0, 2, 1, 4, 3],
+         [0, 2, 4, 1, 3],
+         [0, 1, 2, 4, 3]]
+
+        sage: D = posets.BooleanLattice(3)._hasse_diagram
+        sage: len(list(linear_extension_iterator(D)))
+        48
+
+        sage: D = posets.AntichainPoset(9)._hasse_diagram
+        sage: len(list(linear_extension_iterator(D))) == factorial(9)           # long time
+        True
+    """
+    _le, _a, _b = _linear_extension_prepare(D)
+    _max_pair = len(_a) - 1
+    _is_plus = [True] # this is modified by _linear_extension_switch
+
+    yield _le[:]
+    for e in _linear_extension_gen(D, _le, _a, _b, _is_plus, _max_pair):
+        yield e
+    _linear_extension_switch(_le, _a, _b, _is_plus, _max_pair)
+    if _is_plus[0]:
+        yield _le[:]
+    for e in _linear_extension_gen(D, _le, _a, _b, _is_plus, _max_pair):
+        yield e
diff --git a/src/sage/combinat/posets/linear_extensions.py b/src/sage/combinat/posets/linear_extensions.py
index e14be386228..84471d31e45 100644
--- a/src/sage/combinat/posets/linear_extensions.py
+++ b/src/sage/combinat/posets/linear_extensions.py
@@ -667,7 +667,7 @@ def __iter__(self):
             sage: list(L)
             [[1, 2, 3, 4], [2, 1, 3, 4], [2, 1, 4, 3], [1, 4, 2, 3], [1, 2, 4, 3]]
         """
-        from sage.combinat.combinat_cython import linear_extension_iterator
+        from sage.combinat.posets.linear_extension_iterator import linear_extension_iterator
         vertex_to_element = self._poset._vertex_to_element
         for lin_ext in linear_extension_iterator(self._poset._hasse_diagram):
             yield self._element_constructor_([vertex_to_element(_) for _ in lin_ext])
diff --git a/src/sage/combinat/set_partition.py b/src/sage/combinat/set_partition.py
index 6edf7d963da..f39895daa35 100644
--- a/src/sage/combinat/set_partition.py
+++ b/src/sage/combinat/set_partition.py
@@ -43,8 +43,8 @@
 from sage.rings.infinity import infinity
 from sage.rings.integer import Integer
 from sage.combinat.combinatorial_map import combinatorial_map
-from sage.combinat.combinat_cython import (set_partition_iterator,
-                                           set_partition_iterator_blocks)
+from sage.combinat.set_partition_iterator import (set_partition_iterator,
+                                                  set_partition_iterator_blocks)
 from sage.combinat.partition import Partition, Partitions
 from sage.combinat.combinat import bell_number, stirling_number2 as stirling2
 from sage.combinat.permutation import Permutation
diff --git a/src/sage/combinat/set_partition_iterator.pyx b/src/sage/combinat/set_partition_iterator.pyx
new file mode 100644
index 00000000000..fff6a71fefe
--- /dev/null
+++ b/src/sage/combinat/set_partition_iterator.pyx
@@ -0,0 +1,130 @@
+# cython: binding=True
+r"""
+Fast set partition iterators
+"""
+cimport cython
+
+
+@cython.wraparound(False)
+@cython.boundscheck(False)
+cdef list from_word(list w, list base_set):
+    cdef list sp = []
+    cdef Py_ssize_t i
+    cdef Py_ssize_t b
+    for i in range(len(w)):
+        b = <Py_ssize_t> (w[i])
+        x = base_set[i]
+        if len(sp) <= b:
+            sp.append([x])
+        else:
+            sp[b].append(x)
+    return sp
+
+@cython.wraparound(False)
+@cython.boundscheck(False)
+def set_partition_iterator(base_set):
+    """
+    A fast iterator for the set partitions of the base set, which
+    returns lists of lists instead of set partitions types.
+
+    EXAMPLES::
+
+        sage: from sage.combinat.set_partition_iterator import set_partition_iterator
+        sage: list(set_partition_iterator([1,-1,x]))                                    # optional - sage.symbolic
+        [[[1, -1, x]],
+         [[1, -1], [x]],
+         [[1, x], [-1]],
+         [[1], [-1, x]],
+         [[1], [-1], [x]]]
+    """
+    cdef list base = list(base_set)
+
+    # Knuth, TAOCP 4A 7.2.1.5, Algorithm H
+    cdef Py_ssize_t N = len(base)
+    # H1: initialize
+    cdef list a = [0] * N
+    if N <= 1:
+        yield from_word(a, base)
+        return
+
+    cdef list b = [1] * N
+    cdef Py_ssize_t j
+    cdef Py_ssize_t last = N - 1
+    while True:
+        # H2: visit
+        yield from_word(a, base)
+        if a[last] == b[last]:
+            # H4: find j
+            j = N - 2
+            while a[j] == b[j]:
+                j -= 1
+            # H5: increase a_j
+            if j == 0:
+                break
+            a[j] += 1
+            # H6: zero out a_{j+1},...,a_{n-1}
+            b[last] = b[j] + int(a[j] == b[j])
+            j += 1
+            while j < N - 1:
+                a[j] = 0
+                b[j] = b[last]
+                j += 1
+            a[last] = 0
+        else:
+            # H3: increase a_{n-1}
+            a[last] += 1
+
+@cython.wraparound(False)
+@cython.boundscheck(False)
+def _set_partition_block_gen(Py_ssize_t n, Py_ssize_t k, list a):
+    r"""
+    Recursively generate set partitions of ``n`` with fixed block
+    size ``k`` using Algorithm 4.23 from [Rus2003]_.
+    ``a`` is a list of size ``n``.
+
+    EXAMPLES::
+
+        sage: from sage.combinat.set_partition_iterator import _set_partition_block_gen
+        sage: a = list(range(3))
+        sage: for p in _set_partition_block_gen(3, 2, a):
+        ....:     print(p)
+        [0, 1, 0]
+        [0, 1, 1]
+        [0, 0, 1]
+    """
+    cdef Py_ssize_t i
+    if n == k:
+        yield a
+        return
+
+    for i in range(k):
+        a[n-1] = i
+        for P in _set_partition_block_gen(n-1, k, a):
+            yield P
+        a[n-1] = n-1
+    if k > 1:
+        a[n-1] = k-1
+        for P in _set_partition_block_gen(n-1, k-1, a):
+            yield P
+        a[n-1] = n-1
+
+@cython.wraparound(False)
+@cython.boundscheck(False)
+def set_partition_iterator_blocks(base_set, Py_ssize_t k):
+    """
+    A fast iterator for the set partitions of the base set into the
+    specified number of blocks, which returns lists of lists
+    instead of set partitions types.
+
+    EXAMPLES::
+
+        sage: from sage.combinat.set_partition_iterator import set_partition_iterator_blocks
+        sage: list(set_partition_iterator_blocks([1,-1,x], 2))                          # optional - sage.symbolic
+        [[[1, x], [-1]], [[1], [-1, x]], [[1, -1], [x]]]
+    """
+    cdef list base = list(base_set)
+    cdef Py_ssize_t n = len(base)
+    cdef list a = list(range(n))
+    # TODO: implement _set_partition_block_gen as an iterative algorithm
+    for P in _set_partition_block_gen(n, k, a):
+        yield from_word(<list> P, base)