diff --git a/src/doc/en/reference/references/index.rst b/src/doc/en/reference/references/index.rst
index b7725c7848b..93e142baf93 100644
--- a/src/doc/en/reference/references/index.rst
+++ b/src/doc/en/reference/references/index.rst
@@ -1658,6 +1658,11 @@ REFERENCES:
 .. [CIA] CIA Factbook 09
          https://www.cia.gov/library/publications/the-world-factbook/
 
+.. [CJ2022] \M. Chang, C. Jefferson, *Disjoint direct product decomposition
+            of permutation groups*, Journal of Symbolic Computation (2022),
+            Volume 108, pages 1-16. :doi:`10.1016/j.jsc.2021.04.003`.
+            Preprint: :arxiv:`2004.11618v3`.
+
 .. [CK1986] \R. Calderbank, W.M. Kantor,
             *The geometry of two-weight codes*,
             Bull. London Math. Soc. 18(1986) 97-122.
diff --git a/src/sage/groups/perm_gps/permgroup.py b/src/sage/groups/perm_gps/permgroup.py
index 2d265216519..d97aad93b03 100644
--- a/src/sage/groups/perm_gps/permgroup.py
+++ b/src/sage/groups/perm_gps/permgroup.py
@@ -1581,6 +1581,96 @@ def smallest_moved_point(self):
         p = self._libgap_().SmallestMovedPoint()
         return self._domain_from_gap[Integer(p)]
 
+    @cached_method
+    def disjoint_direct_product_decomposition(self):
+        r"""
+        Return the finest partition of the underlying set such that ``self``
+        is isomorphic to the direct product of the projections of ``self``
+        onto each part of the partition. Each part is a union of orbits
+        of ``self``.
+
+        The algorithm is from [CJ2022]_, which runs in time polynomial in
+        `n \cdot |X|`, where `n` is the degree of the group and `|X|` is
+        the size of a generating set, see Theorem 4.5.
+
+        EXAMPLES:
+
+        The example from the original paper::
+
+            sage: H = PermutationGroup([[(1,2,3),(7,9,8),(10,12,11)],[(4,5,6),(7,8,9),(10,11,12)],[(5,6),(8,9),(11,12)],[(7,8,9),(10,11,12)]])
+            sage: S = H.disjoint_direct_product_decomposition(); S
+            {{1, 2, 3}, {4, 5, 6, 7, 8, 9, 10, 11, 12}}
+            sage: A = libgap.Stabilizer(H, list(S[0]), libgap.OnTuples); A
+            Group([ (7,8,9)(10,11,12), (5,6)(8,9)(11,12), (4,5,6)(7,8,9)(10,11,12) ])
+            sage: B = libgap.Stabilizer(H, list(S[1]), libgap.OnTuples); B
+            Group([ (1,2,3) ])
+            sage: T = PermutationGroup(gap_group=libgap.DirectProduct(A,B))
+            sage: T.is_isomorphic(H)
+            True
+            sage: PermutationGroup(PermutationGroup(gap_group=A).gens(),domain=list(S[1])).disjoint_direct_product_decomposition()
+            {{4, 5, 6, 7, 8, 9, 10, 11, 12}}
+            sage: PermutationGroup(PermutationGroup(gap_group=B).gens(),domain=list(S[0])).disjoint_direct_product_decomposition()
+            {{1, 2, 3}}
+
+        An example with a different domain::
+
+            sage: PermutationGroup([[('a','c','d'),('b','e')]]).disjoint_direct_product_decomposition()
+            {{'a', 'c', 'd'}, {'b', 'e'}}
+            sage: PermutationGroup([[('a','c','d','b','e')]]).disjoint_direct_product_decomposition()
+            {{'a', 'b', 'c', 'd', 'e'}}
+
+        Counting the number of "connected" permutation groups of degree `n`::
+
+            sage: seq = [sum(1 for G in SymmetricGroup(n).conjugacy_classes_subgroups() if len(G.disjoint_direct_product_decomposition()) == 1) for n in range(1,8)]; seq
+            [1, 1, 2, 6, 6, 27, 20]
+            sage: oeis(seq) # optional -- internet
+            0: A005226: Number of atomic species of degree n; also number of connected permutation groups of degree n.
+        """
+        from sage.combinat.set_partition import SetPartition
+        from sage.sets.disjoint_set import DisjointSet
+        H = self._libgap_()
+        # sort each orbit and order list by smallest element of each orbit
+        O = libgap.List([libgap.ShallowCopy(orbit) for orbit in libgap.Orbits(H)])
+        for orbit in O:
+            libgap.Sort(orbit)
+        O.Sort()
+        num_orbits = len(O)
+        OrbitMapping = dict()
+        for i in range(num_orbits):
+            for x in O[i]:
+                OrbitMapping[x] = i
+        C = libgap.StabChain(H, libgap.Concatenation(O))
+        X = libgap.StrongGeneratorsStabChain(C)
+        P = DisjointSet(num_orbits)
+        R = libgap.List([])
+        identity = libgap.Identity(H)
+        for i in range(num_orbits-1):
+            libgap.Append(R, O[i])
+            Xp = libgap.List([])
+            while True:
+                try:
+                    if libgap.IsSubset(O[i], C['orbit']):
+                        C = C['stabilizer']
+                    else:
+                        break
+                except ValueError:
+                    break
+            for x in X:
+                xs = libgap.SiftedPermutation(C, x)
+                if xs != identity:
+                    libgap.Add(Xp, xs)
+                    if libgap.RestrictedPerm(xs, O[i+1]) != identity:
+                        cj = OrbitMapping[libgap.SmallestMovedPoint(libgap.RestrictedPerm(x, R))]
+                        P.union(i+1, cj)
+                else:
+                    libgap.Add(Xp, x)
+            X = Xp
+        return SetPartition([
+            [self._domain_from_gap[Integer(x)]
+             for i in part
+             for x in O[i]] for part in P] +
+             [[x] for x in self.fixed_points()])
+
     def representative_action(self, x, y):
         r"""
         Return an element of ``self`` that maps `x` to `y` if it exists.