Repository: sympy/sympy

Issue:
Simplify Mul of KroneckerDelta
There should be a way to simplify this to zero:
```julia
In [56]: n = Symbol('n', integer=True)                                                                                            

In [57]: expr = KroneckerDelta(0, n) * KroneckerDelta(1, n)                                                                       

In [58]: expr                                                                                                                     
Out[58]: 
δ   ⋅δ   
 0,n  1,n

In [59]: simplify(expr)                                                                                                           
Out[59]: 
δ   ⋅δ   
 0,n  1,n
```

Discussion:

Patch:
diff --git a/doc/src/modules/simplify/simplify.rst b/doc/src/modules/simplify/simplify.rst
index 7f7db11363..5b15a65b33 100644
--- a/doc/src/modules/simplify/simplify.rst
+++ b/doc/src/modules/simplify/simplify.rst
@@ -15,6 +15,10 @@ nthroot
 -------
 .. autofunction:: nthroot
 
+kroneckersimp
+-------------
+.. autofunction:: kroneckersimp
+
 besselsimp
 ----------
 .. autofunction:: besselsimp
diff --git a/sympy/functions/elementary/piecewise.py b/sympy/functions/elementary/piecewise.py
index 9f5bb61203..f62108a332 100644
--- a/sympy/functions/elementary/piecewise.py
+++ b/sympy/functions/elementary/piecewise.py
@@ -1029,6 +1029,65 @@ def _eval_rewrite_as_ITE(self, *args, **kwargs):
             last = ITE(c, a, last)
         return _canonical(last)
 
+    def _eval_rewrite_as_KroneckerDelta(self, *args):
+        from sympy import Ne, Eq, Not, KroneckerDelta
+
+        rules = {
+            And: [False, False],
+            Or: [True, True],
+            Not: [True, False],
+            Eq: [None, None],
+            Ne: [None, None]
+        }
+
+        class UnrecognizedCondition(Exception):
+            pass
+
+        def rewrite(cond):
+            if isinstance(cond, Eq):
+                return KroneckerDelta(*cond.args)
+            if isinstance(cond, Ne):
+                return 1 - KroneckerDelta(*cond.args)
+
+            cls, args = type(cond), cond.args
+            if cls not in rules:
+                raise UnrecognizedCondition(cls)
+
+            b1, b2 = rules[cls]
+            k = 1
+            for c in args:
+                if b1:
+                    k *= 1 - rewrite(c)
+                else:
+                    k *= rewrite(c)
+
+            if b2:
+                return 1 - k
+            return k
+
+        conditions = []
+        true_value = None
+        for value, cond in args:
+            if type(cond) in rules:
+                conditions.append((value, cond))
+            elif cond is S.true:
+                if true_value is None:
+                    true_value = value
+            else:
+                return
+
+        if true_value is not None:
+            result = true_value
+
+            for value, cond in conditions[::-1]:
+                try:
+                    k = rewrite(cond)
+                    result = k * value + (1 - k) * result
+                except UnrecognizedCondition:
+                    return
+
+            return result
+
 
 def piecewise_fold(expr):
     """
diff --git a/sympy/functions/special/tensor_functions.py b/sympy/functions/special/tensor_functions.py
index 66efe8c857..18a056d47d 100644
--- a/sympy/functions/special/tensor_functions.py
+++ b/sympy/functions/special/tensor_functions.py
@@ -444,3 +444,9 @@ def indices(self):
     def _sage_(self):
         import sage.all as sage
         return sage.kronecker_delta(self.args[0]._sage_(), self.args[1]._sage_())
+
+    def _eval_rewrite_as_Piecewise(self, *args, **kwargs):
+        from sympy.functions.elementary.piecewise import Piecewise
+        from sympy.core.relational import Ne
+        i, j = args
+        return Piecewise((0, Ne(i, j)), (1, True))
diff --git a/sympy/simplify/__init__.py b/sympy/simplify/__init__.py
index 1e30271a21..15da10c9aa 100644
--- a/sympy/simplify/__init__.py
+++ b/sympy/simplify/__init__.py
@@ -5,7 +5,7 @@
 the expression (x+x)**2 will be converted into 4*x**2
 """
 from .simplify import (simplify, hypersimp, hypersimilar,
-    logcombine, separatevars, posify, besselsimp,
+    logcombine, separatevars, posify, besselsimp, kroneckersimp,
     signsimp, bottom_up, nsimplify)
 
 from .fu import FU, fu
diff --git a/sympy/simplify/simplify.py b/sympy/simplify/simplify.py
index a64187774d..037f9d8df0 100644
--- a/sympy/simplify/simplify.py
+++ b/sympy/simplify/simplify.py
@@ -4,7 +4,7 @@
 
 from sympy.core import (Basic, S, Add, Mul, Pow, Symbol, sympify, expand_mul,
                         expand_func, Function, Dummy, Expr, factor_terms,
-                        expand_power_exp)
+                        expand_power_exp, Eq)
 from sympy.core.compatibility import iterable, ordered, range, as_int
 from sympy.core.evaluate import global_evaluate
 from sympy.core.function import expand_log, count_ops, _mexpand, _coeff_isneg, nfloat
@@ -17,8 +17,10 @@
 from sympy.functions.elementary.exponential import ExpBase
 from sympy.functions.elementary.hyperbolic import HyperbolicFunction
 from sympy.functions.elementary.integers import ceiling
+from sympy.functions.elementary.piecewise import Piecewise
 from sympy.functions.elementary.trigonometric import TrigonometricFunction
 from sympy.functions.special.bessel import besselj, besseli, besselk, jn, bessely
+from sympy.functions.special.tensor_functions import KroneckerDelta
 from sympy.polys import together, cancel, factor
 from sympy.simplify.combsimp import combsimp
 from sympy.simplify.cse_opts import sub_pre, sub_post
@@ -590,6 +592,9 @@ def shorter(*choices):
 
     expr = piecewise_fold(expr)
 
+    if expr.has(KroneckerDelta):
+        expr = kroneckersimp(expr)
+
     if expr.has(BesselBase):
         expr = besselsimp(expr)
 
@@ -1101,6 +1106,56 @@ def bottom_up(rv, F, atoms=False, nonbasic=False):
     return rv
 
 
+def kroneckersimp(expr):
+    """
+    Simplify expressions with KroneckerDelta.
+
+    The only simplification currently attempted is to identify multiplicative cancellation:
+
+    >>> from sympy import KroneckerDelta, kroneckersimp
+    >>> from sympy.abc import i, j
+    >>> kroneckersimp(1 + KroneckerDelta(0, j) * KroneckerDelta(1, j))
+    1
+    """
+    def args_cancel(args1, args2):
+        for i1 in range(2):
+            for i2 in range(2):
+                a1 = args1[i1]
+                a2 = args2[i2]
+                a3 = args1[(i1 + 1) % 2]
+                a4 = args2[(i2 + 1) % 2]
+                if Eq(a1, a2) is S.true and Eq(a3, a4) is S.false:
+                    return True
+        return False
+
+    def cancel_kronecker_mul(m):
+        from sympy.utilities.iterables import subsets
+
+        args = m.args
+        deltas = [a for a in args if isinstance(a, KroneckerDelta)]
+        for delta1, delta2 in subsets(deltas, 2):
+            args1 = delta1.args
+            args2 = delta2.args
+            if args_cancel(args1, args2):
+                return 0*m
+        return m
+
+    if not expr.has(KroneckerDelta):
+        return expr
+
+    if expr.has(Piecewise):
+        expr = expr.rewrite(KroneckerDelta)
+
+    newexpr = expr
+    expr = None
+
+    while newexpr != expr:
+        expr = newexpr
+        newexpr = expr.replace(lambda e: isinstance(e, Mul), cancel_kronecker_mul)
+
+    return expr
+
+
 def besselsimp(expr):
     """
     Simplify bessel-type functions.


Test patch:
diff --git a/sympy/functions/elementary/tests/test_piecewise.py b/sympy/functions/elementary/tests/test_piecewise.py
index 8596580f32..81b9b39b9a 100644
--- a/sympy/functions/elementary/tests/test_piecewise.py
+++ b/sympy/functions/elementary/tests/test_piecewise.py
@@ -1,7 +1,7 @@
 from sympy import (
     adjoint, And, Basic, conjugate, diff, expand, Eq, Function, I, ITE,
-    Integral, integrate, Interval, lambdify, log, Max, Min, oo, Or, pi,
-    Piecewise, piecewise_fold, Rational, solve, symbols, transpose,
+    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,
     DiracDelta, Heaviside, Add, Mul, factorial, Ge, Contains, Le)
 from sympy.core.expr import unchanged
@@ -1226,3 +1226,79 @@ def test_issue_8458():
     # Test for problem highlighted during review
     p3 = Piecewise((x+1, Eq(x, -1)), (4*x + (y-2)**4, Eq(x, 0) & Eq(x+y, 2)), (sin(x), True))
     assert p3.simplify() == Piecewise((0, Eq(x, -1)), (sin(x), True))
+
+
+def test_eval_rewrite_as_KroneckerDelta():
+    x, y, z, n, t, m = symbols('x y z n t m')
+    K = KroneckerDelta
+    f = lambda p: expand(p.rewrite(K))
+
+    p1 = Piecewise((0, Eq(x, y)), (1, True))
+    assert f(p1) == 1 - K(x, y)
+
+    p2 = Piecewise((x, Eq(y,0)), (z, Eq(t,0)), (n, True))
+    assert f(p2) == n*K(0, t)*K(0, y) - n*K(0, t) - n*K(0, y) + n + \
+           x*K(0, y) - z*K(0, t)*K(0, y) + z*K(0, t)
+
+    p3 = Piecewise((1, Ne(x, y)), (0, True))
+    assert f(p3) == 1 - K(x, y)
+
+    p4 = Piecewise((1, Eq(x, 3)), (4, True), (5, True))
+    assert f(p4) == 4 - 3*K(3, x)
+
+    p5 = Piecewise((3, Ne(x, 2)), (4, Eq(y, 2)), (5, True))
+    assert f(p5) == -K(2, x)*K(2, y) + 2*K(2, x) + 3
+
+    p6 = Piecewise((0, Ne(x, 1) & Ne(y, 4)), (1, True))
+    assert f(p6) == -K(1, x)*K(4, y) + K(1, x) + K(4, y)
+
+    p7 = Piecewise((2, Eq(y, 3) & Ne(x, 2)), (1, True))
+    assert f(p7) == -K(2, x)*K(3, y) + K(3, y) + 1
+
+    p8 = Piecewise((4, Eq(x, 3) & Ne(y, 2)), (1, True))
+    assert f(p8) == -3*K(2, y)*K(3, x) + 3*K(3, x) + 1
+
+    p9 = Piecewise((6, Eq(x, 4) & Eq(y, 1)), (1, True))
+    assert f(p9) == 5 * K(1, y) * K(4, x) + 1
+
+    p10 = Piecewise((4, Ne(x, -4) | Ne(y, 1)), (1, True))
+    assert f(p10) == -3 * K(-4, x) * K(1, y) + 4
+
+    p11 = Piecewise((1, Eq(y, 2) | Ne(x, -3)), (2, True))
+    assert f(p11) == -K(-3, x)*K(2, y) + K(-3, x) + 1
+
+    p12 = Piecewise((-1, Eq(x, 1) | Ne(y, 3)), (1, True))
+    assert f(p12) == -2*K(1, x)*K(3, y) + 2*K(3, y) - 1
+
+    p13 = Piecewise((3, Eq(x, 2) | Eq(y, 4)), (1, True))
+    assert f(p13) == -2*K(2, x)*K(4, y) + 2*K(2, x) + 2*K(4, y) + 1
+
+    p14 = Piecewise((1, Ne(x, 0) | Ne(y, 1)), (3, True))
+    assert f(p14) == 2 * K(0, x) * K(1, y) + 1
+
+    p15 = Piecewise((2, Eq(x, 3) | Ne(y, 2)), (3, Eq(x, 4) & Eq(y, 5)), (1, True))
+    assert f(p15) == -2*K(2, y)*K(3, x)*K(4, x)*K(5, y) + K(2, y)*K(3, x) + \
+           2*K(2, y)*K(4, x)*K(5, y) - K(2, y) + 2
+
+    p16 = Piecewise((0, Ne(m, n)), (1, True))*Piecewise((0, Ne(n, t)), (1, True))\
+          *Piecewise((0, Ne(n, x)), (1, True)) - Piecewise((0, Ne(t, x)), (1, True))
+    assert f(p16) == K(m, n)*K(n, t)*K(n, x) - K(t, x)
+
+    p17 = Piecewise((0, Ne(t, x) & (Ne(m, n) | Ne(n, t) | Ne(n, x))),
+                    (1, Ne(t, x)), (-1, Ne(m, n) | Ne(n, t) | Ne(n, x)), (0, True))
+    assert f(p17) == K(m, n)*K(n, t)*K(n, x) - K(t, x)
+
+    p18 = Piecewise((-4, Eq(y, 1) | (Eq(x, -5) & Eq(x, z))), (4, True))
+    assert f(p18) == 8*K(-5, x)*K(1, y)*K(x, z) - 8*K(-5, x)*K(x, z) - 8*K(1, y) + 4
+
+    p19 = Piecewise((0, x > 2), (1, True))
+    assert f(p19) == p19
+
+    p20 = Piecewise((0, And(x < 2, x > -5)), (1, True))
+    assert f(p20) == p20
+
+    p21 = Piecewise((0, Or(x > 1, x < 0)), (1, True))
+    assert f(p21) == p21
+
+    p22 = Piecewise((0, ~((Eq(y, -1) | Ne(x, 0)) & (Ne(x, 1) | Ne(y, -1)))), (1, True))
+    assert f(p22) == K(-1, y)*K(0, x) - K(-1, y)*K(1, x) - K(0, x) + 1
diff --git a/sympy/functions/special/tests/test_tensor_functions.py b/sympy/functions/special/tests/test_tensor_functions.py
index 831a3d932b..34047c0440 100644
--- a/sympy/functions/special/tests/test_tensor_functions.py
+++ b/sympy/functions/special/tests/test_tensor_functions.py
@@ -1,6 +1,6 @@
 from sympy import (
     adjoint, conjugate, Dummy, Eijk, KroneckerDelta, LeviCivita, Symbol,
-    symbols, transpose,
+    symbols, transpose, Piecewise, Ne
 )
 from sympy.core.compatibility import range, long
 from sympy.physics.secondquant import evaluate_deltas, F
@@ -61,6 +61,8 @@ def test_kronecker_delta():
     # to test if canonical
     assert (KroneckerDelta(i, j) == KroneckerDelta(j, i)) == True
 
+    assert KroneckerDelta(i, j).rewrite(Piecewise) == Piecewise((0, Ne(i, j)), (1, True))
+
 def test_kronecker_delta_secondquant():
     """secondquant-specific methods"""
     D = KroneckerDelta
diff --git a/sympy/simplify/tests/test_simplify.py b/sympy/simplify/tests/test_simplify.py
index bd5a4e4391..d9e1a41b49 100644
--- a/sympy/simplify/tests/test_simplify.py
+++ b/sympy/simplify/tests/test_simplify.py
@@ -3,8 +3,8 @@
     collect,cos, cosh, cot, coth, count_ops, csch, Derivative, diff, E,
     Eq, erf, exp, exp_polar, expand, expand_multinomial, factor,
     factorial, Float, fraction, Function, gamma, GoldenRatio, hyper,
-    hypersimp, I, Integral, integrate, log, logcombine, Lt, Matrix,
-    MatrixSymbol, Mul, nsimplify, O, oo, pi, Piecewise, posify, rad,
+    hypersimp, I, Integral, integrate, KroneckerDelta, log, logcombine, Lt,
+    Matrix, MatrixSymbol, Mul, nsimplify, O, oo, pi, Piecewise, posify, rad,
     Rational, root, S, separatevars, signsimp, simplify, sign, sin,
     sinc, sinh, solve, sqrt, Sum, Symbol, symbols, sympify, tan, tanh,
     zoo)
@@ -837,3 +837,27 @@ def test_issue_17141():
     else:
         raises(RuntimeError, lambda: simplify(2**acos(I+1)**2))
         raises(RuntimeError, lambda: simplify((2**acos(I+1)**2).rewrite('log')))
+
+
+def test_simplify_kroneckerdelta():
+    i, j = symbols("i j")
+    K = KroneckerDelta
+
+    assert simplify(K(i, j)) == K(i, j)
+    assert simplify(K(0, j)) == K(0, j)
+    assert simplify(K(i, 0)) == K(i, 0)
+
+    assert simplify(K(0, j).rewrite(Piecewise) * K(1, j)) == 0
+    assert simplify(K(1, i) + Piecewise((1, Eq(j, 2)), (0, True))) == K(1, i) + K(2, j)
+
+    # issue 17214
+    assert simplify(K(0, j) * K(1, j)) == 0
+
+    n = Symbol('n', integer=True)
+    assert simplify(K(0, n) * K(1, n)) == 0
+
+    M = Matrix(4, 4, lambda i, j: K(j - i, n) if i <= j else 0)
+    assert simplify(M**2) == Matrix([[K(0, n), 0, K(1, n), 0],
+                                     [0, K(0, n), 0, K(1, n)],
+                                     [0, 0, K(0, n), 0],
+                                     [0, 0, 0, K(0, n)]])