Repository: sympy/sympy
Pull request: added _eval_derivative for products
This PR adds the `_eval_derivative` method based on the `Leibniz rule`.

Fixes #20848

Added `_eval_derivative` method based on the `Leibniz rule`.
Also added test cases for the same.

Credits to @Maelstrom6 and @asmeurer 

<!-- BEGIN RELEASE NOTES -->
NO ENTRY
<!-- END RELEASE NOTES -->

Problem statement:
diff(Product(...)) does not evaluate the derivative
Usually Derivative gives an unevaluated derivative and doit is used to evaluate it:
```python
In [4]: Derivative(x**2, x)
Out[4]: 
d ⎛ 2⎞
──⎝x ⎠
dx    

In [5]: Derivative(x**2, x).doit()
Out[5]: 2⋅x

In [6]: diff(x**2, x)
Out[6]: 2⋅x
```
The diff function is used to create the derivative and call doit but it doesn't work with a Product:
```python
In [8]: Derivative(Product(x, (y, 1, z)), x)
Out[8]: 
  ⎛  z    ⎞
∂ ⎜─┬─┬─  ⎟
──⎜ │ │  x⎟
∂x⎜ │ │   ⎟
  ⎝y = 1  ⎠

In [9]: Derivative(Product(x, (y, 1, z)), x).doit()
Out[9]: 
 z  
x ⋅z
────
 x  

In [10]: diff(Product(x, (y, 1, z)), x)
Out[10]: 
  ⎛  z    ⎞
∂ ⎜─┬─┬─  ⎟
──⎜ │ │  x⎟
∂x⎜ │ │   ⎟
  ⎝y = 1  ⎠
```
Here `Out[10]` should be the same as `Out[9]`.

Hints from issue discussion:
Probably because it doesn't actually know how to differentiate a Product. [9] only works because the original Product evaluates first and then the derivative does. We could implement the general Leibniz rule (https://en.wikipedia.org/wiki/Product_rule#A_product_of_more_than_two_factors). I'm not sure if there are concerns with that formula if the product is infinite. 

Outside of this, I'm not sure if Derivative.doit() should recursively call doit() before evaluating. 
> We could implement the general Leibniz rule (https://en.wikipedia.org/wiki/Product_rule#A_product_of_more_than_two_factors). I'm not sure if there are concerns with that formula if the product is infinite.

I think it does work for infinite products if they converge. [Wiki Proof](https://proofwiki.org/wiki/Derivative_of_Infinite_Product_of_Analytic_Functions).
@Maelstrom6, is this issue open, can I take it?
Thanks.
Feel free to work on it👍. #20851 is working on this but the author is fairly busy at the moment. Maybe have a look at the discussion there for some ideas on how to tackle the problem. 
@Maelstrom6 I tried something like with the code you suggested with respect to the discussions:
```
def _eval_derivative(self, x):
        from sympy import Sum
        limits = list(self.limits)
        if any(x in limit[1:].free_symbols for limit in reversed(limits)):
            return None
        return self * Sum(self.function.diff(x) / self.function, self._args[1:], **self.assumptions0)
```
It covers if `x` is in the free symbols of any of the lower or upper limits, it returns unevaluated.
Looks good. You can open up a PR for it and we can do a few finishing touches there. 
Sure, thanks a lot!

Commit message: added test cases

Patch to review:
diff --git a/sympy/concrete/products.py b/sympy/concrete/products.py
index 6ff6a06f66b9..3ea873c8d6e9 100644
--- a/sympy/concrete/products.py
+++ b/sympy/concrete/products.py
@@ -399,6 +399,18 @@ def _eval_product_direct(self, term, limits):
         (k, a, n) = limits
         return Mul(*[term.subs(k, a + i) for i in range(n - a + 1)])
 
+    def _eval_derivative(self, x):
+        if self.is_zero:
+            return S.Zero
+        from sympy.concrete.summations import Sum, Dummy
+        i = Dummy('i')
+        f = self.function
+        limits = list(self.limits)
+        _, a, b = limits[0]
+        if any(x in limit[1:].free_symbols for limit in limits):
+            return None
+        return Sum(f.subs(_, i).diff(x) * Product(f, (_, a, i - 1)) * Product(f, (_, i + 1, b)), (i, a, b), **self.assumptions0)
+
     def is_convergent(self):
         r"""
         See docs of :obj:`.Sum.is_convergent()` for explanation of convergence
diff --git a/sympy/concrete/tests/test_products.py b/sympy/concrete/tests/test_products.py
index 0d2de4a45109..b6d759210509 100644
--- a/sympy/concrete/tests/test_products.py
+++ b/sympy/concrete/tests/test_products.py
@@ -1,6 +1,6 @@
 from sympy import (symbols, Symbol, product, combsimp, factorial, rf, sqrt, cos,
                    Function, Product, Rational, Sum, oo, exp, log, S, pi,
-                   KroneckerDelta)
+                   KroneckerDelta, Derivative, diff, sin, polygamma, gamma, Dummy)
 from sympy.testing.pytest import raises
 from sympy import simplify
 
@@ -382,3 +382,13 @@ def test_rewrite_Sum():
 def test_KroneckerDelta_Product():
     y = Symbol('y')
     assert Product(x*KroneckerDelta(x, y), (x, 0, 1)).doit() == 0
+
+def test_issue_20848():
+    _i = Dummy('_i')
+    t, y, z = symbols('t y z')
+    assert diff(Product(x, (y, 1, z)), x).as_dummy() == Sum(Product(x, (y, 1, _i - 1))*Product(x, (y, _i + 1, z)), (_i, 1, z)).as_dummy()
+    assert diff(Product(x, (y, 1, z)), x).doit() == x**z*z/x
+    assert diff(Product(x, (y, x, z)), x) == Derivative(Product(x, (y, x, z)), x)
+    assert diff(Product(t, (x, 1, z)), x) == S(0)
+    assert Product(sin(n*x), (n, -1, 1)).diff(x).doit() == S(0)
+    assert Product(x*y, (x, 1, z), (y, 1, m)).diff(z).doit() == m*factorial(m)**z*factorial(z)**m*gamma(z + 1)*polygamma(0, z + 1)/factorial(z) + log(factorial(m))*factorial(m)**z*factorial(z)**m


Code review comments:
diff_hunk: @@ -382,3 +382,13 @@ def test_rewrite_Sum():
 def test_KroneckerDelta_Product():
     y = Symbol('y')
     assert Product(x*KroneckerDelta(x, y), (x, 0, 1)).doit() == 0
+
+def test_issue_20848():
+    _i = Dummy('_i')
path: sympy/concrete/tests/test_products.py
text: @user1:
The underscore in `'_i'` is not needed. `Dummy` automatically appends the underscore for the string representation.

@author:
Sure will change.
Thanks.

diff_hunk: @@ -399,6 +399,18 @@ def _eval_product_direct(self, term, limits):
         (k, a, n) = limits
         return Mul(*[term.subs(k, a + i) for i in range(n - a + 1)])
 
+    def _eval_derivative(self, x):
+        if self.is_zero:
+            return S.Zero
+        from sympy.concrete.summations import Sum, Dummy
+        i = Dummy('i')
+        f = self.function
+        limits = list(self.limits)
+        _, a, b = limits[0]
+        if any(x in limit[1:].free_symbols for limit in limits):
+            return None
+        return Sum(f.subs(_, i).diff(x) * Product(f, (_, a, i - 1)) * Product(f, (_, i + 1, b)), (i, a, b), **self.assumptions0)
path: sympy/concrete/products.py
text: @user1:
This would produce wrong answers if `limits` has more than one element. I thought that `Sum` had a similar problem but it does not. Looking at the code of `Sum`: 
- it pops `limits[0]` from the stack
- checks if there are still limits remaining
    * if there are limits remaining, then it replaces `f` with `Sum(f, remaining_limits)`
- it calls `Derivative(f, ...)` instead of `f.diff(x)`

This way, it would keep all the sums as required. 

I think a similar mechanism should be a adopted here. It might even be better to copy the code and change it as you need to. That way, if there is a bug in `Sum.diff` then there will be a similar bug in `Product.diff` and they can be fixed in a similar way. I think it might help with maintainability.

@author:
> This would produce wrong answers if `limits` has more than one element.
>
Can you please provide an example?
I ran it for:
```
In [2]: Product(x*y, (x, 1, z), (y, 1, m)).diff(z).doit()
Out[2]:
    z   m
m⋅m! ⋅z! ⋅Γ(z + 1)⋅polygamma(0, z + 1)             z   m
────────────────────────────────────── + log(m!)⋅m! ⋅z!
                  z!

In [3]: Product(x*y, (x, 1, z), (y, 1, m)).doit().diff(z)
Out[3]:
    z   m
m⋅m! ⋅z! ⋅Γ(z + 1)⋅polygamma(0, z + 1)             z   m
────────────────────────────────────── + log(m!)⋅m! ⋅z!
                  z!
```
This seems to work though.
Thanks.


@user1:
With this, you are differentiating with respect to one of the indexes of the summation. So the `_eval_derivaitve` returns `None` and the value becomes `Derivative(Product(x*y, (x, 1, z), (y, 1, m)), z)` and calling doit will simplify the product and then it takes the derivative. 

I'll quickly come up with an example that shows this issue.

@user1:
```python
from sympy import *
x, y = symbols("x y")
i, j, m, n = symbols("i, j, m, n", integer=True)

expr = Product(sin(x*i/j), (i, 1, m), (j, 1, n))
print(expr.diff(x))
# output:
# Sum(_i*cos(_i*x/j)*Product(sin(i*x/j), (i, 1, _i - 1))*Product(sin(i*x/j), (i, _i + 1, m))/j, (_i, 1, m))
```
As you can see, the `(j, 1, n)` part has completely disappeared from the output.

@author:
@user1 Thank you for correcting, I went ahead as you suggested to start with code in `Sum`.
The implementation below seems to works for the case you gave and for the previous cases:

```
def _eval_derivative(self, x):
        from sympy.concrete.summations import Sum, Dummy, Symbol, Derivative
        if isinstance(x, Symbol) and x not in self.free_symbols:
            return S.Zero
        f, limits = self.function, list(self.limits)
        limit = limits.pop(-1)
        if limits:
            f = self.func(f, *limits)
        _, a, b = limit
        if x in a.free_symbols or x in b.free_symbols:
            return None
        h = Dummy()
        rv = Sum( Product(f, (_, a, h - 1)) * Product(f, (_, h + 1, b)) * Derivative(f, x, evaluate=True).subs(_, h), (h, a, b))
        return rv
```

```
In [1]: i,j=symbols('i j')

In [2]:  expr = Product(sin(x*i/j), (i, 1, n), (j, 1, m))

In [3]: print(expr.diff(x))
Sum(Product(sin(i*x/j), (i, 1, n), (j, 1, _Dummy_23 - 1))*Product(sin(i*x/j), (i, 1, n), (j, _Dummy_23 + 1, m))*Sum(_Dummy_24*cos(_Dummy_24*x/_Dummy_23)*Product(sin(i*x/_Dummy_23), (i, 1, _Dummy_24 - 1))*Product(sin(i*x/_Dummy_23), (i, _Dummy_24 + 1, n))/_Dummy_23, (_Dummy_24, 1, n)), (_Dummy_23, 1, m))
```
Kindly give your opinion.
Thanks.