diff --git a/src/sage/combinat/bijectionist.py b/src/sage/combinat/bijectionist.py
index 42e3c4f1fcc..e5ba051562b 100644
--- a/src/sage/combinat/bijectionist.py
+++ b/src/sage/combinat/bijectionist.py
@@ -614,7 +614,7 @@ def set_constant_blocks(self, P):
         P = sorted(self._sorter["A"](p) for p in P)
         for p in P:
             for a in p:
-                self._P._union(p[0], a)
+                self._P.join(p[0], a)
 
         self._compute_possible_block_values()
 
@@ -1678,7 +1678,7 @@ def merge_until_split():
                         try:
                             solution = different_values(tP[i1], tP[i2])
                         except StopIteration:
-                            tmp_P._union(tP[i1], tP[i2])
+                            tmp_P.join(tP[i1], tP[i2])
                             if len(multiple_preimages[tZ]) == 2:
                                 del multiple_preimages[tZ]
                             else:
@@ -1772,14 +1772,14 @@ def possible_values(self, p=None, optimal=False):
         # convert input to set of block representatives
         blocks = set()
         if p in self._A:
-            blocks.add(self._P._find(p))
+            blocks.add(self._P.find(p))
         elif isinstance(p, list):  # TODO: this looks very brittle
             for p1 in p:
                 if p1 in self._A:
-                    blocks.add(self._P._find(p1))
+                    blocks.add(self._P.find(p1))
                 elif isinstance(p1, list):
                     for p2 in p1:
-                        blocks.add(self._P._find(p2))
+                        blocks.add(self._P.find(p2))
 
         if optimal:
             if self._bmilp is None:
@@ -1941,9 +1941,9 @@ def _find_counterexample(self, P, s0, d, on_blocks):
 
             # try to find a solution which has a different
             # subdistribution on d than s0
-            z_in_d = sum(d[p] * bmilp._x[self._P._find(p), z]
+            z_in_d = sum(d[p] * bmilp._x[self._P.find(p), z]
                          for p in P
-                         if z in self._possible_block_values[self._P._find(p)])
+                         if z in self._possible_block_values[self._P.find(p)])
 
             # it is sufficient to require that z occurs less often as
             # a value among {a | d[a] == 1} than it does in
@@ -2191,14 +2191,14 @@ def _preprocess_intertwining_relations(self):
             # the blocks of the elements of the preimage
             updated_images = defaultdict(set)  # (p_1,...,p_k) to {a_1,....}
             for a_tuple, image_set in images.items():
-                representatives = tuple(P._find(a) for a in a_tuple)
+                representatives = tuple(P.find(a) for a in a_tuple)
                 updated_images[representatives].update(image_set)
 
             # merge blocks
             for a_tuple, image_set in updated_images.items():
                 image = image_set.pop()
                 while image_set:
-                    P._union(image, image_set.pop())
+                    P.join(image, image_set.pop())
                     something_changed = True
                 # we keep a representative
                 image_set.add(image)
@@ -2525,7 +2525,7 @@ def __init__(self, bijectionist: Bijectionist, solutions=None):
         self._solution_cache = []
         if solutions is not None:
             for solution in solutions:
-                self._add_solution({(P._find(a), z): value
+                self._add_solution({(P.find(a), z): value
                                     for (a, z), value in solution.items()})
 
     def show(self, variables=True):
@@ -2813,7 +2813,7 @@ def add_alpha_beta_constraints(self):
         Z_dict = {z: i for i, z in enumerate(Z)}
 
         for a in self._bijectionist._A:
-            p = self._bijectionist._P._find(a)
+            p = self._bijectionist._P.find(a)
             for z in self._bijectionist._possible_block_values[p]:
                 w_index = W_dict[self._bijectionist._alpha(a)]
                 z_index = Z_dict[z]
@@ -2867,7 +2867,7 @@ def add_distribution_constraints(self):
             tA_sum = [zero] * len(Z_dict)
             tZ_sum = [zero] * len(Z_dict)
             for a in tA:
-                p = self._bijectionist._P._find(a)
+                p = self._bijectionist._P.find(a)
                 for z in self._bijectionist._possible_block_values[p]:
                     tA_sum[Z_dict[z]] += self._x[p, z]
             for z in tZ:
@@ -2940,8 +2940,8 @@ def add_intertwining_relation_constraints(self):
                     continue
                 a = pi_rho.pi(*a_tuple)
                 if a in A:
-                    p_tuple = tuple(P._find(a) for a in a_tuple)
-                    p = P._find(a)
+                    p_tuple = tuple(P.find(a) for a in a_tuple)
+                    p = P.find(a)
                     if (p_tuple, p) not in pi_blocks:
                         pi_blocks.add((p_tuple, p))
                         for z_tuple in itertools.product(*[tZ[p] for p in p_tuple]):
@@ -3008,7 +3008,7 @@ def add_quadratic_relation_constraints(self):
                 z0 = phi(p)
                 assert all(phi(a) == z0 for a in block), "phi must be constant on the block %s" % block
                 for z in self._bijectionist._possible_block_values[p]:
-                    p0 = P._find(psi(z))
+                    p0 = P.find(psi(z))
                     if z0 in self._bijectionist._possible_block_values[p0]:
                         c = self._x[p, z] - self._x[p0, z0]
                         if c.is_zero():
@@ -3046,7 +3046,7 @@ def add_homomesic_constraints(self):
         tZ = self._bijectionist._possible_block_values
 
         def sum_q(q):
-            return sum(sum(z * self._x[P._find(a), z] for z in tZ[P._find(a)])
+            return sum(sum(z * self._x[P.find(a), z] for z in tZ[P.find(a)])
                        for a in q)
         q0 = Q[0]
         v0 = sum_q(q0)
@@ -3091,9 +3091,9 @@ def _disjoint_set_roots(d):
 
         sage: from sage.combinat.bijectionist import _disjoint_set_roots
         sage: d = DisjointSet('abcde')
-        sage: d.union("a", "b")
-        sage: d.union("a", "c")
-        sage: d.union("e", "d")
+        sage: d.join("a", "b")
+        sage: d.join("a", "c")
+        sage: d.join("e", "d")
         sage: _disjoint_set_roots(d)
         dict_keys(['a', 'e'])
     """
diff --git a/src/sage/combinat/designs/designs_pyx.pyx b/src/sage/combinat/designs/designs_pyx.pyx
index 9bf3305e92d..155e2d04902 100644
--- a/src/sage/combinat/designs/designs_pyx.pyx
+++ b/src/sage/combinat/designs/designs_pyx.pyx
@@ -492,7 +492,7 @@ def is_group_divisible_design(groups,blocks,v,G=None,K=None,lambd=1,verbose=Fals
         for i in range(n):
             for j in range(i + 1, n):
                 if matrix[i * n + j] == 0:
-                    groups._union(i, j)
+                    groups.join(i, j)
         groups = list(groups.root_to_elements_dict().values())
 
     # Group sizes are element of G
diff --git a/src/sage/combinat/designs/incidence_structures.py b/src/sage/combinat/designs/incidence_structures.py
index e0d0dbd0cbe..cb1cfbcbc2f 100644
--- a/src/sage/combinat/designs/incidence_structures.py
+++ b/src/sage/combinat/designs/incidence_structures.py
@@ -1029,7 +1029,7 @@ def is_connected(self) -> bool:
         for B in self._blocks:
             x = B[0]
             for i in range(1, len(B)):
-                D._union(x, B[i])
+                D.join(x, B[i])
         return D.number_of_subsets() == 1
 
     def is_simple(self) -> bool:
diff --git a/src/sage/combinat/posets/hasse_diagram.py b/src/sage/combinat/posets/hasse_diagram.py
index e0f15d86bf1..2732692074b 100644
--- a/src/sage/combinat/posets/hasse_diagram.py
+++ b/src/sage/combinat/posets/hasse_diagram.py
@@ -3230,7 +3230,7 @@ def fill_to_interval(S):
                 if part:  # Skip empty parts
                     c = part[0]
                     for e in fill_to_interval(part):
-                        cong._union(e, c)
+                        cong.join(e, c)
             t = cong.number_of_subsets()
 
             # Following is needed for cases like
@@ -3238,28 +3238,28 @@ def fill_to_interval(S):
             for c in list(cong):
                 r = c[0]
                 for v in fill_to_interval(c):
-                    cong._union(r, v)
+                    cong.join(r, v)
 
-        todo = {cong._find(e) for part in parts for e in part}
+        todo = {cong.find(e) for part in parts for e in part}
 
         while todo:
 
             # First check if we should stop now.
             for a, b in stop_pairs:
-                if cong._find(a) == cong._find(b):
+                if cong.find(a) == cong.find(b):
                     return None
 
             # We take one block and try to find as big interval
             # as possible to unify as a new block by the quadrilateral
             # argument.
-            block = sorted(cong.root_to_elements_dict()[cong._find(todo.pop())])
+            block = sorted(cong.root_to_elements_dict()[cong.find(todo.pop())])
 
             b = block[-1]
             for a in block:  # Quadrilateral up
                 for c in self.neighbor_out_iterator(a):
                     if c not in block:
                         d = jn[c, b]
-                        if cong._find(d) != cong._find(c):
+                        if cong.find(d) != cong.find(c):
                             break
                 else:
                     continue
@@ -3271,7 +3271,7 @@ def fill_to_interval(S):
                     for d in self.neighbor_in_iterator(b):
                         if d not in block:
                             c = mt[d, a]
-                            if cong._find(c) != cong._find(d):
+                            if cong.find(c) != cong.find(d):
                                 break
                     else:
                         continue
@@ -3289,10 +3289,10 @@ def fill_to_interval(S):
             # recursive process. In particular it may also combine to
             # [a, b] block we just used.
             while c is not None:
-                newblock = cong._find(c)
+                newblock = cong.find(c)
                 for i in self.interval(c, d):
-                    cong._union(newblock, i)
-                C = cong.root_to_elements_dict()[cong._find(newblock)]
+                    cong.join(newblock, i)
+                C = cong.root_to_elements_dict()[cong.find(newblock)]
                 mins = [i for i in C if all(i_ not in C for i_ in self.neighbor_in_iterator(i))]
                 maxs = [i for i in C if all(i_ not in C for i_ in self.neighbor_out_iterator(i))]
                 c = None  # To stop loop, if this is not changed below.
@@ -3305,7 +3305,7 @@ def fill_to_interval(S):
                         d = jn[d, m]
 
             # This removes duplicates from todo.
-            todo = {cong._find(x) for x in todo}
+            todo = {cong.find(x) for x in todo}
 
         return cong
 
diff --git a/src/sage/combinat/set_partition.py b/src/sage/combinat/set_partition.py
index 06b4267dc23..a5c899be44e 100644
--- a/src/sage/combinat/set_partition.py
+++ b/src/sage/combinat/set_partition.py
@@ -2356,7 +2356,7 @@ def from_arcs(self, arcs, n):
         """
         P = DisjointSet(range(1, n + 1))
         for i, j in arcs:
-            P._union(i, j)
+            P.join(i, j)
         return self.element_class(self, P)
 
     def from_rook_placement_gamma(self, rooks, n):
diff --git a/src/sage/combinat/words/finite_word.py b/src/sage/combinat/words/finite_word.py
index b74185f6824..bc863cf125e 100644
--- a/src/sage/combinat/words/finite_word.py
+++ b/src/sage/combinat/words/finite_word.py
@@ -5201,16 +5201,16 @@ def overlap_partition(self, other, delay=0, p=None, involution=None):
         S = zip(islice(self, int(delay), None), other)
         if involution is None:
             for a, b in S:
-                p._union(a, b)
+                p.join(a, b)
         elif isinstance(involution, WordMorphism):
             for a, b in S:
-                p._union(a, b)
+                p.join(a, b)
                 # take the first letter of the word
-                p._union(involution(a)[0], involution(b)[0])
+                p.join(involution(a)[0], involution(b)[0])
         elif callable(involution):
             for a, b in S:
-                p._union(a, b)
-                p._union(involution(a), involution(b))
+                p.join(a, b)
+                p.join(involution(a), involution(b))
         else:
             raise TypeError("involution (=%s) must be callable" % involution)
         return p
diff --git a/src/sage/graphs/connectivity.pyx b/src/sage/graphs/connectivity.pyx
index 4642f72d00c..04452fdda7c 100644
--- a/src/sage/graphs/connectivity.pyx
+++ b/src/sage/graphs/connectivity.pyx
@@ -2540,7 +2540,7 @@ def spqr_tree(G, algorithm="Hopcroft_Tarjan", solver=None, verbose=0,
         if cocycles_count[fe] == 2 and len(virtual_edge_to_cycles[fe]) == 2:
             # This virtual edge is only between 2 cycles
             C1, C2 = virtual_edge_to_cycles[fe]
-            DS._union(C1, C2)
+            DS.join(C1, C2)
             cycles_list[C1].delete_edge(fe)
             cycles_list[C2].delete_edge(fe)
             cocycles_count[fe] -= 2
diff --git a/src/sage/graphs/generators/random.py b/src/sage/graphs/generators/random.py
index 7f31fa9758f..23ac88813bd 100644
--- a/src/sage/graphs/generators/random.py
+++ b/src/sage/graphs/generators/random.py
@@ -576,12 +576,12 @@ def RandomBlockGraph(m, k, kmax=None, incidence_structure=False, seed=None):
         # structure to keep a unique identifier per merged vertices
         DS = DisjointSet([i for u in B for i in B[u]])
         for u, v in T.edges(sort=True, labels=0):
-            DS._union(choice(B[u]), choice(B[v]))
+            DS.join(choice(B[u]), choice(B[v]))
 
         # We relabel vertices in the range [0, m*(k-1)] and build the incidence
         # structure
         new_label = {root: i for i, root in enumerate(DS.root_to_elements_dict())}
-        IS = [[new_label[DS._find(v)] for v in B[u]] for u in B]
+        IS = [[new_label[DS.find(v)] for v in B[u]] for u in B]
 
     if incidence_structure:
         return IS
diff --git a/src/sage/graphs/generic_graph.py b/src/sage/graphs/generic_graph.py
index ea261f6aaf9..a4ca9b04e25 100644
--- a/src/sage/graphs/generic_graph.py
+++ b/src/sage/graphs/generic_graph.py
@@ -12657,11 +12657,11 @@ def contract_edges(self, edges):
         DS = DisjointSet(self.vertex_iterator())
 
         for u, v, label in edge_list:
-            DS._union(u, v)
+            DS.join(u, v)
 
         self.delete_edges(edge_list)
         edges_incident = []
-        vertices = [v for v in vertices if v != DS._find(v)]
+        vertices = [v for v in vertices if v != DS.find(v)]
         if self.is_directed():
             for v in vertices:
                 out_edges = self.edge_boundary([v])
@@ -12674,8 +12674,8 @@ def contract_edges(self, edges):
                 self.delete_vertex(v)
 
         for (u, v, label) in edges_incident:
-            root_u = DS._find(u)
-            root_v = DS._find(v)
+            root_u = DS.find(u)
+            root_v = DS.find(v)
             if root_v != root_u or self.allows_loops():
                 self.add_edge(root_u, root_v, label)
 
diff --git a/src/sage/graphs/graph_decompositions/modular_decomposition.py b/src/sage/graphs/graph_decompositions/modular_decomposition.py
index d80084831b6..412b12ec36c 100644
--- a/src/sage/graphs/graph_decompositions/modular_decomposition.py
+++ b/src/sage/graphs/graph_decompositions/modular_decomposition.py
@@ -456,7 +456,7 @@ def gamma_classes(graph):
             e = frozenset([v1, v])
             for vi in component[1:]:
                 ei = frozenset([vi, v])
-                pieces._union(e, ei)
+                pieces.join(e, ei)
     return {frozenset(chain.from_iterable(loe)): loe for loe in pieces}
 
 
diff --git a/src/sage/graphs/graph_decompositions/tree_decomposition.pyx b/src/sage/graphs/graph_decompositions/tree_decomposition.pyx
index 30f4856f5d6..7910591fea8 100644
--- a/src/sage/graphs/graph_decompositions/tree_decomposition.pyx
+++ b/src/sage/graphs/graph_decompositions/tree_decomposition.pyx
@@ -233,7 +233,7 @@ def is_valid_tree_decomposition(G, T):
         for Xi in X:
             for Xj in T.neighbor_iterator(Xi):
                 if Xj in X:
-                    D._union(Xi, Xj)
+                    D.join(Xi, Xj)
         if D.number_of_subsets() > 1:
             return False
 
diff --git a/src/sage/graphs/partial_cube.py b/src/sage/graphs/partial_cube.py
index 8242775da02..39c21b17d37 100644
--- a/src/sage/graphs/partial_cube.py
+++ b/src/sage/graphs/partial_cube.py
@@ -335,9 +335,9 @@ def is_partial_cube(G, certificate=False):
             if diff not in neighbors:
                 return fail
             neighbor = neighbors[diff]
-            unionfind._union(contracted.edge_label(v, w),
+            unionfind.join(contracted.edge_label(v, w),
                             contracted.edge_label(root, neighbor))
-            unionfind._union(contracted.edge_label(w, v),
+            unionfind.join(contracted.edge_label(w, v),
                             contracted.edge_label(neighbor, root))
             labeled.add_edge(v, w)
 
@@ -356,13 +356,13 @@ def is_partial_cube(G, certificate=False):
                 if vi == wi:
                     return fail
                 if newgraph.has_edge(vi, wi):
-                    unionfind._union(newgraph.edge_label(vi, wi), t)
+                    unionfind.join(newgraph.edge_label(vi, wi), t)
                 else:
                     newgraph.add_edge(vi, wi, t)
         contracted = newgraph
 
     # Make a digraph with edges labeled by the equivalence classes in unionfind
-    g = DiGraph({v: {w: unionfind._find((v, w)) for w in G[v]} for v in G})
+    g = DiGraph({v: {w: unionfind.find((v, w)) for w in G[v]} for v in G})
 
     # Associates to a vertex the token that acts on it, and check that
     # no two edges on a single vertex have the same label
diff --git a/src/sage/graphs/spanning_tree.pyx b/src/sage/graphs/spanning_tree.pyx
index 646fd651cae..fdfb70edd41 100644
--- a/src/sage/graphs/spanning_tree.pyx
+++ b/src/sage/graphs/spanning_tree.pyx
@@ -345,17 +345,17 @@ def kruskal_iterator(G, by_weight=True, weight_function=None, check_weight=False
             yield from G.edge_iterator()
             return
 
-    cdef DisjointSet_of_hashables union_find = DisjointSet_of_hashables(G)
+    cdef DisjointSet_of_hashables disjointset = DisjointSet_of_hashables(G)
     by_weight, weight_function = G._get_weight_function(by_weight=by_weight,
                                                         weight_function=weight_function,
                                                         check_weight=check_weight)
-    yield from kruskal_iterator_from_edges(G.edge_iterator(), union_find,
+    yield from kruskal_iterator_from_edges(G.edge_iterator(), disjointset,
                                            by_weight=by_weight,
                                            weight_function=weight_function,
                                            check_weight=False)
 
 
-def kruskal_iterator_from_edges(edges, union_find, by_weight=True,
+def kruskal_iterator_from_edges(edges, disjointset, by_weight=True,
                                 weight_function=None, check_weight=False):
     """
     Return an iterator implementation of Kruskal algorithm on list of edges.
@@ -364,7 +364,7 @@ def kruskal_iterator_from_edges(edges, union_find, by_weight=True,
 
     - ``edges`` -- list of edges
 
-    - ``union_find`` -- a
+    - ``disjointset`` -- a
       :class:`~sage.sets.disjoint_set.DisjointSet_of_hashables` encoding a
       forest
 
@@ -416,13 +416,13 @@ def kruskal_iterator_from_edges(edges, union_find, by_weight=True,
     # Kruskal's algorithm
     for e in edges:
         # acyclic test via union-find
-        u = union_find._find(e[0])
-        v = union_find._find(e[1])
+        u = disjointset.find(e[0])
+        v = disjointset.find(e[1])
         if u != v:
             yield e
             # merge the trees
-            union_find._union(u, v)
-            if union_find.number_of_subsets() == 1:
+            disjointset.join(u, v)
+            if disjointset.number_of_subsets() == 1:
                 return
 
 
@@ -657,7 +657,7 @@ def filter_kruskal_iterator(G, threshold=10000, by_weight=True, weight_function=
     # Parameter to  equally divide edges with weight equal the to pivot
     cdef bint ch = True
     # Data structure to record the vertices in each tree of the forest
-    cdef DisjointSet_of_hashables union_find = DisjointSet_of_hashables(g)
+    cdef DisjointSet_of_hashables disjointset = DisjointSet_of_hashables(g)
 
     #
     # Iteratively partition the list of edges
@@ -668,12 +668,12 @@ def filter_kruskal_iterator(G, threshold=10000, by_weight=True, weight_function=
         if end - begin < threshold:
             # Filter edges connecting vertices of a same tree
             L = [edges[e_index[i]] for i in range(begin, end + 1)
-                 if union_find._find(edges[e_index[i]][0]) != union_find._find(edges[e_index[i]][1])]
-            yield from kruskal_iterator_from_edges(L, union_find,
+                 if disjointset.find(edges[e_index[i]][0]) != disjointset.find(edges[e_index[i]][1])]
+            yield from kruskal_iterator_from_edges(L, disjointset,
                                                    by_weight=by_weight,
                                                    weight_function=weight_function,
                                                    check_weight=False)
-            if union_find.number_of_subsets() == 1:
+            if disjointset.number_of_subsets() == 1:
                 return
             continue
 
@@ -874,8 +874,8 @@ def boruvka(G, by_weight=True, weight_function=None, check_weight=True, check=Fa
         # each pair of components (trees of the forest), as well as cheapest
         # active edge incident to a component.
         for e, e_weight in edge_list:
-            component1 = partitions._find(e[0])
-            component2 = partitions._find(e[1])
+            component1 = partitions.find(e[0])
+            component2 = partitions.find(e[1])
 
             if component1 != component2:
                 if component1 in cheapest:
@@ -904,11 +904,11 @@ def boruvka(G, by_weight=True, weight_function=None, check_weight=True, check=Fa
         # possible
         for v in cheapest:
             e, e_weight = cheapest[v]
-            component1 = partitions._find(e[0])
-            component2 = partitions._find(e[1])
+            component1 = partitions.find(e[0])
+            component2 = partitions.find(e[1])
 
             if component1 != component2:
-                partitions._union(component1, component2)
+                partitions.join(component1, component2)
                 T.append(e)
                 numConComp = numConComp - 1
 
@@ -1373,7 +1373,7 @@ def edge_disjoint_spanning_trees(G, k, by_weight=False, weight_function=None, ch
     for x, y, _ in G.edges(sort=by_weight, key=weight_function):
         # {x, y} is edge e0 in the algorithm
 
-        if partition[0]._find(x) == partition[0]._find(y):
+        if partition[0].find(x) == partition[0].find(y):
             # x and y are in a same clump. That is x and y are in a same tree
             # in every forest Fi. We proceed with the next edge.
             continue
@@ -1404,7 +1404,7 @@ def edge_disjoint_spanning_trees(G, k, by_weight=False, weight_function=None, ch
             fe = frozenset(e)
             i = (edge_index[fe] % k) + 1
             v, w = e
-            if partition[i]._find(v) != partition[i]._find(w):
+            if partition[i].find(v) != partition[i].find(w):
                 # v and w are in different subtrees of Fi. We have detected an
                 # augmenting sequence since we can join the two subtrees.
                 augmenting_sequence_found = True
@@ -1436,7 +1436,7 @@ def edge_disjoint_spanning_trees(G, k, by_weight=False, weight_function=None, ch
 
         if augmenting_sequence_found:
             # We perform the corresponding augmentation
-            partition[i]._union(v, w)
+            partition[i].join(v, w)
 
             while fe in edge_label:
                 F[edge_index[fe]].delete_edge(fe)
@@ -1450,7 +1450,7 @@ def edge_disjoint_spanning_trees(G, k, by_weight=False, weight_function=None, ch
 
         else:
             # x and y are in a same tree in every Fi, so in a same clump
-            partition[0]._union(x, y)
+            partition[0].join(x, y)
 
     res = [F[i] for i in range(1, k + 1) if F[i].size() == G.order() - 1]
     if len(res) != k:
diff --git a/src/sage/groups/perm_gps/permgroup.py b/src/sage/groups/perm_gps/permgroup.py
index 4b33600db9f..09236f6502d 100644
--- a/src/sage/groups/perm_gps/permgroup.py
+++ b/src/sage/groups/perm_gps/permgroup.py
@@ -5347,7 +5347,7 @@ def __init__(self, gens, action, domain, gap_group=None, category=None, canonica
             for g_orbit in g_orbits:
                 for o in g_orbit:
                     for i in range(1, len(o)):
-                        D.union(o[0], o[i])
+                        D.join(o[0], o[i])
             self._orbits = tuple(tuple(o) for o in D)
 
         PermutationGroup_generic.__init__(self, gens=gens,
diff --git a/src/sage/matroids/graphic_matroid.py b/src/sage/matroids/graphic_matroid.py
index bb31b96d153..5cfd1481d77 100644
--- a/src/sage/matroids/graphic_matroid.py
+++ b/src/sage/matroids/graphic_matroid.py
@@ -286,7 +286,7 @@ def _rank(self, X):
         # This counts components:
         DS_vertices = DisjointSet(vertices)
         for (u, v, l) in edges:
-            DS_vertices._union(u, v)
+            DS_vertices.join(u, v)
         return (len(vertices) - DS_vertices.number_of_subsets())
 
     # Representation:
@@ -665,7 +665,7 @@ def _corank(self, X):
         not_our_edges = self.groundset_to_edges(self._groundset.difference(X))
         DS_vertices = DisjointSet(all_vertices)
         for u, v, l in not_our_edges:
-            DS_vertices._union(u, v)
+            DS_vertices.join(u, v)
         return len(X) - (DS_vertices.number_of_subsets() - Integer(1))
 
     def _is_circuit(self, X):
@@ -776,8 +776,8 @@ def _max_independent(self, X):
         our_set = set()
         DS_vertices = DisjointSet(vertices)
         for (u, v, l) in edges:
-            if DS_vertices._find(u) != DS_vertices._find(v):
-                DS_vertices._union(u, v)
+            if DS_vertices.find(u) != DS_vertices.find(v):
+                DS_vertices.join(u, v)
                 our_set.add(l)
         return frozenset(our_set)
 
@@ -811,13 +811,13 @@ def _max_coindependent(self, X):
         our_set = set()
         DS_vertices = DisjointSet(all_vertices)
         for (u, v, l) in not_our_edges:
-            DS_vertices._union(u, v)
+            DS_vertices.join(u, v)
 
         for (u, v, l) in edges:
-            if DS_vertices._find(u) == DS_vertices._find(v):
+            if DS_vertices.find(u) == DS_vertices.find(v):
                 our_set.add(l)
             else:
-                DS_vertices._union(u, v)
+                DS_vertices.join(u, v)
         return frozenset(our_set)
 
     def _circuit(self, X):
@@ -875,8 +875,8 @@ def _circuit(self, X):
         DS_vertices = DisjointSet(vertices)
         for u, v, l in edges:
             edge_set.add((u, v, l))
-            if DS_vertices._find(u) != DS_vertices._find(v):
-                DS_vertices._union(u, v)
+            if DS_vertices.find(u) != DS_vertices.find(v):
+                DS_vertices.join(u, v)
             else:
                 break
         else:
diff --git a/src/sage/rings/polynomial/multi_polynomial_sequence.py b/src/sage/rings/polynomial/multi_polynomial_sequence.py
index 38e531aa6e5..9c2e098b0fc 100644
--- a/src/sage/rings/polynomial/multi_polynomial_sequence.py
+++ b/src/sage/rings/polynomial/multi_polynomial_sequence.py
@@ -1062,11 +1062,11 @@ def connected_components(self):
         DS = DisjointSet(set().union(*vss))
         for u, *vs in vss:
             for v in vs:
-                DS._union(u, v)
+                DS.join(u, v)
 
         Ps = {}  # map root element -> polynomials in this component
         for f, vs in zip(self, vss):
-            r = DS._find(vs[0])
+            r = DS.find(vs[0])
             if r in Ps:
                 Ps[r].append(f)
             else:
diff --git a/src/sage/sets/disjoint_set.pxd b/src/sage/sets/disjoint_set.pxd
index 5a6a112250e..9be7051dc67 100644
--- a/src/sage/sets/disjoint_set.pxd
+++ b/src/sage/sets/disjoint_set.pxd
@@ -12,16 +12,26 @@
 from sage.groups.perm_gps.partn_ref.data_structures cimport OrbitPartition
 from sage.structure.sage_object cimport SageObject
 
+cpdef DisjointSet(arg)
+
 cdef class DisjointSet_class(SageObject):
     cdef OrbitPartition *_nodes
+    cpdef cardinality(self)
+    cpdef number_of_subsets(self)
 
 cdef class DisjointSet_of_integers(DisjointSet_class):
-    cpdef _find(self, int i)
-    cpdef _union(self, int i, int j)
+    cpdef find(self, int i)
+    cpdef join(self, int i, int j)
+    cpdef root_to_elements_dict(self)
+    cpdef element_to_root_dict(self)
+    cpdef to_digraph(self)
 
 cdef class DisjointSet_of_hashables(DisjointSet_class):
     cdef list _int_to_el
     cdef dict _el_to_int
     cdef DisjointSet_of_integers _d
-    cpdef _find(self, e)
-    cpdef _union(self, e, f)
+    cpdef find(self, e)
+    cpdef join(self, e, f)
+    cpdef root_to_elements_dict(self)
+    cpdef element_to_root_dict(self)
+    cpdef to_digraph(self)
diff --git a/src/sage/sets/disjoint_set.pyx b/src/sage/sets/disjoint_set.pyx
index 246f6f83cfd..53b451bf442 100644
--- a/src/sage/sets/disjoint_set.pyx
+++ b/src/sage/sets/disjoint_set.pyx
@@ -22,9 +22,9 @@ Disjoint set of integers from ``0`` to ``n - 1``::
     sage: s = DisjointSet(6)
     sage: s
     {{0}, {1}, {2}, {3}, {4}, {5}}
-    sage: s.union(2, 4)
-    sage: s.union(1, 3)
-    sage: s.union(5, 1)
+    sage: s.join(2, 4)
+    sage: s.join(1, 3)
+    sage: s.join(5, 1)
     sage: s
     {{0}, {1, 3, 5}, {2, 4}}
     sage: s.find(3)
@@ -39,9 +39,9 @@ Disjoint set of hashables objects::
     sage: d = DisjointSet('abcde')
     sage: d
     {{'a'}, {'b'}, {'c'}, {'d'}, {'e'}}
-    sage: d.union('a','b')
-    sage: d.union('b','c')
-    sage: d.union('c','d')
+    sage: d.join('a','b')
+    sage: d.join('b','c')
+    sage: d.join('c','d')
     sage: d
     {{'a', 'b', 'c', 'd'}, {'e'}}
     sage: d.find('c')
@@ -63,8 +63,7 @@ from sage.structure.sage_object cimport SageObject
 from cpython.object cimport PyObject_RichCompare
 from sage.groups.perm_gps.partn_ref.data_structures cimport *
 
-
-def DisjointSet(arg):
+cpdef DisjointSet(arg):
     r"""
     Constructs a disjoint set where each element of ``arg`` is in its
     own set. If ``arg`` is an integer, then the disjoint set returned is
@@ -77,7 +76,7 @@ def DisjointSet(arg):
     - :meth:`~sage.sets.disjoint_set.DisjointSet_of_hashables.find` --
       Determine which set a particular element is in.
 
-    - :meth:`~sage.sets.disjoint_set.DisjointSet_of_hashables.union` --
+    - :meth:`~sage.sets.disjoint_set.DisjointSet_of_hashables.join` --
       Combine or merge two sets into a single set.
 
     REFERENCES:
@@ -141,7 +140,6 @@ def DisjointSet(arg):
     else:
         return DisjointSet_of_hashables(arg)
 
-
 cdef class DisjointSet_class(SageObject):
     r"""
     Common class and methods for :class:`DisjointSet_of_integers` and
@@ -154,19 +152,19 @@ cdef class DisjointSet_class(SageObject):
         EXAMPLES::
 
             sage: e = DisjointSet(5)
-            sage: e.union(2,4); e._repr_()
+            sage: e.join(2,4); e._repr_()
             '{{0}, {1}, {2, 4}, {3}}'
             sage: e = DisjointSet(5)
-            sage: e.union(4,2); e._repr_()
+            sage: e.join(4,2); e._repr_()
             '{{0}, {1}, {2, 4}, {3}}'
 
         ::
 
             sage: e = DisjointSet(range(5))
-            sage: e.union(2,4); e._repr_()
+            sage: e.join(2,4); e._repr_()
             '{{0}, {1}, {2, 4}, {3}}'
             sage: e = DisjointSet(range(5))
-            sage: e.union(4,2); e._repr_()
+            sage: e.join(4,2); e._repr_()
             '{{0}, {1}, {2, 4}, {3}}'
         """
         res = []
@@ -183,7 +181,7 @@ cdef class DisjointSet_class(SageObject):
         EXAMPLES::
 
             sage: d = DisjointSet(4)
-            sage: d.union(2,0)
+            sage: d.join(2,0)
             sage: sorted(d)
             [[0, 2], [1], [3]]
 
@@ -211,10 +209,10 @@ cdef class DisjointSet_class(SageObject):
 
         ::
 
-            sage: d.union(0,3)
-            sage: d.union(3,4)
-            sage: e.union(4,0)
-            sage: e.union(3,0)
+            sage: d.join(0,3)
+            sage: d.join(3,4)
+            sage: e.join(4,0)
+            sage: e.join(3,0)
             sage: e == d
             True
 
@@ -229,12 +227,12 @@ cdef class DisjointSet_class(SageObject):
 
             sage: d = DisjointSet('abcde')
             sage: e = DisjointSet('abcde')
-            sage: d.union('a','b')
-            sage: d.union('b','c')
-            sage: e.union('c','a')
+            sage: d.join('a','b')
+            sage: d.join('b','c')
+            sage: e.join('c','a')
             sage: e == d
             False
-            sage: e.union('a','b')
+            sage: e.join('a','b')
             sage: e == d
             True
         """
@@ -246,7 +244,7 @@ cdef class DisjointSet_class(SageObject):
             return NotImplemented
         return PyObject_RichCompare(s, t, op)
 
-    def cardinality(self):
+    cpdef cardinality(self):
         r"""
         Return the number of elements in ``self``, *not* the number of subsets.
 
@@ -255,19 +253,19 @@ cdef class DisjointSet_class(SageObject):
             sage: d = DisjointSet(5)
             sage: d.cardinality()
             5
-            sage: d.union(2, 4)
+            sage: d.join(2, 4)
             sage: d.cardinality()
             5
             sage: d = DisjointSet(range(5))
             sage: d.cardinality()
             5
-            sage: d.union(2, 4)
+            sage: d.join(2, 4)
             sage: d.cardinality()
             5
         """
         return self._nodes.degree
 
-    def number_of_subsets(self):
+    cpdef number_of_subsets(self):
         r"""
         Return the number of subsets in ``self``.
 
@@ -276,13 +274,13 @@ cdef class DisjointSet_class(SageObject):
             sage: d = DisjointSet(5)
             sage: d.number_of_subsets()
             5
-            sage: d.union(2, 4)
+            sage: d.join(2, 4)
             sage: d.number_of_subsets()
             4
             sage: d = DisjointSet(range(5))
             sage: d.number_of_subsets()
             5
-            sage: d.union(2, 4)
+            sage: d.join(2, 4)
             sage: d.number_of_subsets()
             4
         """
@@ -297,8 +295,8 @@ cdef class DisjointSet_of_integers(DisjointSet_class):
         sage: d = DisjointSet(5)
         sage: d
         {{0}, {1}, {2}, {3}, {4}}
-        sage: d.union(2,4)
-        sage: d.union(0,2)
+        sage: d.join(2,4)
+        sage: d.join(0,2)
         sage: d
         {{0, 2, 4}, {1}, {3}}
         sage: d.find(2)
@@ -312,7 +310,7 @@ cdef class DisjointSet_of_integers(DisjointSet_class):
 
     ::
 
-        sage: a.union(3,4)
+        sage: a.join(3,4)
         sage: a == loads(dumps(a))
         True
     """
@@ -363,8 +361,8 @@ cdef class DisjointSet_of_integers(DisjointSet_class):
 
         ::
 
-            sage: d.union(2,4)
-            sage: d.union(1,3)
+            sage: d.join(2,4)
+            sage: d.join(1,3)
             sage: d.__reduce__()
             (<built-in function DisjointSet>, (5,), [0, 1, 2, 1, 2])
         """
@@ -379,10 +377,10 @@ cdef class DisjointSet_of_integers(DisjointSet_class):
             sage: d = DisjointSet(5)
             sage: d.__getstate__()
             [0, 1, 2, 3, 4]
-            sage: d.union(2,3)
+            sage: d.join(2,3)
             sage: d.__getstate__()
             [0, 1, 2, 2, 4]
-            sage: d.union(3,0)
+            sage: d.join(3,0)
             sage: d.__getstate__()
             [2, 1, 2, 2, 4]
 
@@ -390,22 +388,23 @@ cdef class DisjointSet_of_integers(DisjointSet_class):
         distinct order::
 
             sage: d = DisjointSet(5)
-            sage: d.union(0,3)
+            sage: d.join(0,3)
             sage: d.__getstate__()
             [0, 1, 2, 0, 4]
-            sage: d.union(2,0)
+            sage: d.join(2,0)
             sage: d.__getstate__()
             [0, 1, 0, 0, 4]
         """
-        cdef list l = []
+        cdef Py_ssize_t card = self.cardinality()
+        cdef list l = [None] * card
         cdef int i
-        for i in range(self.cardinality()):
-            l.append(self._nodes.parent[i])
+        for i in range(card):
+            l[i] = self._nodes.parent[i]
         return l
 
     def __setstate__(self, l):
         r"""
-        Merge the nodes ``i`` and ``l[i]`` (using union) for ``i`` in
+        Merge the nodes ``i`` and ``l[i]`` (using join) for ``i`` in
         ``range(len(l))``.
 
         INPUT:
@@ -441,78 +440,41 @@ cdef class DisjointSet_of_integers(DisjointSet_class):
             {{0, 1, 2, 3}, {4}}
         """
         for i, parent in enumerate(l):
-            self.union(parent, i)
+            self.join(parent, i)
 
-    def find(self, int i):
+    cpdef find(self, int i):
         r"""
         Return the representative of the set that ``i`` currently belongs to.
 
         INPUT:
 
-        - ``i`` -- element in ``self``
+        - ``i`` -- element in ``self`` (no input checking)
 
         EXAMPLES::
 
             sage: e = DisjointSet(5)
-            sage: e.union(4,2); e
+            sage: e.join(4,2); e
             {{0}, {1}, {2, 4}, {3}}
             sage: e.find(2)
             4
             sage: e.find(4)
             4
-            sage: e.union(1,3); e
+            sage: e.join(1,3); e
             {{0}, {1, 3}, {2, 4}}
             sage: e.find(1)
             1
             sage: e.find(3)
             1
-            sage: e.union(3,2); e
-            {{0}, {1, 2, 3, 4}}
-            sage: [e.find(i) for i in range(5)]
-            [0, 1, 1, 1, 1]
-            sage: e.find(5)
-            Traceback (most recent call last):
-            ...
-            ValueError: i(=5) must be between 0 and 4
-        """
-        card = self.cardinality()
-        if i < 0 or i>= card:
-            raise ValueError('i(=%s) must be between 0 and %s' % (i, card - 1))
-        return OP_find(self._nodes, i)
-
-    cpdef _find(self, int i):
-        r"""
-        Return the representative of the set that ``i`` currently belongs to.
-
-        INPUT:
-
-        - ``i`` -- element in ``self`` (no input checking)
-
-        EXAMPLES::
-
-            sage: e = DisjointSet(5)
-            sage: e._union(4,2); e
-            {{0}, {1}, {2, 4}, {3}}
-            sage: e._find(2)
-            4
-            sage: e._find(4)
-            4
-            sage: e._union(1,3); e
-            {{0}, {1, 3}, {2, 4}}
-            sage: e._find(1)
-            1
-            sage: e._find(3)
-            1
-            sage: e._union(3,2); e
+            sage: e.join(3,2); e
             {{0}, {1, 2, 3, 4}}
             sage: [e.find(i) for i in range(5)]
             [0, 1, 1, 1, 1]
-            sage: e._find(5)  # no input checking
+            sage: e.find(5)  # no input checking
             0
         """
         return OP_find(self._nodes, i)
 
-    def union(self, int i, int j):
+    cpdef join(self, int i, int j):
         r"""
         Combine the set of ``i`` and the set of ``j`` into one.
 
@@ -528,52 +490,18 @@ cdef class DisjointSet_of_integers(DisjointSet_class):
 
             sage: d = DisjointSet(5); d
             {{0}, {1}, {2}, {3}, {4}}
-            sage: d.union(0,1); d
-            {{0, 1}, {2}, {3}, {4}}
-            sage: d.union(2,4); d
-            {{0, 1}, {2, 4}, {3}}
-            sage: d.union(1,4); d
-            {{0, 1, 2, 4}, {3}}
-            sage: d.union(1,5)
-            Traceback (most recent call last):
-            ...
-            ValueError: j(=5) must be between 0 and 4
-        """
-        cdef int card = self._nodes.degree
-        if i < 0 or i >= card:
-            raise ValueError('i(=%s) must be between 0 and %s' % (i, card - 1))
-        if j < 0 or j >= card:
-            raise ValueError('j(=%s) must be between 0 and %s' % (j, card - 1))
-        OP_join(self._nodes, i, j)
-
-    cpdef _union(self, int i, int j):
-        r"""
-        Combine the set of ``i`` and the set of ``j`` into one.
-
-        All elements in those two sets will share the same representative
-        that can be gotten using find.
-
-        INPUT:
-
-        - ``i`` -- element in ``self`` (no input checking)
-        - ``j`` -- element in ``self`` (no input checking)
-
-        EXAMPLES::
-
-            sage: d = DisjointSet(5); d
-            {{0}, {1}, {2}, {3}, {4}}
-            sage: d._union(0,1); d
+            sage: d.join(0,1); d
             {{0, 1}, {2}, {3}, {4}}
-            sage: d._union(2,4); d
+            sage: d.join(2,4); d
             {{0, 1}, {2, 4}, {3}}
-            sage: d._union(1,4); d
+            sage: d.join(1,4); d
             {{0, 1, 2, 4}, {3}}
-            sage: d._union(1,5); d
+            sage: d.join(1,5); d  # no input checking
             {{0, 1, 2, 4}, {3}}
         """
         OP_join(self._nodes, i, j)
 
-    def root_to_elements_dict(self):
+    cpdef root_to_elements_dict(self):
         r"""
         Return the dictionary where the keys are the roots of ``self`` and the
         values are the elements in the same set as the root.
@@ -583,17 +511,17 @@ cdef class DisjointSet_of_integers(DisjointSet_class):
             sage: d = DisjointSet(5)
             sage: sorted(d.root_to_elements_dict().items())
             [(0, [0]), (1, [1]), (2, [2]), (3, [3]), (4, [4])]
-            sage: d.union(2,3)
+            sage: d.join(2,3)
             sage: sorted(d.root_to_elements_dict().items())
             [(0, [0]), (1, [1]), (2, [2, 3]), (4, [4])]
-            sage: d.union(3,0)
+            sage: d.join(3,0)
             sage: sorted(d.root_to_elements_dict().items())
             [(1, [1]), (2, [0, 2, 3]), (4, [4])]
             sage: d
             {{0, 2, 3}, {1}, {4}}
         """
         cdef dict s = {}
-        cdef int i
+        cdef int i, o
         for i in range(self.cardinality()):
             o = self.find(i)
             if o not in s:
@@ -601,7 +529,7 @@ cdef class DisjointSet_of_integers(DisjointSet_class):
             s[o].append(i)
         return s
 
-    def element_to_root_dict(self):
+    cpdef element_to_root_dict(self):
         r"""
         Return the dictionary where the keys are the elements of ``self`` and
         the values are their representative inside a list.
@@ -609,8 +537,8 @@ cdef class DisjointSet_of_integers(DisjointSet_class):
         EXAMPLES::
 
             sage: d = DisjointSet(5)
-            sage: d.union(2,3)
-            sage: d.union(4,1)
+            sage: d.join(2,3)
+            sage: d.join(4,1)
             sage: e = d.element_to_root_dict(); e
             {0: 0, 1: 4, 2: 2, 3: 2, 4: 4}
             sage: WordMorphism(e)                                                       # needs sage.combinat
@@ -622,7 +550,7 @@ cdef class DisjointSet_of_integers(DisjointSet_class):
             d[i] = self.find(i)
         return d
 
-    def to_digraph(self):
+    cpdef to_digraph(self):
         r"""
         Return the current digraph of ``self`` where `(a,b)` is an oriented
         edge if `b` is the parent of `a`.
@@ -630,9 +558,9 @@ cdef class DisjointSet_of_integers(DisjointSet_class):
         EXAMPLES::
 
             sage: d = DisjointSet(5)
-            sage: d.union(2,3)
-            sage: d.union(4,1)
-            sage: d.union(3,4)
+            sage: d.join(2,3)
+            sage: d.join(4,1)
+            sage: d.join(3,4)
             sage: d
             {{0}, {1, 2, 3, 4}}
             sage: g = d.to_digraph(); g                                                 # needs sage.graphs
@@ -640,18 +568,18 @@ cdef class DisjointSet_of_integers(DisjointSet_class):
             sage: g.edges(sort=True)                                                    # needs sage.graphs
             [(0, 0, None), (1, 2, None), (2, 2, None), (3, 2, None), (4, 2, None)]
 
-        The result depends on the ordering of the union::
+        The result depends on the ordering of the join::
 
             sage: d = DisjointSet(5)
-            sage: d.union(1,2)
-            sage: d.union(1,3)
-            sage: d.union(1,4)
+            sage: d.join(1,2)
+            sage: d.join(1,3)
+            sage: d.join(1,4)
             sage: d
             {{0}, {1, 2, 3, 4}}
             sage: d.to_digraph().edges(sort=True)                                       # needs sage.graphs
             [(0, 0, None), (1, 1, None), (2, 1, None), (3, 1, None), (4, 1, None)]
         """
-        d = {i: [self._nodes.parent[i]] for i in range(self.cardinality())}
+        cdef dict d = {i: [self._nodes.parent[i]] for i in range(self.cardinality())}
         from sage.graphs.digraph import DiGraph
         return DiGraph(d)
 
@@ -663,7 +591,7 @@ cdef class DisjointSet_of_hashables(DisjointSet_class):
 
         sage: d = DisjointSet('abcde'); d
         {{'a'}, {'b'}, {'c'}, {'d'}, {'e'}}
-        sage: d.union('a', 'c'); d
+        sage: d.join('a', 'c'); d
         {{'a', 'c'}, {'b'}, {'d'}, {'e'}}
         sage: d.find('a')
         'a'
@@ -676,7 +604,7 @@ cdef class DisjointSet_of_hashables(DisjointSet_class):
 
     ::
 
-        sage: a.union('a','c')
+        sage: a.join('a','c')
         sage: a == loads(dumps(a))
         True
     """
@@ -728,8 +656,8 @@ cdef class DisjointSet_of_hashables(DisjointSet_class):
 
          ::
 
-            sage: d.union(2,4)
-            sage: d.union(1,3)
+            sage: d.join(2,4)
+            sage: d.join(1,3)
             sage: d.__reduce__()
             (<built-in function DisjointSet>,
              ([0, 1, 2, 3, 4],),
@@ -746,10 +674,10 @@ cdef class DisjointSet_of_hashables(DisjointSet_class):
             sage: d = DisjointSet('abcde')
             sage: d.__getstate__()
             [('a', 'a'), ('b', 'b'), ('c', 'c'), ('d', 'd'), ('e', 'e')]
-            sage: d.union('c','d')
+            sage: d.join('c','d')
             sage: d.__getstate__()
             [('a', 'a'), ('b', 'b'), ('c', 'c'), ('d', 'c'), ('e', 'e')]
-            sage: d.union('d','a')
+            sage: d.join('d','a')
             sage: d.__getstate__()
             [('a', 'c'), ('b', 'b'), ('c', 'c'), ('d', 'c'), ('e', 'e')]
 
@@ -757,7 +685,7 @@ cdef class DisjointSet_of_hashables(DisjointSet_class):
         different order::
 
             sage: d = DisjointSet('abcde')
-            sage: d.union('d','c')
+            sage: d.join('d','c')
             sage: d.__getstate__()
             [('a', 'a'), ('b', 'b'), ('c', 'd'), ('d', 'd'), ('e', 'e')]
         """
@@ -792,32 +720,32 @@ cdef class DisjointSet_of_hashables(DisjointSet_class):
             {{'a', 'b', 'c', 'd', 'e'}}
         """
         for a, b in l:
-            self.union(a, b)
+            self.join(a, b)
 
-    def find(self, e):
+    cpdef find(self, e):
         r"""
         Return the representative of the set that ``e`` currently belongs to.
 
         INPUT:
 
-        - ``e`` -- element in ``self``
+        - ``e`` -- element in ``self`` (no input checking)
 
         EXAMPLES::
 
             sage: e = DisjointSet(range(5))
-            sage: e.union(4,2); e
+            sage: e.join(4,2); e
             {{0}, {1}, {2, 4}, {3}}
             sage: e.find(2)
             4
             sage: e.find(4)
             4
-            sage: e.union(1,3); e
+            sage: e.join(1,3); e
             {{0}, {1, 3}, {2, 4}}
             sage: e.find(1)
             1
             sage: e.find(3)
             1
-            sage: e.union(3,2); e
+            sage: e.join(3,2); e
             {{0}, {1, 2, 3, 4}}
             sage: [e.find(i) for i in range(5)]
             [0, 1, 1, 1, 1]
@@ -826,74 +754,11 @@ cdef class DisjointSet_of_hashables(DisjointSet_class):
             ...
             KeyError: 5
         """
-        cdef int i = self._el_to_int[e]
-        cdef int r = self._d.find(i)
-        return self._int_to_el[r]
-
-    cpdef _find(self, e):
-        r"""
-        Return the representative of the set that ``e`` currently belongs to.
-
-        INPUT:
-
-        - ``e`` -- element in ``self`` (no input checking)
-
-        EXAMPLES::
-
-            sage: e = DisjointSet(range(5))
-            sage: e._union(4,2); e
-            {{0}, {1}, {2, 4}, {3}}
-            sage: e._find(2)
-            4
-            sage: e._find(4)
-            4
-            sage: e._union(1,3); e
-            {{0}, {1, 3}, {2, 4}}
-            sage: e._find(1)
-            1
-            sage: e._find(3)
-            1
-            sage: e._union(3,2); e
-            {{0}, {1, 2, 3, 4}}
-            sage: [e._find(i) for i in range(5)]
-            [0, 1, 1, 1, 1]
-            sage: e._find(5)
-            Traceback (most recent call last):
-            ...
-            KeyError: 5
-        """
-        cdef int i = self._el_to_int[e]
-        cdef int r = self._d._find(i)
+        cdef int i = <int> self._el_to_int[e]
+        cdef int r = <int> self._d.find(i)
         return self._int_to_el[r]
 
-    def union(self, e, f):
-        r"""
-        Combine the set of ``e`` and the set of ``f`` into one.
-
-        All elements in those two sets will share the same representative
-        that can be gotten using find.
-
-        INPUT:
-
-        - ``e`` -- element in ``self``
-        - ``f`` -- element in ``self``
-
-        EXAMPLES::
-
-            sage: e = DisjointSet('abcde'); e
-            {{'a'}, {'b'}, {'c'}, {'d'}, {'e'}}
-            sage: e.union('a','b'); e
-            {{'a', 'b'}, {'c'}, {'d'}, {'e'}}
-            sage: e.union('c','e'); e
-            {{'a', 'b'}, {'c', 'e'}, {'d'}}
-            sage: e.union('b','e'); e
-            {{'a', 'b', 'c', 'e'}, {'d'}}
-        """
-        cdef int i = self._el_to_int[e]
-        cdef int j = self._el_to_int[f]
-        self._d.union(i, j)
-
-    cpdef _union(self, e, f):
+    cpdef join(self, e, f):
         r"""
         Combine the set of ``e`` and the set of ``f`` into one.
 
@@ -909,18 +774,18 @@ cdef class DisjointSet_of_hashables(DisjointSet_class):
 
             sage: e = DisjointSet('abcde'); e
             {{'a'}, {'b'}, {'c'}, {'d'}, {'e'}}
-            sage: e._union('a','b'); e
+            sage: e.join('a','b'); e
             {{'a', 'b'}, {'c'}, {'d'}, {'e'}}
-            sage: e._union('c','e'); e
+            sage: e.join('c','e'); e
             {{'a', 'b'}, {'c', 'e'}, {'d'}}
-            sage: e._union('b','e'); e
+            sage: e.join('b','e'); e
             {{'a', 'b', 'c', 'e'}, {'d'}}
         """
-        cdef int i = self._el_to_int[e]
-        cdef int j = self._el_to_int[f]
-        self._d._union(i, j)
+        cdef int i = <int> self._el_to_int[e]
+        cdef int j = <int> self._el_to_int[f]
+        self._d.join(i, j)
 
-    def root_to_elements_dict(self):
+    cpdef root_to_elements_dict(self):
         r"""
         Return the dictionary where the keys are the roots of ``self`` and the
         values are the elements in the same set.
@@ -928,13 +793,13 @@ cdef class DisjointSet_of_hashables(DisjointSet_class):
         EXAMPLES::
 
             sage: d = DisjointSet(range(5))
-            sage: d.union(2,3)
-            sage: d.union(4,1)
+            sage: d.join(2,3)
+            sage: d.join(4,1)
             sage: e = d.root_to_elements_dict()
             sage: sorted(e.items())
             [(0, [0]), (2, [2, 3]), (4, [1, 4])]
         """
-        s = {}
+        cdef dict s = {}
         for e in self._int_to_el:
             r = self.find(e)
             if r not in s:
@@ -942,7 +807,7 @@ cdef class DisjointSet_of_hashables(DisjointSet_class):
             s[r].append(e)
         return s
 
-    def element_to_root_dict(self):
+    cpdef element_to_root_dict(self):
         r"""
         Return the dictionary where the keys are the elements of ``self`` and
         the values are their representative inside a list.
@@ -950,20 +815,20 @@ cdef class DisjointSet_of_hashables(DisjointSet_class):
         EXAMPLES::
 
             sage: d = DisjointSet(range(5))
-            sage: d.union(2,3)
-            sage: d.union(4,1)
+            sage: d.join(2,3)
+            sage: d.join(4,1)
             sage: e = d.element_to_root_dict()
             sage: sorted(e.items())
             [(0, 0), (1, 4), (2, 2), (3, 2), (4, 4)]
             sage: WordMorphism(e)                                                       # needs sage.combinat
             WordMorphism: 0->0, 1->4, 2->2, 3->2, 4->4
         """
-        d = {}
+        cdef dict d = {}
         for a in self._int_to_el:
             d[a] = self.find(a)
         return d
 
-    def to_digraph(self):
+    cpdef to_digraph(self):
         r"""
         Return the current digraph of ``self`` where `(a,b)` is an oriented
         edge if `b` is the parent of `a`.
@@ -971,9 +836,9 @@ cdef class DisjointSet_of_hashables(DisjointSet_class):
         EXAMPLES::
 
             sage: d = DisjointSet(range(5))
-            sage: d.union(2,3)
-            sage: d.union(4,1)
-            sage: d.union(3,4)
+            sage: d.join(2,3)
+            sage: d.join(4,1)
+            sage: d.join(3,4)
             sage: d
             {{0}, {1, 2, 3, 4}}
             sage: g = d.to_digraph(); g                                                 # needs sage.graphs
@@ -981,18 +846,19 @@ cdef class DisjointSet_of_hashables(DisjointSet_class):
             sage: g.edges(sort=True)                                                    # needs sage.graphs
             [(0, 0, None), (1, 2, None), (2, 2, None), (3, 2, None), (4, 2, None)]
 
-        The result depends on the ordering of the union::
+        The result depends on the ordering of the join::
 
             sage: d = DisjointSet(range(5))
-            sage: d.union(1,2)
-            sage: d.union(1,3)
-            sage: d.union(1,4)
+            sage: d.join(1,2)
+            sage: d.join(1,3)
+            sage: d.join(1,4)
             sage: d
             {{0}, {1, 2, 3, 4}}
             sage: d.to_digraph().edges(sort=True)                                       # needs sage.graphs
             [(0, 0, None), (1, 1, None), (2, 1, None), (3, 1, None), (4, 1, None)]
         """
-        d = {}
+        cdef dict d = {}
+        cdef int i
         for i in range(self.cardinality()):
             e = self._int_to_el[i]
             p = self._int_to_el[self._nodes.parent[i]]
diff --git a/src/sage/tensor/modules/tensor_free_submodule.py b/src/sage/tensor/modules/tensor_free_submodule.py
index ef6ac449514..b334d7d65cf 100644
--- a/src/sage/tensor/modules/tensor_free_submodule.py
+++ b/src/sage/tensor/modules/tensor_free_submodule.py
@@ -297,11 +297,11 @@ def is_coarsening_of(self_sym_list, other_sym_list):
             for index_set in self_sym_list:
                 i = index_set[0]
                 for j in index_set[1:]:
-                    S._union(i, j)
+                    S.join(i, j)
             for index_set in other_sym_list:
-                i = S._find(index_set[0])
+                i = S.find(index_set[0])
                 for j in index_set[1:]:
-                    if S._find(j) != i:
+                    if S.find(j) != i:
                         return False
             return True