Repository: sympy/sympy

Issue:
Sum((diag(0, 2, -3))**n, (n, 0, 3)).doit() causes an error
```
>>> k  = symbols('k', positive = True, integer = True)
>>> Sum(x**n/(n + 1), (n, 0, k)).doit()
⎧   k + 1                         log(-x + 1)                     
⎪- x     ⋅lerchphi(x, 1, k + 2) - ───────────  for │x│ ≤ 1 ∧ x ≠ 1
⎪                                      x                          
⎪                                                                 
⎪                  k                                              
⎪                 ____                                            
⎪                 ╲                                               
⎨                  ╲      n                                       
⎪                   ╲    x                                        
⎪                   ╱  ─────                        otherwise     
⎪                  ╱   n + 1                                      
⎪                 ╱                                               
⎪                 ‾‾‾‾                                            
⎪                n = 0                                            
⎩                                                                 
```
An evaluation of a finite sum containing a convergence condition

```
>>> from sympy.concrete.gosper import gosper_sum
>>> gosper_sum(x**n, (n, 0, k))
   k    
x⋅x  - 1
────────
 x - 1  
>>> gosper_sum(n*x**n, (n, 0, k))
  ⎛     k      k    k    ⎞
x⋅⎝k⋅x⋅x  - k⋅x  - x  + 1⎠
──────────────────────────
                2         
         (x - 1)          
```
The above have `x*x**k` rather than `x**(k+1)`

```
>>> A = diag(0, 2, -3)
>>> A**0 + A**1 + A**2 + A**3
⎡1  0    0 ⎤
⎢          ⎥
⎢0  15   0 ⎥
⎢          ⎥
⎣0  0   -20⎦
>>> Sum(A**n, (n, 0, 3)).doit() 
Traceback (most recent call last):

  File "<ipython-input-83-bb0cff02949f>", line 1, in <module>
    Sum(A**n, (n, 0, 3))

  File ".\sympy\core\decorators.py", line 132, in binary_op_wrapper
    return func(self, other)

  File ".\sympy\matrices\common.py", line 2052, in __pow__
    return a._matrix_pow_by_jordan_blocks(num)

  File ".\sympy\matrices\matrices.py", line 1888, in _matrix_pow_by_jordan_blocks
    jordan_cell_power(j, num)

  File ".\sympy\matrices\matrices.py", line 1873, in jordan_cell_power
    raise ValueError("Matrix det == 0; not invertible")

ValueError: Matrix det == 0; not invertible
```

```
>>> Sum((n**(-2/n))**n, (n, 1, oo)).doit()
  ∞          
_____        
╲            
 ╲          n
  ╲   ⎛ -2 ⎞ 
   ╲  ⎜ ───⎟ 
   ╱  ⎜  n ⎟ 
  ╱   ⎝n   ⎠ 
 ╱           
╱            
‾‾‾‾‾        
n = 1        
>>> Sum(n**(-2), (n, 1, oo)).doit()
 2
π 
──
6 
```





Discussion:
oscargus: Seems like the problem causing the failed test was fixed elsewhere.
github-actions[bot]: Benchmark results from GitHub Actions

Lower numbers are good, higher numbers are bad. A ratio less than 1
means a speed up and greater than 1 means a slowdown. Green lines
beginning with `+` are slowdowns (the PR is slower then master or
master is slower than the previous release). Red lines beginning
with `-` are speedups.

Significantly changed benchmark results (PR vs master)
```diff

```
Significantly changed benchmark results (master vs previous release)
```diff
       before           after         ratio
     [907895ac]       [b0b66d55]
-         250±4ms          156±5ms     0.62  large_exprs.TimeLargeExpressionOperations.time_subs
-        271±10μs          122±1μs     0.45  matrices.TimeMatrixExpression.time_MatMul
-     16.3±0.08ms      8.90±0.07ms     0.54  matrices.TimeMatrixExpression.time_MatMul_doit
-      4.91±0.01s         370±20ms     0.08  polygon.PolygonArbitraryPoint.time_bench01
+     3.92±0.01ms       6.77±0.1ms     1.73  solve.TimeMatrixOperations.time_det(4, 2)
+     3.94±0.02ms       6.72±0.1ms     1.71  solve.TimeMatrixOperations.time_det_bareiss(4, 2)
+      44.1±0.2ms         83.4±2ms     1.89  solve.TimeMatrixSolvePyDySlow.time_linsolve(1)
+      47.5±0.4ms         82.1±1ms     1.73  solve.TimeMatrixSolvePyDySlow.time_solve(1)

```
Full benchmark results can be found as artifacts in GitHub Actions
(click on checks at the top of the PR).

Patch:
diff --git a/sympy/concrete/products.py b/sympy/concrete/products.py
index 68b0123fb5..b345959e26 100644
--- a/sympy/concrete/products.py
+++ b/sympy/concrete/products.py
@@ -392,7 +392,7 @@ def _eval_product(self, term, limits):
 
     def _eval_simplify(self, **kwargs):
         from sympy.simplify.simplify import product_simplify
-        rv = product_simplify(self)
+        rv = product_simplify(self, **kwargs)
         return rv.doit() if kwargs['doit'] else rv
 
     def _eval_transpose(self):
diff --git a/sympy/concrete/summations.py b/sympy/concrete/summations.py
index 58b73d75f6..bb9faa7dd7 100644
--- a/sympy/concrete/summations.py
+++ b/sympy/concrete/summations.py
@@ -332,8 +332,13 @@ def _eval_difference_delta(self, n, step):
 
     def _eval_simplify(self, **kwargs):
 
+        function = self.function
+
+        if kwargs.get('deep', True):
+            function = function.simplify(**kwargs)
+
         # split the function into adds
-        terms = Add.make_args(expand(self.function))
+        terms = Add.make_args(expand(function))
         s_t = [] # Sum Terms
         o_t = [] # Other Terms
 
@@ -348,7 +353,7 @@ def _eval_simplify(self, **kwargs):
                     # go through each term
                     if isinstance(subterm, Sum):
                         # if it's a sum, simplify it
-                        out_terms.append(subterm._eval_simplify())
+                        out_terms.append(subterm._eval_simplify(**kwargs))
                     else:
                         # otherwise, add it as is
                         out_terms.append(subterm)
diff --git a/sympy/simplify/simplify.py b/sympy/simplify/simplify.py
index 825bb0f784..4d78c277bc 100644
--- a/sympy/simplify/simplify.py
+++ b/sympy/simplify/simplify.py
@@ -714,7 +714,7 @@ def done(e):
             i: factor_terms(i) for i in expr.atoms(Integral)})
 
     if expr.has(Product):
-        expr = product_simplify(expr)
+        expr = product_simplify(expr, **kwargs)
 
     from sympy.physics.units import Quantity
 
@@ -883,15 +883,20 @@ def __refactor(val):
     return Add(self, other)
 
 
-def product_simplify(s):
+def product_simplify(s, **kwargs):
     """Main function for Product simplification"""
     terms = Mul.make_args(s)
     p_t = [] # Product Terms
     o_t = [] # Other Terms
 
+    deep = kwargs.get('deep', True)
     for term in terms:
         if isinstance(term, Product):
-            p_t.append(term)
+            if deep:
+                p_t.append(Product(term.function.simplify(**kwargs),
+                                   *term.limits))
+            else:
+                p_t.append(term)
         else:
             o_t.append(term)
 
@@ -902,8 +907,9 @@ def product_simplify(s):
             if not used[i]:
                 for j, p_term2 in enumerate(p_t):
                     if not used[j] and i != j:
-                        if isinstance(product_mul(p_term1, p_term2, method), Product):
-                            p_t[i] = product_mul(p_term1, p_term2, method)
+                        tmp_prod = product_mul(p_term1, p_term2, method)
+                        if isinstance(tmp_prod, Product):
+                            p_t[i] = tmp_prod
                             used[j] = True
 
     result = Mul(*o_t)


Test patch:
diff --git a/sympy/concrete/tests/test_products.py b/sympy/concrete/tests/test_products.py
index 63c6c41fa4..5195f544aa 100644
--- a/sympy/concrete/tests/test_products.py
+++ b/sympy/concrete/tests/test_products.py
@@ -275,8 +275,9 @@ def test_conjugate_transpose():
     assert p.conjugate().doit() == p.doit().conjugate()
     assert p.transpose().doit() == p.doit().transpose()
 
+
 def test_simplify_prod():
-    y, t, b, c = symbols('y, t, b, c', integer = True)
+    y, t, b, c, v, d = symbols('y, t, b, c, v, d', integer = True)
 
     _simplify = lambda e: simplify(e, doit=False)
     assert _simplify(Product(x*y, (x, n, m), (y, a, k)) * \
@@ -294,6 +295,12 @@ def test_simplify_prod():
     assert _simplify(Product(x, (t, a, b)) * Product(x, (t, b+1, c)) * \
         Product(y, (t, a, b))) == Product(x*y, (t, a, b)) * \
             Product(x, (t, b+1, c))
+    assert _simplify(Product(sin(t)**2 + cos(t)**2 + 1, (t, a, b))) == \
+        Product(2, (t, a, b))
+    assert _simplify(Product(sin(t)**2 + cos(t)**2 - 1, (t, a, b))) == \
+           Product(0, (t, a, b))
+    assert _simplify(Product(v*Product(sin(t)**2 + cos(t)**2, (t, a, b)),
+                             (v, c, d))) == Product(v*Product(1, (t, a, b)), (v, c, d))
 
 
 def test_change_index():
diff --git a/sympy/concrete/tests/test_sums_products.py b/sympy/concrete/tests/test_sums_products.py
index 21172bfd4b..a6d00f60c0 100644
--- a/sympy/concrete/tests/test_sums_products.py
+++ b/sympy/concrete/tests/test_sums_products.py
@@ -1,9 +1,13 @@
 from math import prod
 
+from sympy.concrete.expr_with_intlimits import ReorderError
 from sympy.concrete.products import (Product, product)
-from sympy.concrete.summations import (Sum, summation)
+from sympy.concrete.summations import (Sum, summation, telescopic,
+     eval_sum_residue, _dummy_with_inherited_properties_concrete)
 from sympy.core.function import (Derivative, Function)
 from sympy.core import (Catalan, EulerGamma)
+from sympy.core.facts import InconsistentAssumptions
+from sympy.core.mod import Mod
 from sympy.core.numbers import (E, I, Rational, nan, oo, pi)
 from sympy.core.relational import Eq
 from sympy.core.singleton import S
@@ -25,19 +29,14 @@
 from sympy.logic.boolalg import And, Or
 from sympy.matrices.expressions.matexpr import MatrixSymbol
 from sympy.matrices.expressions.special import Identity
+from sympy.matrices import (Matrix, SparseMatrix,
+    ImmutableDenseMatrix, ImmutableSparseMatrix, diag)
 from sympy.sets.fancysets import Range
 from sympy.sets.sets import Interval
 from sympy.simplify.combsimp import combsimp
 from sympy.simplify.simplify import simplify
 from sympy.tensor.indexed import (Idx, Indexed, IndexedBase)
-from sympy.concrete.summations import (
-    telescopic, _dummy_with_inherited_properties_concrete, eval_sum_residue)
-from sympy.concrete.expr_with_intlimits import ReorderError
-from sympy.core.facts import InconsistentAssumptions
 from sympy.testing.pytest import XFAIL, raises, slow
-from sympy.matrices import (Matrix, SparseMatrix,
-    ImmutableDenseMatrix, ImmutableSparseMatrix)
-from sympy.core.mod import Mod
 from sympy.abc import a, b, c, d, k, m, x, y, z
 
 n = Symbol('n', integer=True)
@@ -863,6 +862,8 @@ def test_simplify_sum():
     assert _simplify(Sum(Sum(d * t, (x, a, b - 1)) + \
                 Sum(d * t, (x, b, c)), (t, a, b))) == \
                     d * Sum(1, (x, a, c)) * Sum(t, (t, a, b))
+    assert _simplify(Sum(sin(t)**2 + cos(t)**2 + 1, (t, a, b))) == \
+        2 * Sum(1, (t, a, b))
 
 
 def test_change_index():
@@ -1164,6 +1165,12 @@ def test_issue_14112():
     assert Sum((-2)**n + (-3)**n, (n, 1, oo)).is_convergent() is S.false
 
 
+def test_issue_14219():
+    A = diag(0, 2, -3)
+    res = diag(1, 15, -20)
+    assert Sum(A**n, (n, 0, 3)).doit() == res
+
+
 def test_sin_times_absolutely_convergent():
     assert Sum(sin(n) / n**3, (n, 1, oo)).is_convergent() is S.true
     assert Sum(sin(n) * log(n) / n**3, (n, 1, oo)).is_convergent() is S.true