Trac #22937: Implement a "distribute" method

Motivation : some "obvious" transformations, sometimes helpful in
routine calculations, cannot be obtained via the current methods. Some
(elementary,
i. e. high-school- or undergrad-level) exemples : given the declarations
:
{{{
var("a,b")
var("j,p", domain="integer")
f,g,X=function("f,g,X")
}}}

the following equalities, true under no or mild conditions, cannot be
deduced currently :
{{{
sum(X(j)+Y(j),j,1,p)==sum(X(j),j,1,p)+sum(Y(j),j,1,p) ## (1)
prod(X(j)*Y(j),j,1,p)==prod(X(j),j,1,p)*prod(Y(j),j,1,p) ## (2)
integrate(f(x)+g(x),x,a,b)==integrate(f(x),x,a,b)+integrate(g(x),x,a,b)
## (3)
}}}

Similarly, if A, B and C are matrices over a suitable ring and of
suitable dimensions, the following equalities cannot be
straightforwardly obtained.
{{{
(A+B).trace()==A.trace()+B.trace() ## (4)
(A*B).det()==A.det()*B.det() ## (5)
}}}

(Curiously, {{{(f(x)+g(x)).diff(x)}}} ''does'' return
{{{diff(f(x),x)+diff(g(x),x)}}} (and a way to return to the original
expression does not seem to exist currently)).

The ability to generate the right-hand form from the left-hand form of
these various examples is sometimes useful in (low-level)
examples. The examples (1) and/or (2) arise naturally when deriving the
maximum-likelihood estimators of the mean and variance of a given
distribution. The third one is often encountered in elementary calculus
exercises.

In all cases, an operator is "distributed" over another one :
(1) : symbolic_sum is distributed over addition.
(2) : symbolic product is distributed over multiplication.
(3) : integration is ditributed over addition.
(4) : trace is distributed over (matrix) addition.
(5) : determinant is distributed over (matrix) multiplication.

Implementing (1) as an extension of the {{{expand()}}} method of
Expression is tempting (see #22915). However, there are situations
where this expansion is '''not''' useful. So, creating a method for such
situations seems useful.

One should note that such "distributions" are not systematically
valid : we have a collection of special cases rather than a general
mechanism. So this method shoud either limit itself to a few
recognized cases or take keyword argument(s) specifying what should be
distributed over what.

Another design choice is to know if these distributions should be
recursive or run only on the top level of a given expression
(possibly controlled by another keyword argument).

The present ticket aims at implementing this method for
{{{sage.symbolic.expression.Expression}}}, at least in cases (1) to (3)
(case (2) needs an implementation of the symbolic products, see #17505).
Further suggestions for other classes are welcome...

URL: https://trac.sagemath.org/22937
Reported by: charpent
Ticket author(s): Emmanuel Charpentier, Ralf Stephan
Reviewer(s): Travis Scrimshaw

Author: Release Manager

src/sage/functions/other.py

@@ -2615,9 +2615,9 @@ def _print_latex_(self, x, var, a, b):
 
             sage: from sage.functions.other import symbolic_sum as ssum
             sage: latex(ssum(x^2, x, 1, 10))
-            \sum_{x=1}^{10} x^2
+            {\sum_{x=1}^{10} x^2}
         """
-        return r"\sum_{{{}={}}}^{{{}}} {}".format(var, a, b, x)
+        return r"{{\sum_{{{}={}}}^{{{}}} {}}}".format(var, a, b, x)
 
 symbolic_sum = Function_sum()
 
@@ -2662,8 +2662,8 @@ def _print_latex_(self, x, var, a, b):
 
             sage: from sage.functions.other import symbolic_product as sprod
             sage: latex(sprod(x^2, x, 1, 10))
-            \prod_{x=1}^{10} x^2
+            {\prod_{x=1}^{10} x^2}
         """
-        return r"\prod_{{{}={}}}^{{{}}} {}".format(var, a, b, x)
+        return r"{{\prod_{{{}={}}}^{{{}}} {}}}".format(var, a, b, x)
 
 symbolic_product = Function_prod()

src/sage/interfaces/maxima_lib.py

@@ -902,17 +902,14 @@ def sr_sum(self,*args):
 
     def sr_prod(self,*args):
         """
-        Helper function to wrap calculus use of Maxima's summation.
+        Helper function to wrap calculus use of Maxima's product.
         """
         try:
             return max_to_sr(maxima_eval([[max_ratsimp],[[max_simplify_prod],([max_prod],[sr_to_max(SR(a)) for a in args])]]));
         except RuntimeError as error:
             s = str(error)
             if "divergent" in s:
-# in pexpect interface, one looks for this;
-# could not find an example where 'Pole encountered' occurred, though
-#            if "divergent" in s or 'Pole encountered' in s:
-                raise ValueError("Sum is divergent.")
+                raise ValueError("Product is divergent.")
             elif "Is" in s: # Maxima asked for a condition
                 self._missing_assumption(s)
             else:

src/sage/symbolic/expression.pyx

@@ -10382,6 +10382,97 @@ cdef class Expression(CommutativeRingElement):
 
     log_expand = expand_log
 
+    def distribute(self, recursive=True):
+        """
+        Distribute some indexed operators over similar operators in
+        order to allow further groupings or simplifications. 
+
+        Implemented cases (so far) :
+
+        - Symbolic sum of a sum ==> sum of symbolic sums
+
+        - Integral (definite or not) of a sum ==> sum of integrals.
+
+        - Symbolic product of a product ==> product of symbolic products.
+
+        INPUT:
+
+        - ``recursive`` -- (default : True) the distribution proceeds
+          along the subtrees of the expression.
+
+        TESTS:
+
+            sage: var("j,k,p,q", domain="integer")
+            (j, k, p, q)
+            sage: X,Y,Z,f,g=function("X,Y,Z,f,g")
+            sage: var("x,a,b")
+            (x, a, b)
+            sage: sum(X(j)+Y(j),j,1,p)
+            sum(X(j) + Y(j), j, 1, p)
+            sage: sum(X(j)+Y(j),j,1,p).distribute()
+            sum(X(j), j, 1, p) + sum(Y(j), j, 1, p)
+            sage: integrate(f(x)+g(x),x)
+            integrate(f(x) + g(x), x)
+            sage: integrate(f(x)+g(x),x).distribute()
+            integrate(f(x), x) + integrate(g(x), x)
+            sage: integrate(f(x)+g(x),x,a,b)
+            integrate(f(x) + g(x), x, a, b)
+            sage: integrate(f(x)+g(x),x,a,b).distribute()
+            integrate(f(x), x, a, b) + integrate(g(x), x, a, b)
+            sage: sum(X(j)+sum(Y(k)+Z(k),k,1,q),j,1,p)
+            sum(X(j) + sum(Y(k) + Z(k), k, 1, q), j, 1, p)
+            sage: sum(X(j)+sum(Y(k)+Z(k),k,1,q),j,1,p).distribute()
+            sum(sum(Y(k), k, 1, q) + sum(Z(k), k, 1, q), j, 1, p) + sum(X(j), j, 1, p)
+            sage: sum(X(j)+sum(Y(k)+Z(k),k,1,q),j,1,p).distribute(recursive=False)
+            sum(X(j), j, 1, p) + sum(sum(Y(k) + Z(k), k, 1, q), j, 1, p)
+            sage: maxima("product(X(j)*Y(j),j,1,p)").sage()
+            product(X(j)*Y(j), j, 1, p)
+            sage: maxima("product(X(j)*Y(j),j,1,p)").sage().distribute()
+            product(X(j), j, 1, p)*product(Y(j), j, 1, p)
+
+
+        AUTHORS:
+
+        - Emmanuel Charpentier, Ralf Stephan (05-2017)
+        """
+        from sage.functions.other import symbolic_sum as opsum, \
+            symbolic_product as opprod
+        from sage.symbolic.integration.integral \
+            import indefinite_integral as opii, definite_integral as opdi
+        from sage.symbolic.operators import add_vararg as opadd, \
+            mul_vararg as opmul
+        from sage.all import prod
+        def treat_term(op, term, args):
+            l=sage.all.copy(args)
+            l.insert(0, term)
+            return(apply(op, l))
+        if self.parent() is not sage.all.SR: return self
+        op = self.operator()
+        if op is None : return self
+        if op in {opsum, opdi, opii}:
+            sa = self.operands()[0].expand()
+            op1 = sa.operator()
+            if op1 is opadd:
+                la = self.operands()[1:]
+                aa = sa.operands()
+                if recursive:
+                    return sum(treat_term(op, t.distribute(), la) for t in aa)
+                return sum(treat_term(op, t, la) for t in aa)
+            return self
+        if op is opprod:
+            sa = self.operands()[0].expand()
+            op1 = sa.operator()
+            if op1 is opmul:
+                la = self.operands()[1:]
+                aa = sa.operands()
+                if recursive:
+                    return prod(treat_term(op, t.distribute(), la) for t in aa)
+                return prod(treat_term(op, t, la) for t in aa)
+            return self
+        if recursive:
+            return apply(op, map(lambda t:t.distribute(), self.operands()))
+        return self
+
 
     def factor(self, dontfactor=[]):
         """
@@ -12421,4 +12512,3 @@ cdef operators compatible_relation(operators lop, operators rop) except <operato
        return greater
     else:
         raise TypeError("incompatible relations")
-