diff --git a/src/sage/algebras/fusion_rings/shm_managers.pyx b/src/sage/algebras/fusion_rings/shm_managers.pyx
index 91aba7ba59f..84c3bb22be4 100644
--- a/src/sage/algebras/fusion_rings/shm_managers.pyx
+++ b/src/sage/algebras/fusion_rings/shm_managers.pyx
@@ -22,7 +22,7 @@ from multiprocessing import shared_memory
 from sage.algebras.fusion_rings.poly_tup_engine cimport poly_to_tup, tup_fixes_sq, _flatten_coeffs
 from sage.rings.integer cimport Integer
 from sage.rings.rational cimport Rational
-from sage.rings.polynomial.multi_polynomial_libsingular cimport MPolynomial_libsingular
+from sage.rings.polynomial.multi_polynomial cimport MPolynomial_libsingular
 from sage.rings.polynomial.polydict cimport ETuple
 
 import numpy as np
diff --git a/src/sage/combinat/kazhdan_lusztig.py b/src/sage/combinat/kazhdan_lusztig.py
index e5a1456b7de..0c213ea9013 100644
--- a/src/sage/combinat/kazhdan_lusztig.py
+++ b/src/sage/combinat/kazhdan_lusztig.py
@@ -20,7 +20,7 @@
 #*****************************************************************************
 
 
-from sage.rings.polynomial.polynomial_element import is_Polynomial
+from sage.rings.polynomial.polynomial_element import Polynomial
 from sage.misc.cachefunc import cached_method
 from sage.rings.polynomial.laurent_polynomial import LaurentPolynomial
 from sage.structure.sage_object import SageObject
@@ -76,7 +76,7 @@ def __init__(self, W, q, trace=False):
         self._trace = trace
         self._one = W.one()
         self._base_ring = q.parent()
-        if is_Polynomial(q):
+        if isinstance(q, Polynomial):
             self._base_ring_type = "polynomial"
         elif isinstance(q, LaurentPolynomial):
             self._base_ring_type = "laurent"
diff --git a/src/sage/combinat/schubert_polynomial.py b/src/sage/combinat/schubert_polynomial.py
index 7e9164db5de..89c857c8efc 100644
--- a/src/sage/combinat/schubert_polynomial.py
+++ b/src/sage/combinat/schubert_polynomial.py
@@ -80,10 +80,9 @@
 from sage.misc.cachefunc import cached_method
 from sage.rings.integer import Integer
 from sage.rings.integer_ring import ZZ
-from sage.rings.polynomial.infinite_polynomial_element import InfinitePolynomial_sparse
+from sage.rings.polynomial.infinite_polynomial_element import InfinitePolynomial
 from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
-from sage.rings.polynomial.multi_polynomial import is_MPolynomial
-
+from sage.rings.polynomial.multi_polynomial import MPolynomial
 import sage.libs.symmetrica.all as symmetrica
 
 
@@ -150,7 +149,7 @@ def expand(self):
             x0
         """
         p = symmetrica.t_SCHUBERT_POLYNOM(self)
-        if not is_MPolynomial(p):
+        if not isinstance(p, MPolynomial):
             R = PolynomialRing(self.parent().base_ring(), 1, 'x0')
             p = R(p)
         return p
@@ -462,9 +461,9 @@ def _element_constructor_(self, x):
         elif isinstance(x, Permutation):
             perm = x.remove_extra_fixed_points()
             return self._from_dict({perm: self.base_ring().one()})
-        elif is_MPolynomial(x):
+        elif isinstance(x, MPolynomial):
             return symmetrica.t_POLYNOM_SCHUBERT(x)
-        elif isinstance(x, InfinitePolynomial_sparse):
+        elif isinstance(x, InfinitePolynomial):
             R = x.polynomial().parent()
             # massage the term order to be what symmetrica expects
             S = PolynomialRing(R.base_ring(),
diff --git a/src/sage/combinat/sf/sfa.py b/src/sage/combinat/sf/sfa.py
index b435764a746..fd86637c9a7 100644
--- a/src/sage/combinat/sf/sfa.py
+++ b/src/sage/combinat/sf/sfa.py
@@ -217,8 +217,8 @@
 from sage.rings.integer import Integer
 from sage.rings.infinity import infinity
 from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
-from sage.rings.polynomial.polynomial_element import is_Polynomial
-from sage.rings.polynomial.multi_polynomial import is_MPolynomial
+from sage.rings.polynomial.polynomial_element import Polynomial
+from sage.rings.polynomial.multi_polynomial import MPolynomial
 from sage.combinat.partition import _Partitions, Partitions, Partitions_n, Partition
 from sage.categories.hopf_algebras import HopfAlgebras
 from sage.categories.hopf_algebras_with_basis import HopfAlgebrasWithBasis
@@ -6409,7 +6409,7 @@ def _nonnegative_coefficients(x):
         sage: _nonnegative_coefficients(x^2-4)
         False
     """
-    if is_Polynomial(x) or is_MPolynomial(x):
+    if isinstance(x, Polynomial) or isinstance(x, MPolynomial):
         return all(c >= 0 for c in x.coefficients(sparse=False))
     else:
         return x >= 0
diff --git a/src/sage/crypto/boolean_function.pyx b/src/sage/crypto/boolean_function.pyx
index d10c0c11546..13b8a1c82a1 100644
--- a/src/sage/crypto/boolean_function.pyx
+++ b/src/sage/crypto/boolean_function.pyx
@@ -40,7 +40,7 @@ from sage.rings.finite_rings.finite_field_constructor import GF
 from sage.rings.polynomial.pbori.pbori import BooleanPolynomial
 from sage.rings.finite_rings.finite_field_constructor import is_FiniteField
 from sage.rings.finite_rings.finite_field_givaro import FiniteField_givaro
-from sage.rings.polynomial.polynomial_element import is_Polynomial
+from sage.rings.polynomial.polynomial_element import Polynomial
 
 from sage.misc.superseded import deprecated_function_alias
 
@@ -327,7 +327,7 @@ cdef class BooleanFunction(SageObject):
             bitset_init(self._truth_table, <mp_bitcnt_t> (1<<self._nvariables))
             bitset_zero(self._truth_table)
 
-        elif is_Polynomial(x):
+        elif isinstance(x, Polynomial):
             K = x.base_ring()
             if is_FiniteField(K) and K.characteristic() == 2:
                 self._nvariables = K.degree()
diff --git a/src/sage/crypto/mq/rijndael_gf.py b/src/sage/crypto/mq/rijndael_gf.py
index 154c84d146c..fc3473d1b06 100644
--- a/src/sage/crypto/mq/rijndael_gf.py
+++ b/src/sage/crypto/mq/rijndael_gf.py
@@ -1549,9 +1549,9 @@ def compose(self, f, g, algorithm='encrypt', f_attr=None, g_attr=None):
         if not isinstance(f, RijndaelGF.Round_Component_Poly_Constr):
             msg = "keyword 'f' must be a Round_Component_Poly_Constr"
             raise TypeError(msg)
-        from sage.rings.polynomial.multi_polynomial import is_MPolynomial
+        from sage.rings.polynomial.multi_polynomial import MPolynomial
         if not isinstance(g, RijndaelGF.Round_Component_Poly_Constr) and \
-           not is_MPolynomial(g):
+           not isinstance(g, MPolynomial):
             msg = ("keyword 'g' must be a Round_Component_Poly_Constr or a "
                    "polynomial over {0}")
             raise TypeError(msg.format(self._F))
diff --git a/src/sage/crypto/sbox.pyx b/src/sage/crypto/sbox.pyx
index 77b7a904cc2..5c63e79bbac 100644
--- a/src/sage/crypto/sbox.pyx
+++ b/src/sage/crypto/sbox.pyx
@@ -173,12 +173,12 @@ cdef class SBox(SageObject):
             sage: S.output_size()
             3
         """
-        from sage.rings.polynomial.polynomial_element import is_Polynomial
+        from sage.rings.polynomial.polynomial_element import Polynomial
 
         if "S" in kwargs:
             args = kwargs["S"]
 
-        if len(args) == 1 and is_Polynomial(args[0]):
+        if len(args) == 1 and isinstance(args[0], Polynomial):
             # SBox defined via Univariate Polynomial, compute lookup table
             # by evaluating the polynomial on every base_ring element
             poly = args[0]
diff --git a/src/sage/crypto/stream.py b/src/sage/crypto/stream.py
index 728e73ca59b..90e4a6f5133 100644
--- a/src/sage/crypto/stream.py
+++ b/src/sage/crypto/stream.py
@@ -20,7 +20,7 @@
 from sage.arith.all import gcd, power_mod
 from sage.rings.finite_rings.finite_field_constructor import FiniteField
 from sage.rings.finite_rings.integer_mod_ring import IntegerModFactory
-from sage.rings.polynomial.polynomial_element import is_Polynomial
+from sage.rings.polynomial.polynomial_element import Polynomial
 
 
 IntegerModRing = IntegerModFactory("IntegerModRing")
@@ -73,7 +73,7 @@ def __call__(self, key):
         if not isinstance(key, (list, tuple)) and len(key) == 2:
             raise TypeError("Argument key (= %s) must be a list of tuple of length 2" % key)
         poly, IS = key
-        if not is_Polynomial(poly):
+        if not isinstance(poly, Polynomial):
             raise TypeError("poly (= %s) must be a polynomial." % poly)
         if not isinstance(IS, (list, tuple)):
             raise TypeError("IS (= %s) must be an initial in the key space." % IS)
diff --git a/src/sage/groups/affine_gps/group_element.py b/src/sage/groups/affine_gps/group_element.py
index cb5fd945239..127c1d33cb5 100644
--- a/src/sage/groups/affine_gps/group_element.py
+++ b/src/sage/groups/affine_gps/group_element.py
@@ -414,13 +414,13 @@ def __call__(self, v):
         if v in parent.vector_space():
             return self._A*v + self._b
 
-        from sage.rings.polynomial.polynomial_element import is_Polynomial
-        if is_Polynomial(v) and parent.degree() == 1:
+        from sage.rings.polynomial.polynomial_element import Polynomial
+        if isinstance(v, Polynomial) and parent.degree() == 1:
             ring = v.parent()
             return ring([self._A[0,0], self._b[0]])
 
-        from sage.rings.polynomial.multi_polynomial import is_MPolynomial
-        if is_MPolynomial(v) and parent.degree() == v.parent().ngens():
+        from sage.rings.polynomial.multi_polynomial import MPolynomial
+        if isinstance(v, MPolynomial) and parent.degree() == v.parent().ngens():
             ring = v.parent()
             from sage.modules.free_module_element import vector
             image_coords = self._A * vector(ring, ring.gens()) + self._b
diff --git a/src/sage/groups/perm_gps/permgroup_element.pyx b/src/sage/groups/perm_gps/permgroup_element.pyx
index c08442d60af..d81fe55f5e7 100644
--- a/src/sage/groups/perm_gps/permgroup_element.pyx
+++ b/src/sage/groups/perm_gps/permgroup_element.pyx
@@ -122,8 +122,8 @@ from cypari2.gen cimport Gen
 
 from sage.ext.stdsage cimport HAS_DICTIONARY
 from sage.rings.all      import ZZ, Integer
-from sage.rings.polynomial.polynomial_element import is_Polynomial
-from sage.rings.polynomial.multi_polynomial import is_MPolynomial
+from sage.rings.polynomial.polynomial_element import Polynomial
+from sage.rings.polynomial.multi_polynomial import MPolynomial
 from sage.structure.element import is_Matrix
 from sage.matrix.all     import MatrixSpace
 from sage.sets.finite_enumerated_set import FiniteEnumeratedSet
@@ -1233,12 +1233,12 @@ cdef class PermutationGroupElement(MultiplicativeGroupElement):
         """
         if not self_on_left:
             left = x
-            if is_Polynomial(left):
+            if isinstance(left, Polynomial):
                 if self != 1:
                     raise ValueError("%s does not act on %s" % (self,
                                                                 left.parent()))
                 return left
-            elif is_MPolynomial(left):
+            elif isinstance(left, MPolynomial):
                 R = left.parent()
                 vars = R.gens()
                 try:
diff --git a/src/sage/libs/symmetrica/sb.pxi b/src/sage/libs/symmetrica/sb.pxi
index 9bde08effa7..b884d33dafd 100644
--- a/src/sage/libs/symmetrica/sb.pxi
+++ b/src/sage/libs/symmetrica/sb.pxi
@@ -111,7 +111,7 @@ def t_POLYNOM_SCHUBERT_symmetrica(a):
 
     cdef OP ca = callocobject(), cres = callocobject()
 
-    if not is_MPolynomial(a):
+    if not isinstance(a, MPolynomial):
         freeall(ca)
         freeall(cres)
         raise TypeError("a (= %s) must be a multivariate polynomial")
diff --git a/src/sage/libs/symmetrica/symmetrica.pxi b/src/sage/libs/symmetrica/symmetrica.pxi
index 7244b7e9655..95f9e52fbda 100644
--- a/src/sage/libs/symmetrica/symmetrica.pxi
+++ b/src/sage/libs/symmetrica/symmetrica.pxi
@@ -401,7 +401,7 @@ cdef void late_import():
            SymmetricFunctions, \
            sqrt, \
            builtinlist, \
-           MPolynomialRing_base, is_MPolynomial,\
+           MPolynomialRing_base, MPolynomial,\
            SchubertPolynomialRing, SchubertPolynomial_class,\
            two, fifteen, thirty, zero, sage_maxint
 
@@ -453,8 +453,8 @@ cdef void late_import():
 
     import sage.rings.polynomial.multi_polynomial_ring
     MPolynomialRing_base = sage.rings.polynomial.multi_polynomial_ring.MPolynomialRing_base
-    import sage.rings.polynomial.multi_polynomial_element
-    is_MPolynomial = sage.rings.polynomial.multi_polynomial_element.is_MPolynomial
+    import sage.rings.polynomial.multi_polynomial
+    MPolynomial = sage.rings.polynomial.multi_polynomial.MPolynomial
 
     import sage.combinat.schubert_polynomial
     SchubertPolynomialRing = sage.combinat.schubert_polynomial.SchubertPolynomialRing
diff --git a/src/sage/manifolds/differentiable/characteristic_cohomology_class.py b/src/sage/manifolds/differentiable/characteristic_cohomology_class.py
index 689760625da..1af1c0c6d74 100644
--- a/src/sage/manifolds/differentiable/characteristic_cohomology_class.py
+++ b/src/sage/manifolds/differentiable/characteristic_cohomology_class.py
@@ -290,7 +290,7 @@
 from .bundle_connection import BundleConnection
 from .levi_civita_connection import LeviCivitaConnection
 from sage.symbolic.expression import Expression
-from sage.rings.polynomial.polynomial_element import is_Polynomial
+from sage.rings.polynomial.polynomial_element import Polynomial
 from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
 
 
@@ -795,7 +795,7 @@ def _element_constructor_(self, x, **kwargs):
             Characteristic cohomology class pontr(TM) of the Tangent bundle
              TM over the 8-dimensional differentiable manifold M
         """
-        if isinstance(x, (str, Expression)) or is_Polynomial(x):
+        if isinstance(x, (str, Expression)) or isinstance(x, Polynomial):
             return self._build_element(x, **kwargs)
 
         R = self.base_ring()
@@ -983,7 +983,7 @@ def _build_element(self, *args, **kwargs):
             val = P(val.taylor(x, 0, pow_range))
 
         # turn polynomial into a characteristic cohomology class via sequences
-        if is_Polynomial(val):
+        if isinstance(val, Polynomial):
             if class_type is None:
                 raise TypeError(f'class_type must be stated if {val} '
                                 f'is a polynomial')
diff --git a/src/sage/modular/modsym/ambient.py b/src/sage/modular/modsym/ambient.py
index 627f1c567e7..e7ca7c669e2 100644
--- a/src/sage/modular/modsym/ambient.py
+++ b/src/sage/modular/modsym/ambient.py
@@ -95,7 +95,7 @@ class ``ModularSymbolsAmbient``, derived from
 from sage.modules.free_module_element import FreeModuleElement
 from sage.rings.integer import Integer
 from sage.rings.integer_ring import ZZ
-from sage.rings.polynomial.multi_polynomial import is_MPolynomial
+from sage.rings.polynomial.multi_polynomial import MPolynomial
 from sage.rings.rational_field import QQ
 from sage.rings.ring import Ring
 from sage.structure.factorization import Factorization
@@ -474,7 +474,7 @@ def _element_constructor_(self, x, computed_with_hecke=False):
             return sum([c*self(y) for c, y in x], self(0))
 
         elif isinstance(x, list):
-            if len(x) == 3 and is_MPolynomial(x[0]):
+            if len(x) == 3 and isinstance(x[0], MPolynomial):
                 return self.modular_symbol_sum(x)
             else:
                 return self.modular_symbol(x)
diff --git a/src/sage/quadratic_forms/binary_qf.py b/src/sage/quadratic_forms/binary_qf.py
index e7cd7086f09..6eebe6f0c46 100755
--- a/src/sage/quadratic_forms/binary_qf.py
+++ b/src/sage/quadratic_forms/binary_qf.py
@@ -130,14 +130,14 @@ def __init__(self, a, b=None, c=None):
             sage: BinaryQF(0)
             0
         """
-        from sage.rings.polynomial.multi_polynomial_element import is_MPolynomial
+        from sage.rings.polynomial.multi_polynomial import MPolynomial
         if b is None and c is None:
             if (isinstance(a, (list, tuple))
                 and len(a) == 3):
                 a, b, c = a
             elif a == 0:
                 a = b = c = 0
-            elif (is_MPolynomial(a) and a.is_homogeneous() and a.base_ring() == ZZ
+            elif (isinstance(a, MPolynomial) and a.is_homogeneous() and a.base_ring() == ZZ
                   and a.degree() == 2 and a.parent().ngens() == 2):
                 x, y = a.parent().gens()
                 a, b, c = [a.monomial_coefficient(mon) for mon in [x**2, x*y, y**2]]
diff --git a/src/sage/quadratic_forms/constructions.py b/src/sage/quadratic_forms/constructions.py
index dd9f83216d3..89fca1bed6e 100644
--- a/src/sage/quadratic_forms/constructions.py
+++ b/src/sage/quadratic_forms/constructions.py
@@ -6,7 +6,7 @@
 ##
 
 from sage.rings.integer_ring import ZZ
-from sage.rings.polynomial.polynomial_element import is_Polynomial
+from sage.rings.polynomial.polynomial_element import Polynomial
 from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
 from sage.quadratic_forms.quadratic_form import QuadraticForm
 
@@ -46,7 +46,7 @@ def BezoutianQuadraticForm(f, g):
 
     """
     # Check that f and g are polynomials with a common base ring
-    if not is_Polynomial(f) or not is_Polynomial(g):
+    if not isinstance(f, Polynomial) or not isinstance(g, Polynomial):
         raise TypeError("one of your inputs is not a polynomial")
     if f.base_ring() != g.base_ring():   # TO DO:  Change this to allow coercion!
         raise TypeError("these polynomials are not defined over the same coefficient ring")
diff --git a/src/sage/quadratic_forms/quadratic_form.py b/src/sage/quadratic_forms/quadratic_form.py
index f53602f49f3..47f3d8be23e 100644
--- a/src/sage/quadratic_forms/quadratic_form.py
+++ b/src/sage/quadratic_forms/quadratic_form.py
@@ -32,8 +32,8 @@
 from sage.rings.ring import is_Ring, PrincipalIdealDomain
 from sage.structure.element import is_Vector
 from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
-from sage.rings.polynomial.polynomial_element import is_Polynomial
-from sage.rings.polynomial.multi_polynomial_element import is_MPolynomial
+from sage.rings.polynomial.polynomial_element import Polynomial
+from sage.rings.polynomial.multi_polynomial import MPolynomial
 from sage.modules.free_module_element import vector
 from sage.quadratic_forms.genera.genus import genera
 from sage.quadratic_forms.quadratic_form__evaluate import QFEvaluateVector, QFEvaluateMatrix
@@ -576,9 +576,9 @@ def __init__(self, R, n=None, entries=None, unsafe_initialization=False, number_
             M_ring = M.base_ring()
             matrix_init_flag = True
 
-        elif is_Polynomial(R) or is_MPolynomial(R):
+        elif isinstance(R, (Polynomial, MPolynomial)):
             p = R
-            
+
             if not p.is_zero() and not (p.is_homogeneous() and p.degree() == 2):
                 raise ValueError("polynomial is neither zero nor homogeneous of degree 2")
 
diff --git a/src/sage/rings/finite_rings/element_ntl_gf2e.pyx b/src/sage/rings/finite_rings/element_ntl_gf2e.pyx
index 109ce513b61..16065f47966 100644
--- a/src/sage/rings/finite_rings/element_ntl_gf2e.pyx
+++ b/src/sage/rings/finite_rings/element_ntl_gf2e.pyx
@@ -1004,8 +1004,8 @@ cdef class FiniteField_ntl_gf2eElement(FinitePolyExtElement):
             sage: e.polynomial()
             a^15 + a^13 + a^11 + a^10 + a^9 + a^8 + a^7 + a^6 + a^4 + a + 1
 
-            sage: from sage.rings.polynomial.polynomial_element import is_Polynomial
-            sage: is_Polynomial(e.polynomial())
+            sage: from sage.rings.polynomial.polynomial_element import Polynomial
+            sage: isinstance(e.polynomial(), Polynomial)
             True
 
             sage: e.polynomial('x')
diff --git a/src/sage/rings/finite_rings/finite_field_base.pyx b/src/sage/rings/finite_rings/finite_field_base.pyx
index 285412aa5bd..8a2d0ce94fa 100644
--- a/src/sage/rings/finite_rings/finite_field_base.pyx
+++ b/src/sage/rings/finite_rings/finite_field_base.pyx
@@ -1521,7 +1521,7 @@ cdef class FiniteField(Field):
             True
         """
         from .finite_field_constructor import GF
-        from sage.rings.polynomial.polynomial_element import is_Polynomial
+        from sage.rings.polynomial.polynomial_element import Polynomial
         from sage.rings.integer import Integer
         if name is None and names is not None:
             name = names
@@ -1532,7 +1532,7 @@ cdef class FiniteField(Field):
                 E = GF((self.characteristic(), modulus), name=name, **kwds)
             elif isinstance(modulus, (list, tuple)):
                 E = GF((self.characteristic(), len(modulus) - 1), name=name, modulus=modulus, **kwds)
-            elif is_Polynomial(modulus):
+            elif isinstance(modulus, Polynomial):
                 if modulus.change_ring(self).is_irreducible():
                     E = GF((self.characteristic(), modulus.degree()), name=name, modulus=modulus, **kwds)
                 else:
diff --git a/src/sage/rings/finite_rings/finite_field_constructor.py b/src/sage/rings/finite_rings/finite_field_constructor.py
index e42b2bed2dd..8b51fd37b2e 100644
--- a/src/sage/rings/finite_rings/finite_field_constructor.py
+++ b/src/sage/rings/finite_rings/finite_field_constructor.py
@@ -173,7 +173,7 @@
 
 from collections import defaultdict
 from sage.structure.category_object import normalize_names
-
+from sage.rings.polynomial.polynomial_element import Polynomial
 from sage.rings.integer import Integer
 
 # the import below is just a redirection
@@ -641,7 +641,7 @@ def create_key_and_extra_args(self, order, name=None, modulus=None, names=None,
                     else:
                         self._modulus_cache[order][modulus] = modulus = R.irreducible_element(n, algorithm=modulus)
                 else:
-                    if sage.rings.polynomial.polynomial_element.is_Polynomial(modulus):
+                    if isinstance(modulus, Polynomial):
                         modulus = modulus.change_variable_name('x')
                     modulus = R(modulus).monic()
 
diff --git a/src/sage/rings/finite_rings/finite_field_givaro.py b/src/sage/rings/finite_rings/finite_field_givaro.py
index b67df86bb30..0a0aa9a8d67 100644
--- a/src/sage/rings/finite_rings/finite_field_givaro.py
+++ b/src/sage/rings/finite_rings/finite_field_givaro.py
@@ -137,8 +137,8 @@ def __init__(self, q, name="a", modulus=None, repr="poly", cache=False):
         from .finite_field_constructor import GF
         FiniteField.__init__(self, GF(p), name, normalize=False)
 
-        from sage.rings.polynomial.polynomial_element import is_Polynomial
-        if not is_Polynomial(modulus):
+        from sage.rings.polynomial.polynomial_element import Polynomial
+        if not isinstance(modulus, Polynomial):
             raise TypeError("modulus must be a polynomial")
 
         self._cache = Cache_givaro(self, p, k, modulus, repr, cache)
diff --git a/src/sage/rings/finite_rings/finite_field_ntl_gf2e.py b/src/sage/rings/finite_rings/finite_field_ntl_gf2e.py
index 82004502674..e750f2fef5e 100644
--- a/src/sage/rings/finite_rings/finite_field_ntl_gf2e.py
+++ b/src/sage/rings/finite_rings/finite_field_ntl_gf2e.py
@@ -138,8 +138,8 @@ def __init__(self, q, names="a", modulus=None, repr="poly"):
             raise ValueError("q must be a 2-power")
         FiniteField.__init__(self, GF2, names, normalize=True)
 
-        from sage.rings.polynomial.polynomial_element import is_Polynomial
-        if not is_Polynomial(modulus):
+        from sage.rings.polynomial.polynomial_element import Polynomial
+        if not isinstance(modulus, Polynomial):
             raise TypeError("modulus must be a polynomial")
 
         self._cache = Cache_ntl_gf2e(self, k, modulus)
diff --git a/src/sage/rings/finite_rings/residue_field.pyx b/src/sage/rings/finite_rings/residue_field.pyx
index 83436bd6582..53ec78e057a 100644
--- a/src/sage/rings/finite_rings/residue_field.pyx
+++ b/src/sage/rings/finite_rings/residue_field.pyx
@@ -170,7 +170,7 @@ from sage.rings.fraction_field import is_FractionField
 
 from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
 from sage.rings.polynomial.polynomial_ring import is_PolynomialRing
-from sage.rings.polynomial.polynomial_element import is_Polynomial
+from sage.rings.polynomial.polynomial_element import Polynomial
 
 from sage.structure.factory import UniqueFactory
 from sage.structure.element cimport parent
@@ -300,7 +300,7 @@ class ResidueFieldFactory(UniqueFactory):
                         p = p.parent().ring_of_integers().ideal(p)
                     else:
                         p = p.parent().ideal(p)
-                elif is_Polynomial(p):
+                elif isinstance(p, Polynomial):
                     p = p.parent().ideal(p)
                 #elif isinstance(p.parent(), FractionField_1poly_field):
                 #    p = p.parent().ring_of_integers().ideal(p)
diff --git a/src/sage/rings/function_field/function_field.py b/src/sage/rings/function_field/function_field.py
index 782aa99cc42..60b834221ba 100644
--- a/src/sage/rings/function_field/function_field.py
+++ b/src/sage/rings/function_field/function_field.py
@@ -1307,8 +1307,8 @@ def __init__(self, polynomial, names, category=None):
             TypeError: unable to evaluate 'x' in Fraction Field of Univariate
             Polynomial Ring in t over Rational Field
         """
-        from sage.rings.polynomial.polynomial_element import is_Polynomial
-        if polynomial.parent().ngens()>1 or not is_Polynomial(polynomial):
+        from sage.rings.polynomial.polynomial_element import Polynomial
+        if polynomial.parent().ngens()>1 or not isinstance(polynomial, Polynomial):
             raise TypeError("polynomial must be univariate a polynomial")
         if names is None:
             names = (polynomial.variable_name(), )
diff --git a/src/sage/rings/laurent_series_ring.py b/src/sage/rings/laurent_series_ring.py
index 9cad038840c..be921f22154 100644
--- a/src/sage/rings/laurent_series_ring.py
+++ b/src/sage/rings/laurent_series_ring.py
@@ -466,8 +466,8 @@ def _element_constructor_(self, x, n=0, prec=infinity):
             x^-3
         """
         from sage.rings.fraction_field_element import is_FractionFieldElement
-        from sage.rings.polynomial.polynomial_element import is_Polynomial
-        from sage.rings.polynomial.multi_polynomial_element import is_MPolynomial
+        from sage.rings.polynomial.polynomial_element import Polynomial
+        from sage.rings.polynomial.multi_polynomial import MPolynomial
         from sage.structure.element import parent
         from sage.libs.pari.all import pari_gen
 
@@ -500,7 +500,7 @@ def _element_constructor_(self, x, n=0, prec=infinity):
                 return (self(self.polynomial_ring()(x)) << n).add_bigoh(prec)
         elif (is_FractionFieldElement(x)
               and (x.base_ring() is self.base_ring() or x.base_ring() == self.base_ring())
-              and (is_Polynomial(x.numerator()) or is_MPolynomial(x.numerator()))):
+              and isinstance(x.numerator(), (Polynomial, MPolynomial))):
             x = self(x.numerator()) / self(x.denominator())
             return (x << n).add_bigoh(prec)
         return self.element_class(self, x, n).add_bigoh(prec)
diff --git a/src/sage/rings/number_field/number_field.py b/src/sage/rings/number_field/number_field.py
index 0d2271844ca..18cf4686c8b 100644
--- a/src/sage/rings/number_field/number_field.py
+++ b/src/sage/rings/number_field/number_field.py
@@ -114,7 +114,7 @@
 import sage.interfaces.gap
 
 import sage.rings.complex_mpfr
-from sage.rings.polynomial.polynomial_element import is_Polynomial
+from sage.rings.polynomial.polynomial_element import Polynomial
 import sage.rings.real_mpfr
 import sage.rings.real_mpfi
 import sage.rings.complex_double
@@ -647,7 +647,7 @@ def create_key_and_extra_args(self, polynomial, name, check, embedding, latex_na
             raise TypeError("You must specify the name of the generator.")
         name = normalize_names(1, name)
 
-        if not is_Polynomial(polynomial):
+        if not isinstance(polynomial, Polynomial):
             try:
                 polynomial = polynomial.polynomial(QQ)
             except (AttributeError, TypeError):
@@ -868,7 +868,7 @@ def NumberFieldTower(polynomials, names, check=True, embeddings=None, latex_name
     f = polynomials[0]
     name = names[0]
     w = NumberFieldTower(polynomials[1:], names=names[1:], check=check, embeddings=embeddings[1:], latex_names=latex_names[1:], assume_disc_small=assume_disc_small, maximize_at_primes=maximize_at_primes, structures=structures[1:])
-    var = f.variable_name() if is_Polynomial(f) else 'x'
+    var = f.variable_name() if isinstance(f, Polynomial) else 'x'
 
     R = w[var]  # polynomial ring
     return w.extension(R(f), name, check=check, embedding=embeddings[0], structure=structures[0], latex_name=latex_names[0])  # currently, extension does not accept assume_disc_small, or maximize_at_primes
diff --git a/src/sage/rings/padics/factory.py b/src/sage/rings/padics/factory.py
index 3e396efa8ac..415aad91712 100644
--- a/src/sage/rings/padics/factory.py
+++ b/src/sage/rings/padics/factory.py
@@ -35,7 +35,7 @@
 from sage.rings.infinity import Infinity
 from sage.structure.factorization import Factorization
 from sage.rings.integer_ring import ZZ
-from sage.rings.polynomial.polynomial_element import is_Polynomial
+from sage.rings.polynomial.polynomial_element import Polynomial
 from sage.structure.element import is_Element
 from .padic_base_leaves import (pAdicRingCappedRelative,
                                 pAdicRingCappedAbsolute,
@@ -2535,7 +2535,7 @@ def Zq(q, prec = None, type = 'capped-rel', modulus = None, names=None,
         if isinstance(names, (list, tuple)):
             names = names[0]
         from sage.structure.element import Expression
-        if not (modulus is None or is_Polynomial(modulus) or isinstance(modulus, Expression)):
+        if not (modulus is None or isinstance(modulus, Polynomial) or isinstance(modulus, Expression)):
             raise TypeError("modulus must be a polynomial")
         if names is not None and not isinstance(names, str):
             names = str(names)
@@ -3278,7 +3278,7 @@ def create_key_and_extra_args(self, base, modulus, prec = None, print_mode = Non
                     raise ValueError("symbolic expression must be in only one variable")
                 exact_modulus = modulus.polynomial(base.exact_field())
                 approx_modulus = modulus.polynomial(base)
-            elif is_Polynomial(modulus):
+            elif isinstance(modulus, Polynomial):
                 if modulus.parent().ngens() != 1:
                     raise ValueError("must use univariate polynomial")
                 exact_modulus = modulus.change_ring(base.exact_field())
@@ -3381,7 +3381,7 @@ def create_object(self, version, key, approx_modulus=None, shift_seed=None):
             from sage.structure.element import Expression
             if isinstance(premodulus, Expression):
                 exact_modulus = premodulus.polynomial(base.exact_field())
-            elif is_Polynomial(premodulus):
+            elif isinstance(premodulus, Polynomial):
                 exact_modulus = premodulus.change_ring(base.exact_field())
             show_prec = None
         else:
diff --git a/src/sage/rings/padics/padic_template_element.pxi b/src/sage/rings/padics/padic_template_element.pxi
index 712f7dc9eeb..ce9e7029002 100644
--- a/src/sage/rings/padics/padic_template_element.pxi
+++ b/src/sage/rings/padics/padic_template_element.pxi
@@ -35,6 +35,7 @@ from sage.rings.infinity import infinity
 from sage.rings.rational import Rational
 from sage.rings.padics.precision_error import PrecisionError
 from sage.rings.padics.misc import trim_zeros
+from sage.rings.polynomial.polynomial_element import Polynomial
 from sage.structure.element import canonical_coercion
 import itertools
 
@@ -155,7 +156,7 @@ cdef class pAdicTemplateElement(pAdicGenericElement):
                 x = x + [k.prime_subfield().zero()] * (k.degree() - len(x))
         elif isinstance(x, (Integer, Rational, list, tuple)):
             pass
-        elif sage.rings.polynomial.polynomial_element.is_Polynomial(x) and x.variable_name() == self.parent().variable_name():
+        elif isinstance(x, Polynomial) and x.variable_name() == self.parent().variable_name():
             x = x.list()
         else:
             x = Rational(x)
diff --git a/src/sage/rings/polynomial/__init__.py b/src/sage/rings/polynomial/__init__.py
deleted file mode 100644
index e69de29bb2d..00000000000
diff --git a/src/sage/rings/polynomial/commutative_polynomial.pxd b/src/sage/rings/polynomial/commutative_polynomial.pxd
new file mode 100644
index 00000000000..c4a8956daa5
--- /dev/null
+++ b/src/sage/rings/polynomial/commutative_polynomial.pxd
@@ -0,0 +1,5 @@
+from sage.structure.element cimport CommutativeAlgebraElement
+
+
+cdef class CommutativePolynomial(CommutativeAlgebraElement):
+    pass
diff --git a/src/sage/rings/polynomial/commutative_polynomial.pyx b/src/sage/rings/polynomial/commutative_polynomial.pyx
new file mode 100644
index 00000000000..dc9f2cab8b7
--- /dev/null
+++ b/src/sage/rings/polynomial/commutative_polynomial.pyx
@@ -0,0 +1,23 @@
+cdef class CommutativePolynomial(CommutativeAlgebraElement):
+    r"""
+    Abstract base class for commutative polynomials in any number of variables
+
+    It is a common base for :class:`~sage.rings.polynomial.polynomial_element.Polynomial`,
+    :class:`~sage.rings.polynomial.multi_polynomial.MPolynomial`, and
+    :class:`~sage.rings.polynomial.infinite_polynomial_element.InfinitePolynomial`.
+
+    EXAMPLES::
+
+        sage: from sage.rings.polynomial.commutative_polynomial import CommutativePolynomial
+        sage: K.<x> = PolynomialRing(QQ)
+        sage: isinstance(x, CommutativePolynomial)
+        True
+        sage: K.<x,y> = PolynomialRing(QQ)
+        sage: isinstance(x, CommutativePolynomial)
+        True
+        sage: X.<x,y> = InfinitePolynomialRing(ZZ, implementation='sparse')
+        sage: isinstance(x[2], CommutativePolynomial)
+        True
+    """
+
+    pass
diff --git a/src/sage/rings/polynomial/infinite_polynomial_element.py b/src/sage/rings/polynomial/infinite_polynomial_element.py
index 5ac9bd00147..c7d43e31d1b 100644
--- a/src/sage/rings/polynomial/infinite_polynomial_element.py
+++ b/src/sage/rings/polynomial/infinite_polynomial_element.py
@@ -103,14 +103,16 @@
 
 from sage.rings.integer_ring import ZZ
 from sage.rings.integer import Integer
-from sage.structure.element import RingElement
 from sage.structure.richcmp import richcmp
 from sage.misc.cachefunc import cached_method
-from sage.rings.polynomial.multi_polynomial_element import is_MPolynomial
+from sage.misc.inherit_comparison import InheritComparisonClasscallMetaclass
+from sage.structure.element import RingElement
+from .commutative_polynomial import CommutativePolynomial
+from .multi_polynomial import MPolynomial
 import copy
 
 
-def InfinitePolynomial(A, p):
+class InfinitePolynomial(CommutativePolynomial, metaclass=InheritComparisonClasscallMetaclass):
     """
     Create an element of a Polynomial Ring with a Countably Infinite Number of Variables.
 
@@ -176,73 +178,72 @@ def InfinitePolynomial(A, p):
         alpha_2^2 + alpha_1^2
 
     """
-    from sage.structure.element import parent
-    if hasattr(A, '_P'):
-        if parent(p) is A._P or (A._P.base_ring().has_coerce_map_from(parent(p))):
-            return InfinitePolynomial_dense(A, p)
-        # MPolynomialRing_polydict is crab. So, in that case, use sage_eval
-        from sage.rings.polynomial.multi_polynomial_ring import MPolynomialRing_polydict
-        if isinstance(A._P, MPolynomialRing_polydict):
-            from sage.rings.polynomial.infinite_polynomial_ring import GenDictWithBasering
-            from sage.misc.sage_eval import sage_eval
-            p = sage_eval(repr(p), GenDictWithBasering(A._P, A._P.gens_dict()))
-            return InfinitePolynomial_dense(A, p)
-        else:
-            # Now there remains to fight the oddities and bugs of libsingular.
-            PP = p.parent()
-            if A._P.has_coerce_map_from(PP):
-                if A._P.ngens() == PP.ngens():  # coercion is sometimes by position!
-                    f = PP.hom(PP.variable_names(), A._P)
-                    try:
-                        return InfinitePolynomial_dense(A, f(p))
-                    except (ValueError, TypeError):
-                        # last desperate attempt: String conversion
-                        from sage.misc.sage_eval import sage_eval
-                        from sage.rings.polynomial.infinite_polynomial_ring import GenDictWithBasering
-                        # the base ring may be a function field, therefore
-                        # we need GenDictWithBasering
-                        return InfinitePolynomial_dense(A, sage_eval(repr(p), GenDictWithBasering(A._P, A._P.gens_dict())))
-                return InfinitePolynomial_dense(A, A._P(p))
-            # there is no coercion, so, we set up a name-preserving map.
-            SV = set(repr(x) for x in p.variables())
-            f = PP.hom([x if x in SV else 0 for x in PP.variable_names()], A._P)
-            try:
-                return InfinitePolynomial_dense(A, f(p))
-            except (ValueError, TypeError):
-                # last desperate attempt: String conversion
-                from sage.misc.sage_eval import sage_eval
-                from sage.rings.polynomial.infinite_polynomial_ring import GenDictWithBasering
-                # the base ring may be a function field, therefore
-                # we need GenDictWithBasering
-                return InfinitePolynomial_dense(A, sage_eval(repr(p), GenDictWithBasering(A._P, A._P.gens_dict())))
-    return InfinitePolynomial_sparse(A, p)
-
-
-class InfinitePolynomial_sparse(RingElement):
-    """
-    Element of a sparse Polynomial Ring with a Countably Infinite Number of Variables.
-
-    INPUT:
-
-    - ``A`` -- an Infinite Polynomial Ring in sparse implementation
-    - ``p`` -- a *classical* polynomial that can be interpreted in ``A``.
-
-    Of course, one should not directly invoke this class, but rather
-    construct elements of ``A`` in the usual way.
 
-    EXAMPLES::
+    @staticmethod
+    def __classcall_private__(cls, A, p):
+        r"""
+        TESTS::
 
-        sage: A.<a> = QQ[]
-        sage: B.<b,c> = InfinitePolynomialRing(A,implementation='sparse')
-        sage: p = a*b[100] + 1/2*c[4]
-        sage: p
-        a*b_100 + 1/2*c_4
-        sage: p.parent()
-        Infinite polynomial ring in b, c over Univariate Polynomial Ring in a over Rational Field
-        sage: p.polynomial().parent()
-        Multivariate Polynomial Ring in b_100, b_0, c_4, c_0 over Univariate Polynomial Ring in a over Rational Field
+            sage: from sage.rings.polynomial.infinite_polynomial_element import InfinitePolynomial
+            sage: X.<x,y> = InfinitePolynomialRing(ZZ, implementation='sparse')
+            sage: xy = (x[0] + y[0]).polynomial()
+            sage: xy.parent()
+            Multivariate Polynomial Ring in x_1, x_0, y_1, y_0 over Integer Ring
+            sage: sparse_xy = InfinitePolynomial(X, xy); sparse_xy
+            x_0 + y_0
+            sage: isinstance(sparse_xy, InfinitePolynomial)
+            True
+            sage: type(sparse_xy)
+            <class 'sage.rings.polynomial.infinite_polynomial_element.InfinitePolynomial_sparse'>
+            sage: X.<x,y> = InfinitePolynomialRing(ZZ, implementation='dense')
+            sage: dense_xy = InfinitePolynomial(X, xy); dense_xy
+            x_0 + y_0
+            sage: isinstance(dense_xy, InfinitePolynomial)
+            True
+            sage: type(dense_xy)
+            <class 'sage.rings.polynomial.infinite_polynomial_element.InfinitePolynomial_dense'>
+        """
+        from sage.structure.element import parent
+        if hasattr(A, '_P'):
+            if parent(p) is A._P or (A._P.base_ring().has_coerce_map_from(parent(p))):
+                return InfinitePolynomial_dense(A, p)
+            # MPolynomialRing_polydict is crab. So, in that case, use sage_eval
+            from sage.rings.polynomial.multi_polynomial_ring import MPolynomialRing_polydict
+            if isinstance(A._P, MPolynomialRing_polydict):
+                from sage.rings.polynomial.infinite_polynomial_ring import GenDictWithBasering
+                from sage.misc.sage_eval import sage_eval
+                p = sage_eval(repr(p), GenDictWithBasering(A._P, A._P.gens_dict()))
+                return InfinitePolynomial_dense(A, p)
+            else:
+                # Now there remains to fight the oddities and bugs of libsingular.
+                PP = p.parent()
+                if A._P.has_coerce_map_from(PP):
+                    if A._P.ngens() == PP.ngens():  # coercion is sometimes by position!
+                        f = PP.hom(PP.variable_names(), A._P)
+                        try:
+                            return InfinitePolynomial_dense(A, f(p))
+                        except (ValueError, TypeError):
+                            # last desperate attempt: String conversion
+                            from sage.misc.sage_eval import sage_eval
+                            from sage.rings.polynomial.infinite_polynomial_ring import GenDictWithBasering
+                            # the base ring may be a function field, therefore
+                            # we need GenDictWithBasering
+                            return InfinitePolynomial_dense(A, sage_eval(repr(p), GenDictWithBasering(A._P, A._P.gens_dict())))
+                    return InfinitePolynomial_dense(A, A._P(p))
+                # there is no coercion, so, we set up a name-preserving map.
+                SV = set(repr(x) for x in p.variables())
+                f = PP.hom([x if x in SV else 0 for x in PP.variable_names()], A._P)
+                try:
+                    return InfinitePolynomial_dense(A, f(p))
+                except (ValueError, TypeError):
+                    # last desperate attempt: String conversion
+                    from sage.misc.sage_eval import sage_eval
+                    from sage.rings.polynomial.infinite_polynomial_ring import GenDictWithBasering
+                    # the base ring may be a function field, therefore
+                    # we need GenDictWithBasering
+                    return InfinitePolynomial_dense(A, sage_eval(repr(p), GenDictWithBasering(A._P, A._P.gens_dict())))
+        return InfinitePolynomial_sparse(A, p)
 
-    """
     # Construction and other basic methods
     # We assume that p is good input. Type checking etc. is now done
     # in the _element_constructor_ of the parent.
@@ -261,7 +262,7 @@ def __init__(self, A, p):
         # the wrong ring and we get here without going through
         # _element_constructor_.  See trac 22514 for examples.
         # So a little extra checking is done here.
-        if not is_MPolynomial(p) or p.base_ring() is not A.base_ring():
+        if not isinstance(p, MPolynomial) or p.base_ring() is not A.base_ring():
             # coerce to a convenient multivariate polynomial ring
             p = A._minP(p)
 
@@ -312,60 +313,6 @@ def polynomial(self):
         """
         return self._p
 
-    def __call__(self, *args, **kwargs):
-        """
-        EXAMPLES::
-
-            sage: X.<x> = InfinitePolynomialRing(QQ,implementation='sparse')
-            sage: a = x[0] + x[1]
-            sage: a(x_0=2,x_1=x[1])
-            x_1 + 2
-            sage: _.parent()
-            Infinite polynomial ring in x over Rational Field
-            sage: a(x_1=3)
-            x_0 + 3
-            sage: _.parent()
-            Infinite polynomial ring in x over Rational Field
-            sage: a(x_1=x[100])
-            x_100 + x_0
-
-            sage: M = matrix([[1,1],[2,0]])
-            sage: a(x_1=M)
-            [x_0 + 1       1]
-            [      2     x_0]
-        """
-        # Replace any InfinitePolynomials by their underlying polynomials
-        if hasattr(self._p, 'variables'):
-            V = [str(x) for x in self._p.variables()]
-        else:
-            V = []
-        for kw in kwargs:
-            value = kwargs[kw]
-            if isinstance(value, InfinitePolynomial_sparse):
-                kwargs[kw] = value._p
-                V.append(kw)
-                if hasattr(value._p, 'variables'):
-                    V.extend([str(x) for x in value._p.variables()])
-        args = list(args)
-        for i, arg in enumerate(args):
-            if isinstance(arg, InfinitePolynomial_sparse):
-                args[i] = arg._p
-                if hasattr(arg._p, 'variables'):
-                    V.extend([str(x) for x in arg._p.variables()])
-        V = list(set(V))
-        V.sort(key=self.parent().varname_key, reverse=True)
-        if V:
-            from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
-            R = PolynomialRing(self._p.base_ring(), V, order=self.parent()._order)
-        else:
-            return self
-        res = R(self._p)(*args, **kwargs)
-        try:
-            from sage.misc.sage_eval import sage_eval
-            return sage_eval(repr(res), self.parent().gens_dict())
-        except Exception:
-            return res
-
     def _getAttributeNames(self):
         """
         This method implements tab completion, see :trac:`6854`.
@@ -649,79 +596,6 @@ def max_index(self):
         """
         return max([Integer(str(X).split('_')[1]) for X in self.variables()]+[-1])
 
-    # Basic arithmetics
-    def _add_(self, x):
-        """
-        EXAMPLES::
-
-            sage: X.<x> = InfinitePolynomialRing(QQ)
-            sage: x[1] + x[2] # indirect doctest
-            x_2 + x_1
-
-        Check adding from a different parent::
-
-            sage: Y.<x_0> = PolynomialRing(QQ)
-            sage: x[0] - x_0
-            0
-        """
-        # One may need a new parent for  self._p and x._p
-        try:
-            return InfinitePolynomial_sparse(self.parent(), self._p + x._p)
-        except Exception:
-            pass
-        # We can now assume that self._p and x._p actually are polynomials,
-        # hence, their parent is not simply the underlying ring.
-        VarList = list(set(self._p.parent().variable_names()).union(set(x._p.parent().variable_names())))
-        VarList.sort(key=self.parent().varname_key, reverse=True)
-        if VarList:
-            from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
-            R = PolynomialRing(self._p.base_ring(), VarList, order=self.parent()._order)
-        else:
-            R = self._p.base_ring()
-        return InfinitePolynomial_sparse(self.parent(), R(self._p) + R(x._p))
-
-    def _mul_(self, x):
-        """
-        EXAMPLES::
-
-            sage: X.<x> = InfinitePolynomialRing(ZZ)
-            sage: x[2]*x[1] # indirect doctest
-            x_2*x_1
-
-        """
-        try:
-            return InfinitePolynomial_sparse(self.parent(), self._p * x._p)
-        except Exception:
-            pass
-        # We can now assume that self._p and x._p actually are polynomials,
-        # hence, their parent is not just the underlying ring.
-        VarList = list(set(self._p.parent().variable_names()).union(set(x._p.parent().variable_names())))
-        VarList.sort(key=self.parent().varname_key, reverse=True)
-        if VarList:
-            from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
-            R = PolynomialRing(self._p.base_ring(), VarList, order=self.parent()._order)
-        else:
-            R = self._p.base_ring()
-        return InfinitePolynomial_sparse(self.parent(), R(self._p) * R(x._p))
-
-    def gcd(self, x):
-        """
-        computes the greatest common divisor
-
-        EXAMPLES::
-
-            sage: R.<x>=InfinitePolynomialRing(QQ)
-            sage: p1=x[0]+x[1]**2
-            sage: gcd(p1,p1+3)
-            1
-            sage: gcd(p1,p1)==p1
-            True
-        """
-        P = self.parent()
-        self._p = P._P(self._p)
-        x._p = P._P(x._p)
-        return InfinitePolynomial_sparse(self.parent(), self._p.gcd(x._p))
-
     def _rmul_(self, left):
         """
         TESTS::
@@ -731,7 +605,7 @@ def _rmul_(self, left):
             4
 
         """
-        return InfinitePolynomial_sparse(self.parent(), left * self._p)
+        return type(self)(self.parent(), left * self._p)
 
     def _lmul_(self, right):
         """
@@ -742,7 +616,7 @@ def _lmul_(self, right):
             4*alpha_3
 
         """
-        return InfinitePolynomial_sparse(self.parent(), self._p * right)
+        return type(self)(self.parent(), self._p * right)
 
     def _div_(self, x):
         r"""
@@ -818,235 +692,72 @@ def _floordiv_(self, x):
             R = self._p.base_ring()
         return InfinitePolynomial_sparse(self.parent(), R(self._p) // R(x._p))
 
-    def _sub_(self, x):
+    @cached_method
+    def lm(self):
         """
+        The leading monomial of ``self``.
+
         EXAMPLES::
 
-            sage: X.<x> = InfinitePolynomialRing(QQ)
-            sage: x[2] - x[1] # indirect doctest
-            x_2 - x_1
+            sage: X.<x,y> = InfinitePolynomialRing(QQ)
+            sage: p = 2*x[10]*y[30]+x[10]*y[1]^3*x[1]^2
+            sage: p.lm()
+            x_10*x_1^2*y_1^3
 
         """
-        try:
-            return InfinitePolynomial_sparse(self.parent(), self._p - x._p)
-        except Exception:
-            pass
-        # We can now assume that self._p and x._p actually are polynomials,
-        # hence, their parent is not just the underlying ring.
-        VarList = list(set(self._p.parent().variable_names()).union(x._p.parent().variable_names()))
-        VarList.sort(key=self.parent().varname_key, reverse=True)
-        if VarList:
-            from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
-            R = PolynomialRing(self._p.base_ring(), VarList, order=self.parent()._order)
-        else:
-            R = self._p.base_ring()
-        return InfinitePolynomial_sparse(self.parent(), R(self._p) - R(x._p))
+        if hasattr(self._p, 'lm'):
+            return InfinitePolynomial(self.parent(), self._p.lm())
+        if self._p == 0:
+            return self
+        if hasattr(self._p, 'variable_name'):  # if it is univariate
+            return InfinitePolynomial(self.parent(),
+                                      self._p.parent().gen() ** max(self._p.exponents()))
+        return self  # if it is scalar
 
-    def __pow__(self, n):
+    @cached_method
+    def lc(self):
         """
-        Exponentiation by an integer, or action by a callable object
+        The coefficient of the leading term of ``self``.
 
-        NOTE:
+        EXAMPLES::
 
-        The callable object must accept non-negative integers as input
-        and return non-negative integers. Typical use case is a
-        permutation, that will result in the corresponding permutation
-        of variables.
+            sage: X.<x,y> = InfinitePolynomialRing(QQ)
+            sage: p = 2*x[10]*y[30]+3*x[10]*y[1]^3*x[1]^2
+            sage: p.lc()
+            3
+
+        """
+        if hasattr(self._p, 'lc'):
+            return self._p.lc()
+        if hasattr(self._p, 'variable_name'):  # univariate case
+            return self._p.leading_coefficient()
+        # scalar case
+        return self._p
+
+    @cached_method
+    def lt(self):
+        """
+        The leading term (= product of coefficient and monomial) of ``self``.
 
         EXAMPLES::
 
-            sage: X.<x,y> = InfinitePolynomialRing(QQ,implementation='sparse')
-            sage: p = x[10]*y[2]+2*x[1]*y[3]
-            sage: P = Permutation(((1,2),(3,4,5)))
-            sage: p^P # indirect doctest
-            x_10*y_1 + 2*x_2*y_4
+            sage: X.<x,y> = InfinitePolynomialRing(QQ)
+            sage: p = 2*x[10]*y[30]+3*x[10]*y[1]^3*x[1]^2
+            sage: p.lt()
+            3*x_10*x_1^2*y_1^3
 
         """
-        P = self.parent()
-        if callable(n):
-            if (self._p.parent() == self._p.base_ring()):
-                return self
-            if not (hasattr(self._p, 'variables') and self._p.variables()):
-                return self
-            if hasattr(n, 'to_cycles') and hasattr(n, '__len__'):  # duck typing Permutation
-                # auxiliary function, necessary since n(m) raises an error if m>len(n)
-                l = len(n)
+        if hasattr(self._p, 'lt'):
+            return InfinitePolynomial(self.parent(), self._p.lt())
+        if self._p == 0:
+            return self
+        if hasattr(self._p, 'variable_name'):  # if it is univariate
+            return InfinitePolynomial(self.parent(), self._p.leading_coefficient()*self._p.parent().gen()**max(self._p.exponents()))
+        return self  # if it is scalar
 
-                def p(m):
-                    return n(m) if 0 < m <= l else m
-            else:  # Permutation group element
-                p = n
-
-            def q(s):
-                return s[0] + '_' + str(p(ZZ(s[1])))
-
-            newVars = [q(X.split('_')) for X in self._p.parent().variable_names()]
-            if not newVars:
-                return self
-            copyVars = copy.copy(newVars)
-            newVars = list(set(list(self._p.parent().variable_names())+newVars))
-            newVars.sort(key=self.parent().varname_key, reverse=True)
-            if newVars == list(self._p.parent().variable_names()):
-                newR = self._p.parent()
-            else:
-                from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
-                newR = PolynomialRing(self._p.base_ring(), newVars, order=P._order)
-            mapR = self._p.parent().hom(copyVars, newR)
-            return InfinitePolynomial_sparse(self.parent(), mapR(self._p))
-        return InfinitePolynomial_sparse(self.parent(), self._p**n)
-
-    # Basic tools for Buchberger algorithm:
-    # order, leading term/monomial, symmetric cancellation order
-
-    def _richcmp_(self, x, op):
-        r"""
-        Comparison of Infinite Polynomials.
-
-        NOTE:
-
-        Let x and y be generators of the parent of self. We only consider
-        monomial orderings in which x[m] > y[n] iff x appears earlier in the
-        list of generators than y, or x==y and m>n
-
-        Under this restriction, the monomial ordering can be 'lex' (default),
-        'degrevlex' or 'deglex'.
-
-        EXAMPLES::
-
-            sage: X.<x,y> = InfinitePolynomialRing(QQ, implementation='sparse')
-            sage: a = x[10]^3
-            sage: b = x[1] + x[2]
-            sage: c = x[1] + x[2]
-            sage: d = y[1] + x[2]
-            sage: a == a # indirect doctest
-            True
-            sage: b == c # indirect doctest
-            True
-            sage: a == b # indirect doctest
-            False
-            sage: c > d # indirect doctest
-            True
-
-        TESTS:
-
-        A classical and an infinite sparse polynomial ring. Note that
-        the Sage coercion system allows comparison only if a common
-        parent for the two rings can be constructed. This is why we
-        have to have the order 'degrevlex'::
-
-            sage: X.<x,y> = InfinitePolynomialRing(ZZ,order='degrevlex', implementation='sparse')
-            sage: Y.<z,x_3,x_1> = QQ[]
-            sage: x[3] == x_3 # indirect doctest
-            True
-
-        Two infinite polynomial rings in different implementation and
-        order::
-
-            sage: Y = InfinitePolynomialRing(QQ,['x','y'],order='deglex',implementation='dense')
-            sage: x[2] == Y(x[2]) # indirect doctest
-            True
-
-        An example in which a previous version had failed::
-
-            sage: X.<x,y> = InfinitePolynomialRing(GF(3), order='degrevlex', implementation='sparse')
-            sage: p = Y('x_3*x_0^2 + x_0*y_3*y_0')
-            sage: q = Y('x_1*x_0^2 + x_0*y_1*y_0')
-            sage: p < q # indirect doctest
-            False
-
-        """
-        # We can assume that self.parent() is x.parent(),
-        # but of course the underlying polynomial rings
-        # may be widely different, and the sage coercion
-        # system can't guess what order we want.
-        from sage.structure.element import parent
-        R1 = parent(self._p)
-        R2 = parent(x._p)
-        if (hasattr(R1, 'has_coerce_map_from') and R1.has_coerce_map_from(R2)) or (hasattr(R2, 'has_coerce_map_from') and R2.has_coerce_map_from(R1)):
-            return richcmp(self._p, x._p, op)
-        VarList = list(set(self._p.parent().variable_names()).union(x._p.parent().variable_names()))
-        VarList.sort(key=self.parent().varname_key, reverse=True)
-        if VarList:
-            from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
-            R = PolynomialRing(self._p.base_ring(), VarList, order=self.parent()._order)
-        else:
-            R = self._p.base_ring()
-        if (self._p.parent() is self._p.base_ring()) or not self._p.parent().gens():
-            fself = self._p.base_ring()
-        else:
-            fself = self._p.parent().hom(self._p.parent().variable_names(), R)
-        if (x._p.parent() is x._p.base_ring()) or not x._p.parent().gens():
-            fx = x._p.base_ring()
-        else:
-            fx = x._p.parent().hom(x._p.parent().variable_names(), R)
-        return richcmp(fself(self._p), fx(x._p), op)
-
-    @cached_method
-    def lm(self):
-        """
-        The leading monomial of ``self``.
-
-        EXAMPLES::
-
-            sage: X.<x,y> = InfinitePolynomialRing(QQ)
-            sage: p = 2*x[10]*y[30]+x[10]*y[1]^3*x[1]^2
-            sage: p.lm()
-            x_10*x_1^2*y_1^3
-
-        """
-        if hasattr(self._p, 'lm'):
-            return InfinitePolynomial(self.parent(), self._p.lm())
-        if self._p == 0:
-            return self
-        if hasattr(self._p, 'variable_name'):  # if it is univariate
-            return InfinitePolynomial(self.parent(),
-                                      self._p.parent().gen() ** max(self._p.exponents()))
-        return self  # if it is scalar
-
-    @cached_method
-    def lc(self):
-        """
-        The coefficient of the leading term of ``self``.
-
-        EXAMPLES::
-
-            sage: X.<x,y> = InfinitePolynomialRing(QQ)
-            sage: p = 2*x[10]*y[30]+3*x[10]*y[1]^3*x[1]^2
-            sage: p.lc()
-            3
-
-        """
-        if hasattr(self._p, 'lc'):
-            return self._p.lc()
-        if hasattr(self._p, 'variable_name'):  # univariate case
-            return self._p.leading_coefficient()
-        # scalar case
-        return self._p
-
-    @cached_method
-    def lt(self):
-        """
-        The leading term (= product of coefficient and monomial) of ``self``.
-
-        EXAMPLES::
-
-            sage: X.<x,y> = InfinitePolynomialRing(QQ)
-            sage: p = 2*x[10]*y[30]+3*x[10]*y[1]^3*x[1]^2
-            sage: p.lt()
-            3*x_10*x_1^2*y_1^3
-
-        """
-        if hasattr(self._p, 'lt'):
-            return InfinitePolynomial(self.parent(), self._p.lt())
-        if self._p == 0:
-            return self
-        if hasattr(self._p, 'variable_name'):  # if it is univariate
-            return InfinitePolynomial(self.parent(), self._p.leading_coefficient()*self._p.parent().gen()**max(self._p.exponents()))
-        return self  # if it is scalar
-
-    def tail(self):
-        """
-        The tail of ``self`` (this is ``self`` minus its leading term).
+    def tail(self):
+        """
+        The tail of ``self`` (this is ``self`` minus its leading term).
 
         EXAMPLES::
 
@@ -1126,7 +837,7 @@ def footprint(self):
             l = len(self.parent()._names)
             # get the pairs (shift,exponent) of the leading monomial, indexed by the variable names
             Vars = self._p.parent().variable_names()
-            from sage.rings.polynomial.multi_polynomial_libsingular import MPolynomial_libsingular
+            from sage.rings.polynomial.multi_polynomial import MPolynomial_libsingular
             if isinstance(self._p, MPolynomial_libsingular):
                 L = [(Vars[i].split('_'), e) for i, e in enumerate(self._p.lm().exponents(as_ETuples=False)[0]) if e]
             elif hasattr(self._p, 'lm'):
@@ -1482,8 +1193,325 @@ def __iter__(self):
                      self.__class__(self.parent(), monomial))
                     for coefficient, monomial in self._p)
 
+    def gcd(self, x):
+        """
+        computes the greatest common divisor
+
+        EXAMPLES::
+
+            sage: R.<x>=InfinitePolynomialRing(QQ)
+            sage: p1=x[0]+x[1]**2
+            sage: gcd(p1,p1+3)
+            1
+            sage: gcd(p1,p1)==p1
+            True
+        """
+        P = self.parent()
+        self._p = P._P(self._p)
+        x._p = P._P(x._p)
+        return self.__class__.__base__(self.parent(), self._p.gcd(x._p))
+
+
+class InfinitePolynomial_sparse(InfinitePolynomial):
+    """
+    Element of a sparse Polynomial Ring with a Countably Infinite Number of Variables.
+
+    INPUT:
+
+    - ``A`` -- an Infinite Polynomial Ring in sparse implementation
+    - ``p`` -- a *classical* polynomial that can be interpreted in ``A``.
+
+    Of course, one should not directly invoke this class, but rather
+    construct elements of ``A`` in the usual way.
+
+    EXAMPLES::
+
+        sage: A.<a> = QQ[]
+        sage: B.<b,c> = InfinitePolynomialRing(A,implementation='sparse')
+        sage: p = a*b[100] + 1/2*c[4]
+        sage: p
+        a*b_100 + 1/2*c_4
+        sage: p.parent()
+        Infinite polynomial ring in b, c over Univariate Polynomial Ring in a over Rational Field
+        sage: p.polynomial().parent()
+        Multivariate Polynomial Ring in b_100, b_0, c_4, c_0 over Univariate Polynomial Ring in a over Rational Field
 
-class InfinitePolynomial_dense(InfinitePolynomial_sparse):
+    """
+
+    def __call__(self, *args, **kwargs):
+        """
+        EXAMPLES::
+
+            sage: X.<x> = InfinitePolynomialRing(QQ,implementation='sparse')
+            sage: a = x[0] + x[1]
+            sage: a(x_0=2,x_1=x[1])
+            x_1 + 2
+            sage: _.parent()
+            Infinite polynomial ring in x over Rational Field
+            sage: a(x_1=3)
+            x_0 + 3
+            sage: _.parent()
+            Infinite polynomial ring in x over Rational Field
+            sage: a(x_1=x[100])
+            x_100 + x_0
+
+            sage: M = matrix([[1,1],[2,0]])
+            sage: a(x_1=M)
+            [x_0 + 1       1]
+            [      2     x_0]
+        """
+        # Replace any InfinitePolynomials by their underlying polynomials
+        if hasattr(self._p, 'variables'):
+            V = [str(x) for x in self._p.variables()]
+        else:
+            V = []
+        for kw in kwargs:
+            value = kwargs[kw]
+            if isinstance(value, InfinitePolynomial):
+                kwargs[kw] = value._p
+                V.append(kw)
+                if hasattr(value._p, 'variables'):
+                    V.extend([str(x) for x in value._p.variables()])
+        args = list(args)
+        for i, arg in enumerate(args):
+            if isinstance(arg, InfinitePolynomial):
+                args[i] = arg._p
+                if hasattr(arg._p, 'variables'):
+                    V.extend([str(x) for x in arg._p.variables()])
+        V = list(set(V))
+        V.sort(key=self.parent().varname_key, reverse=True)
+        if V:
+            from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
+            R = PolynomialRing(self._p.base_ring(), V, order=self.parent()._order)
+        else:
+            return self
+        res = R(self._p)(*args, **kwargs)
+        try:
+            from sage.misc.sage_eval import sage_eval
+            return sage_eval(repr(res), self.parent().gens_dict())
+        except Exception:
+            return res
+
+    # Basic arithmetics
+    def _add_(self, x):
+        """
+        EXAMPLES::
+
+            sage: X.<x> = InfinitePolynomialRing(QQ)
+            sage: x[1] + x[2] # indirect doctest
+            x_2 + x_1
+
+        Check adding from a different parent::
+
+            sage: Y.<x_0> = PolynomialRing(QQ)
+            sage: x[0] - x_0
+            0
+        """
+        # One may need a new parent for  self._p and x._p
+        try:
+            return InfinitePolynomial_sparse(self.parent(), self._p + x._p)
+        except Exception:
+            pass
+        # We can now assume that self._p and x._p actually are polynomials,
+        # hence, their parent is not simply the underlying ring.
+        VarList = list(set(self._p.parent().variable_names()).union(set(x._p.parent().variable_names())))
+        VarList.sort(key=self.parent().varname_key, reverse=True)
+        if VarList:
+            from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
+            R = PolynomialRing(self._p.base_ring(), VarList, order=self.parent()._order)
+        else:
+            R = self._p.base_ring()
+        return InfinitePolynomial_sparse(self.parent(), R(self._p) + R(x._p))
+
+    def _mul_(self, x):
+        """
+        EXAMPLES::
+
+            sage: X.<x> = InfinitePolynomialRing(ZZ)
+            sage: x[2]*x[1] # indirect doctest
+            x_2*x_1
+
+        """
+        try:
+            return InfinitePolynomial_sparse(self.parent(), self._p * x._p)
+        except Exception:
+            pass
+        # We can now assume that self._p and x._p actually are polynomials,
+        # hence, their parent is not just the underlying ring.
+        VarList = list(set(self._p.parent().variable_names()).union(set(x._p.parent().variable_names())))
+        VarList.sort(key=self.parent().varname_key, reverse=True)
+        if VarList:
+            from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
+            R = PolynomialRing(self._p.base_ring(), VarList, order=self.parent()._order)
+        else:
+            R = self._p.base_ring()
+        return InfinitePolynomial_sparse(self.parent(), R(self._p) * R(x._p))
+
+    def _sub_(self, x):
+        """
+        EXAMPLES::
+
+            sage: X.<x> = InfinitePolynomialRing(QQ)
+            sage: x[2] - x[1] # indirect doctest
+            x_2 - x_1
+
+        """
+        try:
+            return InfinitePolynomial_sparse(self.parent(), self._p - x._p)
+        except Exception:
+            pass
+        # We can now assume that self._p and x._p actually are polynomials,
+        # hence, their parent is not just the underlying ring.
+        VarList = list(set(self._p.parent().variable_names()).union(x._p.parent().variable_names()))
+        VarList.sort(key=self.parent().varname_key, reverse=True)
+        if VarList:
+            from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
+            R = PolynomialRing(self._p.base_ring(), VarList, order=self.parent()._order)
+        else:
+            R = self._p.base_ring()
+        return InfinitePolynomial_sparse(self.parent(), R(self._p) - R(x._p))
+
+    def __pow__(self, n):
+        """
+        Exponentiation by an integer, or action by a callable object
+
+        NOTE:
+
+        The callable object must accept non-negative integers as input
+        and return non-negative integers. Typical use case is a
+        permutation, that will result in the corresponding permutation
+        of variables.
+
+        EXAMPLES::
+
+            sage: X.<x,y> = InfinitePolynomialRing(QQ,implementation='sparse')
+            sage: p = x[10]*y[2]+2*x[1]*y[3]
+            sage: P = Permutation(((1,2),(3,4,5)))
+            sage: p^P # indirect doctest
+            x_10*y_1 + 2*x_2*y_4
+
+        """
+        P = self.parent()
+        if callable(n):
+            if (self._p.parent() == self._p.base_ring()):
+                return self
+            if not (hasattr(self._p, 'variables') and self._p.variables()):
+                return self
+            if hasattr(n, 'to_cycles') and hasattr(n, '__len__'):  # duck typing Permutation
+                # auxiliary function, necessary since n(m) raises an error if m>len(n)
+                l = len(n)
+
+                def p(m):
+                    return n(m) if 0 < m <= l else m
+            else:  # Permutation group element
+                p = n
+
+            def q(s):
+                return s[0] + '_' + str(p(ZZ(s[1])))
+
+            newVars = [q(X.split('_')) for X in self._p.parent().variable_names()]
+            if not newVars:
+                return self
+            copyVars = copy.copy(newVars)
+            newVars = list(set(list(self._p.parent().variable_names())+newVars))
+            newVars.sort(key=self.parent().varname_key, reverse=True)
+            if newVars == list(self._p.parent().variable_names()):
+                newR = self._p.parent()
+            else:
+                from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
+                newR = PolynomialRing(self._p.base_ring(), newVars, order=P._order)
+            mapR = self._p.parent().hom(copyVars, newR)
+            return InfinitePolynomial_sparse(self.parent(), mapR(self._p))
+        return InfinitePolynomial_sparse(self.parent(), self._p**n)
+
+    # Basic tools for Buchberger algorithm:
+    # order, leading term/monomial, symmetric cancellation order
+
+    def _richcmp_(self, x, op):
+        r"""
+        Comparison of Infinite Polynomials.
+
+        NOTE:
+
+        Let x and y be generators of the parent of self. We only consider
+        monomial orderings in which x[m] > y[n] iff x appears earlier in the
+        list of generators than y, or x==y and m>n
+
+        Under this restriction, the monomial ordering can be 'lex' (default),
+        'degrevlex' or 'deglex'.
+
+        EXAMPLES::
+
+            sage: X.<x,y> = InfinitePolynomialRing(QQ, implementation='sparse')
+            sage: a = x[10]^3
+            sage: b = x[1] + x[2]
+            sage: c = x[1] + x[2]
+            sage: d = y[1] + x[2]
+            sage: a == a # indirect doctest
+            True
+            sage: b == c # indirect doctest
+            True
+            sage: a == b # indirect doctest
+            False
+            sage: c > d # indirect doctest
+            True
+
+        TESTS:
+
+        A classical and an infinite sparse polynomial ring. Note that
+        the Sage coercion system allows comparison only if a common
+        parent for the two rings can be constructed. This is why we
+        have to have the order 'degrevlex'::
+
+            sage: X.<x,y> = InfinitePolynomialRing(ZZ,order='degrevlex', implementation='sparse')
+            sage: Y.<z,x_3,x_1> = QQ[]
+            sage: x[3] == x_3 # indirect doctest
+            True
+
+        Two infinite polynomial rings in different implementation and
+        order::
+
+            sage: Y = InfinitePolynomialRing(QQ,['x','y'],order='deglex',implementation='dense')
+            sage: x[2] == Y(x[2]) # indirect doctest
+            True
+
+        An example in which a previous version had failed::
+
+            sage: X.<x,y> = InfinitePolynomialRing(GF(3), order='degrevlex', implementation='sparse')
+            sage: p = Y('x_3*x_0^2 + x_0*y_3*y_0')
+            sage: q = Y('x_1*x_0^2 + x_0*y_1*y_0')
+            sage: p < q # indirect doctest
+            False
+
+        """
+        # We can assume that self.parent() is x.parent(),
+        # but of course the underlying polynomial rings
+        # may be widely different, and the sage coercion
+        # system can't guess what order we want.
+        from sage.structure.element import parent
+        R1 = parent(self._p)
+        R2 = parent(x._p)
+        if (hasattr(R1, 'has_coerce_map_from') and R1.has_coerce_map_from(R2)) or (hasattr(R2, 'has_coerce_map_from') and R2.has_coerce_map_from(R1)):
+            return richcmp(self._p, x._p, op)
+        VarList = list(set(self._p.parent().variable_names()).union(x._p.parent().variable_names()))
+        VarList.sort(key=self.parent().varname_key, reverse=True)
+        if VarList:
+            from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
+            R = PolynomialRing(self._p.base_ring(), VarList, order=self.parent()._order)
+        else:
+            R = self._p.base_ring()
+        if (self._p.parent() is self._p.base_ring()) or not self._p.parent().gens():
+            fself = self._p.base_ring()
+        else:
+            fself = self._p.parent().hom(self._p.parent().variable_names(), R)
+        if (x._p.parent() is x._p.base_ring()) or not x._p.parent().gens():
+            fx = x._p.base_ring()
+        else:
+            fx = x._p.parent().hom(x._p.parent().variable_names(), R)
+        return richcmp(fself(self._p), fx(x._p), op)
+
+
+class InfinitePolynomial_dense(InfinitePolynomial):
     """
     Element of a dense Polynomial Ring with a Countably Infinite Number of Variables.
 
@@ -1499,8 +1527,6 @@ class InfinitePolynomial_dense(InfinitePolynomial_sparse):
     :class:`~sage.rings.polynomial.infinite_polynomial_element.InfinitePolynomial_sparse`. See
     there for a description of the methods.
     """
-    # Construction and other basic methods
-##    def __init__(self, A, p): # is inherited from the dense implementation
 
     def __call__(self, *args, **kwargs):
         """
@@ -1524,11 +1550,11 @@ def __call__(self, *args, **kwargs):
         # Replace any InfinitePolynomials by their underlying polynomials
         for kw in kwargs:
             value = kwargs[kw]
-            if isinstance(value, InfinitePolynomial_sparse):
+            if isinstance(value, InfinitePolynomial):
                 kwargs[kw] = value._p
         args = list(args)
         for i, arg in enumerate(args):
-            if isinstance(arg, InfinitePolynomial_sparse):
+            if isinstance(arg, InfinitePolynomial):
                 args[i] = arg._p
         self._p = self.parent().polynomial_ring()(self._p)
         res = self._p(*args, **kwargs)
@@ -1608,28 +1634,6 @@ def _mul_(self, x):
         x._p = P._P(x._p)
         return InfinitePolynomial_dense(self.parent(), self._p * x._p)
 
-    def _rmul_(self, left):
-        """
-        TESTS::
-
-            sage: R.<alpha,beta> = InfinitePolynomialRing(QQ)
-            sage: R.from_base_ring(4)   # indirect doctest
-            4
-
-        """
-        return InfinitePolynomial_dense(self.parent(), left*self._p)
-
-    def _lmul_(self, right):
-        """
-        TESTS::
-
-            sage: R.<alpha,beta> = InfinitePolynomialRing(QQ)
-            sage: alpha[3]*4   # indirect doctest
-            4*alpha_3
-
-        """
-        return InfinitePolynomial_dense(self.parent(), self._p*right)
-
     def _sub_(self, x):
         """
         EXAMPLES::
diff --git a/src/sage/rings/polynomial/infinite_polynomial_ring.py b/src/sage/rings/polynomial/infinite_polynomial_ring.py
index ff237743f4a..c14c142be53 100644
--- a/src/sage/rings/polynomial/infinite_polynomial_ring.py
+++ b/src/sage/rings/polynomial/infinite_polynomial_ring.py
@@ -936,7 +936,7 @@ def _element_constructor_(self, x):
                 raise ValueError("cannot convert %s into an element of %s" % (x, self))
 
         # direct conversion will only be used if the underlying polynomials are libsingular.
-        from sage.rings.polynomial.multi_polynomial_libsingular import MPolynomial_libsingular, MPolynomialRing_libsingular
+        from sage.rings.polynomial.multi_polynomial import MPolynomial_libsingular
         # try interpretation in self._P, if we have a dense implementation
         if hasattr(self, '_P'):
             if x.parent() is self._P:
@@ -945,36 +945,38 @@ def _element_constructor_(self, x):
             # that MPolynomialRing_polydict does not work in complicated settings.
             # So, if self._P is libsingular (and this will be the case in many
             # applications!), we do it "nicely". Otherwise, we have to use sage_eval.
-            if isinstance(x, MPolynomial_libsingular) and isinstance(self._P, MPolynomialRing_libsingular):
-                if xmaxind == -1:  # Otherwise, x has been an InfinitePolynomial
-                    # We infer the correct variable shift.
-                    # Note: Since we are in the "libsingular" case, there are
-                    # no further "variables" hidden in the base ring of x.parent()
+            if isinstance(x, MPolynomial_libsingular):
+                from sage.rings.polynomial.multi_polynomial_libsingular import MPolynomialRing_libsingular
+                if isinstance(self._P, MPolynomialRing_libsingular):
+                    if xmaxind == -1:  # Otherwise, x has been an InfinitePolynomial
+                        # We infer the correct variable shift.
+                        # Note: Since we are in the "libsingular" case, there are
+                        # no further "variables" hidden in the base ring of x.parent()
+                        try:
+                            VarList = [repr(v) for v in x.variables()]
+                            # since interpretation in base ring
+                            # was impossible, it *must* have
+                            # variables
+                            # This tests admissibility on the fly:
+                            VarList.sort(key=self.varname_key, reverse=True)
+                        except ValueError:
+                            raise ValueError("cannot convert %s into an element of %s - variables are not admissible" % (x, self))
+                        xmaxind = max([int(v.split('_')[1]) for v in VarList])
                     try:
-                        VarList = [repr(v) for v in x.variables()]
-                        # since interpretation in base ring
-                        # was impossible, it *must* have
-                        # variables
-                        # This tests admissibility on the fly:
-                        VarList.sort(key=self.varname_key, reverse=True)
-                    except ValueError:
-                        raise ValueError("cannot convert %s into an element of %s - variables are not admissible" % (x, self))
-                    xmaxind = max([int(v.split('_')[1]) for v in VarList])
-                try:
-                    # Apparently, in libsingular, the polynomial conversion is not done by
-                    # name but by position, if the number of variables in the parents coincide.
-                    # So, we shift self._P to achieve xmaxind, and if the number of variables is
-                    # the same then we shift further. We then *must* be
-                    # able to convert x into self._P, or conversion to self is
-                    # impossible (and will be done in InfinitePolynomial(...)
-                    if self._max < xmaxind:
-                        self.gen()[xmaxind]
-                    if self._P.ngens() == x.parent().ngens():
-                        self.gen()[self._max + 1]
-                    # conversion to self._P will be done in InfinitePolynomial.__init__
-                    return InfinitePolynomial(self, x)
-                except (ValueError, TypeError, NameError):
-                    raise ValueError("cannot convert %s (from %s, but variables %s) into an element of %s - no conversion into underlying polynomial ring %s" % (x, x.parent(), x.variables(), self, self._P))
+                        # Apparently, in libsingular, the polynomial conversion is not done by
+                        # name but by position, if the number of variables in the parents coincide.
+                        # So, we shift self._P to achieve xmaxind, and if the number of variables is
+                        # the same then we shift further. We then *must* be
+                        # able to convert x into self._P, or conversion to self is
+                        # impossible (and will be done in InfinitePolynomial(...)
+                        if self._max < xmaxind:
+                            self.gen()[xmaxind]
+                        if self._P.ngens() == x.parent().ngens():
+                            self.gen()[self._max + 1]
+                        # conversion to self._P will be done in InfinitePolynomial.__init__
+                        return InfinitePolynomial(self, x)
+                    except (ValueError, TypeError, NameError):
+                        raise ValueError("cannot convert %s (from %s, but variables %s) into an element of %s - no conversion into underlying polynomial ring %s" % (x, x.parent(), x.variables(), self, self._P))
             # By now, x or self._P are not libsingular. Since MPolynomialRing_polydict
             # is too buggy, we use string evaluation
             try:
@@ -1014,25 +1016,26 @@ def _element_constructor_(self, x):
 
         from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
         R = PolynomialRing(self._base, VarList, order=self._order)
-        if isinstance(R, MPolynomialRing_libsingular) and isinstance(x, MPolynomial_libsingular):  # everything else is so buggy that it's even not worth to try.
-            try:
-                # Problem: If there is only a partial overlap in the variables
-                # of x.parent() and R, then R(x) raises an error (which, I think,
-                # is a bug, since we talk here about conversion, not coercion).
-                # Hence, for being on the safe side, we coerce into a pushout ring:
-                x = R(1) * x
-                return InfinitePolynomial(self, x)
-            except Exception:
-                # OK, last resort, to be on the safe side
+        if isinstance(x, MPolynomial_libsingular):  # everything else is so buggy that it's even not worth to try.
+            from sage.rings.polynomial.multi_polynomial_libsingular import MPolynomialRing_libsingular
+            if isinstance(R, MPolynomialRing_libsingular):
                 try:
-                    return sage_eval(repr(x), self.gens_dict())
-                except (ValueError, TypeError, NameError):
-                    raise ValueError("cannot convert %s into an element of %s; conversion of the underlying polynomial failed" % (x, self))
-        else:
-            try:
-                return sage_eval(repr(x), self.gens_dict())
-            except (ValueError, TypeError, NameError):
-                raise ValueError("cannot convert %s into an element of %s" % (x, self))
+                    # Problem: If there is only a partial overlap in the variables
+                    # of x.parent() and R, then R(x) raises an error (which, I think,
+                    # is a bug, since we talk here about conversion, not coercion).
+                    # Hence, for being on the safe side, we coerce into a pushout ring:
+                    x = R(1) * x
+                    return InfinitePolynomial(self, x)
+                except Exception:
+                    # OK, last resort, to be on the safe side
+                    try:
+                        return sage_eval(repr(x), self.gens_dict())
+                    except (ValueError, TypeError, NameError):
+                        raise ValueError("cannot convert %s into an element of %s; conversion of the underlying polynomial failed" % (x, self))
+        try:
+            return sage_eval(repr(x), self.gens_dict())
+        except (ValueError, TypeError, NameError):
+            raise ValueError("cannot convert %s into an element of %s" % (x, self))
 
     def tensor_with_ring(self, R):
         """
diff --git a/src/sage/rings/polynomial/multi_polynomial.pxd b/src/sage/rings/polynomial/multi_polynomial.pxd
index 44fdf3edfc4..5dc75e6bd3f 100644
--- a/src/sage/rings/polynomial/multi_polynomial.pxd
+++ b/src/sage/rings/polynomial/multi_polynomial.pxd
@@ -1,7 +1,11 @@
-from sage.structure.element cimport CommutativeRingElement
+from .commutative_polynomial cimport CommutativePolynomial
 
-cdef class MPolynomial(CommutativeRingElement):
+
+cdef class MPolynomial(CommutativePolynomial):
     cdef long _hash_c(self) except -1
     cpdef _mod_(self, right)
     cpdef dict _mpoly_dict_recursive(self, tuple vars=*, base_ring=*)
 
+
+cdef class MPolynomial_libsingular(MPolynomial):
+    pass
diff --git a/src/sage/rings/polynomial/multi_polynomial.pyx b/src/sage/rings/polynomial/multi_polynomial.pyx
index 4f8e47049a2..4dc633296cb 100644
--- a/src/sage/rings/polynomial/multi_polynomial.pyx
+++ b/src/sage/rings/polynomial/multi_polynomial.pyx
@@ -20,6 +20,9 @@ from itertools import chain, islice
 from sage.misc.misc_c import prod
 
 def is_MPolynomial(x):
+    from sage.misc.superseded import deprecation
+    deprecation(32709, "the function is_MPolynomial is deprecated; use isinstance(x, sage.structure.element.MPolynomial) instead")
+
     return isinstance(x, MPolynomial)
 
 from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
@@ -32,7 +35,7 @@ from sage.rings.real_mpfr import RealField_class, RealField
 from sage.rings.polynomial.polydict cimport ETuple
 from sage.rings.polynomial.polynomial_element cimport Polynomial
 
-cdef class MPolynomial(CommutativeRingElement):
+cdef class MPolynomial(CommutativePolynomial):
 
     ####################
     # Some standard conversions
@@ -2852,3 +2855,26 @@ cdef remove_from_tuple(e, int ind):
         return w[0]
     else:
         return tuple(w)
+
+
+cdef class MPolynomial_libsingular(MPolynomial):
+    r"""
+    Abstract base class for :class:`~sage.rings.polynomial.multi_polynomial_libsingular.MPolynomial_libsingular`
+
+    This class is defined for the purpose of :func:`isinstance` tests.  It should not be
+    instantiated.
+
+    EXAMPLES::
+
+        sage: R1.<x> = QQ[]
+        sage: isinstance(x, sage.rings.polynomial.multi_polynomial.MPolynomial_libsingular)
+        False
+        sage: R2.<y,z> = QQ[]
+        sage: isinstance(y, sage.rings.polynomial.multi_polynomial.MPolynomial_libsingular)
+        True
+
+    By design, there is a unique direct subclass::
+
+        sage: len(sage.rings.polynomial.multi_polynomial.MPolynomial_libsingular.__subclasses__()) <= 1
+        True
+    """
diff --git a/src/sage/rings/polynomial/multi_polynomial_element.py b/src/sage/rings/polynomial/multi_polynomial_element.py
index 702bf9af7eb..d96e9bfb9cc 100644
--- a/src/sage/rings/polynomial/multi_polynomial_element.py
+++ b/src/sage/rings/polynomial/multi_polynomial_element.py
@@ -61,7 +61,7 @@
 from sage.structure.factorization import Factorization
 from sage.rings.polynomial.polynomial_singular_interface import Polynomial_singular_repr
 from sage.structure.sequence import Sequence
-from .multi_polynomial import MPolynomial
+from .multi_polynomial import MPolynomial, is_MPolynomial
 from sage.categories.morphism import Morphism
 from sage.misc.lazy_attribute import lazy_attribute
 
@@ -71,9 +71,6 @@
 from sage.categories.number_fields import NumberFields
 from sage.rings.real_mpfr import RealField
 
-def is_MPolynomial(x):
-    return isinstance(x, MPolynomial)
-
 
 class MPolynomial_element(MPolynomial):
     def __init__(self, parent, x):
diff --git a/src/sage/rings/polynomial/multi_polynomial_libsingular.pxd b/src/sage/rings/polynomial/multi_polynomial_libsingular.pxd
index f0518c93f9c..c9cec10e2bc 100644
--- a/src/sage/rings/polynomial/multi_polynomial_libsingular.pxd
+++ b/src/sage/rings/polynomial/multi_polynomial_libsingular.pxd
@@ -1,11 +1,11 @@
 from sage.libs.singular.decl cimport poly, ring
 
-from sage.rings.polynomial.multi_polynomial cimport MPolynomial
+from sage.rings.polynomial.multi_polynomial cimport MPolynomial_libsingular as MPolynomial_libsingular_base
 from sage.rings.polynomial.multi_polynomial_ring_base cimport MPolynomialRing_base
 
 cdef class MPolynomialRing_libsingular
 
-cdef class MPolynomial_libsingular(MPolynomial):
+cdef class MPolynomial_libsingular(MPolynomial_libsingular_base):
     cdef poly *_poly
     cdef ring *_parent_ring
     cpdef _add_(self, other)
diff --git a/src/sage/rings/polynomial/multi_polynomial_libsingular.pyx b/src/sage/rings/polynomial/multi_polynomial_libsingular.pyx
index c8d4ad59156..4427921fb70 100644
--- a/src/sage/rings/polynomial/multi_polynomial_libsingular.pyx
+++ b/src/sage/rings/polynomial/multi_polynomial_libsingular.pyx
@@ -1897,7 +1897,7 @@ def unpickle_MPolynomialRing_libsingular(base_ring, names, term_order):
     return _multi_variate(base_ring, tuple(names), None, term_order, None)
 
 
-cdef class MPolynomial_libsingular(MPolynomial):
+cdef class MPolynomial_libsingular(MPolynomial_libsingular_base):
     """
     A multivariate polynomial implemented using libSINGULAR.
     """
diff --git a/src/sage/rings/polynomial/multi_polynomial_ring.py b/src/sage/rings/polynomial/multi_polynomial_ring.py
index f0381318b30..64a49e97ca0 100644
--- a/src/sage/rings/polynomial/multi_polynomial_ring.py
+++ b/src/sage/rings/polynomial/multi_polynomial_ring.py
@@ -420,7 +420,7 @@ def __call__(self, x=0, check=True):
         except TypeError:
             pass
 
-        from .multi_polynomial_libsingular import MPolynomial_libsingular
+        from .multi_polynomial import MPolynomial_libsingular
 
         if isinstance(x, MPolynomial_polydict):
             P = x.parent()
diff --git a/src/sage/rings/polynomial/padics/polynomial_padic_capped_relative_dense.py b/src/sage/rings/polynomial/padics/polynomial_padic_capped_relative_dense.py
index d41eaee7dd1..6aa21db4af7 100644
--- a/src/sage/rings/polynomial/padics/polynomial_padic_capped_relative_dense.py
+++ b/src/sage/rings/polynomial/padics/polynomial_padic_capped_relative_dense.py
@@ -661,7 +661,7 @@ def rshift_coeffs(self, shift, no_list=False):
             return Polynomial_padic_capped_relative_dense(self.parent(), (self._poly // fdiv, 0, [0 if a <= shift else a - shift for a in self._relprecs], False, None, None), construct=True)
 
     # def __floordiv__(self, right):
-    #     if is_Polynomial(right) and right.is_constant() and right[0] in self.base_ring():
+    #     if isinstance(right, Polynomial) and right.is_constant() and right[0] in self.base_ring():
     #         d = self.base_ring()(right[0])
     #     elif (right in self.base_ring()):
     #         d = self.base_ring()(right)
diff --git a/src/sage/rings/polynomial/polynomial_element.pxd b/src/sage/rings/polynomial/polynomial_element.pxd
index 1ba103329c3..083f506b222 100644
--- a/src/sage/rings/polynomial/polynomial_element.pxd
+++ b/src/sage/rings/polynomial/polynomial_element.pxd
@@ -1,11 +1,12 @@
-from sage.structure.element import Element, CommutativeAlgebraElement
+from sage.structure.element import Element
 from sage.structure.element cimport Element, CommutativeAlgebraElement, ModuleElement
 from sage.structure.parent cimport Parent
 from sage.rings.integer cimport Integer
+from .commutative_polynomial cimport CommutativePolynomial
 from .polynomial_compiled cimport CompiledPolynomialFunction
 
 
-cdef class Polynomial(CommutativeAlgebraElement):
+cdef class Polynomial(CommutativePolynomial):
     cdef Polynomial _new_generic(self, list coeffs)
     cdef char _is_gen
     cdef CompiledPolynomialFunction _compiled
diff --git a/src/sage/rings/polynomial/polynomial_element.pyx b/src/sage/rings/polynomial/polynomial_element.pyx
index ceb36976126..3c1ecb7ffa3 100644
--- a/src/sage/rings/polynomial/polynomial_element.pyx
+++ b/src/sage/rings/polynomial/polynomial_element.pyx
@@ -143,6 +143,8 @@ cpdef is_Polynomial(f):
     """
     Return True if f is of type univariate polynomial.
 
+    This function is deprecated.
+
     INPUT:
 
     -  ``f`` -- an object
@@ -152,6 +154,8 @@ cpdef is_Polynomial(f):
         sage: from sage.rings.polynomial.polynomial_element import is_Polynomial
         sage: R.<x> = ZZ[]
         sage: is_Polynomial(x^3 + x + 1)
+        doctest:...: DeprecationWarning: the function is_Polynomial is deprecated; use isinstance(x, sage.structure.element.Polynomial) instead
+        See https://github.com/sagemath/sage/issues/32709 for details.
         True
         sage: S.<y> = R[]
         sage: f = y^3 + x*y -3*x; f
@@ -175,6 +179,9 @@ cpdef is_Polynomial(f):
         sage: is_Polynomial(f)
         False
     """
+    from sage.misc.superseded import deprecation
+    deprecation(32709, "the function is_Polynomial is deprecated; use isinstance(x, sage.structure.element.Polynomial) instead")
+
     return isinstance(f, Polynomial)
 
 from .polynomial_compiled cimport CompiledPolynomialFunction
@@ -182,7 +189,7 @@ from .polynomial_compiled cimport CompiledPolynomialFunction
 from sage.rings.polynomial.polydict cimport ETuple
 
 
-cdef class Polynomial(CommutativeAlgebraElement):
+cdef class Polynomial(CommutativePolynomial):
     """
     A polynomial.
 
diff --git a/src/sage/rings/polynomial/polynomial_integer_dense_ntl.pyx b/src/sage/rings/polynomial/polynomial_integer_dense_ntl.pyx
index c2db35e0bf7..38e2470f1f1 100644
--- a/src/sage/rings/polynomial/polynomial_integer_dense_ntl.pyx
+++ b/src/sage/rings/polynomial/polynomial_integer_dense_ntl.pyx
@@ -51,8 +51,6 @@ from sage.rings.integer_ring import IntegerRing
 from sage.rings.integer_ring cimport IntegerRing_class
 ZZ_sage = IntegerRing()
 
-from sage.rings.polynomial.polynomial_element import is_Polynomial
-
 from sage.libs.ntl.ntl_ZZX cimport ntl_ZZX
 
 from sage.rings.integer_ring import ZZ
@@ -757,7 +755,7 @@ cdef class Polynomial_integer_dense_ntl(Polynomial):
             sage: g // f
             x - 6
         """
-        if is_Polynomial(right) and right.is_constant() and right[0] in ZZ:
+        if isinstance(right, Polynomial) and right.is_constant() and right[0] in ZZ:
             d = ZZ(right[0])
             return self.parent()([c // d for c in self.list()], construct=True)
         elif (right in self.parent().base_ring()):
diff --git a/src/sage/rings/qqbar.py b/src/sage/rings/qqbar.py
index 1e21659d504..da8b8d71b91 100644
--- a/src/sage/rings/qqbar.py
+++ b/src/sage/rings/qqbar.py
@@ -573,7 +573,7 @@
 from sage.rings.complex_interval_field import ComplexIntervalField
 from sage.rings.complex_interval import is_ComplexIntervalFieldElement
 from sage.rings.polynomial.all import PolynomialRing
-from sage.rings.polynomial.polynomial_element import is_Polynomial
+from sage.rings.polynomial.polynomial_element import Polynomial
 from sage.rings.integer_ring import ZZ
 from sage.rings.rational_field import QQ
 from sage.rings.number_field.number_field import NumberField, GaussianField, CyclotomicField
@@ -6724,7 +6724,7 @@ def __init__(self, poly):
             sage: type(P) # indirect doctest
             <class 'sage.rings.qqbar.AlgebraicPolynomialTracker'>
         """
-        if not is_Polynomial(poly):
+        if not isinstance(poly, Polynomial):
             raise ValueError("Trying to create AlgebraicPolynomialTracker on non-Polynomial")
         B = poly.base_ring()
 
diff --git a/src/sage/schemes/berkovich/berkovich_cp_element.py b/src/sage/schemes/berkovich/berkovich_cp_element.py
index b8cc1887957..7c922fefb90 100644
--- a/src/sage/schemes/berkovich/berkovich_cp_element.py
+++ b/src/sage/schemes/berkovich/berkovich_cp_element.py
@@ -83,7 +83,7 @@ def __init__(self, parent, center, radius=None, power=None, prec=20, space_type=
             Type I point centered at 4 + O(5^20)
         """
         from sage.rings.function_field.element import is_FunctionFieldElement
-        from sage.rings.polynomial.polynomial_element import is_Polynomial
+        from sage.rings.polynomial.polynomial_element import Polynomial
         from sage.rings.fraction_field_element import FractionFieldElement_1poly_field
         self._type = None
 
@@ -109,17 +109,17 @@ def __init__(self, parent, center, radius=None, power=None, prec=20, space_type=
 
         # is_FunctionFieldElement calls .parent
         elif hasattr(center, "parent") and hasattr(radius, 'parent'):
-            from sage.rings.polynomial.multi_polynomial_element import is_MPolynomial
-            if is_MPolynomial(center):
+            from sage.rings.polynomial.multi_polynomial import MPolynomial
+            if isinstance(center, MPolynomial):
                 try:
                     center = center.univariate_polynomial()
                 except AttributeError:
                     raise TypeError('center was %s, a multivariable polynomial' % center)
 
             # check if the radius and the center are functions
-            center_func_check = is_FunctionFieldElement(center) or is_Polynomial(center) or\
+            center_func_check = is_FunctionFieldElement(center) or isinstance(center, Polynomial) or\
                 isinstance(center, FractionFieldElement_1poly_field) or isinstance(center, Expression)
-            radius_func_check = is_FunctionFieldElement(radius) or is_Polynomial(radius) or\
+            radius_func_check = is_FunctionFieldElement(radius) or isinstance(radius, Polynomial) or\
                 isinstance(radius, FractionFieldElement_1poly_field) or isinstance(radius, Expression)
 
             if center_func_check:
diff --git a/src/sage/schemes/curves/constructor.py b/src/sage/schemes/curves/constructor.py
index fc3d174de59..c3ff0fa9a42 100644
--- a/src/sage/schemes/curves/constructor.py
+++ b/src/sage/schemes/curves/constructor.py
@@ -36,7 +36,7 @@
 
 from sage.categories.fields import Fields
 
-from sage.rings.polynomial.multi_polynomial_element import is_MPolynomial
+from sage.rings.polynomial.multi_polynomial import MPolynomial
 from sage.rings.polynomial.multi_polynomial_ring import is_MPolynomialRing
 from sage.rings.finite_rings.finite_field_constructor import is_FiniteField
 
@@ -239,7 +239,7 @@ def Curve(F, A=None):
             else:
                 A = ProjectiveSpace(P.ngens()-1, P.base_ring(), names=P.variable_names())
                 A._coordinate_ring = P
-        elif is_MPolynomial(F): # define a plane curve
+        elif isinstance(F, MPolynomial): # define a plane curve
             P = F.parent()
             k = F.base_ring()
 
diff --git a/src/sage/schemes/cyclic_covers/constructor.py b/src/sage/schemes/cyclic_covers/constructor.py
index 0966b0793d5..32bdf239b99 100644
--- a/src/sage/schemes/cyclic_covers/constructor.py
+++ b/src/sage/schemes/cyclic_covers/constructor.py
@@ -8,7 +8,7 @@
 #                  https://www.gnu.org/licenses/
 # *****************************************************************************
 
-from sage.rings.polynomial.polynomial_element import is_Polynomial
+from sage.rings.polynomial.polynomial_element import Polynomial
 from sage.rings.finite_rings.finite_field_constructor import is_FiniteField
 from sage.schemes.affine.affine_space import AffineSpace
 from .cycliccover_generic import CyclicCover_generic
@@ -100,7 +100,7 @@ def CyclicCover(r, f, names=None, check_smooth=True):
 
 
     """
-    if not is_Polynomial(f):
+    if not isinstance(f, Polynomial):
         raise TypeError("Arguments f (= %s) must be a polynomial" % (f,))
     P = f.parent()
     f = P(f)
diff --git a/src/sage/schemes/elliptic_curves/constructor.py b/src/sage/schemes/elliptic_curves/constructor.py
index 1ecdc082cd5..d924bd8f45c 100644
--- a/src/sage/schemes/elliptic_curves/constructor.py
+++ b/src/sage/schemes/elliptic_curves/constructor.py
@@ -29,7 +29,7 @@
 from sage.rings.polynomial.multi_polynomial_ring import is_MPolynomialRing
 from sage.rings.finite_rings.finite_field_constructor import is_FiniteField
 from sage.rings.number_field.number_field import is_NumberField
-from sage.rings.polynomial.multi_polynomial_element import is_MPolynomial
+from sage.rings.polynomial.multi_polynomial import MPolynomial
 from sage.rings.ring import is_Ring
 
 from sage.categories.fields import Fields
@@ -420,7 +420,7 @@ def create_key_and_extra_args(self, x=None, y=None, j=None, minimal_twist=True,
         if isinstance(parent(x), sage.rings.abc.SymbolicRing):
             x = x._polynomial_(rings.QQ['x', 'y'])
 
-        if is_MPolynomial(x):
+        if isinstance(x, MPolynomial):
             if y is None:
                 x = coefficients_from_Weierstrass_polynomial(x)
             else:
diff --git a/src/sage/schemes/elliptic_curves/ell_curve_isogeny.py b/src/sage/schemes/elliptic_curves/ell_curve_isogeny.py
index 441f45c2e12..213ac83d73b 100644
--- a/src/sage/schemes/elliptic_curves/ell_curve_isogeny.py
+++ b/src/sage/schemes/elliptic_curves/ell_curve_isogeny.py
@@ -86,7 +86,7 @@
 from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
 from sage.rings.integer import Integer
 from sage.rings.laurent_series_ring import LaurentSeriesRing
-from sage.rings.polynomial.polynomial_element import is_Polynomial
+from sage.rings.polynomial.polynomial_element import Polynomial
 from sage.rings.fraction_field import FractionField
 
 from sage.schemes.elliptic_curves.all import EllipticCurve
@@ -159,7 +159,7 @@ def _isogeny_determine_algorithm(E, kernel):
         kernel = [kernel]
         kernel_is_list = True
 
-    if is_Polynomial(kernel) or (kernel_is_list and kernel[0] in E.base_ring()):
+    if isinstance(kernel, Polynomial) or (kernel_is_list and kernel[0] in E.base_ring()):
         return "kohel"
 
     if kernel_is_list and kernel[0] in E:
diff --git a/src/sage/schemes/elliptic_curves/jacobian.py b/src/sage/schemes/elliptic_curves/jacobian.py
index b6983ab75f7..6cf48d0760a 100644
--- a/src/sage/schemes/elliptic_curves/jacobian.py
+++ b/src/sage/schemes/elliptic_curves/jacobian.py
@@ -104,8 +104,8 @@ def Jacobian(X, **kwds):
         pass
 
     morphism = kwds.pop('morphism', False)
-    from sage.rings.polynomial.multi_polynomial_element import is_MPolynomial
-    if is_MPolynomial(X):
+    from sage.rings.polynomial.multi_polynomial import MPolynomial
+    if isinstance(X, MPolynomial):
         if morphism:
             from sage.schemes.curves.constructor import Curve
             return Jacobian_of_equation(X, curve=Curve(X), **kwds)
diff --git a/src/sage/schemes/generic/hypersurface.py b/src/sage/schemes/generic/hypersurface.py
index 12138d596af..090ffacf59a 100644
--- a/src/sage/schemes/generic/hypersurface.py
+++ b/src/sage/schemes/generic/hypersurface.py
@@ -18,7 +18,7 @@
 #                  http://www.gnu.org/licenses/
 #*****************************************************************************
 
-from sage.rings.polynomial.multi_polynomial_element import is_MPolynomial
+from sage.rings.polynomial.multi_polynomial import MPolynomial
 from sage.schemes.affine.affine_subscheme import AlgebraicScheme_subscheme_affine
 from sage.schemes.projective.projective_subscheme import AlgebraicScheme_subscheme_projective
 
@@ -92,7 +92,7 @@ def __init__(self, poly, ambient=None):
             sage: H == loads(dumps(H))
             True
         """
-        if not is_MPolynomial(poly):
+        if not isinstance(poly, MPolynomial):
             raise TypeError("Defining polynomial (=%s) must be a multivariate polynomial."%poly)
         if not poly.is_homogeneous():
             raise TypeError("Defining polynomial (=%s) must be homogeneous."%poly)
@@ -177,7 +177,7 @@ def __init__(self, poly, ambient=None):
             sage: H == loads(dumps(H))
             True
         """
-        if not is_MPolynomial(poly):
+        if not isinstance(poly, MPolynomial):
             raise TypeError("Defining polynomial (= %s) must be a multivariate polynomial"%poly)
         if ambient is None:
             R = poly.parent()
diff --git a/src/sage/schemes/hyperelliptic_curves/constructor.py b/src/sage/schemes/hyperelliptic_curves/constructor.py
index 54556e08755..e73efa85558 100644
--- a/src/sage/schemes/hyperelliptic_curves/constructor.py
+++ b/src/sage/schemes/hyperelliptic_curves/constructor.py
@@ -25,7 +25,7 @@
 import sage.rings.abc
 from sage.rings.rational_field import is_RationalField
 from sage.rings.finite_rings.finite_field_constructor import is_FiniteField
-from sage.rings.polynomial.polynomial_element import is_Polynomial
+from sage.rings.polynomial.polynomial_element import Polynomial
 
 from sage.structure.dynamic_class import dynamic_class
 
@@ -197,7 +197,7 @@ def HyperellipticCurve(f, h=0, names=None, PP=None, check_squarefree=True):
     # F is the discriminant; use this for the type check
     # rather than f and h, one of which might be constant.
     F = h**2 + 4*f
-    if not is_Polynomial(F):
+    if not isinstance(F, Polynomial):
         raise TypeError("Arguments f (= %s) and h (= %s) must be polynomials" % (f, h))
     P = F.parent()
     f = P(f)
diff --git a/src/sage/schemes/hyperelliptic_curves/jacobian_homset.py b/src/sage/schemes/hyperelliptic_curves/jacobian_homset.py
index f7b0a1f67f0..74cdccbaa49 100644
--- a/src/sage/schemes/hyperelliptic_curves/jacobian_homset.py
+++ b/src/sage/schemes/hyperelliptic_curves/jacobian_homset.py
@@ -49,7 +49,7 @@
 from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
 from sage.rings.integer_ring import ZZ
 from sage.rings.integer import is_Integer, Integer
-from sage.rings.polynomial.polynomial_element import is_Polynomial
+from sage.rings.polynomial.polynomial_element import Polynomial
 
 from sage.schemes.generic.homset import SchemeHomset_points
 from sage.schemes.generic.morphism import is_SchemeMorphism
@@ -138,15 +138,15 @@ def __call__(self, P):
                     P1 = R(P1)
                     P2 = R(P2)
                     return JacobianMorphism_divisor_class_field(self, (P1, P2))
-                if is_Integer(P1) and is_Polynomial(P2):
+                if is_Integer(P1) and isinstance(P2, Polynomial):
                     R = PolynomialRing(self.value_ring(), 'x')
                     P1 = R(P1)
                     return JacobianMorphism_divisor_class_field(self, (P1, P2))
-                if is_Integer(P2) and is_Polynomial(P1):
+                if is_Integer(P2) and isinstance(P1, Polynomial):
                     R = PolynomialRing(self.value_ring(), 'x')
                     P2 = R(P2)
                     return JacobianMorphism_divisor_class_field(self, (P1, P2))
-                if is_Polynomial(P1) and is_Polynomial(P2):
+                if isinstance(P1, Polynomial) and isinstance(P2, Polynomial):
                     return JacobianMorphism_divisor_class_field(self, tuple(P))
                 if is_SchemeMorphism(P1) and is_SchemeMorphism(P2):
                     return self(P1) - self(P2)
diff --git a/src/sage/schemes/hyperelliptic_curves/monsky_washnitzer.py b/src/sage/schemes/hyperelliptic_curves/monsky_washnitzer.py
index 6d0d2403ba8..687ab8606df 100644
--- a/src/sage/schemes/hyperelliptic_curves/monsky_washnitzer.py
+++ b/src/sage/schemes/hyperelliptic_curves/monsky_washnitzer.py
@@ -71,7 +71,7 @@
 from sage.rings.infinity import Infinity
 from sage.rings.laurent_series_ring import is_LaurentSeriesRing
 from sage.rings.padics.all import pAdicField
-from sage.rings.polynomial.polynomial_element import is_Polynomial
+from sage.rings.polynomial.polynomial_element import Polynomial
 from sage.rings.ring import CommutativeAlgebra
 from sage.schemes.elliptic_curves.constructor import EllipticCurve
 from sage.schemes.elliptic_curves.ell_generic import is_EllipticCurve
@@ -496,7 +496,7 @@ def __init__(self, Q, laurent_series=False):
             ...
             ArithmeticError: 2 and 3 must be invertible in the coefficient ring (=Ring of integers modulo 10) of Q
         """
-        if not is_Polynomial(Q):
+        if not isinstance(Q, Polynomial):
             raise TypeError("Q (=%s) must be a polynomial" % Q)
 
         if Q.degree() != 3 or not Q[2].is_zero():
@@ -2396,7 +2396,7 @@ def __init__(self, Q, R=None, invert_y=True):
             Q = C.hyperelliptic_polynomials()[0].change_ring(R)
             self._curve = C
 
-        if is_Polynomial(Q):
+        if isinstance(Q, Polynomial):
             self._Q = Q.change_ring(R)
             self._coeffs = self._Q.coefficients(sparse=False)
             if self._coeffs.pop() != 1:
diff --git a/src/sage/schemes/plane_conics/constructor.py b/src/sage/schemes/plane_conics/constructor.py
index a2a14190087..abff2935854 100644
--- a/src/sage/schemes/plane_conics/constructor.py
+++ b/src/sage/schemes/plane_conics/constructor.py
@@ -32,7 +32,7 @@
 from sage.rings.ring import IntegralDomain
 from sage.rings.rational_field import is_RationalField
 from sage.rings.finite_rings.finite_field_constructor import is_FiniteField
-from sage.rings.polynomial.multi_polynomial_element import is_MPolynomial
+from sage.rings.polynomial.multi_polynomial import MPolynomial
 from sage.rings.polynomial.polynomial_ring import is_PolynomialRing
 from sage.rings.polynomial.multi_polynomial_ring import is_MPolynomialRing
 from sage.rings.fraction_field import is_FractionField
@@ -204,7 +204,7 @@ def Conic(base_field, F=None, names=None, unique=True):
         temp_ring = PolynomialRing(F.base_ring(), 3, names)
         F = vector(temp_ring.gens()) * F * vector(temp_ring.gens())
 
-    if not is_MPolynomial(F):
+    if not isinstance(F, MPolynomial):
         raise TypeError("F (=%s) must be a three-variable polynomial or "
                         "a sequence of points or coefficients" % F)
 
diff --git a/src/sage/schemes/plane_quartics/quartic_constructor.py b/src/sage/schemes/plane_quartics/quartic_constructor.py
index a0b5ce3c93a..274eaf049b6 100644
--- a/src/sage/schemes/plane_quartics/quartic_constructor.py
+++ b/src/sage/schemes/plane_quartics/quartic_constructor.py
@@ -9,7 +9,7 @@
 #*****************************************************************************
 
 from sage.schemes.projective.projective_space import is_ProjectiveSpace, ProjectiveSpace
-from sage.rings.polynomial.multi_polynomial_element import is_MPolynomial
+from sage.rings.polynomial.multi_polynomial import MPolynomial
 
 from .quartic_generic import QuarticCurve_generic
 
@@ -50,7 +50,7 @@ def QuarticCurve(F, PP=None, check=False):
         ValueError: Argument F (=x^4 + y^4) must be a polynomial in 3 variables
 
     """
-    if not is_MPolynomial(F):
+    if not isinstance(F, MPolynomial):
         raise ValueError(f"Argument F (={F}) must be a multivariate polynomial")
     P = F.parent()
     if not P.ngens() == 3:
diff --git a/src/sage/structure/element.pyx b/src/sage/structure/element.pyx
index ad84ab939ac..4b4ed0a3ab6 100644
--- a/src/sage/structure/element.pyx
+++ b/src/sage/structure/element.pyx
@@ -3361,6 +3361,19 @@ cdef class Expression(CommutativeRingElement):
 
     r"""
     Abstract base class for :class:`~sage.symbolic.expression.Expression`.
+
+    This class is defined for the purpose of :func:`isinstance` tests.  It should not be
+    instantiated.
+
+    EXAMPLES::
+
+        sage: isinstance(SR.var('y'), sage.structure.element.Expression)  # optional - sage.symbolic
+        True
+
+    By design, there is a unique direct subclass::
+
+        sage: len(sage.structure.element.Expression.__subclasses__()) <= 1
+        True
     """
 
     pass
@@ -4301,6 +4314,8 @@ def is_CommutativeAlgebraElement(x):
 cdef class CommutativeAlgebraElement(CommutativeRingElement):
     pass
 
+    ##############################################
+
 def is_InfinityElement(x):
     """
     Return ``True`` if x is of type InfinityElement.