Repository: sympy/sympy

Issue:
Conditions for Justifying Rewrite in `limit`
In `Limit.doit()`, the expression $\lim_{x \to \infty} f(x)$ is handled by rewriting it as $\lim_{x \to +0} f(1/x)$.

https://github.com/sympy/sympy/blob/d5d4d8e4ddf28c42075a3e124752862e229e651a/sympy/series/limits.py#L326-L331

However, this substitution is not always valid. For example, consider the following case:

```python
n = symbols('n', integer=True, positive=True)
limit(floor(n/2), n, oo)
```

Here, `floor(n/2)` is transformed into `floor(1/(2*n))`, but since `n` is a positive integer, `floor(1/(2*n))` is always zero. Previously, due to limitations in `floor`, it wasn’t possible to evaluate `floor(1/(2*n))` as zero, so the final result remained correct. However, with recent improvements (#27769) that allow assumptions to be used to evaluate `floor(1/(2*n))` as zero, this issue has surfaced.

Another example that fails due to the improved evaluation of `floor` is the following test code:

```python
n = Symbol('n')
Sum(S.Half**floor(n/2), (n, 1, oo)).is_convergent()
```

Even though no assumptions are explicitly given for `n`, it is implicitly treated as a positive integer at the following location, leading to the same problem:

https://github.com/sympy/sympy/blob/d5d4d8e4ddf28c42075a3e124752862e229e651a/sympy/concrete/summations.py#L477

A possible workaround might be to ensure that `1/x`, when substituted for `x`, does not inherit the assumptions of the original `x`, or to use a dummy variable instead. However, this feels like a band-aid solution. What are the conditions under which this substitution can be safely applied? I believe it's necessary to clarify this and improve the `limit` implementation accordingly.


Discussion:
sympy-bot: :white_check_mark:

Hi, I am the [SymPy bot](https://github.com/sympy/sympy-bot). I'm here to help you write a release notes entry. Please read the [guide on how to write release notes](https://github.com/sympy/sympy/wiki/Writing-Release-Notes).



* No release notes entry will be added for this pull request.

<details><summary>Click here to see the pull request description that was parsed.</summary>


    <!-- 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 #27773.

    #### Brief description of what is fixed or changed
    Previously, `Limit.doit()` would rewrite limits of the form  

    $$
    \lim_{x \to \infty} f(x) \quad \text{as} \quad \lim_{x \to 0^+} f(1/x)
    $$

    without accounting for cases where $x$ is a discrete variable (e.g., an integer).  

    This caused incorrect evaluations when the function involved was discontinuous, such as `floor`, `ceiling`, or `mod`.  

    Before this fix, evaluating the following limit:

    $$
    \lim_{n \to \infty} \lfloor \frac{n}{2} \rfloor
    $$

    resulted in:

    $$
    \lim_{n \to 0^+} \lfloor \frac{1}{2n} \rfloor = 0
    $$

    which is incorrect.

    To resolve this issue, substitution $x$ to $1/x$ is avoided when the limit variable is an integer and the function being evaluated contains discontinuous operations (`floor`, `ceiling`, `mod`).  

    In such cases, the limit is now computed using the Gruntz algorithm directly

    #### Other comments
    This change allows the removal of a previously XFAILED test, as the issue it covered is now resolved.

    #### 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 -->
    NO ENTRY
    <!-- END RELEASE NOTES -->


</details><p>

voaidesr: So, at this point, I think I have an answer to the question in #27773.

> What are the conditions under which this substitution can be safely applied? I believe it's necessary to clarify this and improve the limit implementation accordingly.

As I mentioned initially, the problem arises with discrete variables whose limits go to infinity.
Therefore, in the last commit, I sought to make the substitution without preserving the integer assumption, while still maintaining the sign assumptions.

Initially, I considered using a dummy variable with no assumptions at all, but this leads to issues - such as in `test_issue_22334`:
```python
assert limit((n+1)**k / ((n+1)**(k+1) - n**(k+1)).expand(), n, oo) == 1/(k + 1)
```
In this test, if `k` does not preserve its sign, an error will occur: `Not sure of sign of 1 - 2*k`
oscarbenjamin: Rather than preserving the sign of the symbol would it not make sense to just set the assumption that is implied by the given limit?

For example if the limit is `z -> oo` then in the substitution `z -> 1/z` we should have that both `z` and `1/z` are positive.
voaidesr: That was my initial approach as well, but it led to further complications.

For instance:
```python
def test_issue_15202():
    e = (2**x*(2 + 2**(-x)*(-2*2**x + x + 2))/(x + 1))**(x + 1)
    assert limit(e, x, oo) == exp(1)
```
This test runs indefinitely when the dummy variable is forced to be strictly positive as `z -> oo`.

I haven’t investigated the root of this issue in detail, because I think that attempting to patch such a narrow case by introducing additional constraints would make the implementation overly specific. It risks introducing subtle bugs elsewhere.

Preserving the sign assumptions of the original variable seems to be the most natural and elegant solution. So far, I haven’t encountered any cases where maintaining these assumptions has caused issues, which further supports keeping them intact.


oscarbenjamin: I think this looks good but can you squash the commits?
voaidesr: Yes, I'll do this now. 

Patch:
diff --git a/.mailmap b/.mailmap
index 655f82501b..b63d7d43f0 100644
--- a/.mailmap
+++ b/.mailmap
@@ -1303,7 +1303,7 @@ Roman Inflianskas <infroma@gmail.com>
 Ronan Lamy <ronan.lamy@gmail.com> <Ronan.Lamy@normalesup.org>
 Rudr Tiwari <rudrtiwari@gmail.com>
 Rupesh Harode <rupeshharode@gmail.com>
-Rushabh Mehta <mehtarushabh2005@gmail.com> Rushabh Mehta <139112780+RushabhMehta2005@users.noreply.github.com>
+Rushabh Mehta <mehtarushabh2005@gmail.com> RushabhMehta2005 <139112780+RushabhMehta2005@users.noreply.github.com>
 Ruslan Pisarev <rpisarev@cloudlinux.com> <ruslan@rpisarev.org.ua>
 Ryan Krauss <ryanlists@gmail.com>
 Rémy Léone <rleone@online.net>
@@ -1576,6 +1576,7 @@ Vladimir Lagunov <werehuman@gmail.com>
 Vladimir Perić <vlada.peric@gmail.com>
 Vladimir Poluhsin <vovapolu@gmail.com>
 Vladimir Sereda <voffch@gmail.com>
+Voaides Negustor Robert <134785947+voaidesr@users.noreply.github.com> <voaides.r@yahoo.com>
 Waldir Pimenta <waldyrious@gmail.com>
 Wang Ran (汪然) <wangr@smail.nju.edu.cn>
 Warren Jacinto <warrenjacinto@gmail.com>
diff --git a/sympy/algebras/quaternion.py b/sympy/algebras/quaternion.py
index 1bf4b36354..3251d92af3 100644
--- a/sympy/algebras/quaternion.py
+++ b/sympy/algebras/quaternion.py
@@ -1095,6 +1095,7 @@ def integrate(self, *args):
         4 + 8*i + 12*j + 16*k
 
         """
+        # TODO: is this expression correct?
         return Quaternion(integrate(self.a, *args), integrate(self.b, *args),
                           integrate(self.c, *args), integrate(self.d, *args))
 
diff --git a/sympy/core/mod.py b/sympy/core/mod.py
index 8be0c56e49..a5d7449c95 100644
--- a/sympy/core/mod.py
+++ b/sympy/core/mod.py
@@ -5,7 +5,6 @@
 from .logic import fuzzy_and, fuzzy_not
 from .mul import Mul
 from .numbers import equal_valued
-from .relational import is_le, is_lt, is_ge, is_gt
 from .singleton import S
 
 
@@ -97,20 +96,21 @@ def number_eval(p, q):
 
             # by difference
             # -2|q| < p < 2|q|
-            if q.is_positive:
-                comp1, comp2 = is_le, is_lt
-            elif q.is_negative:
-                comp1, comp2 = is_ge, is_gt
-            else:
-                return
-            ls = -2*q
-            r = p - q
-            for _ in range(4):
-                if not comp1(ls, p):
-                    return
-                if comp2(r, ls):
-                    return p - ls
-                ls += q
+            d = abs(p)
+            for _ in range(2):
+                d -= abs(q)
+                if d.is_negative:
+                    if q.is_positive:
+                        if p.is_positive:
+                            return d + q
+                        elif p.is_negative:
+                            return -d
+                    elif q.is_negative:
+                        if p.is_positive:
+                            return d
+                        elif p.is_negative:
+                            return -d + q
+                    break
 
         rv = number_eval(p, q)
         if rv is not None:
diff --git a/sympy/physics/mechanics/pathway.py b/sympy/physics/mechanics/pathway.py
index b86ba85b1d..0e5e2bd955 100644
--- a/sympy/physics/mechanics/pathway.py
+++ b/sympy/physics/mechanics/pathway.py
@@ -322,7 +322,7 @@ class ObstacleSetPathway(PathwayBase):
     Parameters
     ==========
 
-    attachments : tuple[Point, ...]
+    attachments : tuple[Point, Point]
         The set of ``Point`` objects that define the segmented obstacle-set
         pathway.
 
diff --git a/sympy/polys/monomials.py b/sympy/polys/monomials.py
index 43a1722386..e5897a0998 100644
--- a/sympy/polys/monomials.py
+++ b/sympy/polys/monomials.py
@@ -85,8 +85,8 @@ def itermonomials(variables, max_degrees, min_degrees=None):
         >>> sorted(itermonomials([x, y], [2, 4], [1, 2]), reverse=True, key=monomial_key('lex', [x, y]))
         [x**2*y**4, x**2*y**3, x**2*y**2, x*y**4, x*y**3, x*y**2]
     """
+    n = len(variables)
     if is_sequence(max_degrees):
-        n = len(variables)
         if len(max_degrees) != n:
             raise ValueError('Argument sizes do not match')
         if min_degrees is None:
@@ -98,13 +98,7 @@ def itermonomials(variables, max_degrees, min_degrees=None):
                 raise ValueError('Argument sizes do not match')
             if any(i < 0 for i in min_degrees):
                 raise ValueError("min_degrees cannot contain negative numbers")
-        if any(min_degrees[i] > max_degrees[i] for i in range(n)):
-            raise ValueError('min_degrees[i] must be <= max_degrees[i] for all i')
-        power_lists = []
-        for var, min_d, max_d in zip(variables, min_degrees, max_degrees):
-            power_lists.append([var**i for i in range(min_d, max_d + 1)])
-        for powers in product(*power_lists):
-            yield Mul(*powers)
+        total_degree = False
     else:
         max_degree = max_degrees
         if max_degree < 0:
@@ -115,6 +109,8 @@ def itermonomials(variables, max_degrees, min_degrees=None):
             if min_degrees < 0:
                 raise ValueError("min_degrees cannot be negative")
             min_degree = min_degrees
+        total_degree = True
+    if total_degree:
         if min_degree > max_degree:
             return
         if not variables or max_degree == 0:
@@ -123,21 +119,33 @@ def itermonomials(variables, max_degrees, min_degrees=None):
         # Force to list in case of passed tuple or other incompatible collection
         variables = list(variables) + [S.One]
         if all(variable.is_commutative for variable in variables):
-            it = combinations_with_replacement(variables, max_degree)
+            monomials_list_comm = []
+            for item in combinations_with_replacement(variables, max_degree):
+                powers = dict.fromkeys(variables, 0)
+                for variable in item:
+                    if variable != 1:
+                        powers[variable] += 1
+                if sum(powers.values()) >= min_degree:
+                    monomials_list_comm.append(Mul(*item))
+            yield from set(monomials_list_comm)
         else:
-            it = product(variables, repeat=max_degree)
-        monomials_set = set()
-        d = max_degree - min_degree
-        for item in it:
-            count = 0
-            for variable in item:
-                if variable == 1:
-                    count += 1
-                    if d < count:
-                        break
-            else:
-                monomials_set.add(Mul(*item))
-        yield from monomials_set
+            monomials_list_non_comm = []
+            for item in product(variables, repeat=max_degree):
+                powers = dict.fromkeys(variables, 0)
+                for variable in item:
+                    if variable != 1:
+                        powers[variable] += 1
+                if sum(powers.values()) >= min_degree:
+                    monomials_list_non_comm.append(Mul(*item))
+            yield from set(monomials_list_non_comm)
+    else:
+        if any(min_degrees[i] > max_degrees[i] for i in range(n)):
+            raise ValueError('min_degrees[i] must be <= max_degrees[i] for all i')
+        power_lists = []
+        for var, min_d, max_d in zip(variables, min_degrees, max_degrees):
+            power_lists.append([var**i for i in range(min_d, max_d + 1)])
+        for powers in product(*power_lists):
+            yield Mul(*powers)
 
 def monomial_count(V, N):
     r"""
diff --git a/sympy/polys/numberfields/subfield.py b/sympy/polys/numberfields/subfield.py
index c56d0662e4..b959ddeb27 100644
--- a/sympy/polys/numberfields/subfield.py
+++ b/sympy/polys/numberfields/subfield.py
@@ -347,12 +347,6 @@ def primitive_element(extension, x=None, *, ex=False, polys=False):
     else:
         x, cls = Dummy('x'), PurePoly
 
-    def _canonicalize(f):
-        _, f = f.primitive()
-        if f.LC() < 0:
-            f = -f
-        return f
-
     if not ex:
         gen, coeffs = extension[0], [1]
         g = minimal_polynomial(gen, x, polys=True)
@@ -366,7 +360,6 @@ def _canonicalize(f):
             gen += s*ext
             coeffs.append(s)
 
-        g = _canonicalize(g)
         if not polys:
             return g.as_expr(), coeffs
         else:
@@ -400,8 +393,6 @@ def _canonicalize(f):
         f = cls(x, domain=QQ)
     else:
         H = [_.to_list() for _ in reps]
-
-    f = _canonicalize(f)
     if not polys:
         return f.as_expr(), coeffs, H
     else:
diff --git a/sympy/series/limits.py b/sympy/series/limits.py
index 9100624963..e15f7a1243 100644
--- a/sympy/series/limits.py
+++ b/sympy/series/limits.py
@@ -3,6 +3,7 @@
 from sympy.core.exprtools import factor_terms
 from sympy.core.numbers import Float, _illegal
 from sympy.core.function import AppliedUndef
+from sympy.core.symbol import Dummy
 from sympy.functions.combinatorial.factorials import factorial
 from sympy.functions.elementary.complexes import (Abs, sign, arg, re)
 from sympy.functions.elementary.exponential import (exp, log)
@@ -326,13 +327,16 @@ def set_signs(expr):
         if z0 is S.Infinity:
             if e.is_Mul:
                 e = factor_terms(e)
-            newe = e.subs(z, 1/z)
+            dummy = Dummy('z', positive=z.is_positive, negative=z.is_negative, real=z.is_real)
+            newe = e.subs(z, 1/dummy)
             # cdir changes sign as oo- should become 0+
             cdir = -cdir
+            newz = dummy
         else:
             newe = e.subs(z, z + z0)
+            newz = z
         try:
-            coeff, ex = newe.leadterm(z, cdir=cdir)
+            coeff, ex = newe.leadterm(newz, cdir=cdir)
         except (ValueError, NotImplementedError, PoleError):
             # The NotImplementedError catching is for custom functions
             from sympy.simplify.powsimp import powsimp
@@ -342,9 +346,9 @@ def set_signs(expr):
                 if r is not None:
                     return r
             try:
-                coeff = newe.as_leading_term(z, cdir=cdir)
+                coeff = newe.as_leading_term(newz, cdir=cdir)
                 if coeff != newe and (coeff.has(exp) or coeff.has(S.Exp1)):
-                    return gruntz(coeff, z, 0, "-" if re(cdir).is_negative else "+")
+                    return gruntz(coeff, newz, 0, "-" if re(cdir).is_negative else "+")
             except (ValueError, NotImplementedError, PoleError):
                 pass
         else:
@@ -352,7 +356,7 @@ def set_signs(expr):
                 return coeff
             if coeff.has(S.Infinity, S.NegativeInfinity, S.ComplexInfinity, S.NaN):
                 return self
-            if not coeff.has(z):
+            if not coeff.has(newz):
                 if ex.is_positive:
                     return S.Zero
                 elif ex == 0:


Test patch:
diff --git a/sympy/algebras/tests/test_quaternion.py b/sympy/algebras/tests/test_quaternion.py
index a4331cd6af..87b68f54b5 100644
--- a/sympy/algebras/tests/test_quaternion.py
+++ b/sympy/algebras/tests/test_quaternion.py
@@ -189,15 +189,6 @@ def test_quaternion_functions():
     assert integrate(Quaternion(x, x, x, x), x) == \
     Quaternion(x**2 / 2, x**2 / 2, x**2 / 2, x**2 / 2)
 
-    assert Quaternion(1, x, x**2, x**3).integrate(x) == \
-    Quaternion(x, x**2/2, x**3/3, x**4/4)
-
-    assert Quaternion(sin(x), cos(x), sin(2*x), cos(2*x)).integrate(x) == \
-    Quaternion(-cos(x), sin(x), -cos(2*x)/2, sin(2*x)/2)
-
-    assert Quaternion(x**2, y**2, z**2, x*y*z).integrate(x, y) == \
-    Quaternion(x**3*y/3, x*y**3/3, x*y*z**2, x**2*y**2*z/4)
-
     assert Quaternion.rotate_point((1, 1, 1), q1) == (S.One / 5, 1, S(7) / 5)
     n = Symbol('n')
     raises(TypeError, lambda: q1**n)
diff --git a/sympy/core/tests/test_arit.py b/sympy/core/tests/test_arit.py
index 251fc4c423..b989651a94 100644
--- a/sympy/core/tests/test_arit.py
+++ b/sympy/core/tests/test_arit.py
@@ -1908,12 +1908,6 @@ def test_Mod():
     assert Mod(-p - 5, -p - 3) == -2
     assert Mod(p + 1, p - 1).func is Mod
 
-    # issue 27749
-    n = symbols('n', integer=True, positive=True)
-    assert unchanged(Mod, 1, n)
-    n = symbols('n', prime=True)
-    assert Mod(1, n) == 1
-
     # handling sums
     assert (x + 3) % 1 == Mod(x, 1)
     assert (x + 3.0) % 1 == Mod(1.*x, 1)
diff --git a/sympy/core/tests/test_priority.py b/sympy/core/tests/test_priority.py
index 276e653100..7be30f2403 100644
--- a/sympy/core/tests/test_priority.py
+++ b/sympy/core/tests/test_priority.py
@@ -2,7 +2,6 @@
 from sympy.core.expr import Expr
 from sympy.core.mod import Mod
 from sympy.core.numbers import Integer
-from sympy.core.singleton import S
 from sympy.core.symbol import Symbol
 from sympy.functions.elementary.integers import floor
 
@@ -15,7 +14,7 @@ class Higher(Integer):
     '''
 
     _op_priority = 20.0
-    result: Expr = S.One
+    result = 1
 
     def __new__(cls):
         obj = Expr.__new__(cls)
@@ -87,7 +86,7 @@ class Lower(Higher):
     '''
 
     _op_priority = 5.0
-    result: Expr = S.NegativeOne
+    result = -1
 
     def __new__(cls):
         obj = Expr.__new__(cls)
diff --git a/sympy/polys/numberfields/tests/test_subfield.py b/sympy/polys/numberfields/tests/test_subfield.py
index b152dd684a..34b93229ce 100644
--- a/sympy/polys/numberfields/tests/test_subfield.py
+++ b/sympy/polys/numberfields/tests/test_subfield.py
@@ -12,7 +12,6 @@
     primitive_element,
     to_number_field,
 )
-from sympy.polys.domains import QQ
 from sympy.polys.polyerrors import IsomorphismFailed
 from sympy.polys.polytools import Poly
 from sympy.polys.rootoftools import CRootOf
@@ -303,15 +302,3 @@ def test_issue_22736():
     a._reset()
     b = exp(2*I*pi/5)
     assert field_isomorphism(a, b) == [1, 0]
-
-
-def test_issue_27798():
-    # https://github.com/sympy/sympy/issues/27798
-    a, b = CRootOf(49*x**3 - 49*x**2 + 14*x - 1, 2), CRootOf(49*x**3 - 49*x**2 + 14*x - 1, 0)
-    assert primitive_element([a, b], polys=True)[0].primitive()[0] == 1
-    assert primitive_element([a, b], polys=True, ex=True)[0].primitive()[0] == 1
-
-    f1, f2 = QQ.algebraic_field(a), QQ.algebraic_field(b)
-    f3 = f1.unify(f2)
-    assert f3.mod.primitive()[0] == 1
-    assert Poly(x, x, domain=f1) + Poly(x, x, domain=f2) == Poly(2*x, x, domain=f3)
diff --git a/sympy/series/tests/test_limits.py b/sympy/series/tests/test_limits.py
index c5d97cbbce..aa3ab7683f 100644
--- a/sympy/series/tests/test_limits.py
+++ b/sympy/series/tests/test_limits.py
@@ -954,7 +954,6 @@ def test_issue_10978():
     assert LambertW(x).limit(x, 0) == 0
 
 
-@XFAIL # issue 27773
 def test_issue_14313_comment():
     assert limit(floor(n/2), n, oo) is oo