Repository: sympy/sympy
Pull request: Avoid infinite integration by parts in manualintegrate
#### References to other Issues or PRs

Fixes #9858 

#### Brief description of what is fixed or changed

With an integral like `exp(x) * cos(exp(x))`, manualintegrate attempts integration by parts with `u = cos(exp(x))`, which "succeeds", but leads to an infinite chain of different integrals such as `exp(n*x) * cos(exp(x))` and `exp(n*x) * sin(exp(x))`. While repeated integration with the same `dv` term can be useful (like `x**9 * exp(x)`), the same is not true for repeated use of `u` term, i.e., differentiating a factor and getting the same factor back. This is now prevented by keeping track of u-terms used by integration by parts.
 
#### Release Notes
 
<!-- BEGIN RELEASE NOTES -->
* integrals
  * prevented creation of infinite chains of integration by parts 
<!-- END RELEASE NOTES -->


Problem statement:
manualintegrate(exp(x)*cos(exp(x)), x) infinite loop
`manualintegrate(exp(x)*cos(exp(x)), x)` goes into an infinite loop when it tries to do integration by parts. 

It picks `cos(exp(x))` as u and `exp(x)` as dv by the LIATE rule, and then tries to evaluate the integral v*du. This causes it to try integration by parts again and causes the infinite recursion.

Hints from issue discussion:

Commit message: Avoid infinite integration by parts in manualintegrate

With an integral like `exp(x) * cos(exp(x))`, manualintegrate attempts
integration by parts with `u = cos(exp(x))`, which "succeeds", but
leads to an infinite chain of different integrals of the form `exp(n*x)
* cos(exp(x))` and `exp(n*x) * sin(exp(x))`. While repeated integration
with the same `dv` term can be useful (like `x**9 * exp(x)`), the same
is not true for repeated use of `u` term, i.e., differentiating a
factor and getting the same factor back. This is now prevented by
keeping track of u-terms used by integration by parts.

Patch to review:
diff --git a/sympy/integrals/manualintegrate.py b/sympy/integrals/manualintegrate.py
index 4938ebc88fe1..65e995a325c6 100644
--- a/sympy/integrals/manualintegrate.py
+++ b/sympy/integrals/manualintegrate.py
@@ -18,7 +18,7 @@
 """
 from __future__ import print_function, division
 
-from collections import namedtuple
+from collections import namedtuple, defaultdict
 
 import sympy
 
@@ -496,6 +496,12 @@ def parts_rule(integral):
         if isinstance(v, sympy.Integral):
             return
 
+        # Set a limit on the number of times u can be used
+        cachekey = u.xreplace({symbol: _cache_dummy})
+        if _parts_u_cache[cachekey] > 5:
+            return
+        _parts_u_cache[cachekey] += 1
+
         while True:
             if (integrand / (v * du)).cancel() == 1:
                 break
@@ -541,6 +547,7 @@ def make_second_step(steps, integrand):
         if (constant != 1) and rule:
             rule = ConstantTimesRule(constant, integrand, rule,
                                      constant * integrand, symbol)
+        cachekey = u.xreplace({symbol: _cache_dummy})
         return rule
 
 
@@ -1008,7 +1015,9 @@ def fallback_rule(integral):
 
 # Cache is used to break cyclic integrals.
 # Need to use the same dummy variable in cached expressions for them to match.
+# Also record "u" of integration by parts, to avoid infinite repetition.
 _integral_cache = {}
+_parts_u_cache = defaultdict(int)
 _cache_dummy = sympy.Dummy("z")
 
 def integral_steps(integrand, symbol, **options):
@@ -1390,6 +1399,8 @@ def manualintegrate(f, var):
     sympy.integrals.integrals.Integral
     """
     result = _manualintegrate(integral_steps(f, var))
+    # Clear the cache of u-parts
+    _parts_u_cache.clear()
     # If we got Piecewise with two parts, put generic first
     if isinstance(result, Piecewise) and len(result.args) == 2:
         cond = result.args[0][1]
diff --git a/sympy/integrals/tests/test_manual.py b/sympy/integrals/tests/test_manual.py
index a6d0622d9a06..b431c5ad5a9d 100644
--- a/sympy/integrals/tests/test_manual.py
+++ b/sympy/integrals/tests/test_manual.py
@@ -387,3 +387,12 @@ def test_issue_13297():
 def test_issue_14470():
     assert manualintegrate(1/(x*sqrt(x + 1)), x) == \
         log(-1 + 1/sqrt(x + 1)) - log(1 + 1/sqrt(x + 1))
+
+
+def test_issue_9858():
+    assert manualintegrate(exp(x)*cos(exp(x)), x) == sin(exp(x))
+    assert manualintegrate(exp(2*x)*cos(exp(x)), x) == \
+        exp(x)*sin(exp(x)) + cos(exp(x))
+    res = manualintegrate(exp(10*x)*sin(exp(x)), x)
+    assert not res.has(Integral)
+    assert res.diff(x) == exp(10*x)*sin(exp(x))


Code review comments:
diff_hunk: @@ -496,6 +496,12 @@ def parts_rule(integral):
         if isinstance(v, sympy.Integral):
             return
 
+        # Set a limit on the number of times u can be used
+        cachekey = u.xreplace({symbol: _cache_dummy})
+        if _parts_u_cache[cachekey] > 5:
path: sympy/integrals/manualintegrate.py
text: @user1:
I tried this with bound 2. It seems to work as well.

@author:
In the sense that the existing tests pass... but manualintegrate tests are not very extensive, and there are more complex integrals out there. An integral may have many parts that require picking the same `u` term. In fact, I just found a reasonable example that fails in this branch due to the limit imposed on integration by parts with `u=x`. 
```
>>> manualintegrate(sum([x*exp(k*x) for k in range(1, 8)]), x)
x*exp(7*x)/7 + x*exp(6*x)/6 + x*exp(5*x)/5 + x*exp(4*x)/4 + x*exp(3*x)/3 + x*exp(2*x)/2 - exp(7*x)/49 - exp(6*x)/36 - exp(5*x)/25 - exp(4*x)/16 - exp(3*x)/9 - exp(2*x)/4 + Integral(x*exp(x), x)
```
I guess this should be fixed. 

One way is to limit caching to `u` being an instance of `sin, cos, exp` -- those have the potential for infinite loops when differentiated, not other things like the powers of x.

Another is to try to clear the cache at a more appropriate place (not yet sure where), so that when multiple terms of a sum are being integrated, the choices made in one do not affect others. 

@user2:
This only catches infinite recursion from the same u term. I wonder if it can happen even with this logic, from slightly different terms each time. My guess is not currently because it only handles LIATE terms but would that continue to be the case if we extend it (e.g., https://github.com/sympy/sympy/issues/15431). 

@author:
Differentiating functions other than trigonometric / exponential is unlikely to return the integral to its original form. This check is meant for these types only. 

With the more targeted check, I reduced the bound to 2 as suggested by @user1, and added a test to make sure that `manualintegrate` is still able to do a bunch of similar integrations by parts.  

I don't know how to prevent manualintegrate from creating an infinite sequence of integrands (which is also possible with other methods, like substitution) without imposing an arbitrary limit on the recursion depth (which would hurt with integrands like `x**n * exp(x)`).