Repository: sympy/sympy

Issue:
RecursionError with Piecewise.xreplace containing Indexed objects
Minimal example to reproduce the issue:

```python
from sympy import Piecewise, Symbol, ITE, MatrixSymbol

x = Symbol('x')
y = Symbol('y')

inp = MatrixSymbol('inp', 1, 1)
rep_dict = {y: inp[0, 0], x: inp[0, 0]}

p = Piecewise((1, ITE(y > 0, x < 0, True)))
z = p.xreplace(rep_dict)
```

```
---------------------------------------------------------------------------
RecursionError                            Traceback (most recent call last)
<ipython-input-44-46c498a48320> in <module>()
----> 1 z = p.xreplace(rep_dict)
      2 print(z)

/home/rotis/anaconda/envs/calphadpy3/lib/python3.5/site-packages/sympy/core/basic.py in xreplace(self, rule)
   1113 
   1114         """
-> 1115         value, _ = self._xreplace(rule)
   1116         return value
   1117 

/home/rotis/anaconda/envs/calphadpy3/lib/python3.5/site-packages/sympy/core/basic.py in _xreplace(self, rule)
   1134             args = tuple(args)
   1135             if changed:
-> 1136                 return self.func(*args), True
   1137         return self, False
   1138 

/home/rotis/anaconda/envs/calphadpy3/lib/python3.5/site-packages/sympy/functions/elementary/piecewise.py in __new__(cls, *args, **options)
    134 
    135         if options.pop('evaluate', True):
--> 136             r = cls.eval(*newargs)
    137         else:
    138             r = None

/home/rotis/anaconda/envs/calphadpy3/lib/python3.5/site-packages/sympy/functions/elementary/piecewise.py in eval(cls, *_args)
    183                     x = free.pop()
    184                     try:
--> 185                         c = c.as_set().as_relational(x)
    186                     except NotImplementedError:
    187                         pass

/home/rotis/anaconda/envs/calphadpy3/lib/python3.5/site-packages/sympy/logic/boolalg.py in as_set(self)
    144             reps = {}
    145             for r in self.atoms(Relational):
--> 146                 if periodicity(r, x) not in (0, None):
    147                     s = r._eval_as_set()
    148                     if s in (S.EmptySet, S.UniversalSet, S.Reals):

/home/rotis/anaconda/envs/calphadpy3/lib/python3.5/site-packages/sympy/calculus/util.py in periodicity(f, symbol, check)
    471     elif period is None:
    472         from sympy.solvers.decompogen import compogen
--> 473         g_s = decompogen(f, symbol)
    474         num_of_gs = len(g_s)
    475         if num_of_gs > 1:

/home/rotis/anaconda/envs/calphadpy3/lib/python3.5/site-packages/sympy/solvers/decompogen.py in decompogen(f, symbol)
     55         f1 = f.subs(gens[0], symbol)
     56         f2 = gens[0]
---> 57         result += [f1] + decompogen(f2, symbol)
     58         return result
     59 

/home/rotis/anaconda/envs/calphadpy3/lib/python3.5/site-packages/sympy/solvers/decompogen.py in decompogen(f, symbol)
     55         f1 = f.subs(gens[0], symbol)
     56         f2 = gens[0]
---> 57         result += [f1] + decompogen(f2, symbol)
     58         return result
     59 

/home/rotis/anaconda/envs/calphadpy3/lib/python3.5/site-packages/sympy/solvers/decompogen.py in decompogen(f, symbol)
     55         f1 = f.subs(gens[0], symbol)
     56         f2 = gens[0]
---> 57         result += [f1] + decompogen(f2, symbol)
     58         return result
     59 

... last 3 frames repeated, from the frame below ...

/home/rotis/anaconda/envs/calphadpy3/lib/python3.5/site-packages/sympy/solvers/decompogen.py in decompogen(f, symbol)
     55         f1 = f.subs(gens[0], symbol)
     56         f2 = gens[0]
---> 57         result += [f1] + decompogen(f2, symbol)
     58         return result
     59 

RecursionError: maximum recursion depth exceeded
```

Note that this will occur even if you pass `evaluate=False` to `Piecewise`.
This may be related to #14899

Discussion:

Patch:
diff --git a/sympy/calculus/util.py b/sympy/calculus/util.py
index 1f36402efe..e604985c60 100644
--- a/sympy/calculus/util.py
+++ b/sympy/calculus/util.py
@@ -4,8 +4,9 @@
 from sympy.core.compatibility import iterable
 from sympy.core.expr import AtomicExpr, Expr
 from sympy.core.function import expand_mul
+from sympy.core.kind import NumberKind
 from sympy.core.numbers import _sympifyit, oo, zoo
-from sympy.core.relational import is_le, is_lt, is_ge, is_gt
+from sympy.core.relational import is_le, is_lt, is_ge, is_gt, Relational
 from sympy.core.sympify import _sympify
 from sympy.functions.elementary.miscellaneous import Min, Max
 from sympy.logic.boolalg import And
@@ -392,7 +393,6 @@ def periodicity(f, symbol, check=False):
     >>> periodicity(exp(x), x)
     """
     from sympy.core.mod import Mod
-    from sympy.core.relational import Relational
     from sympy.functions.elementary.exponential import exp
     from sympy.functions.elementary.complexes import Abs
     from sympy.functions.elementary.trigonometric import (
@@ -401,6 +401,8 @@ def periodicity(f, symbol, check=False):
     from sympy.solvers.decompogen import decompogen
     from sympy.polys.polytools import degree
 
+    if symbol.kind is not NumberKind:
+        raise NotImplementedError("Cannot use symbol of kind %s" % symbol.kind)
     temp = Dummy('x', real=True)
     f = f.subs(symbol, temp)
     symbol = temp
diff --git a/sympy/logic/boolalg.py b/sympy/logic/boolalg.py
index 404acbc784..bbe46ca356 100644
--- a/sympy/logic/boolalg.py
+++ b/sympy/logic/boolalg.py
@@ -10,6 +10,7 @@
 from sympy.core.compatibility import ordered, as_int
 from sympy.core.decorators import sympify_method_args, sympify_return
 from sympy.core.function import Application, Derivative
+from sympy.core.kind import NumberKind
 from sympy.core.numbers import Number
 from sympy.core.operations import LatticeOp
 from sympy.core.singleton import Singleton, S
@@ -155,18 +156,20 @@ def as_set(self):
         free = self.free_symbols
         if len(free) == 1:
             x = free.pop()
-            reps = {}
-            for r in self.atoms(Relational):
-                if periodicity(r, x) not in (0, None):
-                    s = r._eval_as_set()
-                    if s in (S.EmptySet, S.UniversalSet, S.Reals):
-                        reps[r] = s.as_relational(x)
-                        continue
-                    raise NotImplementedError(filldedent('''
-                        as_set is not implemented for relationals
-                        with periodic solutions
-                        '''))
-            return self.subs(reps)._eval_as_set()
+            if x.kind is NumberKind:
+                reps = {}
+                for r in self.atoms(Relational):
+                    if periodicity(r, x) not in (0, None):
+                        s = r._eval_as_set()
+                        if s in (S.EmptySet, S.UniversalSet, S.Reals):
+                            reps[r] = s.as_relational(x)
+                            continue
+                        raise NotImplementedError(filldedent('''
+                            as_set is not implemented for relationals
+                            with periodic solutions
+                            '''))
+                return self.subs(reps)._eval_as_set()
+            return self._eval_as_set()
         else:
             raise NotImplementedError("Sorry, as_set has not yet been"
                                       " implemented for multivariate"


Test patch:
diff --git a/sympy/calculus/tests/test_util.py b/sympy/calculus/tests/test_util.py
index e7add2c894..758ca6f5d6 100644
--- a/sympy/calculus/tests/test_util.py
+++ b/sympy/calculus/tests/test_util.py
@@ -1,6 +1,6 @@
 from sympy import (Symbol, S, exp, log, sqrt, oo, E, zoo, pi, tan, sin, cos,
                    cot, sec, csc, Abs, symbols, I, re, simplify,
-                   expint, Rational, Piecewise)
+                   expint, Rational, Piecewise, MatrixSymbol)
 from sympy.calculus.util import (function_range, continuous_domain, not_empty_in,
                                  periodicity, lcim, AccumBounds, is_convex,
                                  stationary_points, minimum, maximum)
@@ -178,6 +178,12 @@ def test_periodicity():
     p = Piecewise((0, x < -1), (x**2, x <= 1), (log(x), True))
     assert periodicity(p, x) is None
 
+    m = MatrixSymbol('m', 3, 3)
+    raises(NotImplementedError, lambda: periodicity(sin(m), m))
+    raises(NotImplementedError, lambda: periodicity(sin(m[0, 0]), m))
+    raises(NotImplementedError, lambda: periodicity(sin(m), m[0, 0]))
+    raises(NotImplementedError, lambda: periodicity(sin(m[0, 0]), m[0, 0]))
+
 
 def test_periodicity_check():
     x = Symbol('x')
@@ -635,6 +641,7 @@ def test_issue_16469():
     f = abs(x)
     assert function_range(f, x, S.Reals) == Interval(0, oo, False, True)
 
+
 @_both_exp_pow
 def test_issue_18747():
     assert periodicity(exp(pi*I*(x/4+S.Half/2)), x) == 8
diff --git a/sympy/functions/elementary/tests/test_piecewise.py b/sympy/functions/elementary/tests/test_piecewise.py
index 0701bf2752..85b68689a7 100644
--- a/sympy/functions/elementary/tests/test_piecewise.py
+++ b/sympy/functions/elementary/tests/test_piecewise.py
@@ -3,7 +3,7 @@
     Integral, integrate, Interval, KroneckerDelta, lambdify, log, Max, Min,
     oo, Or, pi, Piecewise, piecewise_fold, Rational, solve, symbols, transpose,
     cos, sin, exp, Abs, Ne, Not, Symbol, S, sqrt, Sum, Tuple, zoo, Float,
-    DiracDelta, Heaviside, Add, Mul, factorial, Ge, Contains)
+    DiracDelta, Heaviside, Add, Mul, factorial, Ge, Contains, MatrixSymbol)
 from sympy.core.expr import unchanged
 from sympy.functions.elementary.piecewise import Undefined, ExprCondPair
 from sympy.printing import srepr
@@ -1353,9 +1353,22 @@ def test_issue_7370():
     v = integrate(f, (x, 0, Float("252.4", 30)))
     assert str(v) == '252.400000000000000000000000000'
 
+
+def test_issue_14933():
+    x = Symbol('x')
+    y = Symbol('y')
+
+    inp = MatrixSymbol('inp', 1, 1)
+    rep_dict = {y: inp[0, 0], x: inp[0, 0]}
+
+    p = Piecewise((1, ITE(y > 0, x < 0, True)))
+    assert p.xreplace(rep_dict) == Piecewise((1, ITE(inp[0, 0] > 0, inp[0, 0] < 0, True)))
+
+
 def test_issue_16715():
     raises(NotImplementedError, lambda: Piecewise((x, x<0), (0, y>1)).as_expr_set_pairs())
 
+
 def test_issue_20360():
     t, tau = symbols("t tau", real=True)
     n = symbols("n", integer=True)
diff --git a/sympy/series/tests/test_gruntz.py b/sympy/series/tests/test_gruntz.py
index d1d9aa5b61..f9be75352c 100644
--- a/sympy/series/tests/test_gruntz.py
+++ b/sympy/series/tests/test_gruntz.py
@@ -129,12 +129,12 @@ def test_grunts_eval_special_slow_sometimes_fail():
     assert gruntz(exp(gamma(x - exp(-x))*exp(1/x)) - exp(gamma(x)), x, oo) is oo
 
 
-def test_gruntz_Ei():
-    assert gruntz((Ei(x - exp(-exp(x))) - Ei(x)) *exp(-x)*exp(exp(x))*x, x, oo) == -1
-
-
 @XFAIL
 def test_gruntz_eval_special_fail():
+    # TODO exponential integral Ei
+    assert gruntz(
+        (Ei(x - exp(-exp(x))) - Ei(x)) *exp(-x)*exp(exp(x))*x, x, oo) == -1
+
     # TODO zeta function series
     assert gruntz(
         exp((log(2) + 1)*x) * (zeta(x + exp(-x)) - zeta(x)), x, oo) == -log(2)
diff --git a/sympy/sets/tests/test_setexpr.py b/sympy/sets/tests/test_setexpr.py
index ab8077a4fe..e6490e7757 100644
--- a/sympy/sets/tests/test_setexpr.py
+++ b/sympy/sets/tests/test_setexpr.py
@@ -2,7 +2,7 @@
 from sympy.sets import Interval, FiniteSet, Intersection, ImageSet, Union
 
 from sympy import (Expr, Set, exp, log, cos, Symbol, Min, Max, S, oo, I,
-        symbols, Lambda, Dummy, Rational, sqrt)
+        symbols, Lambda, Dummy, Rational)
 
 
 a, x = symbols("a, x")
@@ -216,7 +216,7 @@ def test_SetExpr_Interval_div():
 
     assert 1/SetExpr(Interval(0, 2)) == SetExpr(Interval(S.Half, oo))
     assert (-1)/SetExpr(Interval(0, 2)) == SetExpr(Interval(-oo, Rational(-1, 2)))
-    assert 1/SetExpr(Interval(-oo, 0)) == SetExpr(Interval.open(-oo, 0))
+    # assert 1/SetExpr(Interval(-oo, 0)) == SetExpr(Interval.open(-oo, 0))
     assert 1/SetExpr(Interval(-1, 0)) == SetExpr(Interval(-oo, -1))
     # assert (-2)/SetExpr(Interval(-oo, 0)) == SetExpr(Interval(0, oo))
     # assert 1/SetExpr(Interval(-oo, -1)) == SetExpr(Interval(-1, 0))
@@ -244,7 +244,7 @@ def test_SetExpr_Interval_pow():
     assert SetExpr(Interval(-1, 1))**0 == SetExpr(FiniteSet(1))
 
 
-    assert SetExpr(Interval(1, 2))**Rational(5, 2) == SetExpr(Interval(1, 4*sqrt(2)))
+    #assert SetExpr(Interval(1, 2))**Rational(5, 2) == SetExpr(Interval(1, 4*sqrt(2)))
     #assert SetExpr(Interval(-1, 2))**Rational(1, 3) == SetExpr(Interval(-1, 2**Rational(1, 3)))
     #assert SetExpr(Interval(0, 2))**S.Half == SetExpr(Interval(0, sqrt(2)))