diff --git a/build/pkgs/configure/checksums.ini b/build/pkgs/configure/checksums.ini
index 0fb544d529f..f54312eb0f5 100644
--- a/build/pkgs/configure/checksums.ini
+++ b/build/pkgs/configure/checksums.ini
@@ -1,4 +1,4 @@
 tarball=configure-VERSION.tar.gz
-sha1=b2c2a31975aecc02df78e38082e915419c3354e6
-md5=2441f31fbe9e5a3ca1f6070fcee5e42c
-cksum=4212643350
+sha1=b6d662b1301c2ef342d28ffa22c6f09be4ead821
+md5=e7ed474499cb651b218ef6ab4f1b3285
+cksum=1633939417
diff --git a/build/pkgs/configure/package-version.txt b/build/pkgs/configure/package-version.txt
index 867f16b0119..dc5cb53e6d8 100644
--- a/build/pkgs/configure/package-version.txt
+++ b/build/pkgs/configure/package-version.txt
@@ -1 +1 @@
-f32860b2015ad42b5cb9c30a0af1351e597bf314
+821c5f5d7e9ab1476094459252ae9b42de71cf98
diff --git a/src/sage/combinat/cluster_algebra_quiver/quiver.py b/src/sage/combinat/cluster_algebra_quiver/quiver.py
index 7956cd48e7d..23d4e5c3a62 100644
--- a/src/sage/combinat/cluster_algebra_quiver/quiver.py
+++ b/src/sage/combinat/cluster_algebra_quiver/quiver.py
@@ -36,8 +36,8 @@
 #  Distributed under the terms of the GNU General Public License (GPL)
 #                  https://www.gnu.org/licenses/
 # ****************************************************************************
-
 from copy import copy
+from itertools import product
 
 from sage.structure.sage_object import SageObject
 from sage.rings.integer_ring import ZZ
@@ -46,7 +46,6 @@
 from sage.graphs.graph import Graph
 from sage.graphs.views import EdgesView
 from sage.arith.misc import gcd
-from sage.categories.cartesian_product import cartesian_product
 from sage.misc.misc_c import prod
 from sage.misc.lazy_import import lazy_import
 from sage.rings.rational_field import QQ
@@ -2109,8 +2108,8 @@ def poincare_semistable(self, theta, d):
         mu_d = theta.dot_product(d) / sum(d)
 
         Li = [0 * d]
-        it = (vector(e) for e in cartesian_product([range(d_i + 1)
-                                                    for d_i in d]))
+        it = (vector(e) for e in product(*[range(d_i + 1)
+                                           for d_i in d]))
         Li += [e for e in it if e.dot_product(theta) > mu_d * sum(e)]
         Li.append(d)
         N = len(Li) - 1
diff --git a/src/sage/combinat/constellation.py b/src/sage/combinat/constellation.py
index 8f7274cc256..cbefe06054d 100644
--- a/src/sage/combinat/constellation.py
+++ b/src/sage/combinat/constellation.py
@@ -48,7 +48,8 @@
 # (at your option) any later version.
 #                  https://www.gnu.org/licenses/
 # ****************************************************************************
-from itertools import repeat
+from itertools import repeat, product
+
 from sage.structure.element import parent
 from sage.structure.parent import Parent
 from sage.structure.element import Element
@@ -1408,8 +1409,6 @@ def __iter__(self):
             g1 ('a','d','b')('c')
             g2 ('a','b')('c','d')
         """
-        from sage.misc.mrange import cartesian_product_iterator
-
         if self._cd._length == 1:
             if self._cd._degree == 1:
                 yield self([[0]])
@@ -1417,8 +1416,7 @@ def __iter__(self):
 
         S = self._cd._sym
         profile = list(self._profile)[:-1]
-        for p in cartesian_product_iterator([S.conjugacy_class(pi)
-                                             for pi in profile]):
+        for p in product(*[S.conjugacy_class(pi) for pi in profile]):
             if self._cd._connected and not perms_are_connected(p, self._cd._degree):
                 continue
             c = self._cd(list(p) + [None], check=False)
diff --git a/src/sage/combinat/diagram.py b/src/sage/combinat/diagram.py
index fdbb6912ee1..00b94920726 100644
--- a/src/sage/combinat/diagram.py
+++ b/src/sage/combinat/diagram.py
@@ -23,6 +23,7 @@
 # (at your option) any later version.
 #                  https://www.gnu.org/licenses/
 # ****************************************************************************
+from itertools import product
 
 from sage.categories.infinite_enumerated_sets import InfiniteEnumeratedSets
 from sage.combinat.composition import Composition
@@ -1523,8 +1524,7 @@ def RotheDiagram(w):
 
     N = w.size()
     winv = w.inverse()
-    from sage.misc.mrange import cartesian_product_iterator
-    cells = [c for c in cartesian_product_iterator((range(N), range(N)))
+    cells = [c for c in product(range(N), range(N))
              if c[0] + 1 < winv(c[1] + 1) and c[1] + 1 < w(c[0] + 1)]
 
     return NorthwestDiagram(cells, n_rows=N, n_cols=N, check=False)
diff --git a/src/sage/combinat/posets/posets.py b/src/sage/combinat/posets/posets.py
index a832e6a673f..c31597c72e9 100644
--- a/src/sage/combinat/posets/posets.py
+++ b/src/sage/combinat/posets/posets.py
@@ -288,6 +288,7 @@
 from __future__ import annotations
 from collections import defaultdict
 from copy import copy, deepcopy
+from itertools import product
 from typing import List
 
 from sage.misc.cachefunc import cached_method
@@ -5296,7 +5297,6 @@ def factor(self):
            Discrete Applied Mathematics 15 (1986) 105-110
            :doi:`10.1016/0166-218X(86)90023-5`
         """
-        from sage.misc.mrange import cartesian_product_iterator
         from sage.graphs.graph import Graph
         from sage.misc.flatten import flatten
         dg = self._hasse_diagram
@@ -5327,7 +5327,7 @@ def edge_color(va, vb):
                           if edge_color(x, y) == i0]
                 neigh1 = [z for z in prod_dg.neighbor_iterator(x)
                           if edge_color(x, z) == i1]
-                for x0, x1 in cartesian_product_iterator([neigh0, neigh1]):
+                for x0, x1 in product(neigh0, neigh1):
                     x2 = list(x0)
                     x2[i1] = x1[i1]
                     x2 = tuple(x2)
diff --git a/src/sage/groups/abelian_gps/abelian_group.py b/src/sage/groups/abelian_gps/abelian_group.py
index a72ad3e180c..898ca9cebb5 100644
--- a/src/sage/groups/abelian_gps/abelian_group.py
+++ b/src/sage/groups/abelian_gps/abelian_group.py
@@ -201,6 +201,7 @@
 # (at your option) any later version.
 #                  https://www.gnu.org/licenses/
 # ****************************************************************************
+from itertools import product
 
 from sage.arith.functions import lcm
 from sage.arith.misc import divisors, gcd
@@ -211,7 +212,7 @@
 from sage.matrix.special import diagonal_matrix
 from sage.misc.cachefunc import cached_method
 from sage.misc.misc_c import prod
-from sage.misc.mrange import cartesian_product_iterator, mrange
+from sage.misc.mrange import mrange
 from sage.modules.free_module_element import vector
 from sage.rings.infinity import infinity
 from sage.rings.integer import Integer
@@ -1498,7 +1499,7 @@ def subgroups(self, check=False):
                 # H = the subgroup of *index* H.
                 its = [range(0, H, H // gcd(H, G.gen(i).order()))
                        for i in range(ngens)]
-                for f in cartesian_product_iterator(its):
+                for f in product(*its):
                     verbose("using hom from G to C_%s sending gens to %s" % (H, f))
                     new_sub = []
                     for a in range(ngens):
diff --git a/src/sage/modular/multiple_zeta.py b/src/sage/modular/multiple_zeta.py
index 471b39325c3..22cabe23344 100644
--- a/src/sage/modular/multiple_zeta.py
+++ b/src/sage/modular/multiple_zeta.py
@@ -168,11 +168,11 @@
 from __future__ import annotations
 import numbers
 from typing import Iterator
+from itertools import product
 
 from sage.misc.fast_methods import Singleton
 from sage.structure.richcmp import op_EQ, op_NE
 from sage.structure.element import parent
-from sage.categories.cartesian_product import cartesian_product
 from sage.categories.graded_algebras_with_basis import GradedAlgebrasWithBasis
 from sage.categories.rings import Rings
 from sage.categories.domains import Domains
@@ -578,7 +578,7 @@ def extend_multiplicative_basis(B, n) -> Iterator:
         [((7,),), ((5,), (2,)), ((3,), (2,), (2,))]
     """
     for pi in Partitions(n, min_part=2):
-        yield from cartesian_product([B[i] for i in pi])
+        yield from product(*[B[i] for i in pi])
 
 
 # several classes for the algebra of MZV
diff --git a/src/sage/modules/fg_pid/fgp_module.py b/src/sage/modules/fg_pid/fgp_module.py
index 1ca35212682..f2b1063ce64 100644
--- a/src/sage/modules/fg_pid/fgp_module.py
+++ b/src/sage/modules/fg_pid/fgp_module.py
@@ -211,6 +211,7 @@
 # (at your option) any later version.
 #                  https://www.gnu.org/licenses/
 # ****************************************************************************
+from itertools import product
 
 from sage.modules.module import Module
 from sage.modules.free_module import is_FreeModule
@@ -1835,8 +1836,7 @@ def __iter__(self):
             raise NotImplementedError("currently self must be finite to iterate over")
         B = self.optimized()[0].V().basis_matrix()
         V = self.base_ring()**B.nrows()
-        from sage.misc.mrange import cartesian_product_iterator
-        for a in cartesian_product_iterator([range(k) for k in v]):
+        for a in product(*[range(k) for k in v]):
             b = V(a) * B
             yield self(b)
 
diff --git a/src/sage/schemes/affine/affine_rational_point.py b/src/sage/schemes/affine/affine_rational_point.py
index 739d1795086..4ac17a60d79 100644
--- a/src/sage/schemes/affine/affine_rational_point.py
+++ b/src/sage/schemes/affine/affine_rational_point.py
@@ -50,10 +50,10 @@
 # (at your option) any later version.
 #                  https://www.gnu.org/licenses/
 # ****************************************************************************
+from itertools import product
 
 from sage.rings.integer_ring import ZZ
 from sage.rings.rational_field import QQ
-from sage.misc.mrange import cartesian_product_iterator
 from sage.schemes.generic.scheme import is_Scheme
 
 
@@ -300,7 +300,7 @@ def enum_affine_finite_field(X):
     n = X.codomain().ambient_space().ngens()
     F = X.value_ring()
     pts = []
-    for c in cartesian_product_iterator([F] * n):
+    for c in product(*([F] * n)):
         try:
             pts.append(X(c))
         except Exception:
diff --git a/src/sage/schemes/affine/affine_space.py b/src/sage/schemes/affine/affine_space.py
index 927a065cf60..181221488b2 100644
--- a/src/sage/schemes/affine/affine_space.py
+++ b/src/sage/schemes/affine/affine_space.py
@@ -9,6 +9,7 @@
 #
 #                  https://www.gnu.org/licenses/
 # ****************************************************************************
+from itertools import product
 
 from sage.functions.orthogonal_polys import chebyshev_T, chebyshev_U
 from sage.rings.integer import Integer
@@ -22,7 +23,6 @@
 from sage.categories.fields import Fields
 from sage.categories.number_fields import NumberFields
 from sage.misc.latex import latex
-from sage.misc.mrange import cartesian_product_iterator
 from sage.matrix.constructor import matrix
 from sage.structure.category_object import normalize_names
 from sage.schemes.generic.scheme import AffineScheme
@@ -225,7 +225,7 @@ def __iter__(self):
         AHom = self.point_homset()
         C = AHom.codomain()
 
-        for v in cartesian_product_iterator([R for _ in range(n)]):
+        for v in product(*[R for _ in range(n)]):
             yield C._point(AHom, v, check=False)
 
     def ngens(self):
diff --git a/src/sage/schemes/elliptic_curves/ell_rational_field.py b/src/sage/schemes/elliptic_curves/ell_rational_field.py
index 0ab4d47d1e1..47dc2fc2a64 100644
--- a/src/sage/schemes/elliptic_curves/ell_rational_field.py
+++ b/src/sage/schemes/elliptic_curves/ell_rational_field.py
@@ -49,6 +49,7 @@
 #
 #                  https://www.gnu.org/licenses/
 ##############################################################################
+from itertools import product
 
 from . import constructor
 from . import BSD
@@ -6814,7 +6815,7 @@ def S_integral_x_coords_with_abs_bounded_by(abs_bound):
                 # denom_maxpa is a list of pairs (d,q) where d runs
                 # through possible denominators, and q=p^a is the
                 # maximum prime power divisor of d:
-                denom_maxpa = [(prod(tmp),max(tmp)) for tmp in cartesian_product_iterator(p_pow_alpha)]
+                denom_maxpa = [(prod(tmp), max(tmp)) for tmp in product(*p_pow_alpha)]
 #               The maximum denominator is this (not used):
 #                denom = [prod([pp[-1] for pp in p_pow_alpha],1)]
                 for de,maxpa in denom_maxpa:
@@ -6840,12 +6841,11 @@ def S_integral_x_coords_with_abs_bounded_by(abs_bound):
                 return set(xs)
         # -------------------------------------------------------------------
         # End internal functions ############################################
-        from sage.misc.mrange import cartesian_product_iterator
 
         E = self
         tors_points = E.torsion_points()
 
-        if (r == 0):#only Torsionpoints to consider
+        if (r == 0):  # only Torsionpoints to consider
             int_points = [P for P in tors_points if not P.is_zero()]
             int_points = [P for P in int_points if P[0].is_S_integral(S)]
             if not both_signs:
diff --git a/src/sage/schemes/elliptic_curves/heegner.py b/src/sage/schemes/elliptic_curves/heegner.py
index 793c1652e53..c021dc616e3 100644
--- a/src/sage/schemes/elliptic_curves/heegner.py
+++ b/src/sage/schemes/elliptic_curves/heegner.py
@@ -95,6 +95,7 @@
 # (at your option) any later version.
 #                  https://www.gnu.org/licenses/
 # ****************************************************************************
+from itertools import product
 
 import sage.rings.abc
 import sage.rings.number_field.number_field_element
@@ -5513,11 +5514,10 @@ def kolyvagin_sigma_operator(self, D, c, r, bound=None):
             X = I.cyclic_right_subideals(p, alpha_quaternions[i])
             J_lists.append(dict(enumerate(X)))
 
-        ans = [0]*B.dimension()
-        from sage.misc.mrange import cartesian_product_iterator
-        for v in cartesian_product_iterator([range(1,p+1) for p,_ in F]):
+        ans = [0] * B.dimension()
+        for v in product(*[range(1, p + 1) for p, _ in F]):
             J = J_lists[0][v[0]]
-            for i in range(1,len(J_lists)):
+            for i in range(1, len(J_lists)):
                 J = J.intersection(J_lists[i][v[i]])
             J_theta = tuple(J.theta_series_vector(bound))
             d = theta_dict[J_theta]
diff --git a/src/sage/schemes/elliptic_curves/height.py b/src/sage/schemes/elliptic_curves/height.py
index 97020957574..64d025a6481 100644
--- a/src/sage/schemes/elliptic_curves/height.py
+++ b/src/sage/schemes/elliptic_curves/height.py
@@ -29,6 +29,7 @@
 import numpy
 import math
 import bisect
+from itertools import product
 
 from sage.rings.integer_ring import ZZ
 from sage.rings.rational_field import QQ
@@ -41,7 +42,6 @@
 from sage.rings.real_mpfr import RR
 
 from sage.misc.cachefunc import cached_method
-from sage.misc.mrange import cartesian_product_iterator
 from sage.arith.functions import lcm
 from sage.arith.misc import factorial
 from sage.ext.fast_callable import fast_callable
@@ -685,7 +685,7 @@ def rat_term_CIF(z, try_strict=True):
 
         corner_reals = []
         corner_imags = []
-        for a, b in cartesian_product_iterator([z.real().endpoints(), z.imag().endpoints()]):
+        for a, b in product(z.real().endpoints(), z.imag().endpoints()):
             zz = CDF(a,b)
             u = (two_pi_i_CDF*zz).exp()
             f = u/(1-u)**2
diff --git a/src/sage/schemes/projective/projective_rational_point.py b/src/sage/schemes/projective/projective_rational_point.py
index 31dc0703b5a..0525883886d 100644
--- a/src/sage/schemes/projective/projective_rational_point.py
+++ b/src/sage/schemes/projective/projective_rational_point.py
@@ -45,23 +45,22 @@
 - Raghukul Raman <raghukul.raman01@gmail.com> (2018): Added sieve algorithm
 """
 
-#*****************************************************************************
+# ****************************************************************************
 #       Copyright (C) 2010 William Stein, David Kohel, John Cremona, Charlie Turner
 #
 # 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.
-#                  http://www.gnu.org/licenses/
-#*****************************************************************************
-
+#                  https://www.gnu.org/licenses/
+# ****************************************************************************
+from itertools import product
 
 from sage.arith.misc import gcd, next_prime, previous_prime, crt
 from sage.arith.srange import srange
 from sage.rings.integer_ring import ZZ
 from sage.rings.real_mpfr import RR
 from sage.rings.finite_rings.finite_field_constructor import FiniteField as GF
-from sage.misc.mrange import cartesian_product_iterator
 from sage.misc.misc_c import prod
 from sage.misc.mrange import xmrange
 from sage.schemes.generic.scheme import is_Scheme
@@ -140,7 +139,7 @@ def enum_projective_rational_field(X, B):
     n = X.codomain().ambient_space().ngens()
     zero = (0,) * n
     pts = []
-    for c in cartesian_product_iterator([srange(-B,B+1) for _ in range(n)]):
+    for c in product(*[srange(-B, B + 1) for _ in range(n)]):
         if gcd(c) == 1 and c > zero:
             try:
                 pts.append(X(c))
@@ -296,13 +295,13 @@ def enum_projective_finite_field(X):
     elif not is_ProjectiveSpace(X.codomain().ambient_space()):
         raise TypeError("codomain must be projective space over a finite field")
 
-    n = X.codomain().ambient_space().ngens()-1
+    n = X.codomain().ambient_space().ngens() - 1
     F = X.value_ring()
     pts = []
-    for k in range(n+1):
-        for c in cartesian_product_iterator([F for _ in range(k)]):
+    for k in range(n + 1):
+        for c in product(*[F for _ in range(k)]):
             try:
-                pts.append(X(list(c)+[1]+[0]*(n-k)))
+                pts.append(X(list(c) + [1] + [0] * (n - k)))
             except TypeError:
                 pass
     pts.sort()
diff --git a/src/sage/schemes/projective/projective_space.py b/src/sage/schemes/projective/projective_space.py
index ecb34f8230a..7a06401f39f 100644
--- a/src/sage/schemes/projective/projective_space.py
+++ b/src/sage/schemes/projective/projective_space.py
@@ -76,8 +76,9 @@
 #
 #  Distributed under the terms of the GNU General Public License (GPL)
 #
-#                  http://www.gnu.org/licenses/
+#                  https://www.gnu.org/licenses/
 # ****************************************************************************
+from itertools import product
 
 from sage.arith.misc import gcd, binomial
 
@@ -98,7 +99,6 @@
 from sage.categories.map import Map
 from sage.misc.latex import latex
 from sage.misc.misc_c import prod
-from sage.misc.mrange import cartesian_product_iterator
 from sage.misc.persist import register_unpickle_override
 
 from sage.structure.category_object import normalize_names
@@ -2352,7 +2352,7 @@ def __iter__(self):
         C = PHom.codomain()
 
         for k in range(n + 1): # position of last 1 before the 0's
-            for v in cartesian_product_iterator([R for _ in range(n - k)]):
+            for v in product(*[R for _ in range(n - k)]):
                 yield C._point(PHom, v + one + zero * k, check=False)
 
     def rational_points(self, F=None):
diff --git a/src/sage/symbolic/relation.py b/src/sage/symbolic/relation.py
index 66e88706e84..ea3f5b170ea 100644
--- a/src/sage/symbolic/relation.py
+++ b/src/sage/symbolic/relation.py
@@ -357,7 +357,7 @@
 - William Stein (2007-07-16): added arithmetic with symbolic equations
 
 """
-
+from itertools import product
 import operator
 
 
@@ -1570,7 +1570,6 @@ def solve_mod(eqns, modulus, solution_dict=False):
     from sage.rings.integer import Integer
     from sage.rings.integer_ring import crt_basis
     from sage.structure.element import Expression
-    from sage.misc.mrange import cartesian_product_iterator
     from sage.modules.free_module_element import vector
     from sage.matrix.constructor import matrix
 
@@ -1602,7 +1601,7 @@ def solve_mod(eqns, modulus, solution_dict=False):
 
     ans = []
     if has_solution:
-        for solution in cartesian_product_iterator(solutions):
+        for solution in product(*solutions):
             solution_mat = matrix(Integers(modulus), solution)
             ans.append(tuple(c.dot_product(crt_basis) for c in solution_mat.columns()))
 
@@ -1680,12 +1679,10 @@ def _solve_mod_prime_power(eqns, p, m, vars):
         sage: [sorted(_solve_mod_prime_power([x^2==41], 10, i, [x]))[0][0] for i in [1..13]]
         [1, 21, 71, 1179, 2429, 47571, 1296179, 8703821, 26452429, 526452429,
         13241296179, 19473547571, 2263241296179]
-
     """
     from sage.rings.finite_rings.integer_mod_ring import Integers
     from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
     from sage.modules.free_module_element import vector
-    from sage.misc.mrange import cartesian_product_iterator
 
     mrunning = 1
     ans = []
@@ -1695,10 +1692,10 @@ def _solve_mod_prime_power(eqns, p, m, vars):
         S = PolynomialRing(R, len(vars), vars)
         eqns_mod = [S(eq) for eq in eqns]
         if mi == 0:
-            possibles = cartesian_product_iterator([range(len(R)) for _ in range(len(vars))])
+            possibles = product(*[range(len(R)) for _ in range(len(vars))])
         else:
-            shifts = cartesian_product_iterator([range(p) for _ in range(len(vars))])
-            pairs = cartesian_product_iterator([shifts, ans])
+            shifts = product(*[range(p) for _ in range(len(vars))])
+            pairs = product(shifts, ans)
             possibles = (tuple(vector(t) + vector(shift) * (mrunning // p))
                          for shift, t in pairs)
         ans = list(t for t in possibles if all(e(*t) == 0 for e in eqns_mod))