diff --git a/src/sage/combinat/posets/hasse_cython_flint.pyx b/src/sage/combinat/posets/hasse_cython_flint.pyx
index e2b4589b238..28f3f026430 100644
--- a/src/sage/combinat/posets/hasse_cython_flint.pyx
+++ b/src/sage/combinat/posets/hasse_cython_flint.pyx
@@ -20,6 +20,53 @@ 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
+from sage.libs.flint.fmpz_poly_sage cimport Fmpz_poly
+
+
+cpdef Fmpz_poly chain_poly(list positions):
+    r"""
+    Return the chain polynomial of a poset.
+
+    INPUT:
+
+    - ``positions`` -- a list of sets of integers describing the poset, as
+      given by the lazy attribute ``_leq_storage`` of Hasse diagrams
+
+    OUTPUT: a Flint polynomial in one variable over `\ZZ`.
+
+    EXAMPLES::
+
+        sage: from sage.combinat.posets.hasse_cython_flint import chain_poly
+        sage: D = [{0, 1}, {1}]
+        sage: chain_poly(D)
+        3  1 2 1
+        sage: P = posets.TamariLattice(5)
+        sage: H = P._hasse_diagram
+        sage: D = H._leq_storage
+        sage: chain_poly(D)
+        12  1 42 357 1385 3133 4635 4758 3468 1778 612 127 12
+    """
+    cdef Py_ssize_t n = len(positions)
+    cdef Py_ssize_t i, j
+
+    q = Fmpz_poly([0, 1])
+    zero = Fmpz_poly(0)
+    one = Fmpz_poly(1)
+
+    cdef list chain_polys = [zero] * n
+
+    # chain_polys[i] will be the generating function for the
+    # chains with lowest vertex i (in the labelling of the
+    # Hasse diagram).
+    for i in range(n - 1, -1, -1):
+        cpi = q
+        for j in positions[i]:
+            cpi += q * chain_polys[j]
+        chain_polys[i] = cpi
+    total = one
+    for i in range(n):
+        total += chain_polys[i]
+    return total
 
 
 cpdef Matrix_integer_dense moebius_matrix_fast(list positions):
diff --git a/src/sage/combinat/posets/hasse_diagram.py b/src/sage/combinat/posets/hasse_diagram.py
index 0056ad37eed..56b51e96371 100644
--- a/src/sage/combinat/posets/hasse_diagram.py
+++ b/src/sage/combinat/posets/hasse_diagram.py
@@ -27,7 +27,8 @@
 from sage.rings.integer_ring import ZZ
 
 lazy_import('sage.combinat.posets.hasse_cython_flint',
-            ['moebius_matrix_fast', 'coxeter_matrix_fast'])
+            ['moebius_matrix_fast', 'coxeter_matrix_fast',
+             'chain_poly'])
 lazy_import('sage.matrix.constructor', 'matrix')
 lazy_import('sage.rings.finite_rings.finite_field_constructor', 'GF')
 
@@ -2413,6 +2414,23 @@ def chains(self, element_class=list, exclude=None, conversion=None):
         """
         return IncreasingChains(self._leq_storage, element_class, exclude, conversion)
 
+    def chain_polynomial(self):
+        """
+        Return the chain polynomial of the poset.
+
+        The coefficient of `q^k` is the number of chains of `k`
+        elements in the poset. List of coefficients of this polynomial
+        is also called a *f-vector* of the poset.
+
+        EXAMPLES::
+
+            sage: P = posets.ChainPoset(3)
+            sage: H = P._hasse_diagram
+            sage: t = H.chain_polynomial(); t
+            q^3 + 3*q^2 + 3*q + 1
+        """
+        return chain_poly(self._leq_storage)._sage_('q')  # noqa: F821
+
     def is_linear_interval(self, t_min, t_max) -> bool:
         """
         Return whether the interval ``[t_min, t_max]`` is linear.
diff --git a/src/sage/combinat/posets/posets.py b/src/sage/combinat/posets/posets.py
index a111609b4b0..6471163fe00 100644
--- a/src/sage/combinat/posets/posets.py
+++ b/src/sage/combinat/posets/posets.py
@@ -7682,18 +7682,7 @@ def chain_polynomial(self):
             sage: R.chain_polynomial()
             q + 1
         """
-        hasse = self._hasse_diagram
-        q = polygen(ZZ, 'q')
-        one = q.parent().one()
-        hasse_size = hasse.cardinality()
-        chain_polys = [0] * hasse_size
-        # chain_polys[i] will be the generating function for the
-        # chains with topmost vertex i (in the labelling of the
-        # Hasse diagram).
-        for i in range(hasse_size):
-            chain_polys[i] = q + sum(q * chain_polys[j]
-                                     for j in hasse.principal_order_ideal(i))
-        return one + sum(chain_polys)
+        return self._hasse_diagram.chain_polynomial()
 
     def order_polynomial(self):
         r"""