Repository: sympy/sympy

Issue:
limit ValueError: The argument 'nan' is not comparable.
I got this exception when calculate `limit((cos((k + 1)/2) - 1)/(cos(k/2) - 1)*x, k, oo)`. Possibly, should give some interval.


Discussion:

Patch:
diff --git a/doc/src/modules/solvers/diophantine.rst b/doc/src/modules/solvers/diophantine.rst
index 06edbfc73f..5afe1f8d0f 100644
--- a/doc/src/modules/solvers/diophantine.rst
+++ b/doc/src/modules/solvers/diophantine.rst
@@ -324,7 +324,7 @@ References
 .. [3] Methods to solve Ax^2 + Bxy + Cy^2 + Dx + Ey + F = 0,[online],
     Available: http://www.alpertron.com.ar/METHODS.HTM
 .. [4] Solving the equation ax^2+ bxy + cy^2 + dx + ey + f= 0, [online],
-    Available: https://web.archive.org/web/20160323033111/http://www.jpr2718.org/ax2p.pdf
+    Available: http://www.jpr2718.org/ax2p.pdf
 
 User Functions
 --------------
diff --git a/sympy/calculus/util.py b/sympy/calculus/util.py
index 15b8042aa1..9f427d9633 100644
--- a/sympy/calculus/util.py
+++ b/sympy/calculus/util.py
@@ -1159,16 +1159,18 @@ def __rsub__(self, other):
 
     @_sympifyit('other', NotImplemented)
     def __mul__(self, other):
+        if self.args == (-oo, oo):
+            return self
         if isinstance(other, Expr):
             if isinstance(other, AccumBounds):
-                return AccumBounds(Min(Mul(self.min, other.min),
-                                       Mul(self.min, other.max),
-                                       Mul(self.max, other.min),
-                                       Mul(self.max, other.max)),
-                                   Max(Mul(self.min, other.min),
-                                       Mul(self.min, other.max),
-                                       Mul(self.max, other.min),
-                                       Mul(self.max, other.max)))
+                if other.args == (-oo, oo):
+                    return other
+                v = set()
+                for i in self.args:
+                    vi = other*i
+                    for i in vi.args or (vi,):
+                        v.add(i)
+                return AccumBounds(Min(*v), Max(*v))
             if other is S.Infinity:
                 if self.min.is_zero:
                     return AccumBounds(0, oo)
@@ -1181,8 +1183,6 @@ def __mul__(self, other):
                     return AccumBounds(0, oo)
             if other.is_extended_real:
                 if other.is_zero:
-                    if self == AccumBounds(-oo, oo):
-                        return AccumBounds(-oo, oo)
                     if self.max is S.Infinity:
                         return AccumBounds(0, oo)
                     if self.min is S.NegativeInfinity:
diff --git a/sympy/physics/vector/frame.py b/sympy/physics/vector/frame.py
index e181a7d367..aa9e21e4b5 100644
--- a/sympy/physics/vector/frame.py
+++ b/sympy/physics/vector/frame.py
@@ -4,8 +4,6 @@
 from sympy.physics.vector.vector import Vector, _check_vector
 from sympy.utilities.misc import translate
 
-from warnings import warn
-
 __all__ = ['CoordinateSym', 'ReferenceFrame']
 
 
@@ -556,25 +554,6 @@ def _dcm(self, parent, parent_orient):
             self._dcm_dict = self._dlist[0] = {}
         # Reset the _dcm_cache
             self._dcm_cache = {}
-
-        else:
-        #Check for loops and raise warning accordingly.
-            visited = []
-            queue = list(frames)
-            cont = True #Flag to control queue loop.
-            while queue and cont:
-                node = queue.pop(0)
-                if node not in visited:
-                    visited.append(node)
-                    neighbors = node._dcm_dict.keys()
-                    for neighbor in neighbors:
-                        if neighbor == parent:
-                            warn('Loops are defined among the orientation of frames.' + \
-                                ' This is likely not desired and may cause errors in your calculations.')
-                            cont = False
-                            break
-                        queue.append(neighbor)
-
         # Add the dcm relationship to _dcm_dict
         self._dcm_dict.update({parent: parent_orient.T})
         parent._dcm_dict.update({self: parent_orient})
@@ -600,12 +579,6 @@ def orient_axis(self, parent, axis, angle):
         angle : sympifiable
             Angle in radians by which it the frame is to be rotated.
 
-        Warns
-        ======
-
-        UserWarning
-            If the orientation creates a kinematic loop.
-
         Examples
         ========
 
@@ -684,12 +657,6 @@ def orient_explicit(self, parent, dcm):
             Direction cosine matrix that specifies the relative rotation
             between the two reference frames.
 
-        Warns
-        ======
-
-        UserWarning
-            If the orientation creates a kinematic loop.
-
         Examples
         ========
 
@@ -794,12 +761,6 @@ def orient_body_fixed(self, parent, angles, rotation_order):
             Tait-Bryan): zxz, xyx, yzy, zyz, xzx, yxy, xyz, yzx, zxy, xzy, zyx,
             and yxz.
 
-        Warns
-        ======
-
-        UserWarning
-            If the orientation creates a kinematic loop.
-
         Examples
         ========
 
@@ -812,7 +773,6 @@ def orient_body_fixed(self, parent, angles, rotation_order):
         >>> B = ReferenceFrame('B')
         >>> B1 = ReferenceFrame('B1')
         >>> B2 = ReferenceFrame('B2')
-        >>> B3 = ReferenceFrame('B3')
 
         For example, a classic Euler Angle rotation can be done by:
 
@@ -830,8 +790,8 @@ def orient_body_fixed(self, parent, angles, rotation_order):
 
         >>> B1.orient_axis(N, N.x, q1)
         >>> B2.orient_axis(B1, B1.y, q2)
-        >>> B3.orient_axis(B2, B2.x, q3)
-        >>> B3.dcm(N)
+        >>> B.orient_axis(B2, B2.x, q3)
+        >>> B.dcm(N)
         Matrix([
         [        cos(q2),                            sin(q1)*sin(q2),                           -sin(q2)*cos(q1)],
         [sin(q2)*sin(q3), -sin(q1)*sin(q3)*cos(q2) + cos(q1)*cos(q3),  sin(q1)*cos(q3) + sin(q3)*cos(q1)*cos(q2)],
@@ -911,12 +871,6 @@ def orient_space_fixed(self, parent, angles, rotation_order):
             ``'131'``, or the integer ``131``. There are 12 unique valid
             rotation orders.
 
-        Warns
-        ======
-
-        UserWarning
-            If the orientation creates a kinematic loop.
-
         Examples
         ========
 
@@ -929,7 +883,6 @@ def orient_space_fixed(self, parent, angles, rotation_order):
         >>> B = ReferenceFrame('B')
         >>> B1 = ReferenceFrame('B1')
         >>> B2 = ReferenceFrame('B2')
-        >>> B3 = ReferenceFrame('B3')
 
         >>> B.orient_space_fixed(N, (q1, q2, q3), '312')
         >>> B.dcm(N)
@@ -942,8 +895,8 @@ def orient_space_fixed(self, parent, angles, rotation_order):
 
         >>> B1.orient_axis(N, N.z, q1)
         >>> B2.orient_axis(B1, N.x, q2)
-        >>> B3.orient_axis(B2, N.y, q3)
-        >>> B3.dcm(N).simplify()
+        >>> B.orient_axis(B2, N.y, q3)
+        >>> B.dcm(N).simplify()
         Matrix([
         [ sin(q1)*sin(q2)*sin(q3) + cos(q1)*cos(q3), sin(q1)*cos(q2), sin(q1)*sin(q2)*cos(q3) - sin(q3)*cos(q1)],
         [-sin(q1)*cos(q3) + sin(q2)*sin(q3)*cos(q1), cos(q1)*cos(q2), sin(q1)*sin(q3) + sin(q2)*cos(q1)*cos(q3)],
@@ -1039,12 +992,6 @@ def orient_quaternion(self, parent, numbers):
             The four quaternion scalar numbers as defined above: ``q0``,
             ``q1``, ``q2``, ``q3``.
 
-        Warns
-        ======
-
-        UserWarning
-            If the orientation creates a kinematic loop.
-
         Examples
         ========
 
@@ -1151,12 +1098,6 @@ def orient(self, parent, rot_type, amounts, rot_order=''):
             ``'123'`` and integer ``123`` are equivalent, for example. Required
             for ``'Body'`` and ``'Space'``.
 
-        Warns
-        ======
-
-        UserWarning
-            If the orientation creates a kinematic loop.
-
         """
 
         _check_frame(parent)
diff --git a/sympy/solvers/diophantine/diophantine.py b/sympy/solvers/diophantine/diophantine.py
index 22967b709f..b9a47ede31 100644
--- a/sympy/solvers/diophantine/diophantine.py
+++ b/sympy/solvers/diophantine/diophantine.py
@@ -470,7 +470,7 @@ class BinaryQuadratic(DiophantineEquationType):
     .. [1] Methods to solve Ax^2 + Bxy + Cy^2 + Dx + Ey + F = 0, [online],
           Available: http://www.alpertron.com.ar/METHODS.HTM
     .. [2] Solving the equation ax^2+ bxy + cy^2 + dx + ey + f= 0, [online],
-          Available: https://web.archive.org/web/20160323033111/http://www.jpr2718.org/ax2p.pdf
+          Available: http://www.jpr2718.org/ax2p.pdf
 
     """
 
@@ -570,7 +570,7 @@ def solve(self, parameters=None, limit=None) -> DiophantineSolutionSet:
 
         # (3) Method used when B**2 - 4*A*C is a square, is described in p. 6 of the below paper
         # by John P. Robertson.
-        # https://web.archive.org/web/20160323033111/http://www.jpr2718.org/ax2p.pdf
+        # http://www.jpr2718.org/ax2p.pdf
 
         elif is_square(discr):
             if A != 0:
@@ -1923,7 +1923,7 @@ def diop_quadratic(eq, param=symbols("t", integer=True)):
     .. [1] Methods to solve Ax^2 + Bxy + Cy^2 + Dx + Ey + F = 0, [online],
           Available: http://www.alpertron.com.ar/METHODS.HTM
     .. [2] Solving the equation ax^2+ bxy + cy^2 + dx + ey + f= 0, [online],
-          Available: https://web.archive.org/web/20160323033111/http://www.jpr2718.org/ax2p.pdf
+          Available: http://www.jpr2718.org/ax2p.pdf
 
     See Also
     ========
@@ -2006,7 +2006,7 @@ def diop_DN(D, N, t=symbols("t", integer=True)):
 
     .. [1] Solving the generalized Pell equation x**2 - D*y**2 = N, John P.
         Robertson, July 31, 2004, Pages 16 - 17. [online], Available:
-        https://web.archive.org/web/20160323033128/http://www.jpr2718.org/pell.pdf
+        http://www.jpr2718.org/pell.pdf
     """
     if D < 0:
         if N == 0:
@@ -2349,7 +2349,7 @@ def PQa(P_0, Q_0, D):
     ==========
 
     .. [1] Solving the generalized Pell equation x^2 - Dy^2 = N, John P.
-        Robertson, July 31, 2004, Pages 4 - 8. https://web.archive.org/web/20160323033128/http://www.jpr2718.org/pell.pdf
+        Robertson, July 31, 2004, Pages 4 - 8. http://www.jpr2718.org/pell.pdf
     """
     A_i_2 = B_i_1 = 0
     A_i_1 = B_i_2 = 1
@@ -2420,7 +2420,7 @@ def diop_bf_DN(D, N, t=symbols("t", integer=True)):
     ==========
 
     .. [1] Solving the generalized Pell equation x**2 - D*y**2 = N, John P.
-        Robertson, July 31, 2004, Page 15. https://web.archive.org/web/20160323033128/http://www.jpr2718.org/pell.pdf
+        Robertson, July 31, 2004, Page 15. http://www.jpr2718.org/pell.pdf
     """
     D = as_int(D)
     N = as_int(N)
@@ -2501,7 +2501,7 @@ def equivalent(u, v, r, s, D, N):
     ==========
 
     .. [1] Solving the generalized Pell equation x**2 - D*y**2 = N, John P.
-        Robertson, July 31, 2004, Page 12. https://web.archive.org/web/20160323033128/http://www.jpr2718.org/pell.pdf
+        Robertson, July 31, 2004, Page 12. http://www.jpr2718.org/pell.pdf
 
     """
     return divisible(u*r - D*v*s, N) and divisible(u*s - v*r, N)
@@ -2627,7 +2627,7 @@ def transformation_to_DN(eq):
 
     .. [1] Solving the equation ax^2 + bxy + cy^2 + dx + ey + f = 0,
            John P.Robertson, May 8, 2003, Page 7 - 11.
-           https://web.archive.org/web/20160323033111/http://www.jpr2718.org/ax2p.pdf
+           http://www.jpr2718.org/ax2p.pdf
     """
 
     var, coeff, diop_type = classify_diop(eq, _dict=False)
@@ -2723,7 +2723,7 @@ def find_DN(eq):
 
     .. [1] Solving the equation ax^2 + bxy + cy^2 + dx + ey + f = 0,
            John P.Robertson, May 8, 2003, Page 7 - 11.
-           https://web.archive.org/web/20160323033111/http://www.jpr2718.org/ax2p.pdf
+           http://www.jpr2718.org/ax2p.pdf
     """
     var, coeff, diop_type = classify_diop(eq, _dict=False)
     if diop_type == BinaryQuadratic.name:


Test patch:
diff --git a/sympy/calculus/tests/test_util.py b/sympy/calculus/tests/test_util.py
index 41514e7210..36d660df5d 100644
--- a/sympy/calculus/tests/test_util.py
+++ b/sympy/calculus/tests/test_util.py
@@ -359,7 +359,13 @@ def test_AccumBounds_mul():
     assert AccumBounds(1, 2)*2 == AccumBounds(2, 4)
     assert 2*AccumBounds(1, 2) == AccumBounds(2, 4)
     assert AccumBounds(1, 2)*AccumBounds(2, 3) == AccumBounds(2, 6)
-
+    assert AccumBounds(0, 2)*AccumBounds(2, oo) == AccumBounds(0, oo)
+    l, r = AccumBounds(-oo, oo), AccumBounds(-a, a)
+    assert l*r == AccumBounds(-oo, oo)
+    assert r*l == AccumBounds(-oo, oo)
+    l, r = AccumBounds(1, oo), AccumBounds(-3, -2)
+    assert l*r == AccumBounds(-oo, -2)
+    assert r*l == AccumBounds(-oo, -2)
     assert AccumBounds(1, 2)*0 == 0
     assert AccumBounds(1, oo)*0 == AccumBounds(0, oo)
     assert AccumBounds(-oo, 1)*0 == AccumBounds(-oo, 0)
diff --git a/sympy/physics/vector/tests/test_frame.py b/sympy/physics/vector/tests/test_frame.py
index 819ba0cb54..61fe2030bb 100644
--- a/sympy/physics/vector/tests/test_frame.py
+++ b/sympy/physics/vector/tests/test_frame.py
@@ -6,7 +6,6 @@
 from sympy.physics.vector.frame import _check_frame
 from sympy.physics.vector.vector import VectorTypeError
 from sympy.testing.pytest import raises
-import warnings
 
 Vector.simp = True
 
@@ -473,22 +472,6 @@ def test_orient_quaternion():
     B.orient_quaternion(A, (0,0,0,0))
     assert B.dcm(A) == Matrix([[0, 0, 0], [0, 0, 0], [0, 0, 0]])
 
-def test_looped_frame_warning():
-    A = ReferenceFrame('A')
-    B = ReferenceFrame('B')
-    C = ReferenceFrame('C')
-
-    a, b, c = symbols('a b c')
-    B.orient_axis(A, A.x, a)
-    C.orient_axis(B, B.x, b)
-
-    with warnings.catch_warnings(record = True) as w:
-        warnings.simplefilter("always")
-        A.orient_axis(C, C.x, c)
-        assert issubclass(w[-1].category, UserWarning)
-        assert 'Loops are defined among the orientation of frames. ' + \
-            'This is likely not desired and may cause errors in your calculations.' in str(w[-1].message)
-
 def test_frame_dict():
     A = ReferenceFrame('A')
     B = ReferenceFrame('B')
diff --git a/sympy/solvers/diophantine/tests/test_diophantine.py b/sympy/solvers/diophantine/tests/test_diophantine.py
index 26a9280710..1f97568928 100644
--- a/sympy/solvers/diophantine/tests/test_diophantine.py
+++ b/sympy/solvers/diophantine/tests/test_diophantine.py
@@ -195,7 +195,7 @@ def test_quadratic_non_perfect_slow():
 def test_DN():
     # Most of the test cases were adapted from,
     # Solving the generalized Pell equation x**2 - D*y**2 = N, John P. Robertson, July 31, 2004.
-    # https://web.archive.org/web/20160323033128/http://www.jpr2718.org/pell.pdf
+    # http://www.jpr2718.org/pell.pdf
     # others are verified using Wolfram Alpha.
 
     # Covers cases where D <= 0 or D > 0 and D is a square or N = 0