Repository: sympy/sympy
Pull request: polys: Add PuiseuxRing and remove the ring cache
<!-- Your title above should be a short description of what
was changed. Do not include the issue number in the title. -->

#### References to other Issues or PRs
<!-- If this pull request fixes an issue, write "Fixes #NNNN" in that exact
format, e.g. "Fixes #1234" (see
https://tinyurl.com/auto-closing for more information). Also, please
write a comment on that issue linking back to this pull request once it is
open. -->

Fixes gh-24581.

Includes the changes from gh-24585 to remove the PolyRing ring cache and the dynamically created PolyElement classes.

This is analogous to gh-25691 which made changes to DMP in preparation for being able to make it so that DMP can use python-flint internally in gh-25722.

#### Brief description of what is fixed or changed

Remove the ring cache for polynomial rings and the dynamically created PolyElement subclasses. Some other changes are needed to avoid performance regressions as a result of not having all rings cached.

Add a new PuiseuxRing type that can represent "polynomials" that have negative or fractional exponents like `x**(1/2)*y**-1`. This is needed by the `ring_series` module but currently `PolyElement` is abused for this just by inserting invalid exponents into the data structure. That abuse makes it impossible to use python-flint's multivariate polynomials as the internal representation of PolyElement.

The ring_series module tests and docs are changed to use PuiseuxRing when necessary instead of PolyElement to represent series having negative or fractional exponents.

I have also been through the whole test suite verifying that the types of coefficients in PolyElement are always correct. This required a few changes in some places.

Various places that create a polynomial by mutating the zero polynomial have been changed to build up a dict instead.

True division (`a / b`) with polynomials is no longer equivalent to floor division (`a // b`) and now uses `exquo` instead. This is a Python 2 hangover where division like `a / b` with integers would give floor division. Since we now have two separate division operators it does not make sense for true division to return floor division. Various other methods like `__rdiv__` etc where changed to be more consistent with true division.

#### Other comments

In future it would be good to expand on PuiseuxRing to make a proper PuiseuxSeries type that can keep track of the "precision" i.e. the number of terms in a series. This would require changes in the `ring_series` module functions though so for now the PuiseuxRing type just provides the functionality that the ring_series module currently expects.

#### Release Notes

<!-- Write the release notes for this release below between the BEGIN and END
statements. The basic format is a bulleted list with the name of the subpackage
and the release note for this PR. For example:

* solvers
  * Added a new solver for logarithmic equations.

* functions
  * Fixed a bug with log of integers. Formerly, `log(-x)` incorrectly gave `-log(x)`.

* physics.units
  * Corrected a semantical error in the conversion between volt and statvolt which
    reported the volt as being larger than the statvolt.

or if no release note(s) should be included use:

NO ENTRY

See https://github.com/sympy/sympy/wiki/Writing-Release-Notes for more
information on how to write release notes. The bot will check your release
notes automatically to see if they are formatted correctly. -->

<!-- BEGIN RELEASE NOTES -->
* polys
   * A memory leak caused by using many polynomial rings is fixed by not "caching" all rings permanently in memory. Elements of polynomial rings no longer use dynamically created classes.
   * A new PuiseuxRing type is added that can represent "polynomials" with negative or fractional exponents. These are used in the ring_series module to represent truncated Puiseux series.
   * **BREAKING CHANGE**: Using polynomial rings with the ring_series module now only works for series that have nonnegative integer exponents. For series with negative or fractional exponents the PuiseuxRing type must be used instead.
   * **BREAKING CHANGE**: True division like `a / b` with elements of polynomial rings now computes an exact division (`exquo`) rather than floor division. For floor division use `a // b`.
<!-- END RELEASE NOTES -->


Problem statement:
memory leak, not a cache issue
Hello

I have the following code that uses the latest sympy (1.11.1) tested with python 3.8.10 on Windows 10 and with python 3.11.1 on macOS 13.1, the memory is constantly growing. Any ideas why?

```
from sympy import Point, Ray, Circle
while True:
    circle = Circle(Point(0.0, 0.0), 0.5)
    ray = Ray(Point(0.2, 0.3), Point(0.3, 0.2))
    res = circle.intersection(ray)
```

NOTE: the growth in geometry has been addressed, but the scope of the issue is greater than that module so the issue is still open.

Hints from issue discussion:
I can reproduce this but no idea why it happens. Clearing the cache doesn't help. Calling `gc.collect` doesn't help.
With this diff the memory usage stays constant:
```diff
diff --git a/sympy/geometry/ellipse.py b/sympy/geometry/ellipse.py
index 0c1c5d0..35d92af 100644
--- a/sympy/geometry/ellipse.py
+++ b/sympy/geometry/ellipse.py
@@ -14,7 +14,7 @@
 from sympy.core.logic import fuzzy_bool
 from sympy.core.numbers import Rational, oo
 from sympy.core.sorting import ordered
-from sympy.core.symbol import Dummy, uniquely_named_symbol, _symbol
+from sympy.core.symbol import Symbol, Dummy, uniquely_named_symbol, _symbol
 from sympy.simplify import simplify, trigsimp
 from sympy.functions.elementary.miscellaneous import sqrt, Max
 from sympy.functions.elementary.trigonometric import cos, sin
@@ -664,8 +664,8 @@ def intersection(self, o):
         [Point2D(-17/5, -12/5), Point2D(-17/5, 12/5), Point2D(7/5, -12/5), Point2D(7/5, 12/5)]
         """
         # TODO: Replace solve with nonlinsolve, when nonlinsolve will be able to solve in real domain
-        x = Dummy('x', real=True)
-        y = Dummy('y', real=True)
+        x = Symbol('x', real=True)
+        y = Symbol('y', real=True)
 
         if isinstance(o, Point):
             if o in self:
```
So what is happening is that every time `intersection` is called two new `Dummy` symbols are created. Somewhere there is some sort of cache that grows without bound if more and more symbols are created.
It's here (and a similar one in fields.py):
https://github.com/sympy/sympy/blob/f8e33851e174bd686ac8dc88d75f955cf5dc8eeb/sympy/polys/rings.py#L194
This cache grows without bound every time a polynomial ring with new symbols is used.

A fix could be to use `lru_cache` instead of a global cache.

I am not entirely sure that this is only a cache though because it also associates with a particular PolyElement class and it is possible that that class needs to be globally unique.
Can't that be solved by using a WeakDictionary? 

I've never been a huge fan of this class design in the polys FWIW.
I think that a WeakValue dictionary could work but it would probably not achieve the performance benefits intended by the cache in the first place. The point is that this second call is faster:
```python
In [1]: %time K = QQ[x,y]
CPU times: user 4 ms, sys: 0 ns, total: 4 ms
Wall time: 2.38 ms

In [2]: %time K = QQ[x,y]
CPU times: user 0 ns, sys: 0 ns, total: 0 ns
Wall time: 89.6 µs
```
Usually these rings are created transiently as part of an operation like `factor`, `cancel` etc. Since the ring is only a temporary object there will probably not be many in existence at any one time so with a WeakValue dictionary the ring will always disappear. In that case we would be better off just not having any cache. On the other hand of many operations like `factor`, `cancel` etc are used as part of a complex operation it is likely that the same ring will be reconstructed many times so keeping a cache can save time on that.

I think that probably a LRU cache is good here but given the size of these rings it should probably be kept quite small:
```diff
diff --git a/sympy/polys/fields.py b/sympy/polys/fields.py
index a3f239c..cdb079a 100644
--- a/sympy/polys/fields.py
+++ b/sympy/polys/fields.py
@@ -6,6 +6,7 @@
 
 from operator import add, mul, lt, le, gt, ge
 
+from sympy.core.cache import __cacheit as _cacheit
 from sympy.core.expr import Expr
 from sympy.core.mod import Mod
 from sympy.core.numbers import Exp1
@@ -99,11 +100,11 @@ def sfield(exprs, *symbols, **options):
     else:
         return (_field, fracs)
 
-_field_cache: dict[Any, Any] = {}
 
 class FracField(DefaultPrinting):
     """Multivariate distributed rational function field. """
 
+    @_cacheit(10)
     def __new__(cls, symbols, domain, order=lex):
         from sympy.polys.rings import PolyRing
         ring = PolyRing(symbols, domain, order)
@@ -113,7 +114,7 @@ def __new__(cls, symbols, domain, order=lex):
         order = ring.order
 
         _hash_tuple = (cls.__name__, symbols, ngens, domain, order)
-        obj = _field_cache.get(_hash_tuple)
+        obj = None
 
         if obj is None:
             obj = object.__new__(cls)
@@ -138,8 +139,6 @@ def __new__(cls, symbols, domain, order=lex):
                     if not hasattr(obj, name):
                         setattr(obj, name, generator)
 
-            _field_cache[_hash_tuple] = obj
-
         return obj
 
     def _gens(self):
diff --git a/sympy/polys/rings.py b/sympy/polys/rings.py
index 0db1897..3a56486 100644
--- a/sympy/polys/rings.py
+++ b/sympy/polys/rings.py
@@ -7,6 +7,7 @@
 from functools import reduce
 from types import GeneratorType
 
+from sympy.core.cache import __cacheit as _cacheit
 from sympy.core.expr import Expr
 from sympy.core.numbers import igcd, oo
 from sympy.core.symbol import Symbol, symbols as _symbols
@@ -191,11 +192,11 @@ def _parse_symbols(symbols):
 
     raise GeneratorsError("expected a string, Symbol or expression or a non-empty sequence of strings, Symbols or expressions")
 
-_ring_cache: dict[Any, Any] = {}
 
 class PolyRing(DefaultPrinting, IPolys):
     """Multivariate distributed polynomial ring. """
 
+    @_cacheit(10)
     def __new__(cls, symbols, domain, order=lex):
         symbols = tuple(_parse_symbols(symbols))
         ngens = len(symbols)
@@ -203,7 +204,7 @@ def __new__(cls, symbols, domain, order=lex):
         order = OrderOpt.preprocess(order)
 
         _hash_tuple = (cls.__name__, symbols, ngens, domain, order)
-        obj = _ring_cache.get(_hash_tuple)
+        obj = None
 
         if obj is None:
             if domain.is_Composite and set(symbols) & set(domain.symbols):
@@ -257,8 +258,6 @@ def __new__(cls, symbols, domain, order=lex):
                     if not hasattr(obj, name):
                         setattr(obj, name, generator)
 
-            _ring_cache[_hash_tuple] = obj
-
         return obj
 
     def _gens(self):
```
In the long run the best solution would be to make the ring objects lighter weight and faster to construct.
> With this diff the memory usage stays constant:

Those two Dummy symbols should be created at the top of ellipse.py and then used the 4 times when needed in the various classes -- there is no need to keep creating new Dummy symbols.
To be clear plenty of other operations can cause this cache to grow and even to grow more quickly:
```
while True:
    cancel(Dummy('x') + 1)
```
I thought the point was that these objects break if they aren't singletonized? I seem to remember that sort of thing being the case for these classes in the polys, but correct me if I am wrong. If that's the case, you'd need a version of lru_cache that also acts like a weak dictionary (i.e., least recently used items aren't removed from the cache if they are also referenced somewhere). 
> Those two Dummy symbols should be created at the top of ellipse.py and then used the 4 times when needed in the various classes -- there is no need to keep creating new Dummy symbols.

Wouldn't that change the semantics? A dummy symbol is supposed to be unequal to everything except for itself. Two separate intersections should have dummies that are unequal, but reusing the same Dummy would make them equal. 
> I thought the point was that these objects break if they aren't singletonized?

Yes, actually they do. I just tried removing the cache and a bunch of poly tests failed. The problem is e.g. this:
https://github.com/sympy/sympy/blob/f8e33851e174bd686ac8dc88d75f955cf5dc8eeb/sympy/polys/rings.py#L419
Each ring dynamically creates a class for its elements and distinct copies of the same ring would have different classes.

I've opened #24585 which for now removes the ring cache altogether so we can see if there's a noticeable impact on performance and if the need for the ring cache is gone.

Longer term I would like to remove the dynamic class generation altogether. Also we should just make the construction of rings faster.

Commit message: fix(integrals): Fix PRDE inexact division

In the risch integrator Poly.quo_ground was used when the ground domain
was a ring leading to an invalid PolyElement with negative exponent.
This commit converts the ground domain to a field so that the division
can be exact.

Patch to review:
diff --git a/doc/src/modules/polys/ringseries.rst b/doc/src/modules/polys/ringseries.rst
index 43d5a9e6ce7e..0a57ede91961 100644
--- a/doc/src/modules/polys/ringseries.rst
+++ b/doc/src/modules/polys/ringseries.rst
@@ -35,47 +35,22 @@ Taylor series, we extend it to allow Laurent and even Puiseux series (with
 fractional exponents)::
 
     >>> from sympy.polys.ring_series import rs_cos, rs_tan
-    >>> R, x, y = ring('x, y', QQ)
+    >>> from sympy.polys.puiseux import puiseux_ring
+    >>> R, x, y = puiseux_ring('x, y', QQ)
 
     >>> rs_cos(x + x*y, x, 3)/x**3
-    -1/2*x**(-1)*y**2 - x**(-1)*y - 1/2*x**(-1) + x**(-3)
+    x**(-3) + -1/2*x**(-1) + -1*x**(-1)*y + -1/2*x**(-1)*y**2
 
     >>> rs_tan(x**QQ(2, 5)*y**QQ(1, 2), x, 2)
-    1/3*x**(6/5)*y**(3/2) + x**(2/5)*y**(1/2)
-
-By default, ``PolyElement`` did not allow non-natural numbers as exponents. It
-converted a fraction to an integer and raised an error on getting negative
-exponents. The goal of the ``ring series`` module is fast series expansion, and
-not to use the ``polys`` module. The reason we use it as our backend is simply
-because it implements a sparse representation and most of the basic functions
-that we need. However, this default behaviour of ``polys`` was limiting for
-``ring series``.
-
-Note that there is no such constraint (in having rational exponents) in the
-data-structure used by ``polys``- ``dict``. Sparse polynomials
-(``PolyElement``) use the Python dict to store a polynomial term by term, where
-a tuple of exponents is the key and the coefficient of that term is the value.
-There is no reason why we can't have rational values in the ``dict`` so as to
-support rational exponents.
-
-So the approach we took was to modify sparse ``polys`` to allow non-natural
-exponents. And it turned out to be quite simple. We only had to delete the
-conversion to ``int`` of exponents in the ``__pow__`` method of
-``PolyElement``. So::
-
-    >>> x**QQ(3, 4)
-    x**(3/4)
-
-and not ``1`` as was the case earlier.
-
-Though this change violates the definition of a polynomial, it doesn't break
-anything yet.  Ideally, we shouldn't modify ``polys`` in any way. But to have
-all the ``series`` capabilities we want, no other simple way was found. If need
-be, we can separate the modified part of ``polys`` from core ``polys``. It
-would be great if any other elegant solution is found.
-
-All series returned by the functions of this module are instances of the
-``PolyElement`` class. To use them with other SymPy types, convert them  to
+    x**(2/5)*y**(1/2) + 1/3*x**(6/5)*y**(3/2)
+
+Since polynomial rings cannot handle negative or fractional exponents, we use
+the :func:`sympy.polys.puiseux.puiseux_ring` function to create a ring that can
+represent such series.
+
+All series returned by the functions of this module are instances of
+``PolyElement`` or ``PuiseuxPoly``. To use them with other SymPy types, convert
+them  to
 ``Expr``::
 
     >>> from sympy.polys.ring_series import rs_exp
@@ -213,6 +188,7 @@ by ``polys.ring.ring``.
 
 **Utility functions**
 
+.. autofunction:: rs_series
 .. autofunction:: rs_is_puiseux
 .. autofunction:: rs_puiseux
 .. autofunction:: rs_puiseux2
@@ -220,3 +196,15 @@ by ``polys.ring.ring``.
 .. autofunction:: rs_fun
 .. autofunction:: mul_xin
 .. autofunction:: pow_xin
+
+**Puiseux rings**
+
+.. currentmodule:: sympy.polys.puiseux
+
+.. autofunction:: puiseux_ring
+
+.. autoclass:: PuiseuxRing
+    :members:
+
+.. autoclass:: PuiseuxPoly
+    :members:
diff --git a/pyproject.toml b/pyproject.toml
index 3430924e523c..67391201b5f9 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -30,6 +30,14 @@ markers = [
     "tooslow",
 ]
 
+[tool.coverage.report]
+
+exclude_lines = [
+    "pragma: no cover",
+    "if TYPE_CHECKING:",
+    "assert False",
+]
+
 [tool.ruff]
 # Enable Pyflakes `E` and `F` codes by default.
 lint.select = [
diff --git a/sympy/integrals/prde.py b/sympy/integrals/prde.py
index 4488cbfc4000..28e91ea0ff3a 100644
--- a/sympy/integrals/prde.py
+++ b/sympy/integrals/prde.py
@@ -533,7 +533,7 @@ def param_poly_rischDE(a, b, q, n, DE):
     if a.is_ground:
         # Normalization: a = 1.
         a = a.LC()
-        b, q = b.quo_ground(a), [qi.quo_ground(a) for qi in q]
+        b, q = b.to_field().exquo_ground(a), [qi.to_field().exquo_ground(a) for qi in q]
 
         if not b.is_zero and (DE.case == 'base' or
                 b.degree() > max(0, DE.d.degree() - 1)):
diff --git a/sympy/integrals/tests/test_integrals.py b/sympy/integrals/tests/test_integrals.py
index 51b81775bb14..41e1ef3aa363 100644
--- a/sympy/integrals/tests/test_integrals.py
+++ b/sympy/integrals/tests/test_integrals.py
@@ -1813,8 +1813,8 @@ def test_issue_15810():
 def test_issue_21024():
     x = Symbol('x', real=True, nonzero=True)
     f = log(x)*log(4*x) + log(3*x + exp(2))
-    F = x*log(x)**2 + x*(1 - 2*log(2)) + (-2*x + 2*x*log(2))*log(x) + \
-        (x + exp(2)/6)*log(3*x + exp(2)) + exp(2)*log(3*x + exp(2))/6
+    F = x*log(x)**2 + x*log(3*x + exp(2)) + x*(1 - 2*log(2)) + \
+        (-2*x + 2*x*log(2))*log(x) + exp(2)*log(3*x + exp(2))/3
     assert F == integrate(f, x)
 
     f = (x + exp(3))/x**2
diff --git a/sympy/polys/domains/complexfield.py b/sympy/polys/domains/complexfield.py
index de02e46d190b..69f0bff2c1b3 100644
--- a/sympy/polys/domains/complexfield.py
+++ b/sympy/polys/domains/complexfield.py
@@ -142,10 +142,7 @@ def from_RealField(self, element, base):
         return self.dtype(element)
 
     def from_ComplexField(self, element, base):
-        if self == base:
-            return element
-        else:
-            return self.dtype(element)
+        return self.dtype(element)
 
     def get_ring(self):
         """Returns a ring associated with ``self``. """
diff --git a/sympy/polys/domains/domain.py b/sympy/polys/domains/domain.py
index 1c2b0d3171d6..1d7fc1eac618 100644
--- a/sympy/polys/domains/domain.py
+++ b/sympy/polys/domains/domain.py
@@ -116,7 +116,7 @@ class (``dtype``) for the elements of the domain. For example the
     ZZ[x]
     >>> type(K)             # class of the domain
     <class 'sympy.polys.domains.polynomialring.PolynomialRing'>
-    >>> K.dtype             # class of the elements
+    >>> K.dtype             # doctest: +SKIP
     <class 'sympy.polys.rings.PolyElement'>
     >>> p_expr = x**2 + 1   # Expr
     >>> p_expr
@@ -469,7 +469,7 @@ def convert(self, element, base=None):
 
     def of_type(self, element):
         """Check if ``a`` is of type ``dtype``. """
-        return isinstance(element, self.tp) # XXX: this isn't correct, e.g. PolyElement
+        return isinstance(element, self.tp)
 
     def __contains__(self, a):
         """Check if ``a`` belongs to this domain. """
diff --git a/sympy/polys/domains/fractionfield.py b/sympy/polys/domains/fractionfield.py
index 47bc25436b8e..78f5054ddd54 100644
--- a/sympy/polys/domains/fractionfield.py
+++ b/sympy/polys/domains/fractionfield.py
@@ -37,6 +37,10 @@ def __init__(self, domain_or_field, symbols=None, order=None):
     def new(self, element):
         return self.field.field_new(element)
 
+    def of_type(self, element):
+        """Check if ``a`` is of type ``dtype``. """
+        return self.field.is_element(element)
+
     @property
     def zero(self):
         return self.field.zero
@@ -53,13 +57,13 @@ def __str__(self):
         return str(self.domain) + '(' + ','.join(map(str, self.symbols)) + ')'
 
     def __hash__(self):
-        return hash((self.__class__.__name__, self.dtype.field, self.domain, self.symbols))
+        return hash((self.__class__.__name__, self.field, self.domain, self.symbols))
 
     def __eq__(self, other):
         """Returns ``True`` if two domains are equivalent. """
-        return isinstance(other, FractionField) and \
-            (self.dtype.field, self.domain, self.symbols) ==\
-            (other.dtype.field, other.domain, other.symbols)
+        if not isinstance(other, FractionField):
+            return NotImplemented
+        return self.field == other.field
 
     def to_sympy(self, a):
         """Convert ``a`` to a SymPy object. """
diff --git a/sympy/polys/domains/polynomialring.py b/sympy/polys/domains/polynomialring.py
index bad73208f866..daccdcdede4d 100644
--- a/sympy/polys/domains/polynomialring.py
+++ b/sympy/polys/domains/polynomialring.py
@@ -43,6 +43,10 @@ def __init__(self, domain_or_ring, symbols=None, order=None):
     def new(self, element):
         return self.ring.ring_new(element)
 
+    def of_type(self, element):
+        """Check if ``a`` is of type ``dtype``. """
+        return self.ring.is_element(element)
+
     @property
     def zero(self):
         return self.ring.zero
@@ -59,13 +63,13 @@ def __str__(self):
         return str(self.domain) + '[' + ','.join(map(str, self.symbols)) + ']'
 
     def __hash__(self):
-        return hash((self.__class__.__name__, self.dtype.ring, self.domain, self.symbols))
+        return hash((self.__class__.__name__, self.ring, self.domain, self.symbols))
 
     def __eq__(self, other):
         """Returns `True` if two domains are equivalent. """
-        return isinstance(other, PolynomialRing) and \
-            (self.dtype.ring, self.domain, self.symbols) == \
-            (other.dtype.ring, other.domain, other.symbols)
+        if not isinstance(other, PolynomialRing):
+            return NotImplemented
+        return self.ring == other.ring
 
     def is_unit(self, a):
         """Returns ``True`` if ``a`` is a unit of ``self``"""
diff --git a/sympy/polys/domains/realfield.py b/sympy/polys/domains/realfield.py
index 79ada6f70737..cb7fac2218c1 100644
--- a/sympy/polys/domains/realfield.py
+++ b/sympy/polys/domains/realfield.py
@@ -171,10 +171,7 @@ def from_AlgebraicField(self, element, base):
         return self.from_sympy(base.to_sympy(element).evalf(self.dps))
 
     def from_RealField(self, element, base):
-        if self == base:
-            return element
-        else:
-            return self.dtype(element)
+        return self.dtype(element)
 
     def from_ComplexField(self, element, base):
         if not element.imag:
diff --git a/sympy/polys/domains/tests/test_domains.py b/sympy/polys/domains/tests/test_domains.py
index 5cdcb9f0403a..405665726b49 100644
--- a/sympy/polys/domains/tests/test_domains.py
+++ b/sympy/polys/domains/tests/test_domains.py
@@ -19,9 +19,9 @@
 from sympy.polys.domains.realfield import RealField
 
 from sympy.polys.numberfields.subfield import field_isomorphism
-from sympy.polys.rings import ring
+from sympy.polys.rings import ring, PolyElement
 from sympy.polys.specialpolys import cyclotomic_poly
-from sympy.polys.fields import field
+from sympy.polys.fields import field, FracElement
 
 from sympy.polys.agca.extensions import FiniteExtension
 
@@ -657,7 +657,12 @@ def test_Domain_is_unit():
 def test_Domain_convert():
 
     def check_element(e1, e2, K1, K2, K3):
-        assert type(e1) is type(e2), '%s, %s: %s %s -> %s' % (e1, e2, K1, K2, K3)
+        if isinstance(e1, PolyElement):
+            assert isinstance(e2, PolyElement) and e1.ring == e2.ring
+        elif isinstance(e1, FracElement):
+            assert isinstance(e2, FracElement) and e1.field == e2.field
+        else:
+            assert type(e1) is type(e2), '%s, %s: %s %s -> %s' % (e1, e2, K1, K2, K3)
         assert e1 == e2, '%s, %s: %s %s -> %s' % (e1, e2, K1, K2, K3)
 
     def check_domains(K1, K2):
diff --git a/sympy/polys/fields.py b/sympy/polys/fields.py
index e45063b5f7ad..ee844df55690 100644
--- a/sympy/polys/fields.py
+++ b/sympy/polys/fields.py
@@ -1,7 +1,6 @@
 """Sparse rational function fields. """
 
 from __future__ import annotations
-from typing import Any
 from functools import reduce
 
 from operator import add, mul, lt, le, gt, ge
@@ -100,7 +99,6 @@ def sfield(exprs, *symbols, **options):
     else:
         return (_field, fracs)
 
-_field_cache: dict[Any, Any] = {}
 
 class FracField(DefaultPrinting):
     """Multivariate distributed rational function field. """
@@ -120,32 +118,29 @@ def __new__(cls, symbols, domain, order=lex):
         order = ring.order
 
         _hash_tuple = (cls.__name__, symbols, ngens, domain, order)
-        obj = _field_cache.get(_hash_tuple)
 
-        if obj is None:
-            obj = object.__new__(cls)
-            obj._hash_tuple = _hash_tuple
-            obj._hash = hash(_hash_tuple)
-            obj.ring = ring
-            obj.dtype = type("FracElement", (FracElement,), {"field": obj})
-            obj.symbols = symbols
-            obj.ngens = ngens
-            obj.domain = domain
-            obj.order = order
+        obj = object.__new__(cls)
+        obj._hash_tuple = _hash_tuple
+        obj._hash = hash(_hash_tuple)
+        obj.ring = ring
+        obj.symbols = symbols
+        obj.ngens = ngens
+        obj.domain = domain
+        obj.order = order
 
-            obj.zero = obj.dtype(ring.zero)
-            obj.one = obj.dtype(ring.one)
+        obj.dtype = FracElement(obj, ring.zero).raw_new
 
-            obj.gens = obj._gens()
+        obj.zero = obj.dtype(ring.zero)
+        obj.one = obj.dtype(ring.one)
 
-            for symbol, generator in zip(obj.symbols, obj.gens):
-                if isinstance(symbol, Symbol):
-                    name = symbol.name
+        obj.gens = obj._gens()
 
-                    if not hasattr(obj, name):
-                        setattr(obj, name, generator)
+        for symbol, generator in zip(obj.symbols, obj.gens):
+            if isinstance(symbol, Symbol):
+                name = symbol.name
 
-            _field_cache[_hash_tuple] = obj
+                if not hasattr(obj, name):
+                    setattr(obj, name, generator)
 
         return obj
 
@@ -160,7 +155,7 @@ def __hash__(self):
         return self._hash
 
     def index(self, gen):
-        if isinstance(gen, self.dtype):
+        if self.is_element(gen):
             return self.ring.index(gen.to_poly())
         else:
             raise ValueError("expected a %s, got %s instead" % (self.dtype,gen))
@@ -173,8 +168,13 @@ def __eq__(self, other):
     def __ne__(self, other):
         return not self == other
 
+    def is_element(self, element):
+        """True if ``element`` is an element of this field. False otherwise. """
+        return isinstance(element, FracElement) and element.field == self
+
     def raw_new(self, numer, denom=None):
         return self.dtype(numer, denom)
+
     def new(self, numer, denom=None):
         if denom is None: denom = self.ring.one
         numer, denom = numer.cancel(denom)
@@ -292,17 +292,19 @@ def to_ring(self):
 class FracElement(DomainElement, DefaultPrinting, CantSympify):
     """Element of multivariate distributed rational function field. """
 
-    def __init__(self, numer, denom=None):
+    def __init__(self, field, numer, denom=None):
         if denom is None:
-            denom = self.field.ring.one
+            denom = field.ring.one
         elif not denom:
             raise ZeroDivisionError("zero denominator")
 
+        self.field = field
         self.numer = numer
         self.denom = denom
 
-    def raw_new(f, numer, denom):
-        return f.__class__(numer, denom)
+    def raw_new(f, numer, denom=None):
+        return f.__class__(f.field, numer, denom)
+
     def new(f, numer, denom):
         return f.raw_new(*numer.cancel(denom))
 
@@ -356,7 +358,7 @@ def sort_key(self):
         return (self.denom.sort_key(), self.numer.sort_key())
 
     def _cmp(f1, f2, op):
-        if isinstance(f2, f1.field.dtype):
+        if f1.field.is_element(f2):
             return op(f1.sort_key(), f2.sort_key())
         else:
             return NotImplemented
@@ -406,12 +408,12 @@ def __add__(f, g):
             return f
         elif not f:
             return g
-        elif isinstance(g, field.dtype):
+        elif field.is_element(g):
             if f.denom == g.denom:
                 return f.new(f.numer + g.numer, f.denom)
             else:
                 return f.new(f.numer*g.denom + f.denom*g.numer, f.denom*g.denom)
-        elif isinstance(g, field.ring.dtype):
+        elif field.ring.is_element(g):
             return f.new(f.numer + f.denom*g, f.denom)
         else:
             if isinstance(g, FracElement):
@@ -430,7 +432,7 @@ def __add__(f, g):
         return f.__radd__(g)
 
     def __radd__(f, c):
-        if isinstance(c, f.field.ring.dtype):
+        if f.field.ring.is_element(c):
             return f.new(f.numer + f.denom*c, f.denom)
 
         op, g_numer, g_denom = f._extract_ground(c)
@@ -450,12 +452,12 @@ def __sub__(f, g):
             return f
         elif not f:
             return -g
-        elif isinstance(g, field.dtype):
+        elif field.is_element(g):
             if f.denom == g.denom:
                 return f.new(f.numer - g.numer, f.denom)
             else:
                 return f.new(f.numer*g.denom - f.denom*g.numer, f.denom*g.denom)
-        elif isinstance(g, field.ring.dtype):
+        elif field.ring.is_element(g):
             return f.new(f.numer - f.denom*g, f.denom)
         else:
             if isinstance(g, FracElement):
@@ -481,7 +483,7 @@ def __sub__(f, g):
             return f.new(f.numer*g_denom - f.denom*g_numer, f.denom*g_denom)
 
     def __rsub__(f, c):
-        if isinstance(c, f.field.ring.dtype):
+        if f.field.ring.is_element(c):
             return f.new(-f.numer + f.denom*c, f.denom)
 
         op, g_numer, g_denom = f._extract_ground(c)
@@ -499,9 +501,9 @@ def __mul__(f, g):
 
         if not f or not g:
             return field.zero
-        elif isinstance(g, field.dtype):
+        elif field.is_element(g):
             return f.new(f.numer*g.numer, f.denom*g.denom)
-        elif isinstance(g, field.ring.dtype):
+        elif field.ring.is_element(g):
             return f.new(f.numer*g, f.denom)
         else:
             if isinstance(g, FracElement):
@@ -520,7 +522,7 @@ def __mul__(f, g):
         return f.__rmul__(g)
 
     def __rmul__(f, c):
-        if isinstance(c, f.field.ring.dtype):
+        if f.field.ring.is_element(c):
             return f.new(f.numer*c, f.denom)
 
         op, g_numer, g_denom = f._extract_ground(c)
@@ -538,9 +540,9 @@ def __truediv__(f, g):
 
         if not g:
             raise ZeroDivisionError
-        elif isinstance(g, field.dtype):
+        elif field.is_element(g):
             return f.new(f.numer*g.denom, f.denom*g.numer)
-        elif isinstance(g, field.ring.dtype):
+        elif field.ring.is_element(g):
             return f.new(f.numer, f.denom*g)
         else:
             if isinstance(g, FracElement):
@@ -568,7 +570,7 @@ def __truediv__(f, g):
     def __rtruediv__(f, c):
         if not f:
             raise ZeroDivisionError
-        elif isinstance(c, f.field.ring.dtype):
+        elif f.field.ring.is_element(c):
             return f.new(f.denom*c, f.numer)
 
         op, g_numer, g_denom = f._extract_ground(c)
diff --git a/sympy/polys/modulargcd.py b/sympy/polys/modulargcd.py
index 20dfd33d9197..00d1920f69fe 100644
--- a/sympy/polys/modulargcd.py
+++ b/sympy/polys/modulargcd.py
@@ -609,7 +609,8 @@ def _chinese_remainder_reconstruction_multivariate(hp, hq, p, q):
     hpmonoms.difference_update(monoms)
     hqmonoms.difference_update(monoms)
 
-    zero = hp.ring.domain.zero
+    domain = hp.ring.domain
+    zero = domain.zero
 
     hpq = hp.ring.zero
 
@@ -617,7 +618,7 @@ def _chinese_remainder_reconstruction_multivariate(hp, hq, p, q):
         crt_ = _chinese_remainder_reconstruction_multivariate
     else:
         def crt_(cp, cq, p, q):
-            return crt([p, q], [cp, cq], symmetric=True)[0]
+            return domain(crt([p, q], [cp, cq], symmetric=True)[0])
 
     for monom in monoms:
         hpq[monom] = crt_(hp[monom], hq[monom], p, q)
diff --git a/sympy/polys/monomials.py b/sympy/polys/monomials.py
index f464ba97f137..e5897a09986d 100644
--- a/sympy/polys/monomials.py
+++ b/sympy/polys/monomials.py
@@ -4,6 +4,7 @@
 from itertools import combinations_with_replacement, product
 from textwrap import dedent
 
+from sympy.core.cache import cacheit
 from sympy.core import Mul, S, Tuple, sympify
 from sympy.polys.polyerrors import ExactQuotientFailed
 from sympy.polys.polyutils import PicklableWithSlots, dict_from_expr
@@ -394,8 +395,14 @@ def term_div(a, b, domain):
 class MonomialOps:
     """Code generator of fast monomial arithmetic functions. """
 
-    def __init__(self, ngens):
-        self.ngens = ngens
+    @cacheit
+    def __new__(cls, ngens):
+        obj = super().__new__(cls)
+        obj.ngens = ngens
+        return obj
+
+    def __getnewargs__(self):
+        return (self.ngens,)
 
     def _build(self, code, name):
         ns = {}
@@ -405,6 +412,7 @@ def _build(self, code, name):
     def _vars(self, name):
         return [ "%s%s" % (name, i) for i in range(self.ngens) ]
 
+    @cacheit
     def mul(self):
         name = "monomial_mul"
         template = dedent("""\
@@ -419,6 +427,7 @@ def %(name)s(A, B):
         code = template % {"name": name, "A": ", ".join(A), "B": ", ".join(B), "AB": ", ".join(AB)}
         return self._build(code, name)
 
+    @cacheit
     def pow(self):
         name = "monomial_pow"
         template = dedent("""\
@@ -431,6 +440,7 @@ def %(name)s(A, k):
         code = template % {"name": name, "A": ", ".join(A), "Ak": ", ".join(Ak)}
         return self._build(code, name)
 
+    @cacheit
     def mulpow(self):
         name = "monomial_mulpow"
         template = dedent("""\
@@ -445,6 +455,7 @@ def %(name)s(A, B, k):
         code = template % {"name": name, "A": ", ".join(A), "B": ", ".join(B), "ABk": ", ".join(ABk)}
         return self._build(code, name)
 
+    @cacheit
     def ldiv(self):
         name = "monomial_ldiv"
         template = dedent("""\
@@ -459,6 +470,7 @@ def %(name)s(A, B):
         code = template % {"name": name, "A": ", ".join(A), "B": ", ".join(B), "AB": ", ".join(AB)}
         return self._build(code, name)
 
+    @cacheit
     def div(self):
         name = "monomial_div"
         template = dedent("""\
@@ -475,6 +487,7 @@ def %(name)s(A, B):
         code = template % {"name": name, "A": ", ".join(A), "B": ", ".join(B), "RAB": "\n    ".join(RAB), "R": ", ".join(R)}
         return self._build(code, name)
 
+    @cacheit
     def lcm(self):
         name = "monomial_lcm"
         template = dedent("""\
@@ -489,6 +502,7 @@ def %(name)s(A, B):
         code = template % {"name": name, "A": ", ".join(A), "B": ", ".join(B), "AB": ", ".join(AB)}
         return self._build(code, name)
 
+    @cacheit
     def gcd(self):
         name = "monomial_gcd"
         template = dedent("""\
diff --git a/sympy/polys/puiseux.py b/sympy/polys/puiseux.py
new file mode 100644
index 000000000000..f02c86ddc1ae
--- /dev/null
+++ b/sympy/polys/puiseux.py
@@ -0,0 +1,795 @@
+"""
+Puiseux rings. These are used by the ring_series module to represented
+truncated Puiseux series. Elements of a Puiseux ring are like polynomials
+except that the exponents can be negative or rational rather than just
+non-negative integers.
+"""
+
+# Previously the ring_series module used PolyElement to represent Puiseux
+# series. This is problematic because it means that PolyElement has to support
+# negative and non-integer exponents which most polynomial representations do
+# not support. This module provides an implementation of a ring for Puiseux
+# series that can be used by ring_series without breaking the basic invariants
+# of polynomial rings.
+#
+# Ideally there would be more of a proper series type that can keep track of
+# not not just the leading terms of a truncated series but also the precision
+# of the series. For now the rings here are just introduced to keep the
+# interface that ring_series was using before.
+
+from __future__ import annotations
+
+from sympy.polys.domains import QQ
+from sympy.polys.rings import PolyRing, PolyElement
+from sympy.core.add import Add
+from sympy.core.mul import Mul
+from sympy.external.gmpy import gcd, lcm
+
+
+from typing import TYPE_CHECKING
+
+
+if TYPE_CHECKING:
+    from typing import Any, Unpack
+    from sympy.core.expr import Expr
+    from sympy.polys.domains import Domain
+    from collections.abc import Iterable, Iterator
+
+
+def puiseux_ring(
+    symbols: str | list[Expr], domain: Domain
+) -> tuple[PuiseuxRing, Unpack[tuple[PuiseuxPoly, ...]]]:
+    """Construct a Puiseux ring.
+
+    This function constructs a Puiseux ring with the given symbols and domain.
+
+    >>> from sympy.polys.domains import QQ
+    >>> from sympy.polys.puiseux import puiseux_ring
+    >>> R, x, y = puiseux_ring('x y', QQ)
+    >>> R
+    PuiseuxRing((x, y), QQ)
+    >>> p = 5*x**QQ(1,2) + 7/y
+    >>> p
+    7*y**(-1) + 5*x**(1/2)
+    """
+    ring = PuiseuxRing(symbols, domain)
+    return (ring,) + ring.gens # type: ignore
+
+
+class PuiseuxRing:
+    """Ring of Puiseux polynomials.
+
+    A Puiseux polynomial is a truncated Puiseux series. The exponents of the
+    monomials can be negative or rational numbers. This ring is used by the
+    ring_series module:
+
+    >>> from sympy.polys.domains import QQ
+    >>> from sympy.polys.puiseux import puiseux_ring
+    >>> from sympy.polys.ring_series import rs_exp, rs_nth_root
+    >>> ring, x, y = puiseux_ring('x y', QQ)
+    >>> f = x**2 + y**3
+    >>> f
+    y**3 + x**2
+    >>> f.diff(x)
+    2*x
+    >>> rs_exp(x, x, 5)
+    1 + x + 1/2*x**2 + 1/6*x**3 + 1/24*x**4
+
+    Importantly the Puiseux ring can represent truncated series with negative
+    and fractional exponents:
+
+    >>> f = 1/x + 1/y**2
+    >>> f
+    x**(-1) + y**(-2)
+    >>> f.diff(x)
+    -1*x**(-2)
+
+    >>> rs_nth_root(8*x + x**2 + x**3, 3, x, 5)
+    2*x**(1/3) + 1/12*x**(4/3) + 23/288*x**(7/3) + -139/20736*x**(10/3)
+
+    See Also
+    ========
+
+    sympy.polys.ring_series.rs_series
+    PuiseuxPoly
+    """
+    def __init__(self, symbols: str | list[Expr], domain: Domain):
+
+        poly_ring = PolyRing(symbols, domain)
+
+        domain = poly_ring.domain
+        ngens = poly_ring.ngens
+
+        self.poly_ring = poly_ring
+        self.domain = domain
+
+        self.symbols = poly_ring.symbols
+        self.gens = tuple([self.from_poly(g) for g in poly_ring.gens])
+        self.ngens = ngens
+
+        self.zero = self.from_poly(poly_ring.zero)
+        self.one = self.from_poly(poly_ring.one)
+
+        self.zero_monom = poly_ring.zero_monom # type: ignore
+        self.monomial_mul = poly_ring.monomial_mul # type: ignore
+
+    def __repr__(self) -> str:
+        return f"PuiseuxRing({self.symbols}, {self.domain})"
+
+    def __eq__(self, other: Any) -> bool:
+        if not isinstance(other, PuiseuxRing):
+            return NotImplemented
+        return self.symbols == other.symbols and self.domain == other.domain
+
+    def from_poly(self, poly: PolyElement) -> PuiseuxPoly:
+        """Create a Puiseux polynomial from a polynomial.
+
+        >>> from sympy.polys.domains import QQ
+        >>> from sympy.polys.rings import ring
+        >>> from sympy.polys.puiseux import puiseux_ring
+        >>> R1, x1 = ring('x', QQ)
+        >>> R2, x2 = puiseux_ring('x', QQ)
+        >>> R2.from_poly(x1**2)
+        x**2
+        """
+        return PuiseuxPoly(poly, self)
+
+    def from_dict(self, terms: dict[tuple[int, ...], Any]) -> PuiseuxPoly:
+        """Create a Puiseux polynomial from a dictionary of terms.
+
+        >>> from sympy.polys.domains import QQ
+        >>> from sympy.polys.puiseux import puiseux_ring
+        >>> R, x = puiseux_ring('x', QQ)
+        >>> R.from_dict({(QQ(1,2),): QQ(3)})
+        3*x**(1/2)
+        """
+        return PuiseuxPoly.from_dict(terms, self)
+
+    def from_int(self, n: int) -> PuiseuxPoly:
+        """Create a Puiseux polynomial from an integer.
+
+        >>> from sympy.polys.domains import QQ
+        >>> from sympy.polys.puiseux import puiseux_ring
+        >>> R, x = puiseux_ring('x', QQ)
+        >>> R.from_int(3)
+        3
+        """
+        return self.from_poly(self.poly_ring(n))
+
+    def domain_new(self, arg: Any) -> Any:
+        """Create a new element of the domain.
+
+        >>> from sympy.polys.domains import QQ
+        >>> from sympy.polys.puiseux import puiseux_ring
+        >>> R, x = puiseux_ring('x', QQ)
+        >>> R.domain_new(3)
+        3
+        >>> QQ.of_type(_)
+        True
+        """
+        return self.poly_ring.domain_new(arg)
+
+    def ground_new(self, arg: Any) -> PuiseuxPoly:
+        """Create a new element from a ground element.
+
+        >>> from sympy.polys.domains import QQ
+        >>> from sympy.polys.puiseux import puiseux_ring, PuiseuxPoly
+        >>> R, x = puiseux_ring('x', QQ)
+        >>> R.ground_new(3)
+        3
+        >>> isinstance(_, PuiseuxPoly)
+        True
+        """
+        return self.from_poly(self.poly_ring.ground_new(arg))
+
+    def __call__(self, arg: Any) -> PuiseuxPoly:
+        """Coerce an element into the ring.
+
+        >>> from sympy.polys.domains import QQ
+        >>> from sympy.polys.puiseux import puiseux_ring
+        >>> R, x = puiseux_ring('x', QQ)
+        >>> R(3)
+        3
+        >>> R({(QQ(1,2),): QQ(3)})
+        3*x**(1/2)
+        """
+        if isinstance(arg, dict):
+            return self.from_dict(arg)
+        else:
+            return self.from_poly(self.poly_ring(arg))
+
+    def index(self, x: PuiseuxPoly) -> int:
+        """Return the index of a generator.
+
+        >>> from sympy.polys.domains import QQ
+        >>> from sympy.polys.puiseux import puiseux_ring
+        >>> R, x, y = puiseux_ring('x y', QQ)
+        >>> R.index(x)
+        0
+        >>> R.index(y)
+        1
+        """
+        return self.gens.index(x)
+
+
+def _div_poly_monom(poly: PolyElement, monom: Iterable[int]) -> PolyElement:
+    ring = poly.ring
+    div = ring.monomial_div
+    return ring.from_dict({div(m, monom): c for m, c in poly.terms()})
+
+
+def _mul_poly_monom(poly: PolyElement, monom: Iterable[int]) -> PolyElement:
+    ring = poly.ring
+    mul = ring.monomial_mul
+    return ring.from_dict({mul(m, monom): c for m, c in poly.terms()})
+
+
+def _div_monom(monom: Iterable[int], div: Iterable[int]) -> tuple[int, ...]:
+    return tuple(mi - di for mi, di in zip(monom, div))
+
+
+class PuiseuxPoly:
+    """Puiseux polynomial. Represents a truncated Puiseux series.
+
+    See the :class:`PuiseuxRing` class for more information.
+
+    >>> from sympy import QQ
+    >>> from sympy.polys.puiseux import puiseux_ring
+    >>> R, x, y = puiseux_ring('x, y', QQ)
+    >>> p = 5*x**2 + 7*y**3
+    >>> p
+    7*y**3 + 5*x**2
+
+    The internal representation of a Puiseux polynomial wraps a normal
+    polynomial. To support negative powers the polynomial is considered to be
+    divided by a monomial.
+
+    >>> p2 = 1/x + 1/y**2
+    >>> p2.monom # x*y**2
+    (1, 2)
+    >>> p2.poly
+    x + y**2
+    >>> (y**2 + x) / (x*y**2) == p2
+    True
+
+    To support fractional powers the polynomial is considered to be a function
+    of ``x**(1/nx) * y**(1/ny) * ...``. The representation keeps track of a
+    monomial and a list of exponent denominators so that the polynomial can be
+    used to represent both negative and fractional powers.
+
+    >>> p3 = x**QQ(1,2) + y**QQ(2,3)
+    >>> p3.ns
+    (2, 3)
+    >>> p3.poly
+    x + y**2
+
+    See Also
+    ========
+
+    sympy.polys.puiseux.PuiseuxRing
+    sympy.polys.rings.PolyElement
+    """
+
+    ring: PuiseuxRing
+    poly: PolyElement
+    monom: tuple[int, ...] | None
+    ns: tuple[int, ...] | None
+
+    def __new__(cls, poly: PolyElement, ring: PuiseuxRing) -> PuiseuxPoly:
+        return cls._new(ring, poly, None, None)
+
+    @classmethod
+    def _new(
+        cls,
+        ring: PuiseuxRing,
+        poly: PolyElement,
+        monom: tuple[int, ...] | None,
+        ns: tuple[int, ...] | None,
+    ) -> PuiseuxPoly:
+        poly, monom, ns = cls._normalize(poly, monom, ns)
+        return cls._new_raw(ring, poly, monom, ns)
+
+    @classmethod
+    def _new_raw(
+        cls,
+        ring: PuiseuxRing,
+        poly: PolyElement,
+        monom: tuple[int, ...] | None,
+        ns: tuple[int, ...] | None,
+    ) -> PuiseuxPoly:
+        obj = object.__new__(cls)
+        obj.ring = ring
+        obj.poly = poly
+        obj.monom = monom
+        obj.ns = ns
+        return obj
+
+    def __eq__(self, other: Any) -> bool:
+        if isinstance(other, PuiseuxPoly):
+            return (
+                self.poly == other.poly
+                and self.monom == other.monom
+                and self.ns == other.ns
+            )
+        elif self.monom is None and self.ns is None:
+            return self.poly.__eq__(other)
+        else:
+            return NotImplemented
+
+    @classmethod
+    def _normalize(
+        cls,
+        poly: PolyElement,
+        monom: tuple[int, ...] | None,
+        ns: tuple[int, ...] | None,
+    ) -> tuple[PolyElement, tuple[int, ...] | None, tuple[int, ...] | None]:
+        if monom is None and ns is None:
+            return poly, None, None
+
+        if monom is not None:
+            degs = [max(d, 0) for d in poly.tail_degrees()]
+            if all(di >= mi for di, mi in zip(degs, monom)):
+                poly = _div_poly_monom(poly, monom)
+                monom = None
+            elif any(degs):
+                poly = _div_poly_monom(poly, degs)
+                monom = _div_monom(monom, degs)
+
+        if ns is not None:
+            factors_d, [poly_d] = poly.deflate()
+            degrees = poly.degrees()
+            monom_d = monom if monom is not None else [0] * len(degrees)
+            ns_new = []
+            monom_new = []
+            inflations = []
+            for fi, ni, di, mi in zip(factors_d, ns, degrees, monom_d):
+                if di == 0:
+                    g = gcd(ni, mi)
+                else:
+                    g = gcd(fi, ni, mi)
+                ns_new.append(ni // g)
+                monom_new.append(mi // g)
+                inflations.append(fi // g)
+
+            if any(infl > 1 for infl in inflations):
+                poly_d = poly_d.inflate(inflations)
+
+            poly = poly_d
+
+            if monom is not None:
+                monom = tuple(monom_new)
+
+            if all(n == 1 for n in ns_new):
+                ns = None
+            else:
+                ns = tuple(ns_new)
+
+        return poly, monom, ns
+
+    @classmethod
+    def _monom_fromint(
+        cls,
+        monom: tuple[int, ...],
+        dmonom: tuple[int, ...] | None,
+        ns: tuple[int, ...] | None,
+    ) -> tuple[Any, ...]:
+        if dmonom is not None and ns is not None:
+            return tuple(QQ(mi - di, ni) for mi, di, ni in zip(monom, dmonom, ns))
+        elif dmonom is not None:
+            return tuple(QQ(mi - di) for mi, di in zip(monom, dmonom))
+        elif ns is not None:
+            return tuple(QQ(mi, ni) for mi, ni in zip(monom, ns))
+        else:
+            return tuple(QQ(mi) for mi in monom)
+
+    @classmethod
+    def _monom_toint(
+        cls,
+        monom: tuple[Any, ...],
+        dmonom: tuple[int, ...] | None,
+        ns: tuple[int, ...] | None,
+    ) -> tuple[int, ...]:
+        if dmonom is not None and ns is not None:
+            return tuple(
+                int((mi * ni).numerator + di) for mi, di, ni in zip(monom, dmonom, ns)
+            )
+        elif dmonom is not None:
+            return tuple(int(mi.numerator + di) for mi, di in zip(monom, dmonom))
+        elif ns is not None:
+            return tuple(int((mi * ni).numerator) for mi, ni in zip(monom, ns))
+        else:
+            return tuple(int(mi.numerator) for mi in monom)
+
+    def itermonoms(self) -> Iterator[tuple[Any, ...]]:
+        """Iterate over the monomials of a Puiseux polynomial.
+
+        >>> from sympy import QQ
+        >>> from sympy.polys.puiseux import puiseux_ring
+        >>> R, x, y = puiseux_ring('x, y', QQ)
+        >>> p = 5*x**2 + 7*y**3
+        >>> list(p.itermonoms())
+        [(2, 0), (0, 3)]
+        >>> p[(2, 0)]
+        5
+        """
+        monom, ns = self.monom, self.ns
+        for m in self.poly.itermonoms():
+            yield self._monom_fromint(m, monom, ns)
+
+    def monoms(self) -> list[tuple[Any, ...]]:
+        """Return a list of the monomials of a Puiseux polynomial."""
+        return list(self.itermonoms())
+
+    def __iter__(self) -> Iterator[tuple[tuple[Any, ...], Any]]:
+        return self.itermonoms()
+
+    def __getitem__(self, monom: tuple[int, ...]) -> Any:
+        monom = self._monom_toint(monom, self.monom, self.ns)
+        return self.poly[monom]
+
+    def __len__(self) -> int:
+        return len(self.poly)
+
+    def iterterms(self) -> Iterator[tuple[tuple[Any, ...], Any]]:
+        """Iterate over the terms of a Puiseux polynomial.
+
+        >>> from sympy import QQ
+        >>> from sympy.polys.puiseux import puiseux_ring
+        >>> R, x, y = puiseux_ring('x, y', QQ)
+        >>> p = 5*x**2 + 7*y**3
+        >>> list(p.iterterms())
+        [((2, 0), 5), ((0, 3), 7)]
+        """
+        monom, ns = self.monom, self.ns
+        for m, coeff in self.poly.iterterms():
+            mq = self._monom_fromint(m, monom, ns)
+            yield mq, coeff
+
+    def terms(self) -> list[tuple[tuple[Any, ...], Any]]:
+        """Return a list of the terms of a Puiseux polynomial."""
+        return list(self.iterterms())
+
+    @property
+    def is_term(self) -> bool:
+        """Return True if the Puiseux polynomial is a single term."""
+        return self.poly.is_term
+
+    def to_dict(self) -> dict[tuple[int, ...], Any]:
+        """Return a dictionary representation of a Puiseux polynomial."""
+        return dict(self.iterterms())
+
+    @classmethod
+    def from_dict(
+        cls, terms: dict[tuple[Any, ...], Any], ring: PuiseuxRing
+    ) -> PuiseuxPoly:
+        """Create a Puiseux polynomial from a dictionary of terms.
+
+        >>> from sympy import QQ
+        >>> from sympy.polys.puiseux import puiseux_ring, PuiseuxPoly
+        >>> R, x = puiseux_ring('x', QQ)
+        >>> PuiseuxPoly.from_dict({(QQ(1,2),): QQ(3)}, R)
+        3*x**(1/2)
+        >>> R.from_dict({(QQ(1,2),): QQ(3)})
+        3*x**(1/2)
+        """
+        ns = [1] * ring.ngens
+        mon = [0] * ring.ngens
+        for mo in terms:
+            ns = [lcm(n, m.denominator) for n, m in zip(ns, mo)]
+            mon = [min(m, n) for m, n in zip(mo, mon)]
+
+        if not any(mon):
+            monom = None
+        else:
+            monom = tuple(-int((m * n).numerator) for m, n in zip(mon, ns))
+
+        if all(n == 1 for n in ns):
+            ns_final = None
+        else:
+            ns_final = tuple(ns)
+
+        terms_p = {cls._monom_toint(m, monom, ns_final): coeff for m, coeff in terms.items()}
+
+        poly = ring.poly_ring.from_dict(terms_p)
+
+        return cls._new(ring, poly, monom, ns_final)
+
+    def as_expr(self) -> Expr:
+        """Convert a Puiseux polynomial to :class:`~sympy.core.expr.Expr`.
+
+        >>> from sympy import QQ, Expr
+        >>> from sympy.polys.puiseux import puiseux_ring
+        >>> R, x = puiseux_ring('x', QQ)
+        >>> p = 5*x**2 + 7*x**3
+        >>> p.as_expr()
+        7*x**3 + 5*x**2
+        >>> isinstance(_, Expr)
+        True
+        """
+        ring = self.ring
+        dom = ring.domain
+        symbols = ring.symbols
+        terms = []
+        for monom, coeff in self.iterterms():
+            coeff_expr = dom.to_sympy(coeff)
+            monoms_expr = []
+            for i, m in enumerate(monom):
+                monoms_expr.append(symbols[i] ** m)
+            terms.append(Mul(coeff_expr, *monoms_expr))
+        return Add(*terms)
+
+    def __repr__(self) -> str:
+
+        def format_power(base: str, exp: int) -> str:
+            if exp == 1:
+                return base
+            elif exp >= 0 and int(exp) == exp:
+                return f"{base}**{exp}"
+            else:
+                return f"{base}**({exp})"
+
+        ring = self.ring
+        dom = ring.domain
+
+        syms = [str(s) for s in ring.symbols]
+        terms_str = []
+        for monom, coeff in sorted(self.terms()):
+            monom_str = "*".join(format_power(s, e) for s, e in zip(syms, monom) if e)
+            if coeff == dom.one:
+                if monom_str:
+                    terms_str.append(monom_str)
+                else:
+                    terms_str.append("1")
+            elif not monom_str:
+                terms_str.append(str(coeff))
+            else:
+                terms_str.append(f"{coeff}*{monom_str}")
+
+        return " + ".join(terms_str)
+
+    def _unify(
+        self, other: PuiseuxPoly
+    ) -> tuple[
+        PolyElement, PolyElement, tuple[int, ...] | None, tuple[int, ...] | None
+    ]:
+        """Bring two Puiseux polynomials to a common monom and ns."""
+        poly1, monom1, ns1 = self.poly, self.monom, self.ns
+        poly2, monom2, ns2 = other.poly, other.monom, other.ns
+
+        if monom1 == monom2 and ns1 == ns2:
+            return poly1, poly2, monom1, ns1
+
+        if ns1 == ns2:
+            ns = ns1
+        elif ns1 is not None and ns2 is not None:
+            ns = tuple(lcm(n1, n2) for n1, n2 in zip(ns1, ns2))
+            f1 = [n // n1 for n, n1 in zip(ns, ns1)]
+            f2 = [n // n2 for n, n2 in zip(ns, ns2)]
+            poly1 = poly1.inflate(f1)
+            poly2 = poly2.inflate(f2)
+            if monom1 is not None:
+                monom1 = tuple(m * f for m, f in zip(monom1, f1))
+            if monom2 is not None:
+                monom2 = tuple(m * f for m, f in zip(monom2, f2))
+        elif ns2 is not None:
+            ns = ns2
+            poly1 = poly1.inflate(ns)
+            if monom1 is not None:
+                monom1 = tuple(m * n for m, n in zip(monom1, ns))
+        elif ns1 is not None:
+            ns = ns1
+            poly2 = poly2.inflate(ns)
+            if monom2 is not None:
+                monom2 = tuple(m * n for m, n in zip(monom2, ns))
+        else:
+            assert False
+
+        if monom1 == monom2:
+            monom = monom1
+        elif monom1 is not None and monom2 is not None:
+            monom = tuple(max(m1, m2) for m1, m2 in zip(monom1, monom2))
+            poly1 = _mul_poly_monom(poly1, _div_monom(monom, monom1))
+            poly2 = _mul_poly_monom(poly2, _div_monom(monom, monom2))
+        elif monom2 is not None:
+            monom = monom2
+            poly1 = _mul_poly_monom(poly1, monom2)
+        elif monom1 is not None:
+            monom = monom1
+            poly2 = _mul_poly_monom(poly2, monom1)
+        else:
+            assert False
+
+        return poly1, poly2, monom, ns
+
+    def __pos__(self) -> PuiseuxPoly:
+        return self
+
+    def __neg__(self) -> PuiseuxPoly:
+        return self._new_raw(self.ring, -self.poly, self.monom, self.ns)
+
+    def __add__(self, other: Any) -> PuiseuxPoly:
+        if isinstance(other, PuiseuxPoly):
+            if self.ring != other.ring:
+                raise ValueError("Cannot add Puiseux polynomials from different rings")
+            return self._add(other)
+        domain = self.ring.domain
+        if isinstance(other, int):
+            return self._add_ground(domain.convert_from(QQ(other), QQ))
+        elif domain.of_type(other):
+            return self._add_ground(other)
+        else:
+            return NotImplemented
+
+    def __radd__(self, other: Any) -> PuiseuxPoly:
+        domain = self.ring.domain
+        if isinstance(other, int):
+            return self._add_ground(domain.convert_from(QQ(other), QQ))
+        elif domain.of_type(other):
+            return self._add_ground(other)
+        else:
+            return NotImplemented
+
+    def __sub__(self, other: Any) -> PuiseuxPoly:
+        if isinstance(other, PuiseuxPoly):
+            if self.ring != other.ring:
+                raise ValueError(
+                    "Cannot subtract Puiseux polynomials from different rings"
+                )
+            return self._sub(other)
+        domain = self.ring.domain
+        if isinstance(other, int):
+            return self._sub_ground(domain.convert_from(QQ(other), QQ))
+        elif domain.of_type(other):
+            return self._sub_ground(other)
+        else:
+            return NotImplemented
+
+    def __rsub__(self, other: Any) -> PuiseuxPoly:
+        domain = self.ring.domain
+        if isinstance(other, int):
+            return self._rsub_ground(domain.convert_from(QQ(other), QQ))
+        elif domain.of_type(other):
+            return self._rsub_ground(other)
+        else:
+            return NotImplemented
+
+    def __mul__(self, other: Any) -> PuiseuxPoly:
+        if isinstance(other, PuiseuxPoly):
+            if self.ring != other.ring:
+                raise ValueError(
+                    "Cannot multiply Puiseux polynomials from different rings"
+                )
+            return self._mul(other)
+        domain = self.ring.domain
+        if isinstance(other, int):
+            return self._mul_ground(domain.convert_from(QQ(other), QQ))
+        elif domain.of_type(other):
+            return self._mul_ground(other)
+        else:
+            return NotImplemented
+
+    def __rmul__(self, other: Any) -> PuiseuxPoly:
+        domain = self.ring.domain
+        if isinstance(other, int):
+            return self._mul_ground(domain.convert_from(QQ(other), QQ))
+        elif domain.of_type(other):
+            return self._mul_ground(other)
+        else:
+            return NotImplemented
+
+    def __pow__(self, other: Any) -> PuiseuxPoly:
+        if isinstance(other, int):
+            if other >= 0:
+                return self._pow_pint(other)
+            else:
+                return self._pow_nint(-other)
+        elif QQ.of_type(other):
+            return self._pow_rational(other)
+        else:
+            return NotImplemented
+
+    def __truediv__(self, other: Any) -> PuiseuxPoly:
+        if isinstance(other, PuiseuxPoly):
+            if self.ring != other.ring:
+                raise ValueError(
+                    "Cannot divide Puiseux polynomials from different rings"
+                )
+            return self._mul(other._inv())
+        domain = self.ring.domain
+        if isinstance(other, int):
+            return self._mul_ground(domain.convert_from(QQ(1, other), QQ))
+        elif domain.of_type(other):
+            return self._div_ground(other)
+        else:
+            return NotImplemented
+
+    def __rtruediv__(self, other: Any) -> PuiseuxPoly:
+        if isinstance(other, int):
+            return self._inv()._mul_ground(self.ring.domain.convert_from(QQ(other), QQ))
+        elif self.ring.domain.of_type(other):
+            return self._inv()._mul_ground(other)
+        else:
+            return NotImplemented
+
+    def _add(self, other: PuiseuxPoly) -> PuiseuxPoly:
+        poly1, poly2, monom, ns = self._unify(other)
+        return self._new(self.ring, poly1 + poly2, monom, ns)
+
+    def _add_ground(self, ground: Any) -> PuiseuxPoly:
+        return self._add(self.ring.ground_new(ground))
+
+    def _sub(self, other: PuiseuxPoly) -> PuiseuxPoly:
+        poly1, poly2, monom, ns = self._unify(other)
+        return self._new(self.ring, poly1 - poly2, monom, ns)
+
+    def _sub_ground(self, ground: Any) -> PuiseuxPoly:
+        return self._sub(self.ring.ground_new(ground))
+
+    def _rsub_ground(self, ground: Any) -> PuiseuxPoly:
+        return self.ring.ground_new(ground)._sub(self)
+
+    def _mul(self, other: PuiseuxPoly) -> PuiseuxPoly:
+        poly1, poly2, monom, ns = self._unify(other)
+        if monom is not None:
+            monom = tuple(2 * e for e in monom)
+        return self._new(self.ring, poly1 * poly2, monom, ns)
+
+    def _mul_ground(self, ground: Any) -> PuiseuxPoly:
+        return self._new_raw(self.ring, self.poly * ground, self.monom, self.ns)
+
+    def _div_ground(self, ground: Any) -> PuiseuxPoly:
+        return self._new_raw(self.ring, self.poly / ground, self.monom, self.ns)
+
+    def _pow_pint(self, n: int) -> PuiseuxPoly:
+        assert n >= 0
+        monom = self.monom
+        if monom is not None:
+            monom = tuple(m * n for m in monom)
+        return self._new(self.ring, self.poly**n, monom, self.ns)
+
+    def _pow_nint(self, n: int) -> PuiseuxPoly:
+        return self._inv()._pow_pint(n)
+
+    def _pow_rational(self, n: Any) -> PuiseuxPoly:
+        if not self.is_term:
+            raise ValueError("Only monomials can be raised to a rational power")
+        [(monom, coeff)] = self.terms()
+        domain = self.ring.domain
+        if not domain.is_one(coeff):
+            raise ValueError("Only monomials can be raised to a rational power")
+        monom = tuple(m * n for m in monom)
+        return self.ring.from_dict({monom: domain.one})
+
+    def _inv(self) -> PuiseuxPoly:
+        if not self.is_term:
+            raise ValueError("Only terms can be inverted")
+        [(monom, coeff)] = self.terms()
+        domain = self.ring.domain
+        if not domain.is_Field and not domain.is_one(coeff):
+            raise ValueError("Cannot invert non-unit coefficient")
+        monom = tuple(-m for m in monom)
+        coeff = 1 / coeff
+        return self.ring.from_dict({monom: coeff})
+
+    def diff(self, x: PuiseuxPoly) -> PuiseuxPoly:
+        """Differentiate a Puiseux polynomial with respect to a variable.
+
+        >>> from sympy import QQ
+        >>> from sympy.polys.puiseux import puiseux_ring
+        >>> R, x, y = puiseux_ring('x, y', QQ)
+        >>> p = 5*x**2 + 7*y**3
+        >>> p.diff(x)
+        10*x
+        >>> p.diff(y)
+        21*y**2
+        """
+        ring = self.ring
+        i = ring.index(x)
+        g = {}
+        for expv, coeff in self.iterterms():
+            n = expv[i]
+            if n:
+                e = list(expv)
+                e[i] -= 1
+                g[tuple(e)] = coeff * n
+        return ring(g)
diff --git a/sympy/polys/ring_series.py b/sympy/polys/ring_series.py
index d08b0c0507d1..4afcba37627c 100644
--- a/sympy/polys/ring_series.py
+++ b/sympy/polys/ring_series.py
@@ -43,6 +43,7 @@
 
 from sympy.polys.domains import QQ, EX
 from sympy.polys.rings import PolyElement, ring, sring
+from sympy.polys.puiseux import PuiseuxPoly
 from sympy.polys.polyerrors import DomainError
 from sympy.polys.monomials import (monomial_min, monomial_mul, monomial_div,
                                    monomial_ldiv)
@@ -89,7 +90,8 @@ def _invert_monoms(p1):
 
 def _giant_steps(target):
     """Return a list of precision steps for the Newton's method"""
-    res = giant_steps(2, target)
+    # We use ceil here because giant_steps cannot handle flint.fmpq
+    res = giant_steps(2, math.ceil(target))
     if res[0] != 2:
         res = [2] + res
     return res
@@ -113,13 +115,13 @@ def rs_trunc(p1, x, prec):
     x**5 + x + 1
     """
     R = p1.ring
-    p = R.zero
+    p = {}
     i = R.gens.index(x)
     for exp1 in p1:
         if exp1[i] >= prec:
             continue
         p[exp1] = p1[exp1]
-    return p
+    return R(p)
 
 def rs_is_puiseux(p, x):
     """
@@ -131,15 +133,15 @@ def rs_is_puiseux(p, x):
     ========
 
     >>> from sympy.polys.domains import QQ
-    >>> from sympy.polys.rings import ring
+    >>> from sympy.polys.puiseux import puiseux_ring
     >>> from sympy.polys.ring_series import rs_is_puiseux
-    >>> R, x = ring('x', QQ)
+    >>> R, x = puiseux_ring('x', QQ)
     >>> p = x**QQ(2,5) + x**QQ(2,3) + x
     >>> rs_is_puiseux(p, x)
     True
     """
     index = p.ring.gens.index(x)
-    for k in p:
+    for k in p.itermonoms():
         if k[index] != int(k[index]):
             return True
         if k[index] < 0:
@@ -156,12 +158,12 @@ def rs_puiseux(f, p, x, prec):
     ========
 
     >>> from sympy.polys.domains import QQ
-    >>> from sympy.polys.rings import ring
+    >>> from sympy.polys.puiseux import puiseux_ring
     >>> from sympy.polys.ring_series import rs_puiseux, rs_exp
-    >>> R, x = ring('x', QQ)
+    >>> R, x = puiseux_ring('x', QQ)
     >>> p = x**QQ(2,5) + x**QQ(2,3) + x
     >>> rs_puiseux(rs_exp,p, x, 1)
-    1/2*x**(4/5) + x**(2/3) + x**(2/5) + 1
+    1 + x**(2/5) + x**(2/3) + 1/2*x**(4/5)
     """
     index = p.ring.gens.index(x)
     n = 1
@@ -229,18 +231,18 @@ def rs_mul(p1, p2, x, prec):
     3*x**2 + 3*x + 1
     """
     R = p1.ring
-    p = R.zero
+    p = {}
     if R.__class__ != p2.ring.__class__ or R != p2.ring:
         raise ValueError('p1 and p2 must have the same ring')
     iv = R.gens.index(x)
-    if not isinstance(p2, PolyElement):
+    if not isinstance(p2, (PolyElement, PuiseuxPoly)):
         raise ValueError('p2 must be a polynomial')
     if R == p2.ring:
         get = p.get
-        items2 = list(p2.items())
+        items2 = p2.terms()
         items2.sort(key=lambda e: e[0][iv])
         if R.ngens == 1:
-            for exp1, v1 in p1.items():
+            for exp1, v1 in p1.iterterms():
                 for exp2, v2 in items2:
                     exp = exp1[0] + exp2[0]
                     if exp < prec:
@@ -250,7 +252,7 @@ def rs_mul(p1, p2, x, prec):
                         break
         else:
             monomial_mul = R.monomial_mul
-            for exp1, v1 in p1.items():
+            for exp1, v1 in p1.iterterms():
                 for exp2, v2 in items2:
                     if exp1[iv] + exp2[iv] < prec:
                         exp = monomial_mul(exp1, exp2)
@@ -258,8 +260,7 @@ def rs_mul(p1, p2, x, prec):
                     else:
                         break
 
-    p.strip_zero()
-    return p
+    return R(p)
 
 def rs_square(p1, x, prec):
     """
@@ -277,10 +278,10 @@ def rs_square(p1, x, prec):
     6*x**2 + 4*x + 1
     """
     R = p1.ring
-    p = R.zero
+    p = {}
     iv = R.gens.index(x)
     get = p.get
-    items = list(p1.items())
+    items = p1.terms()
     items.sort(key=lambda e: e[0][iv])
     monomial_mul = R.monomial_mul
     for i in range(len(items)):
@@ -292,14 +293,13 @@ def rs_square(p1, x, prec):
                 p[exp] = get(exp, 0) + v1*v2
             else:
                 break
-    p = p.imul_num(2)
+    p = {m: 2*v for m, v in p.items()}
     get = p.get
-    for expv, v in p1.items():
+    for expv, v in p1.iterterms():
         if 2*expv[iv] < prec:
             e2 = monomial_mul(expv, expv)
             p[e2] = get(e2, 0) + v**2
-    p.strip_zero()
-    return p
+    return R(p)
 
 def rs_pow(p1, n, x, prec):
     """
@@ -753,7 +753,7 @@ def rs_diff(p, x):
     """
     R = p.ring
     n = R.gens.index(x)
-    p1 = R.zero
+    p1 = {}
     mn = [0]*R.ngens
     mn[n] = 1
     mn = tuple(mn)
@@ -761,7 +761,7 @@ def rs_diff(p, x):
         if expv[n]:
             e = monomial_ldiv(expv, mn)
             p1[e] = R.domain_new(p[expv]*expv[n])
-    return p1
+    return R(p1)
 
 def rs_integrate(p, x):
     """
@@ -784,7 +784,7 @@ def rs_integrate(p, x):
     1/3*x**3*y**3 + 1/2*x**2
     """
     R = p.ring
-    p1 = R.zero
+    p1 = {}
     n = R.gens.index(x)
     mn = [0]*R.ngens
     mn[n] = 1
@@ -793,7 +793,7 @@ def rs_integrate(p, x):
     for expv in p:
         e = monomial_mul(expv, mn)
         p1[e] = R.domain_new(p[expv]/(expv[n] + 1))
-    return p1
+    return R(p1)
 
 def rs_fun(p, f, *args):
     r"""
@@ -858,31 +858,31 @@ def mul_xin(p, i, n):
     `x\_i` is the ith variable in ``p``.
     """
     R = p.ring
-    q = R(0)
-    for k, v in p.items():
+    q = {}
+    for k, v in p.terms():
         k1 = list(k)
         k1[i] += n
         q[tuple(k1)] = v
-    return q
+    return R(q)
 
 def pow_xin(p, i, n):
     """
     >>> from sympy.polys.domains import QQ
-    >>> from sympy.polys.rings import ring
+    >>> from sympy.polys.puiseux import puiseux_ring
     >>> from sympy.polys.ring_series import pow_xin
-    >>> R, x, y = ring('x, y', QQ)
+    >>> R, x, y = puiseux_ring('x, y', QQ)
     >>> p = x**QQ(2,5) + x + x**QQ(2,3)
     >>> index = p.ring.gens.index(x)
     >>> pow_xin(p, index, 15)
-    x**15 + x**10 + x**6
+    x**6 + x**10 + x**15
     """
     R = p.ring
-    q = R(0)
-    for k, v in p.items():
+    q = {}
+    for k, v in p.terms():
         k1 = list(k)
         k1[i] *= n
         q[tuple(k1)] = v
-    return q
+    return R(q)
 
 def _nth_root1(p, n, x, prec):
     """
@@ -973,7 +973,7 @@ def rs_nth_root(p, n, x, prec):
         c = p[zm]
         if R.domain is EX:
             c_expr = c.as_expr()
-            const = c_expr**QQ(1, n)
+            const = EX(c_expr**QQ(1, n))
         elif isinstance(c, PolyElement):
             try:
                 c_expr = c.as_expr()
@@ -991,7 +991,7 @@ def rs_nth_root(p, n, x, prec):
     else:
         res = _nth_root1(p, n, x, prec)
     if m:
-        m = QQ(m, n)
+        m = QQ(m) / n
         res = mul_xin(res, index, m)
     return res
 
@@ -1008,13 +1008,13 @@ def rs_log(p, x, prec):
     ========
 
     >>> from sympy.polys.domains import QQ
-    >>> from sympy.polys.rings import ring
+    >>> from sympy.polys.puiseux import puiseux_ring
     >>> from sympy.polys.ring_series import rs_log
-    >>> R, x = ring('x', QQ)
+    >>> R, x = puiseux_ring('x', QQ)
     >>> rs_log(1 + x, x, 8)
-    1/7*x**7 - 1/6*x**6 + 1/5*x**5 - 1/4*x**4 + 1/3*x**3 - 1/2*x**2 + x
+    x + -1/2*x**2 + 1/3*x**3 + -1/4*x**4 + 1/5*x**5 + -1/6*x**6 + 1/7*x**7
     >>> rs_log(x**QQ(3, 2) + 1, x, 5)
-    1/3*x**(9/2) - 1/2*x**3 + x**(3/2)
+    x**(3/2) + -1/2*x**3 + 1/3*x**(9/2)
     """
     if rs_is_puiseux(p, x):
         return rs_puiseux(rs_log, p, x, prec)
@@ -1030,7 +1030,7 @@ def rs_log(p, x, prec):
             c_expr = c.as_expr()
             if R.domain is EX:
                 const = log(c_expr)
-            elif isinstance(c, PolyElement):
+            elif isinstance(c, (PolyElement, PuiseuxPoly)):
                 try:
                     const = R(log(c_expr))
                 except ValueError:
@@ -1400,13 +1400,13 @@ def rs_sin(p, x, prec):
     ========
 
     >>> from sympy.polys.domains import QQ
-    >>> from sympy.polys.rings import ring
+    >>> from sympy.polys.puiseux import puiseux_ring
     >>> from sympy.polys.ring_series import rs_sin
-    >>> R, x, y = ring('x, y', QQ)
+    >>> R, x, y = puiseux_ring('x, y', QQ)
     >>> rs_sin(x + x*y, x, 4)
-    -1/6*x**3*y**3 - 1/2*x**3*y**2 - 1/2*x**3*y - 1/6*x**3 + x*y + x
+    x + x*y + -1/6*x**3 + -1/2*x**3*y + -1/2*x**3*y**2 + -1/6*x**3*y**3
     >>> rs_sin(x**QQ(3, 2) + x*y**QQ(7, 5), x, 4)
-    -1/2*x**(7/2)*y**(14/5) - 1/6*x**3*y**(21/5) + x**(3/2) + x*y**(7/5)
+    x*y**(7/5) + x**(3/2) + -1/6*x**3*y**(21/5) + -1/2*x**(7/2)*y**(14/5)
 
     See Also
     ========
@@ -1470,13 +1470,13 @@ def rs_cos(p, x, prec):
     ========
 
     >>> from sympy.polys.domains import QQ
-    >>> from sympy.polys.rings import ring
+    >>> from sympy.polys.puiseux import puiseux_ring
     >>> from sympy.polys.ring_series import rs_cos
-    >>> R, x, y = ring('x, y', QQ)
+    >>> R, x, y = puiseux_ring('x, y', QQ)
     >>> rs_cos(x + x*y, x, 4)
-    -1/2*x**2*y**2 - x**2*y - 1/2*x**2 + 1
+    1 + -1/2*x**2 + -1*x**2*y + -1/2*x**2*y**2
     >>> rs_cos(x + x*y, x, 4)/x**QQ(7, 5)
-    -1/2*x**(3/5)*y**2 - x**(3/5)*y - 1/2*x**(3/5) + x**(-7/5)
+    x**(-7/5) + -1/2*x**(3/5) + -1*x**(3/5)*y + -1/2*x**(3/5)*y**2
 
     See Also
     ========
@@ -1830,7 +1830,7 @@ def rs_compose_add(p1, p2):
     np2e = rs_hadamard_exp(np2)
     np3e = rs_mul(np1e, np2e, x, prec)
     np3 = rs_hadamard_exp(np3e, True)
-    np3a = (np3[(0,)] - np3)/x
+    np3a = (np3[(0,)] - np3) / x
     q = rs_integrate(np3a, x)
     q = rs_exp(q, x, prec)
     q = _invert_monoms(q)
@@ -1960,8 +1960,8 @@ def rs_series(expr, a, prec):
     Parameters
     ==========
 
-    expr : :class:`Expr`
-    a : :class:`Symbol` with respect to which expr is to be expanded
+    expr : :class:`~.Expr`
+    a : :class:`~.Symbol` with respect to which expr is to be expanded
     prec : order of the series expansion
 
     Currently supports multivariate Taylor series expansion. This is much
diff --git a/sympy/polys/rings.py b/sympy/polys/rings.py
index 9103b1737af1..2e902f30f809 100644
--- a/sympy/polys/rings.py
+++ b/sympy/polys/rings.py
@@ -1,12 +1,12 @@
 """Sparse polynomial rings. """
 
 from __future__ import annotations
-from typing import Any
 
 from operator import add, mul, lt, le, gt, ge
 from functools import reduce
 from types import GeneratorType
 
+from sympy.core.cache import cacheit
 from sympy.core.expr import Expr
 from sympy.core.intfunc import igcd
 from sympy.core.symbol import Symbol, symbols as _symbols
@@ -192,7 +192,6 @@ def _parse_symbols(symbols):
 
     raise GeneratorsError("expected a string, Symbol or expression or a non-empty sequence of strings, Symbols or expressions")
 
-_ring_cache: dict[Any, Any] = {}
 
 class PolyRing(DefaultPrinting, IPolys):
     """Multivariate distributed polynomial ring. """
@@ -210,61 +209,58 @@ def __new__(cls, symbols, domain, order=lex):
         order = OrderOpt.preprocess(order)
 
         _hash_tuple = (cls.__name__, symbols, ngens, domain, order)
-        obj = _ring_cache.get(_hash_tuple)
-
-        if obj is None:
-            if domain.is_Composite and set(symbols) & set(domain.symbols):
-                raise GeneratorsError("polynomial ring and it's ground domain share generators")
-
-            obj = object.__new__(cls)
-            obj._hash_tuple = _hash_tuple
-            obj._hash = hash(_hash_tuple)
-            obj.dtype = type("PolyElement", (PolyElement,), {"ring": obj})
-            obj.symbols = symbols
-            obj.ngens = ngens
-            obj.domain = domain
-            obj.order = order
-
-            obj.zero_monom = (0,)*ngens
-            obj.gens = obj._gens()
-            obj._gens_set = set(obj.gens)
-
-            obj._one = [(obj.zero_monom, domain.one)]
-
-            if ngens:
-                # These expect monomials in at least one variable
-                codegen = MonomialOps(ngens)
-                obj.monomial_mul = codegen.mul()
-                obj.monomial_pow = codegen.pow()
-                obj.monomial_mulpow = codegen.mulpow()
-                obj.monomial_ldiv = codegen.ldiv()
-                obj.monomial_div = codegen.div()
-                obj.monomial_lcm = codegen.lcm()
-                obj.monomial_gcd = codegen.gcd()
-            else:
-                monunit = lambda a, b: ()
-                obj.monomial_mul = monunit
-                obj.monomial_pow = monunit
-                obj.monomial_mulpow = lambda a, b, c: ()
-                obj.monomial_ldiv = monunit
-                obj.monomial_div = monunit
-                obj.monomial_lcm = monunit
-                obj.monomial_gcd = monunit
-
-
-            if order is lex:
-                obj.leading_expv = max
-            else:
-                obj.leading_expv = lambda f: max(f, key=order)
 
-            for symbol, generator in zip(obj.symbols, obj.gens):
-                if isinstance(symbol, Symbol):
-                    name = symbol.name
+        if domain.is_Composite and set(symbols) & set(domain.symbols):
+            raise GeneratorsError("polynomial ring and it's ground domain share generators")
+
+        obj = object.__new__(cls)
+        obj._hash_tuple = _hash_tuple
+        obj._hash = hash(_hash_tuple)
+        obj.symbols = symbols
+        obj.ngens = ngens
+        obj.domain = domain
+        obj.order = order
+
+        obj.dtype = PolyElement(obj, ()).new
+
+        obj.zero_monom = (0,)*ngens
+        obj.gens = obj._gens()
+        obj._gens_set = set(obj.gens)
+
+        obj._one = [(obj.zero_monom, domain.one)]
+
+        if ngens:
+            # These expect monomials in at least one variable
+            codegen = MonomialOps(ngens)
+            obj.monomial_mul = codegen.mul()
+            obj.monomial_pow = codegen.pow()
+            obj.monomial_mulpow = codegen.mulpow()
+            obj.monomial_ldiv = codegen.ldiv()
+            obj.monomial_div = codegen.div()
+            obj.monomial_lcm = codegen.lcm()
+            obj.monomial_gcd = codegen.gcd()
+        else:
+            monunit = lambda a, b: ()
+            obj.monomial_mul = monunit
+            obj.monomial_pow = monunit
+            obj.monomial_mulpow = lambda a, b, c: ()
+            obj.monomial_ldiv = monunit
+            obj.monomial_div = monunit
+            obj.monomial_lcm = monunit
+            obj.monomial_gcd = monunit
+
+
+        if order is lex:
+            obj.leading_expv = max
+        else:
+            obj.leading_expv = lambda f: max(f, key=order)
 
-                    if not hasattr(obj, name):
-                        setattr(obj, name, generator)
+        for symbol, generator in zip(obj.symbols, obj.gens):
+            if isinstance(symbol, Symbol):
+                name = symbol.name
 
-            _ring_cache[_hash_tuple] = obj
+                if not hasattr(obj, name):
+                    setattr(obj, name, generator)
 
         return obj
 
@@ -304,6 +300,13 @@ def __ne__(self, other):
         return not self == other
 
     def clone(self, symbols=None, domain=None, order=None):
+        # Need a hashable tuple for cacheit to work
+        if symbols is not None and isinstance(symbols, list):
+            symbols = tuple(symbols)
+        return self._clone(symbols, domain, order)
+
+    @cacheit
+    def _clone(self, symbols, domain, order):
         return self.__class__(symbols or self.symbols, domain or self.domain, order or self.order)
 
     def monomial_basis(self, i):
@@ -314,12 +317,16 @@ def monomial_basis(self, i):
 
     @property
     def zero(self):
-        return self.dtype()
+        return self.dtype([])
 
     @property
     def one(self):
         return self.dtype(self._one)
 
+    def is_element(self, element):
+        """True if ``element`` is an element of this ring. False otherwise. """
+        return isinstance(element, PolyElement) and element.ring == self
+
     def domain_new(self, element, orig_domain=None):
         return self.domain.convert(element, orig_domain)
 
@@ -423,7 +430,7 @@ def index(self, gen):
                 i = -i - 1
             else:
                 raise ValueError("invalid generator index: %s" % gen)
-        elif isinstance(gen, self.dtype):
+        elif self.is_element(gen):
             try:
                 i = self.gens.index(gen)
             except ValueError:
@@ -579,8 +586,24 @@ def symmetric_poly(self, n):
 class PolyElement(DomainElement, DefaultPrinting, CantSympify, dict):
     """Element of multivariate distributed polynomial ring. """
 
+    def __init__(self, ring, init):
+        super().__init__(init)
+        self.ring = ring
+        # This check would be too slow to run every time:
+        # self._check()
+
+    def _check(self):
+        assert isinstance(self, PolyElement)
+        assert isinstance(self.ring, PolyRing)
+        dom = self.ring.domain
+        assert isinstance(dom, Domain)
+        for monom, coeff in self.items():
+            assert dom.of_type(coeff)
+            assert len(monom) == self.ring.ngens
+            assert all(isinstance(exp, int) and exp >= 0 for exp in monom)
+
     def new(self, init):
-        return self.__class__(init)
+        return self.__class__(self.ring, init)
 
     def parent(self):
         return self.ring.to_domain()
@@ -695,7 +718,7 @@ def __eq__(p1, p2):
         """
         if not p2:
             return not p1
-        elif isinstance(p2, PolyElement) and p2.ring == p1.ring:
+        elif p1.ring.is_element(p2):
             return dict.__eq__(p1, p2)
         elif len(p1) > 1:
             return False
@@ -709,7 +732,7 @@ def almosteq(p1, p2, tolerance=None):
         """Approximate equality test for polynomials. """
         ring = p1.ring
 
-        if isinstance(p2, ring.dtype):
+        if ring.is_element(p2):
             if set(p1.keys()) != set(p2.keys()):
                 return False
 
@@ -733,7 +756,7 @@ def sort_key(self):
         return (len(self), self.terms())
 
     def _cmp(p1, p2, op):
-        if isinstance(p2, p1.ring.dtype):
+        if p1.ring.is_element(p2):
             return op(p1.sort_key(), p2.sort_key())
         else:
             return NotImplemented
@@ -956,7 +979,7 @@ def __add__(p1, p2):
         if not p2:
             return p1.copy()
         ring = p1.ring
-        if isinstance(p2, ring.dtype):
+        if ring.is_element(p2):
             p = p1.copy()
             get = p.get
             zero = ring.domain.zero
@@ -1032,7 +1055,7 @@ def __sub__(p1, p2):
         if not p2:
             return p1.copy()
         ring = p1.ring
-        if isinstance(p2, ring.dtype):
+        if ring.is_element(p2):
             p = p1.copy()
             get = p.get
             zero = ring.domain.zero
@@ -1092,6 +1115,7 @@ def __rsub__(p1, n):
             for expv in p1:
                 p[expv] = -p1[expv]
             p += n
+            # p._check()
             return p
 
     def __mul__(p1, p2):
@@ -1114,7 +1138,7 @@ def __mul__(p1, p2):
         p = ring.zero
         if not p1 or not p2:
             return p
-        elif isinstance(p2, ring.dtype):
+        elif ring.is_element(p2):
             get = p.get
             zero = ring.domain.zero
             monomial_mul = ring.monomial_mul
@@ -1124,6 +1148,7 @@ def __mul__(p1, p2):
                     exp = monomial_mul(exp1, exp2)
                     p[exp] = get(exp, zero) + v1*v2
             p.strip_zero()
+            # p._check()
             return p
         elif isinstance(p2, PolyElement):
             if isinstance(ring.domain, PolynomialRing) and ring.domain.ring == p2.ring:
@@ -1142,6 +1167,7 @@ def __mul__(p1, p2):
                 v = v1*p2
                 if v:
                     p[exp1] = v
+            # p._check()
             return p
 
     def __rmul__(p1, p2):
@@ -1188,6 +1214,11 @@ def __pow__(self, n):
         x**3 + 3*x**2*y**2 + 3*x*y**4 + y**6
 
         """
+        if not isinstance(n, int):
+            raise TypeError("exponent must be an integer, got %s" % n)
+        elif n < 0:
+            raise ValueError("exponent must be a non-negative integer, got %s" % n)
+
         ring = self.ring
 
         if not n:
@@ -1202,6 +1233,7 @@ def __pow__(self, n):
                 p[ring.monomial_pow(monom, n)] = coeff
             else:
                 p[ring.monomial_pow(monom, n)] = coeff**n
+            # p._check()
             return p
 
         # For ring series, we need negative and rational exponent support only
@@ -1300,6 +1332,7 @@ def square(self):
             k2 = monomial_mul(k, k)
             p[k2] = get(k2, zero) + v**2
         p.strip_zero()
+        # p._check()
         return p
 
     def __divmod__(p1, p2):
@@ -1307,7 +1340,7 @@ def __divmod__(p1, p2):
 
         if not p2:
             raise ZeroDivisionError("polynomial division")
-        elif isinstance(p2, ring.dtype):
+        elif ring.is_element(p2):
             return p1.div(p2)
         elif isinstance(p2, PolyElement):
             if isinstance(ring.domain, PolynomialRing) and ring.domain.ring == p2.ring:
@@ -1325,14 +1358,20 @@ def __divmod__(p1, p2):
             return (p1.quo_ground(p2), p1.rem_ground(p2))
 
     def __rdivmod__(p1, p2):
-        return NotImplemented
+        ring = p1.ring
+        try:
+            p2 = ring.ground_new(p2)
+        except CoercionFailed:
+            return NotImplemented
+        else:
+            return p2.div(p1)
 
     def __mod__(p1, p2):
         ring = p1.ring
 
         if not p2:
             raise ZeroDivisionError("polynomial division")
-        elif isinstance(p2, ring.dtype):
+        elif ring.is_element(p2):
             return p1.rem(p2)
         elif isinstance(p2, PolyElement):
             if isinstance(ring.domain, PolynomialRing) and ring.domain.ring == p2.ring:
@@ -1350,18 +1389,21 @@ def __mod__(p1, p2):
             return p1.rem_ground(p2)
 
     def __rmod__(p1, p2):
-        return NotImplemented
+        ring = p1.ring
+        try:
+            p2 = ring.ground_new(p2)
+        except CoercionFailed:
+            return NotImplemented
+        else:
+            return p2.rem(p1)
 
-    def __truediv__(p1, p2):
+    def __floordiv__(p1, p2):
         ring = p1.ring
 
         if not p2:
             raise ZeroDivisionError("polynomial division")
-        elif isinstance(p2, ring.dtype):
-            if p2.is_monomial:
-                return p1*(p2**(-1))
-            else:
-                return p1.quo(p2)
+        elif ring.is_element(p2):
+            return p1.quo(p2)
         elif isinstance(p2, PolyElement):
             if isinstance(ring.domain, PolynomialRing) and ring.domain.ring == p2.ring:
                 pass
@@ -1377,13 +1419,45 @@ def __truediv__(p1, p2):
         else:
             return p1.quo_ground(p2)
 
-    def __rtruediv__(p1, p2):
-        return NotImplemented
+    def __rfloordiv__(p1, p2):
+        ring = p1.ring
+        try:
+            p2 = ring.ground_new(p2)
+        except CoercionFailed:
+            return NotImplemented
+        else:
+            return p2.quo(p1)
 
-    __floordiv__ = __truediv__
-    __rfloordiv__ = __rtruediv__
+    def __truediv__(p1, p2):
+        ring = p1.ring
+
+        if not p2:
+            raise ZeroDivisionError("polynomial division")
+        elif ring.is_element(p2):
+            return p1.exquo(p2)
+        elif isinstance(p2, PolyElement):
+            if isinstance(ring.domain, PolynomialRing) and ring.domain.ring == p2.ring:
+                pass
+            elif isinstance(p2.ring.domain, PolynomialRing) and p2.ring.domain.ring == ring:
+                return p2.__rtruediv__(p1)
+            else:
+                return NotImplemented
 
-    # TODO: use // (__floordiv__) for exquo()?
+        try:
+            p2 = ring.domain_new(p2)
+        except CoercionFailed:
+            return NotImplemented
+        else:
+            return p1.quo_ground(p2)
+
+    def __rtruediv__(p1, p2):
+        ring = p1.ring
+        try:
+            p2 = ring.ground_new(p2)
+        except CoercionFailed:
+            return NotImplemented
+        else:
+            return p2.exquo(p1)
 
     def _term_div(self):
         zm = self.ring.zero_monom
@@ -1738,7 +1812,7 @@ def coeff(self, element):
         """
         if element == 1:
             return self._get_coeff(self.ring.zero_monom)
-        elif isinstance(element, self.ring.dtype):
+        elif self.ring.is_element(element):
             terms = list(element.iterterms())
             if len(terms) == 1:
                 monom, coeff = terms[0]
diff --git a/sympy/polys/tests/test_fields.py b/sympy/polys/tests/test_fields.py
index da9f39101599..4f85a00d75dc 100644
--- a/sympy/polys/tests/test_fields.py
+++ b/sympy/polys/tests/test_fields.py
@@ -29,19 +29,10 @@ def test_FracField___hash__():
 
 def test_FracField___eq__():
     assert field("x,y,z", QQ)[0] == field("x,y,z", QQ)[0]
-    assert field("x,y,z", QQ)[0] is field("x,y,z", QQ)[0]
-
     assert field("x,y,z", QQ)[0] != field("x,y,z", ZZ)[0]
-    assert field("x,y,z", QQ)[0] is not field("x,y,z", ZZ)[0]
-
     assert field("x,y,z", ZZ)[0] != field("x,y,z", QQ)[0]
-    assert field("x,y,z", ZZ)[0] is not field("x,y,z", QQ)[0]
-
     assert field("x,y,z", QQ)[0] != field("x,y", QQ)[0]
-    assert field("x,y,z", QQ)[0] is not field("x,y", QQ)[0]
-
     assert field("x,y", QQ)[0] != field("x,y,z", QQ)[0]
-    assert field("x,y", QQ)[0] is not field("x,y,z", QQ)[0]
 
 def test_sfield():
     x = symbols("x")
@@ -99,34 +90,34 @@ def test_FracElement_from_expr():
     F, X, Y, Z = field((x, y, z), ZZ)
 
     f = F.from_expr(1)
-    assert f == 1 and isinstance(f, F.dtype)
+    assert f == 1 and F.is_element(f)
 
     f = F.from_expr(Rational(3, 7))
-    assert f == F(3)/7 and isinstance(f, F.dtype)
+    assert f == F(3)/7 and F.is_element(f)
 
     f = F.from_expr(x)
-    assert f == X and isinstance(f, F.dtype)
+    assert f == X and F.is_element(f)
 
     f = F.from_expr(Rational(3,7)*x)
-    assert f == X*Rational(3, 7) and isinstance(f, F.dtype)
+    assert f == X*Rational(3, 7) and F.is_element(f)
 
     f = F.from_expr(1/x)
-    assert f == 1/X and isinstance(f, F.dtype)
+    assert f == 1/X and F.is_element(f)
 
     f = F.from_expr(x*y*z)
-    assert f == X*Y*Z and isinstance(f, F.dtype)
+    assert f == X*Y*Z and F.is_element(f)
 
     f = F.from_expr(x*y/z)
-    assert f == X*Y/Z and isinstance(f, F.dtype)
+    assert f == X*Y/Z and F.is_element(f)
 
     f = F.from_expr(x*y*z + x*y + x)
-    assert f == X*Y*Z + X*Y + X and isinstance(f, F.dtype)
+    assert f == X*Y*Z + X*Y + X and F.is_element(f)
 
     f = F.from_expr((x*y*z + x*y + x)/(x*y + 7))
-    assert f == (X*Y*Z + X*Y + X)/(X*Y + 7) and isinstance(f, F.dtype)
+    assert f == (X*Y*Z + X*Y + X)/(X*Y + 7) and F.is_element(f)
 
     f = F.from_expr(x**3*y*z + x**2*y**7 + 1)
-    assert f == X**3*Y*Z + X**2*Y**7 + 1 and isinstance(f, F.dtype)
+    assert f == X**3*Y*Z + X**2*Y**7 + 1 and F.is_element(f)
 
     raises(ValueError, lambda: F.from_expr(2**x))
     raises(ValueError, lambda: F.from_expr(7*x + sqrt(2)))
diff --git a/sympy/polys/tests/test_puiseux.py b/sympy/polys/tests/test_puiseux.py
new file mode 100644
index 000000000000..031881e9d12c
--- /dev/null
+++ b/sympy/polys/tests/test_puiseux.py
@@ -0,0 +1,204 @@
+#
+# Tests for PuiseuxRing and PuiseuxPoly
+#
+
+from sympy.testing.pytest import raises
+
+from sympy import ZZ, QQ, ring
+from sympy.polys.puiseux import PuiseuxRing, PuiseuxPoly, puiseux_ring
+
+from sympy.abc import x, y
+
+
+def test_puiseux_ring():
+    R, px = puiseux_ring('x', QQ)
+    R2, px2 = puiseux_ring([x], QQ)
+    assert isinstance(R, PuiseuxRing)
+    assert isinstance(px, PuiseuxPoly)
+    assert R == R2
+    assert px == px2
+    assert R == PuiseuxRing('x', QQ)
+    assert R == PuiseuxRing([x], QQ)
+    assert R != PuiseuxRing('y', QQ)
+    assert R != PuiseuxRing('x', ZZ)
+    assert R != PuiseuxRing('x, y', QQ)
+    assert R != QQ
+    assert str(R) == 'PuiseuxRing((x,), QQ)'
+
+
+def test_puiseux_ring_attributes():
+    R1, px1, py1 = ring('x, y', QQ)
+    R2, px2, py2 = puiseux_ring('x, y', QQ)
+    assert R2.domain == QQ
+    assert R2.symbols == (x, y)
+    assert R2.gens == (px2, py2)
+    assert R2.ngens == 2
+    assert R2.poly_ring == R1
+    assert R2.zero == PuiseuxPoly(R1.zero, R2)
+    assert R2.one == PuiseuxPoly(R1.one, R2)
+    assert R2.zero_monom == R1.zero_monom == (0, 0) # type: ignore
+    assert R2.monomial_mul((1, 2), (3, 4)) == (4, 6)
+
+
+def test_puiseux_ring_methods():
+    R1, px1, py1 = ring('x, y', QQ)
+    R2, px2, py2 = puiseux_ring('x, y', QQ)
+    assert R2({(1, 2): 3}) == 3*px2*py2**2
+    assert R2(px1) == px2
+    assert R2(1) == R2.one
+    assert R2(QQ(1,2)) == QQ(1,2)*R2.one
+    assert R2.from_poly(px1) == px2
+    assert R2.from_poly(px1) != py2
+    assert R2.from_dict({(1, 2): QQ(3)}) == 3*px2*py2**2
+    assert R2.from_dict({(QQ(1,2), 2): QQ(3)}) == 3*px2**QQ(1,2)*py2**2
+    assert R2.from_int(3) == 3*R2.one
+    assert R2.domain_new(3) == QQ(3)
+    assert QQ.of_type(R2.domain_new(3))
+    assert R2.ground_new(3) == 3*R2.one
+    assert isinstance(R2.ground_new(3), PuiseuxPoly)
+    assert R2.index(px2) == 0
+    assert R2.index(py2) == 1
+
+
+def test_puiseux_poly():
+    R1, px1 = ring('x', QQ)
+    R2, px2 = puiseux_ring('x', QQ)
+    assert PuiseuxPoly(px1, R2) == px2
+    assert px2.ring == R2
+    assert px2.as_expr() == px1.as_expr() == x
+    assert px1 != px2
+    assert R2.one == px2**0 == 1
+    assert px2 == px1
+    assert px2 != 2.0
+    assert px2**QQ(1,2) != px1
+
+
+def test_puiseux_poly_normalization():
+    R, x = puiseux_ring('x', QQ)
+    assert (x**2 + 1) / x == x + 1/x == R({(1,): 1, (-1,): 1})
+    assert (x**QQ(1,6))**2 == x**QQ(1,3) == R({(QQ(1,3),): 1})
+    assert (x**QQ(1,6))**(-2) == x**(-QQ(1,3)) == R({(-QQ(1,3),): 1})
+    assert (x**QQ(1,6))**QQ(1,2) == x**QQ(1,12) == R({(QQ(1,12),): 1})
+    assert (x**QQ(1,6))**6 == x == R({(1,): 1})
+    assert x**QQ(1,6) * x**QQ(1,3) == x**QQ(1,2) == R({(QQ(1,2),): 1})
+    assert 1/x * x**2 == x == R({(1,): 1})
+    assert 1/x**QQ(1,3) * x**QQ(1,3) == 1 == R({(0,): 1})
+
+
+def test_puiseux_poly_monoms():
+    R, x = puiseux_ring('x', QQ)
+    assert x.monoms() == [(1,)]
+    assert list(x) == [(1,)]
+    assert (x**2 + 1).monoms() == [(2,), (0,)]
+    assert R({(1,): 1, (-1,): 1}).monoms() == [(1,), (-1,)]
+    assert R({(QQ(1,3),): 1}).monoms() == [(QQ(1,3),)]
+    assert R({(-QQ(1,3),): 1}).monoms() == [(-QQ(1,3),)]
+    p = x**QQ(1,6)
+    assert p[(QQ(1,6),)] == 1
+    raises(KeyError, lambda: p[(1,)])
+    assert p.to_dict() == {(QQ(1,6),): 1}
+    assert R(p.to_dict()) == p
+    assert PuiseuxPoly.from_dict({(QQ(1,6),): 1}, R) == p
+
+
+def test_puiseux_poly_repr():
+    R, x = puiseux_ring('x', QQ)
+    assert repr(x) == 'x'
+    assert repr(x**QQ(1,2)) == 'x**(1/2)'
+    assert repr(1/x) == 'x**(-1)'
+    assert repr(2*x**2 + 1) == '1 + 2*x**2'
+    assert repr(R.one) == '1'
+    assert repr(2*R.one) == '2'
+
+
+def test_puiseux_poly_unify():
+    R, x = puiseux_ring('x', QQ)
+    assert 1/x + x == x + 1/x == R({(1,): 1, (-1,): 1})
+    assert repr(1/x + x) == 'x**(-1) + x'
+    assert 1/x + 1/x == 2/x == R({(-1,): 2})
+    assert repr(1/x + 1/x) == '2*x**(-1)'
+    assert x**QQ(1,2) + x**QQ(1,2) == 2*x**QQ(1,2) == R({(QQ(1,2),): 2})
+    assert repr(x**QQ(1,2) + x**QQ(1,2)) == '2*x**(1/2)'
+    assert x**QQ(1,2) + x**QQ(1,3) == R({(QQ(1,2),): 1, (QQ(1,3),): 1})
+    assert repr(x**QQ(1,2) + x**QQ(1,3)) == 'x**(1/3) + x**(1/2)'
+    assert x + x**QQ(1,2) == R({(1,): 1, (QQ(1,2),): 1})
+    assert repr(x + x**QQ(1,2)) == 'x**(1/2) + x'
+    assert 1/x**QQ(1,2) + 1/x**QQ(1,3) == R({(-QQ(1,2),): 1, (-QQ(1,3),): 1})
+    assert repr(1/x**QQ(1,2) + 1/x**QQ(1,3)) == 'x**(-1/2) + x**(-1/3)'
+    assert 1/x + x**QQ(1,2) == x**QQ(1,2) + 1/x == R({(-1,): 1, (QQ(1,2),): 1})
+    assert repr(1/x + x**QQ(1,2)) == 'x**(-1) + x**(1/2)'
+
+
+def test_puiseux_poly_arit():
+    R, x = puiseux_ring('x', QQ)
+    R2, y = puiseux_ring('y', QQ)
+    p = x**2 + 1
+    assert +p == p
+    assert -p == -1 - x**2
+    assert p + p == 2*p == 2*x**2 + 2
+    assert p + 1 == 1 + p == x**2 + 2
+    assert p + QQ(1,2) == QQ(1,2) + p == x**2 + QQ(3,2)
+    assert p - p == 0
+    assert p - 1 == -1 + p == x**2
+    assert p - QQ(1,2) == -QQ(1,2) + p == x**2 + QQ(1,2)
+    assert 1 - p == -p + 1 == -x**2
+    assert QQ(1,2) - p == -p + QQ(1,2) == -x**2 - QQ(1,2)
+    assert p * p == x**4 + 2*x**2 + 1
+    assert p * 1 == 1 * p == p
+    assert 2 * p == p * 2 == 2*x**2 + 2
+    assert p * QQ(1,2) == QQ(1,2) * p == QQ(1,2)*x**2 + QQ(1,2)
+    assert x**QQ(1,2) * x**QQ(1,2) == x
+    raises(ValueError, lambda: x + y)
+    raises(ValueError, lambda: x - y)
+    raises(ValueError, lambda: x * y)
+    raises(TypeError, lambda: x + None)
+    raises(TypeError, lambda: x - None)
+    raises(TypeError, lambda: x * None)
+    raises(TypeError, lambda: None + x)
+    raises(TypeError, lambda: None - x)
+    raises(TypeError, lambda: None * x)
+
+
+def test_puiseux_poly_div():
+    R, x = puiseux_ring('x', QQ)
+    R2, y = puiseux_ring('y', QQ)
+    p = x**2 - 1
+    assert p / 1 == p
+    assert p / QQ(1,2) == 2*p == 2*x**2 - 2
+    assert p / x == x - 1/x == R({(1,): 1, (-1,): -1})
+    assert 2 / x == 2*x**-1 == R({(-1,): 2})
+    assert QQ(1,2) / x == QQ(1,2)*x**-1 == 1/(2*x) == 1/x/2 == R({(-1,): QQ(1,2)})
+    raises(ZeroDivisionError, lambda: p / 0)
+    raises(ValueError, lambda: (x + 1) / (x + 2))
+    raises(ValueError, lambda: (x + 1) / (x + 1))
+    raises(ValueError, lambda: x / y)
+    raises(TypeError, lambda: x / None)
+    raises(TypeError, lambda: None / x)
+
+
+def test_puiseux_poly_pow():
+    R, x = puiseux_ring('x', QQ)
+    Rz, xz = puiseux_ring('x', ZZ)
+    assert x**0 == 1 == R({(0,): 1})
+    assert x**1 == x == R({(1,): 1})
+    assert x**2 == x*x == R({(2,): 1})
+    assert x**QQ(1,2) == R({(QQ(1,2),): 1})
+    assert x**-1 == 1/x == R({(-1,): 1})
+    assert x**-QQ(1,2) == 1/x**QQ(1,2) == R({(-QQ(1,2),): 1})
+    assert (2*x)**-1 == 1/(2*x) == QQ(1,2)/x == QQ(1,2)*x**-1 == R({(-1,): QQ(1,2)})
+    assert 2/x**2 == 2*x**-2 == R({(-2,): 2})
+    assert 2/xz**2 == 2*xz**-2 == Rz({(-2,): 2})
+    raises(TypeError, lambda: x**None)
+    raises(ValueError, lambda: (x + 1)**-1)
+    raises(ValueError, lambda: (x + 1)**QQ(1,2))
+    raises(ValueError, lambda: (2*x)**QQ(1,2))
+    raises(ValueError, lambda: (2*xz)**-1)
+
+
+def test_puiseux_poly_diff():
+    R, x, y = puiseux_ring('x, y', QQ)
+    assert (x**2 + 1).diff(x) == 2*x
+    assert (x**2 + 1).diff(y) == 0
+    assert (x**2 + y**2).diff(x) == 2*x
+    assert (x**QQ(1,2) + y**QQ(1,2)).diff(x) == QQ(1,2)*x**-QQ(1,2)
+    assert ((x*y)**QQ(1,2)).diff(x) == QQ(1,2)*y**QQ(1,2)*x**-QQ(1,2)
diff --git a/sympy/polys/tests/test_ring_series.py b/sympy/polys/tests/test_ring_series.py
index 0f70c05d3888..b19156fbaceb 100644
--- a/sympy/polys/tests/test_ring_series.py
+++ b/sympy/polys/tests/test_ring_series.py
@@ -1,5 +1,6 @@
 from sympy.polys.domains import QQ, EX, RR
 from sympy.polys.rings import ring
+from sympy.polys.puiseux import puiseux_ring
 from sympy.polys.ring_series import (_invert_monoms, rs_integrate,
     rs_trunc, rs_mul, rs_square, rs_pow, _has_constant_term, rs_hadamard_exp,
     rs_series_from_list, rs_exp, rs_log, rs_newton, rs_series_inversion,
@@ -141,11 +142,11 @@ def test_series_from_list():
         p2 += cx*rs_pow(p, i, x, h)
     assert p1 == p2
 
+
 def test_log():
     R, x = ring('x', QQ)
     p = 1 + x
-    p1 = rs_log(p, x, 4)/x**2
-    assert p1 == Rational(1, 3)*x - S.Half + x**(-1)
+    assert rs_log(p, x, 4) == x - x**2/2 + x**3/3
     p = 1 + x +2*x**2/3
     p1 = rs_log(p, x, 9)
     assert p1 == -17*x**8/648 + 13*x**7/189 - 11*x**6/162 - x**5/45 + \
@@ -172,6 +173,7 @@ def test_log():
     p = x + x**2 + 3
     assert rs_log(p, x, 10).compose(x, 5) == EX(log(3) + Rational(19281291595, 9920232))
 
+
 def test_exp():
     R, x = ring('x', QQ)
     p = x + x**4
@@ -222,8 +224,9 @@ def test_fun():
     assert rs_fun(p, rs_tan, x, 10) == rs_tan(p, x, 10)
     assert rs_fun(p, _tan1, x, 10) == _tan1(p, x, 10)
 
+
 def test_nth_root():
-    R, x, y = ring('x, y', QQ)
+    R, x, y = puiseux_ring('x, y', QQ)
     assert rs_nth_root(1 + x**2*y, 4, x, 10) == -77*x**8*y**4/2048 + \
         7*x**6*y**3/128 - 3*x**4*y**2/32 + x**2*y/4 + 1
     assert rs_nth_root(1 + x*y + x**2*y**3, 3, x, 5) == -x**4*y**6/9 + \
@@ -236,14 +239,15 @@ def test_nth_root():
 
     # Constant term in series
     a = symbols('a')
-    R, x, y = ring('x, y', EX)
-    assert rs_nth_root(x + a, 3, x, 4) == EX(5/(81*a**QQ(8, 3)))*x**3 - \
+    R, x, y = puiseux_ring('x, y', EX)
+    assert rs_nth_root(x + EX(a), 3, x, 4) == EX(5/(81*a**QQ(8, 3)))*x**3 - \
         EX(1/(9*a**QQ(5, 3)))*x**2 + EX(1/(3*a**QQ(2, 3)))*x + EX(a**QQ(1, 3))
     assert rs_nth_root(x**QQ(2, 3) + x**2*y + 5, 2, x, 3) == -EX(sqrt(5)/100)*\
         x**QQ(8, 3)*y - EX(sqrt(5)/16000)*x**QQ(8, 3) + EX(sqrt(5)/10)*x**2*y + \
         EX(sqrt(5)/2000)*x**2 - EX(sqrt(5)/200)*x**QQ(4, 3) + \
         EX(sqrt(5)/10)*x**QQ(2, 3) + EX(sqrt(5))
 
+
 def test_atan():
     R, x, y = ring('x, y', QQ)
     assert rs_atan(x, x, 9) == -x**7/7 + x**5/5 - x**3/3 + x
@@ -272,8 +276,7 @@ def test_asin():
 
 def test_tan():
     R, x, y = ring('x, y', QQ)
-    assert rs_tan(x, x, 9)/x**5 == \
-        Rational(17, 315)*x**2 + Rational(2, 15) + Rational(1, 3)*x**(-2) + x**(-4)
+    assert rs_tan(x, x, 9) == x + x**3/3 + QQ(2,15)*x**5 + QQ(17,315)*x**7
     assert rs_tan(x*y + x**2*y**3, x, 9) == 4*x**8*y**11/3 + 17*x**8*y**9/45 + \
         4*x**7*y**9/3 + 17*x**7*y**7/315 + x**6*y**9/3 + 2*x**6*y**7/3 + \
         x**5*y**7 + 2*x**5*y**5/15 + x**4*y**5 + x**3*y**3/3 + x**2*y**3 + x*y
@@ -301,18 +304,19 @@ def test_tan():
     assert rs_atan(p, x, 10).compose(x, 10) == EX(atan(5) + S(67701870330562640) / \
         668083460499)
 
+
 def test_cot():
-    R, x, y = ring('x, y', QQ)
+    R, x, y = puiseux_ring('x, y', QQ)
     assert rs_cot(x**6 + x**7, x, 8) == x**(-6) - x**(-5) + x**(-4) - \
         x**(-3) + x**(-2) - x**(-1) + 1 - x + x**2 - x**3 + x**4 - x**5 + \
         2*x**6/3 - 4*x**7/3
     assert rs_cot(x + x**2*y, x, 5) == -x**4*y**5 - x**4*y/15 + x**3*y**4 - \
         x**3/45 - x**2*y**3 - x**2*y/3 + x*y**2 - x/3 - y + x**(-1)
 
+
 def test_sin():
     R, x, y = ring('x, y', QQ)
-    assert rs_sin(x, x, 9)/x**5 == \
-        Rational(-1, 5040)*x**2 + Rational(1, 120) - Rational(1, 6)*x**(-2) + x**(-4)
+    assert rs_sin(x, x, 9) == x - x**3/6 + x**5/120 - x**7/5040
     assert rs_sin(x*y + x**2*y**3, x, 9) == x**8*y**11/12 - \
         x**8*y**9/720 + x**7*y**9/12 - x**7*y**7/5040 - x**6*y**9/6 + \
         x**6*y**7/24 - x**5*y**7/2 + x**5*y**5/120 - x**4*y**5/2 - \
@@ -337,8 +341,7 @@ def test_sin():
 
 def test_cos():
     R, x, y = ring('x, y', QQ)
-    assert rs_cos(x, x, 9)/x**5 == \
-        Rational(1, 40320)*x**3 - Rational(1, 720)*x + Rational(1, 24)*x**(-1) - S.Half*x**(-3) + x**(-5)
+    assert rs_cos(x, x, 9) == 1 - x**2/2 + x**4/24 - x**6/720 + x**8/40320
     assert rs_cos(x*y + x**2*y**3, x, 9) == x**8*y**12/24 - \
         x**8*y**10/48 + x**8*y**8/40320 + x**7*y**10/6 - \
         x**7*y**8/120 + x**6*y**8/4 - x**6*y**6/720 + x**5*y**6/6 - \
@@ -372,7 +375,7 @@ def test_cos_sin():
 
 def test_atanh():
     R, x, y = ring('x, y', QQ)
-    assert rs_atanh(x, x, 9)/x**5 == Rational(1, 7)*x**2 + Rational(1, 5) + Rational(1, 3)*x**(-2) + x**(-4)
+    assert rs_atanh(x, x, 9) == x + x**3/3 + x**5/5 + x**7/7
     assert rs_atanh(x*y + x**2*y**3, x, 9) == 2*x**8*y**11 + x**8*y**9 + \
         2*x**7*y**9 + x**7*y**7/7 + x**6*y**9/3 + x**6*y**7 + x**5*y**7 + \
         x**5*y**5/5 + x**4*y**5 + x**3*y**3/3 + x**2*y**3 + x*y
@@ -395,7 +398,7 @@ def test_atanh():
 
 def test_sinh():
     R, x, y = ring('x, y', QQ)
-    assert rs_sinh(x, x, 9)/x**5 == Rational(1, 5040)*x**2 + Rational(1, 120) + Rational(1, 6)*x**(-2) + x**(-4)
+    assert rs_sinh(x, x, 9) == x + x**3/6 + x**5/120 + x**7/5040
     assert rs_sinh(x*y + x**2*y**3, x, 9) == x**8*y**11/12 + \
         x**8*y**9/720 + x**7*y**9/12 + x**7*y**7/5040 + x**6*y**9/6 + \
         x**6*y**7/24 + x**5*y**7/2 + x**5*y**5/120 + x**4*y**5/2 + \
@@ -403,8 +406,7 @@ def test_sinh():
 
 def test_cosh():
     R, x, y = ring('x, y', QQ)
-    assert rs_cosh(x, x, 9)/x**5 == Rational(1, 40320)*x**3 + Rational(1, 720)*x + Rational(1, 24)*x**(-1) + \
-        S.Half*x**(-3) + x**(-5)
+    assert rs_cosh(x, x, 9) == 1 + x**2/2 + x**4/24 + x**6/720 + x**8/40320
     assert rs_cosh(x*y + x**2*y**3, x, 9) == x**8*y**12/24 + \
         x**8*y**10/48 + x**8*y**8/40320 + x**7*y**10/6 + \
         x**7*y**8/120 + x**6*y**8/4 + x**6*y**6/720 + x**5*y**6/6 + \
@@ -412,7 +414,7 @@ def test_cosh():
 
 def test_tanh():
     R, x, y = ring('x, y', QQ)
-    assert rs_tanh(x, x, 9)/x**5 == Rational(-17, 315)*x**2 + Rational(2, 15) - Rational(1, 3)*x**(-2) + x**(-4)
+    assert rs_tanh(x, x, 9) == x - QQ(1,3)*x**3 + QQ(2,15)*x**5 - QQ(17,315)*x**7
     assert rs_tanh(x*y + x**2*y**3, x, 9) == 4*x**8*y**11/3 - \
         17*x**8*y**9/45 + 4*x**7*y**9/3 - 17*x**7*y**7/315 - x**6*y**9/3 + \
         2*x**6*y**7/3 - x**5*y**7 + 2*x**5*y**5/15 - x**4*y**5 - \
@@ -443,8 +445,9 @@ def test_RR():
     q = ((2 + a)**QQ(1, 5)).series(a, 0, 5).removeO()
     is_close(p.as_expr(), q.subs(a, 5).n())
 
+
 def test_is_regular():
-    R, x, y = ring('x, y', QQ)
+    R, x, y = puiseux_ring('x, y', QQ)
     p = 1 + 2*x + x**2 + 3*x**3
     assert not rs_is_puiseux(p, x)
 
@@ -455,8 +458,9 @@ def test_is_regular():
     p = x + x**2*y**QQ(1,5)*y
     assert not rs_is_puiseux(p, x)
 
+
 def test_puiseux():
-    R, x, y = ring('x, y', QQ)
+    R, x, y = puiseux_ring('x, y', QQ)
     p = x**QQ(2,5) + x**QQ(2,3) + x
 
     r = rs_series_inversion(p, x, 1)
@@ -518,20 +522,21 @@ def test_puiseux():
     assert r == -x**QQ(9,5) - x**QQ(26,15) - x**QQ(22,15) - x**QQ(6,5)/3 + \
         x + x**QQ(2,3) + x**QQ(2,5)
 
+
 def test_puiseux_algebraic(): # https://github.com/sympy/sympy/issues/24395
 
     K = QQ.algebraic_field(sqrt(2))
     sqrt2 = K.from_sympy(sqrt(2))
     x, y = symbols('x, y')
-    R, xr, yr = ring([x, y], K)
+    R, xr, yr = puiseux_ring([x, y], K)
     p = (1+sqrt2)*xr**QQ(1,2) + (1-sqrt2)*yr**QQ(2,3)
 
-    assert dict(p) == {(QQ(1,2),QQ(0)):1+sqrt2, (QQ(0),QQ(2,3)):1-sqrt2}
+    assert p.to_dict() == {(QQ(1,2),QQ(0)):1+sqrt2, (QQ(0),QQ(2,3)):1-sqrt2}
     assert p.as_expr() == (1 + sqrt(2))*x**(S(1)/2) + (1 - sqrt(2))*y**(S(2)/3)
 
 
 def test1():
-    R, x = ring('x', QQ)
+    R, x = puiseux_ring('x', QQ)
     r = rs_sin(x, x, 15)*x**(-5)
     assert r == x**8/6227020800 - x**6/39916800 + x**4/362880 - x**2/5040 + \
         QQ(1,120) - x**-2/6 + x**-4
@@ -556,9 +561,10 @@ def test1():
         x**3/720 + x**QQ(5,2)/120 + x**2/24 + x**QQ(3,2)/6 + x/2 + \
         x**QQ(1,2) + 1
 
+
 def test_puiseux2():
     R, y = ring('y', QQ)
-    S, x = ring('x', R)
+    S, x = puiseux_ring('x', R.to_domain())
 
     p = x + x**QQ(1,5)*y
     r = rs_atan(p, x, 3)
diff --git a/sympy/polys/tests/test_rings.py b/sympy/polys/tests/test_rings.py
index 3a48d45a6f15..4fd3c4c05eff 100644
--- a/sympy/polys/tests/test_rings.py
+++ b/sympy/polys/tests/test_rings.py
@@ -58,19 +58,10 @@ def test_PolyRing___hash__():
 
 def test_PolyRing___eq__():
     assert ring("x,y,z", QQ)[0] == ring("x,y,z", QQ)[0]
-    assert ring("x,y,z", QQ)[0] is ring("x,y,z", QQ)[0]
-
     assert ring("x,y,z", QQ)[0] != ring("x,y,z", ZZ)[0]
-    assert ring("x,y,z", QQ)[0] is not ring("x,y,z", ZZ)[0]
-
     assert ring("x,y,z", ZZ)[0] != ring("x,y,z", QQ)[0]
-    assert ring("x,y,z", ZZ)[0] is not ring("x,y,z", QQ)[0]
-
     assert ring("x,y,z", QQ)[0] != ring("x,y", QQ)[0]
-    assert ring("x,y,z", QQ)[0] is not ring("x,y", QQ)[0]
-
     assert ring("x,y", QQ)[0] != ring("x,y,z", QQ)[0]
-    assert ring("x,y", QQ)[0] is not ring("x,y,z", QQ)[0]
 
 def test_PolyRing_ring_new():
     R, x, y, z = ring("x,y,z", QQ)
@@ -288,23 +279,23 @@ def test_PolyElement_from_expr():
     R, X, Y, Z = ring((x, y, z), ZZ)
 
     f = R.from_expr(1)
-    assert f == 1 and isinstance(f, R.dtype)
+    assert f == 1 and R.is_element(f)
 
     f = R.from_expr(x)
-    assert f == X and isinstance(f, R.dtype)
+    assert f == X and R.is_element(f)
 
     f = R.from_expr(x*y*z)
-    assert f == X*Y*Z and isinstance(f, R.dtype)
+    assert f == X*Y*Z and R.is_element(f)
 
     f = R.from_expr(x*y*z + x*y + x)
-    assert f == X*Y*Z + X*Y + X and isinstance(f, R.dtype)
+    assert f == X*Y*Z + X*Y + X and R.is_element(f)
 
     f = R.from_expr(x**3*y*z + x**2*y**7 + 1)
-    assert f == X**3*Y*Z + X**2*Y**7 + 1 and isinstance(f, R.dtype)
+    assert f == X**3*Y*Z + X**2*Y**7 + 1 and R.is_element(f)
 
     r, F = sring([exp(2)])
     f = r.from_expr(exp(2))
-    assert f == F[0] and isinstance(f, r.dtype)
+    assert f == F[0] and r.is_element(f)
 
     raises(ValueError, lambda: R.from_expr(1/x))
     raises(ValueError, lambda: R.from_expr(2**x))
@@ -312,7 +303,7 @@ def test_PolyElement_from_expr():
 
     R, = ring("", ZZ)
     f = R.from_expr(1)
-    assert f == 1 and isinstance(f, R.dtype)
+    assert f == 1 and R.is_element(f)
 
 def test_PolyElement_degree():
     R, x,y,z = ring("x,y,z", ZZ)
@@ -615,9 +606,9 @@ def test_PolyElement___truediv__():
     assert (x**2 - 1).quo(x) == x
     assert (x**2 - x).quo(x) == x - 1
 
-    assert (x**2 - 1)/x == x - x**(-1)
+    raises(ExactQuotientFailed, lambda: (x**2 - 1)/x)
     assert (x**2 - x)/x == x - 1
-    assert (x**2 - 1)/(2*x) == x/2 - x**(-1)/2
+    raises(ExactQuotientFailed, lambda: (x**2 - 1)/(2*x))
 
     assert (x**2 - 1).quo(2*x) == 0
     assert (x**2 - x)/(x - 1) == (x**2 - x).quo(x - 1) == x
@@ -634,7 +625,7 @@ def test_PolyElement___truediv__():
     Rxyz, x,y,z = ring("x,y,z", Ruv)
 
     assert dict((u**2*x + u)/u) == {(1, 0, 0): u, (0, 0, 0): 1}
-    raises(TypeError, lambda: u/(u**2*x + u))
+    raises(ExactQuotientFailed, lambda: u/(u**2*x + u))
 
     raises(TypeError, lambda: t/x)
     raises(TypeError, lambda: x/t)
@@ -670,7 +661,8 @@ def test_PolyElement___truediv__():
 
     assert f.div(g) == divmod(f, g) == (q, r)
     assert f.rem(g) == f % g == r
-    assert f.quo(g) == f / g == q
+    assert f.quo(g) == q
+    raises(ExactQuotientFailed, lambda: f / g)
     raises(ExactQuotientFailed, lambda: f.exquo(g))
 
     f, g = 3*x**3 + x**2 + x + 5, 5*x**2 - 3*x + 1
@@ -678,7 +670,8 @@ def test_PolyElement___truediv__():
 
     assert f.div(g) == divmod(f, g) == (q, r)
     assert f.rem(g) == f % g == r
-    assert f.quo(g) == f / g == q
+    assert f.quo(g) == q
+    raises(ExactQuotientFailed, lambda: f / g)
     raises(ExactQuotientFailed, lambda: f.exquo(g))
 
     f, g = 5*x**4 + 4*x**3 + 3*x**2 + 2*x + 1, x**2 + 2*x + 3
@@ -686,7 +679,8 @@ def test_PolyElement___truediv__():
 
     assert f.div(g) == divmod(f, g) == (q, r)
     assert f.rem(g) == f % g == r
-    assert f.quo(g) == f / g == q
+    assert f.quo(g) == q
+    raises(ExactQuotientFailed, lambda: f / g)
     raises(ExactQuotientFailed, lambda: f.exquo(g))
 
     f, g = 5*x**5 + 4*x**4 + 3*x**3 + 2*x**2 + x, x**4 + 2*x**3 + 9
@@ -694,7 +688,8 @@ def test_PolyElement___truediv__():
 
     assert f.div(g) == divmod(f, g) == (q, r)
     assert f.rem(g) == f % g == r
-    assert f.quo(g) == f / g == q
+    assert f.quo(g) == q
+    raises(ExactQuotientFailed, lambda: f / g)
     raises(ExactQuotientFailed, lambda: f.exquo(g))
 
     R, x = ring("x", QQ)
@@ -704,7 +699,8 @@ def test_PolyElement___truediv__():
 
     assert f.div(g) == divmod(f, g) == (q, r)
     assert f.rem(g) == f % g == r
-    assert f.quo(g) == f / g == q
+    assert f.quo(g) == q
+    raises(ExactQuotientFailed, lambda: f / g)
     raises(ExactQuotientFailed, lambda: f.exquo(g))
 
     f, g = 3*x**3 + x**2 + x + 5, 5*x**2 - 3*x + 1
@@ -712,7 +708,8 @@ def test_PolyElement___truediv__():
 
     assert f.div(g) == divmod(f, g) == (q, r)
     assert f.rem(g) == f % g == r
-    assert f.quo(g) == f / g == q
+    assert f.quo(g) == q
+    raises(ExactQuotientFailed, lambda: f / g)
     raises(ExactQuotientFailed, lambda: f.exquo(g))
 
     R, x,y = ring("x,y", ZZ)
@@ -722,15 +719,16 @@ def test_PolyElement___truediv__():
 
     assert f.div(g) == divmod(f, g) == (q, r)
     assert f.rem(g) == f % g == r
-    assert f.quo(g) == f / g == q
-    assert f.exquo(g) == q
+    assert f.quo(g) == q
+    assert f.exquo(g) == f / g == q
 
     f, g = x**2 + y**2, x - y
     q, r = x + y, 2*y**2
 
     assert f.div(g) == divmod(f, g) == (q, r)
     assert f.rem(g) == f % g == r
-    assert f.quo(g) == f / g == q
+    assert f.quo(g) == q
+    raises(ExactQuotientFailed, lambda: f / g)
     raises(ExactQuotientFailed, lambda: f.exquo(g))
 
     f, g = x**2 + y**2, -x + y
@@ -738,7 +736,8 @@ def test_PolyElement___truediv__():
 
     assert f.div(g) == divmod(f, g) == (q, r)
     assert f.rem(g) == f % g == r
-    assert f.quo(g) == f / g == q
+    assert f.quo(g) == q
+    raises(ExactQuotientFailed, lambda: f / g)
     raises(ExactQuotientFailed, lambda: f.exquo(g))
 
     f, g = x**2 + y**2, 2*x - 2*y
@@ -746,7 +745,8 @@ def test_PolyElement___truediv__():
 
     assert f.div(g) == divmod(f, g) == (q, r)
     assert f.rem(g) == f % g == r
-    assert f.quo(g) == f / g == q
+    assert f.quo(g) == q
+    raises(ExactQuotientFailed, lambda: f / g)
     raises(ExactQuotientFailed, lambda: f.exquo(g))
 
     R, x,y = ring("x,y", QQ)
@@ -756,15 +756,16 @@ def test_PolyElement___truediv__():
 
     assert f.div(g) == divmod(f, g) == (q, r)
     assert f.rem(g) == f % g == r
-    assert f.quo(g) == f / g == q
-    assert f.exquo(g) == q
+    assert f.quo(g) == q
+    assert f.exquo(g) == f / g == q
 
     f, g = x**2 + y**2, x - y
     q, r = x + y, 2*y**2
 
     assert f.div(g) == divmod(f, g) == (q, r)
     assert f.rem(g) == f % g == r
-    assert f.quo(g) == f / g == q
+    assert f.quo(g) == q
+    raises(ExactQuotientFailed, lambda: f / g)
     raises(ExactQuotientFailed, lambda: f.exquo(g))
 
     f, g = x**2 + y**2, -x + y
@@ -772,7 +773,8 @@ def test_PolyElement___truediv__():
 
     assert f.div(g) == divmod(f, g) == (q, r)
     assert f.rem(g) == f % g == r
-    assert f.quo(g) == f / g == q
+    assert f.quo(g) == q
+    raises(ExactQuotientFailed, lambda: f / g)
     raises(ExactQuotientFailed, lambda: f.exquo(g))
 
     f, g = x**2 + y**2, 2*x - 2*y
@@ -780,7 +782,8 @@ def test_PolyElement___truediv__():
 
     assert f.div(g) == divmod(f, g) == (q, r)
     assert f.rem(g) == f % g == r
-    assert f.quo(g) == f / g == q
+    assert f.quo(g) == q
+    raises(ExactQuotientFailed, lambda: f / g)
     raises(ExactQuotientFailed, lambda: f.exquo(g))
 
 def test_PolyElement___pow__():
@@ -791,8 +794,6 @@ def test_PolyElement___pow__():
     assert f**1 == f
     raises(ValueError, lambda: f**(-1))
 
-    assert x**(-1) == x**(-1)
-
     assert f**2 == f._pow_generic(2) == f._pow_multinomial(2) == 4*x**2 + 12*x + 9
     assert f**3 == f._pow_generic(3) == f._pow_multinomial(3) == 8*x**3 + 36*x**2 + 54*x + 27
     assert f**4 == f._pow_generic(4) == f._pow_multinomial(4) == 16*x**4 + 96*x**3 + 216*x**2 + 216*x + 81
@@ -1172,13 +1173,13 @@ def test_PolyElement_evaluate():
     f = (x*y)**3 + 4*(x*y)**2 + 2*x*y + 3
 
     r = f.evaluate(x, 0)
-    assert r == 3 and isinstance(r, R.drop(x).dtype)
+    assert r == 3 and R.drop(x).is_element(r)
     r = f.evaluate([(x, 0), (y, 0)])
-    assert r == 3 and isinstance(r, R.drop(x, y).dtype)
+    assert r == 3 and R.drop(x, y).is_element(r)
     r = f.evaluate(y, 0)
-    assert r == 3 and isinstance(r, R.drop(y).dtype)
+    assert r == 3 and R.drop(y).is_element(r)
     r = f.evaluate([(y, 0), (x, 0)])
-    assert r == 3 and isinstance(r, R.drop(y, x).dtype)
+    assert r == 3 and R.drop(y, x).is_element(r)
 
     r = f.evaluate([(x, 0), (y, 0), (z, 0)])
     assert r == 3 and not isinstance(r, PolyElement)
@@ -1192,7 +1193,7 @@ def test_PolyElement_subs():
     f = x**3 + 4*x**2 + 2*x + 3
 
     r = f.subs(x, 0)
-    assert r == 3 and isinstance(r, R.dtype)
+    assert r == 3 and R.is_element(r)
 
     raises(CoercionFailed, lambda: f.subs(x, QQ(1,7)))
 
@@ -1200,9 +1201,9 @@ def test_PolyElement_subs():
     f = x**3 + 4*x**2 + 2*x + 3
 
     r = f.subs(x, 0)
-    assert r == 3 and isinstance(r, R.dtype)
+    assert r == 3 and R.is_element(r)
     r = f.subs([(x, 0), (y, 0)])
-    assert r == 3 and isinstance(r, R.dtype)
+    assert r == 3 and R.is_element(r)
 
     raises(CoercionFailed, lambda: f.subs([(x, 1), (y, QQ(1,7))]))
     raises(CoercionFailed, lambda: f.subs([(x, QQ(1,7)), (y, 1)]))
@@ -1252,7 +1253,7 @@ def test_PolyElement_compose():
     f = x**3 + 4*x**2 + 2*x + 3
 
     r = f.compose(x, 0)
-    assert r == 3 and isinstance(r, R.dtype)
+    assert r == 3 and R.is_element(r)
 
     assert f.compose(x, x) == f
     assert f.compose(x, x**2) == x**6 + 4*x**4 + 2*x**2 + 3
@@ -1263,13 +1264,13 @@ def test_PolyElement_compose():
     f = x**3 + 4*x**2 + 2*x + 3
 
     r = f.compose(x, 0)
-    assert r == 3 and isinstance(r, R.dtype)
+    assert r == 3 and R.is_element(r)
     r = f.compose([(x, 0), (y, 0)])
-    assert r == 3 and isinstance(r, R.dtype)
+    assert r == 3 and R.is_element(r)
 
     r = (x**3 + 4*x**2 + 2*x*y*z + 3).compose(x, y*z**2 - 1)
     q = (y*z**2 - 1)**3 + 4*(y*z**2 - 1)**2 + 2*(y*z**2 - 1)*y*z + 3
-    assert r == q and isinstance(r, R.dtype)
+    assert r == q and R.is_element(r)
 
 def test_PolyElement_is_():
     R, x,y,z = ring("x,y,z", QQ)
@@ -1350,7 +1351,7 @@ def test_PolyElement_drop():
 
     assert R(1).drop(0).ring == PolyRing("y,z", ZZ, lex)
     assert R(1).drop(0).drop(0).ring == PolyRing("z", ZZ, lex)
-    assert isinstance(R(1).drop(0).drop(0).drop(0), R.dtype) is False
+    assert R.is_element(R(1).drop(0).drop(0).drop(0)) is False
 
     raises(ValueError, lambda: z.drop(0).drop(0).drop(0))
     raises(ValueError, lambda: x.drop(0))
diff --git a/sympy/polys/tests/test_solvers.py b/sympy/polys/tests/test_solvers.py
index 9b7c2b3c9f74..bf8708314466 100644
--- a/sympy/polys/tests/test_solvers.py
+++ b/sympy/polys/tests/test_solvers.py
@@ -13,7 +13,7 @@ def test_solve_lin_sys_2x2_one():
            2*x1 - x2]
     sol = {x1: QQ(5, 3), x2: QQ(10, 3)}
     _sol = solve_lin_sys(eqs, domain)
-    assert _sol == sol and all(isinstance(s, domain.dtype) for s in _sol)
+    assert _sol == sol and all(s.ring == domain for s in _sol)
 
 def test_solve_lin_sys_2x4_none():
     domain, x1,x2 = ring("x1,x2", QQ)
diff --git a/sympy/printing/tests/test_str.py b/sympy/printing/tests/test_str.py
index 5140136b2465..675212964b03 100644
--- a/sympy/printing/tests/test_str.py
+++ b/sympy/printing/tests/test_str.py
@@ -463,8 +463,6 @@ def test_PolyElement():
     assert str(x - 1) == "x - 1"
     assert str(x + 1) == "x + 1"
     assert str(x**2) == "x**2"
-    assert str(x**(-2)) == "x**(-2)"
-    assert str(x**QQ(1, 2)) == "x**(1/2)"
 
     assert str((u**2 + 3*u*v + 1)*x**2*y + u + 1) == "(u**2 + 3*u*v + 1)*x**2*y + u + 1"
     assert str((u**2 + 3*u*v + 1)*x**2*y + (u + 1)*x) == "(u**2 + 3*u*v + 1)*x**2*y + (u + 1)*x"


Code review comments:
diff_hunk: @@ -0,0 +1,795 @@
+"""
+Puiseux rings. These are used by the ring_series module to represented
+truncated Puiseux series. Elements of a Puiseux ring are like polynomials
+except that the exponents can be negative or rational rather than just
+non-negative integers.
+"""
+
+# Previously the ring_series module used PolyElement to represent Puiseux
+# series. This is problematic because it means that PolyElement has to support
+# negative and non-integer exponents which most polynomial representations do
+# not support. This module provides an implementation of a ring for Puiseux
+# series that can be used by ring_series without breaking the basic invariants
+# of polynomial rings.
+#
+# Ideally there would be more of a proper series type that can keep track of
+# not not just the leading terms of a truncated series but also the precision
+# of the series. For now the rings here are just introduced to keep the
+# interface that ring_series was using before.
+
+from __future__ import annotations
+
+from sympy.polys.domains import QQ
+from sympy.polys.rings import PolyRing, PolyElement
+from sympy.core.add import Add
+from sympy.core.mul import Mul
+from sympy.external.gmpy import gcd, lcm
+
+
+from typing import TYPE_CHECKING
+
+
+if TYPE_CHECKING:
+    from typing import Any, Unpack
+    from sympy.core.expr import Expr
+    from sympy.polys.domains import Domain
+    from collections.abc import Iterable, Iterator
+
+
+def puiseux_ring(
+    symbols: str | list[Expr], domain: Domain
+) -> tuple[PuiseuxRing, Unpack[tuple[PuiseuxPoly, ...]]]:
+    """Construct a Puiseux ring.
+
+    This function constructs a Puiseux ring with the given symbols and domain.
+
+    >>> from sympy.polys.domains import QQ
+    >>> from sympy.polys.puiseux import puiseux_ring
+    >>> R, x, y = puiseux_ring('x y', QQ)
+    >>> R
+    PuiseuxRing((x, y), QQ)
+    >>> p = 5*x**QQ(1,2) + 7/y
+    >>> p
+    7*y**(-1) + 5*x**(1/2)
+    """
+    ring = PuiseuxRing(symbols, domain)
+    return (ring,) + ring.gens # type: ignore
+
+
+class PuiseuxRing:
+    """Ring of Puiseux polynomials.
+
+    A Puiseux polynomial is a truncated Puiseux series. The exponents of the
+    monomials can be negative or rational numbers. This ring is used by the
+    ring_series module:
+
+    >>> from sympy.polys.domains import QQ
+    >>> from sympy.polys.puiseux import puiseux_ring
+    >>> from sympy.polys.ring_series import rs_exp, rs_nth_root
+    >>> ring, x, y = puiseux_ring('x y', QQ)
+    >>> f = x**2 + y**3
+    >>> f
+    y**3 + x**2
+    >>> f.diff(x)
+    2*x
+    >>> rs_exp(x, x, 5)
+    1 + x + 1/2*x**2 + 1/6*x**3 + 1/24*x**4
+
+    Importantly the Puiseux ring can represent truncated series with negative
+    and fractional exponents:
+
+    >>> f = 1/x + 1/y**2
+    >>> f
+    x**(-1) + y**(-2)
+    >>> f.diff(x)
+    -1*x**(-2)
+
+    >>> rs_nth_root(8*x + x**2 + x**3, 3, x, 5)
+    2*x**(1/3) + 1/12*x**(4/3) + 23/288*x**(7/3) + -139/20736*x**(10/3)
+
+    See Also
+    ========
+
+    sympy.polys.ring_series.rs_series
+    PuiseuxPoly
+    """
+    def __init__(self, symbols: str | list[Expr], domain: Domain):
+
+        poly_ring = PolyRing(symbols, domain)
+
+        domain = poly_ring.domain
+        ngens = poly_ring.ngens
+
+        self.poly_ring = poly_ring
+        self.domain = domain
+
+        self.symbols = poly_ring.symbols
+        self.gens = tuple([self.from_poly(g) for g in poly_ring.gens])
+        self.ngens = ngens
+
+        self.zero = self.from_poly(poly_ring.zero)
+        self.one = self.from_poly(poly_ring.one)
+
+        self.zero_monom = poly_ring.zero_monom # type: ignore
+        self.monomial_mul = poly_ring.monomial_mul # type: ignore
+
+    def __repr__(self) -> str:
+        return f"PuiseuxRing({self.symbols}, {self.domain})"
+
+    def __eq__(self, other: Any) -> bool:
+        if not isinstance(other, PuiseuxRing):
+            return NotImplemented
+        return self.symbols == other.symbols and self.domain == other.domain
+
+    def from_poly(self, poly: PolyElement) -> PuiseuxPoly:
+        """Create a Puiseux polynomial from a polynomial.
+
+        >>> from sympy.polys.domains import QQ
+        >>> from sympy.polys.rings import ring
+        >>> from sympy.polys.puiseux import puiseux_ring
+        >>> R1, x1 = ring('x', QQ)
+        >>> R2, x2 = puiseux_ring('x', QQ)
+        >>> R2.from_poly(x1**2)
+        x**2
+        """
+        return PuiseuxPoly(poly, self)
+
+    def from_dict(self, terms: dict[tuple[int, ...], Any]) -> PuiseuxPoly:
+        """Create a Puiseux polynomial from a dictionary of terms.
+
+        >>> from sympy.polys.domains import QQ
+        >>> from sympy.polys.puiseux import puiseux_ring
+        >>> R, x = puiseux_ring('x', QQ)
+        >>> R.from_dict({(QQ(1,2),): QQ(3)})
+        3*x**(1/2)
+        """
+        return PuiseuxPoly.from_dict(terms, self)
+
+    def from_int(self, n: int) -> PuiseuxPoly:
+        """Create a Puiseux polynomial from an integer.
+
+        >>> from sympy.polys.domains import QQ
+        >>> from sympy.polys.puiseux import puiseux_ring
+        >>> R, x = puiseux_ring('x', QQ)
+        >>> R.from_int(3)
+        3
+        """
+        return self.from_poly(self.poly_ring(n))
+
+    def domain_new(self, arg: Any) -> Any:
+        """Create a new element of the domain.
+
+        >>> from sympy.polys.domains import QQ
+        >>> from sympy.polys.puiseux import puiseux_ring
+        >>> R, x = puiseux_ring('x', QQ)
+        >>> R.domain_new(3)
+        3
+        >>> QQ.of_type(_)
+        True
+        """
+        return self.poly_ring.domain_new(arg)
+
+    def ground_new(self, arg: Any) -> PuiseuxPoly:
+        """Create a new element from a ground element.
+
+        >>> from sympy.polys.domains import QQ
+        >>> from sympy.polys.puiseux import puiseux_ring, PuiseuxPoly
+        >>> R, x = puiseux_ring('x', QQ)
+        >>> R.ground_new(3)
+        3
+        >>> isinstance(_, PuiseuxPoly)
+        True
+        """
+        return self.from_poly(self.poly_ring.ground_new(arg))
+
+    def __call__(self, arg: Any) -> PuiseuxPoly:
+        """Coerce an element into the ring.
+
+        >>> from sympy.polys.domains import QQ
+        >>> from sympy.polys.puiseux import puiseux_ring
+        >>> R, x = puiseux_ring('x', QQ)
+        >>> R(3)
+        3
+        >>> R({(QQ(1,2),): QQ(3)})
+        3*x**(1/2)
+        """
+        if isinstance(arg, dict):
+            return self.from_dict(arg)
+        else:
+            return self.from_poly(self.poly_ring(arg))
+
+    def index(self, x: PuiseuxPoly) -> int:
+        """Return the index of a generator.
+
+        >>> from sympy.polys.domains import QQ
+        >>> from sympy.polys.puiseux import puiseux_ring
+        >>> R, x, y = puiseux_ring('x y', QQ)
+        >>> R.index(x)
+        0
+        >>> R.index(y)
+        1
+        """
+        return self.gens.index(x)
+
+
+def _div_poly_monom(poly: PolyElement, monom: Iterable[int]) -> PolyElement:
+    ring = poly.ring
+    div = ring.monomial_div
+    return ring.from_dict({div(m, monom): c for m, c in poly.terms()})
+
+
+def _mul_poly_monom(poly: PolyElement, monom: Iterable[int]) -> PolyElement:
+    ring = poly.ring
+    mul = ring.monomial_mul
+    return ring.from_dict({mul(m, monom): c for m, c in poly.terms()})
+
+
+def _div_monom(monom: Iterable[int], div: Iterable[int]) -> tuple[int, ...]:
+    return tuple(mi - di for mi, di in zip(monom, div))
+
+
+class PuiseuxPoly:
+    """Puiseux polynomial. Represents a truncated Puiseux series.
+
+    See the :class:`PuiseuxRing` class for more information.
+
+    >>> from sympy import QQ
+    >>> from sympy.polys.puiseux import puiseux_ring
+    >>> R, x, y = puiseux_ring('x, y', QQ)
+    >>> p = 5*x**2 + 7*y**3
+    >>> p
+    7*y**3 + 5*x**2
+
+    The internal representation of a Puiseux polynomial wraps a normal
+    polynomial. To support negative powers the polynomial is considered to be
+    divided by a monomial.
+
+    >>> p2 = 1/x + 1/y**2
+    >>> p2.monom # x*y**2
+    (1, 2)
+    >>> p2.poly
+    x + y**2
+    >>> (y**2 + x) / (x*y**2) == p2
+    True
+
+    To support fractional powers the polynomial is considered to be a function
+    of ``x**(1/nx) * y**(1/ny) * ...``. The representation keeps track of a
path: sympy/polys/puiseux.py
text: @user1:
Should this be
```suggestion
    of ``x**(1/nx), y**(1/ny), ...``. The representation keeps track of a
```
?