Repository: sympy/sympy
Issue category: Performance Issue

Issue:
Potential hang in monomials during integration
In the example given below, the `intermonomials` function is called from the `heurisch` function of the integration process. The problem I encounter is that before calling `internomomials`, the script subs some expression with polynomials. The expressions look something like this, `(27000.0 - 6902.92767267848*exp(-1.66666666666667*s) - 24100.2086229868*exp(-s/120))*exp(k10*s)` and the mapping it uses to sub over this example would be, `[(exp(k10*s), _x1), (exp(-s/120), _x0), (exp(-1.66666666666667*s), _x3), (s, _x2)] `. As you can see, the third exponential is exactly the second one at the 200th power. So the subbed expression would come out as, `_x1*(-6902.92767267848*_x0**200 - 24100.2086229868*_x0 + 27000.0)`, which isn't false, but it causes the `intermonomials` to go on for hours (maybe days, I have never seen it end) in the `for variable in item:` loop at line 85.
```
import sympy as sy

Z0 = sy.Function('Z0')
Z1 = sy.Function('Z1')
k01, k10, t, s= sy.symbols('k01 k10 t s', real=True, positive=True)

m = sy.Matrix([[sy.exp(-k10*t)]])

[a, b, c, d, e, f, g] = [100, 0.5, 0.83, 50, 0.6, 2, 120]
AIF_btf = a*(d*e*(1-sy.exp(-(t-b)/e))+f*g*(1-sy.exp(-(t-b)/g)))
AIF_atf = a*(d*e*sy.exp(-(t-b)/e)*(sy.exp((c-b)/e)-1)+f*g*sy.exp(-(t-b)/g)*(sy.exp((c-b)/g)-1))
AIF_sym = sy.Piecewise((0, t < b), (AIF_btf, sy.And(b<=t, t<c)), (AIF_atf, c<=t))

aif_eq = sy.Eq(Z0(t), AIF_sym)
c_vec  = sy.Matrix([[Z1(t)]])
f_vec  = sy.Matrix([[k01*Z0(t)]])

integrand = m*(m.subs(t, s)**-1)*f_vec.subs(aif_eq.lhs, aif_eq.rhs).subs(t, s)
solution  = sy.integrate(integrand[0], (s, 0, t))
```
This behavior can be "hotfixed" in the `_subterms` function of `heurisch` by replacing the `return expr.subs(mapping)` of line 460 by `return expr.xreplace(dict(mapping))` but I know this probably would be far from an ideal fix. Maybe this doesn't actually need to be fixed and I'm just doing something wrong, too.

Now here is the interesting part, I first talked about this issue and #15494 in the gitter chat thinking they weren't related at all. But for some reason I don't understand, when applying the fix I just talked about, the behavior from #15494 disappeared. Also, since @normalhuman talked about how SymPy prefers working with known rationals, I tried replacing `[a, b, c, d, e, f, g] = [100, 0.5, 0.83, 50, 0.6, 2, 120]` with `[a, b, c, d, e, f, g] = [100, sy.Rational('0.5'), 0.83, 50, sy.Rational('0.6'), 2, 120]` and got NaN as the result of the integral. I then needed to use one of the fixes I proposed in #15494 to get the good result, but it all worked out in the end. Seeing as both issues seem related only in my code, I'm assuming my script is the source of the problem, but I've been searching for a while without success.


Hints:
> [a, b, c, d, e, f, g] = [100, sy.Rational('0.5'), 0.83, 50, sy.Rational('0.6'), 2, 120]

If there is a single floating point number `0.83`, then all numbers are considered inexact. That will be much slower than computing with exact numbers as the arithmetic laws are not valid, and therefore coefficients are treated as expressions in `heurisch`.
Actually, it seems that `nan` is created by `manualintegrate` as a result of division by zero. In principle, `heurisch` should be able to compute the integral but it does not succeed in this case, probably because of some unlucky choices in its cache. It looks like following patch would fix `manualintegrate`:
```
--- a/sympy/integrals/manualintegrate.py
+++ b/sympy/integrals/manualintegrate.py
@@ -130,6 +130,8 @@ def find_substitutions(integrand, symbol, u_var):
     results = []
 
     def test_subterm(u, u_diff):
+        if u_diff == 0:
+            return False
         substituted = integrand / u_diff
         if symbol not in substituted.free_symbols:
             # replaced everything already
```
> Actually, it seems that `nan` is created by `manualintegrate` as a result of division by zero.

Yes, I referenced this issue already in #15494. I had proposed two initial fixes, one being quite close to what you just wrote. I think I like your format better though.

Should we go ahead and make a PR for it? 

Patch:
diff --git a/sympy/integrals/heurisch.py b/sympy/integrals/heurisch.py
index 344edf250a2e..2a1b61c27da3 100644
--- a/sympy/integrals/heurisch.py
+++ b/sympy/integrals/heurisch.py
@@ -1,7 +1,8 @@
 from __future__ import annotations
 
-from itertools import permutations
+from collections import defaultdict
 from functools import reduce
+from itertools import permutations
 
 from sympy.core.add import Add
 from sympy.core.basic import Basic
@@ -503,7 +504,16 @@ def heurisch(f, x, rewrite=False, hints=None, mappings=None, retries=3,
         # optimizing the number of permutations of mapping              #
         assert mapping[-1][0] == x # if not, find it and correct this comment
         unnecessary_permutations = [mapping.pop(-1)]
-        mappings = permutations(mapping)
+        # only permute types of objects and let the ordering
+        # of types take care of the order of replacement
+        types = defaultdict(list)
+        for i in mapping:
+            types[type(i)].append(i)
+        mapping = [types[i] for i in types]
+        def _iter_mappings():
+            for i in permutations(mapping):
+                yield [j for i in i for j in i]
+        mappings = _iter_mappings()
     else:
         unnecessary_permutations = unnecessary_permutations or []
 


Test patch:
diff --git a/sympy/integrals/tests/test_heurisch.py b/sympy/integrals/tests/test_heurisch.py
index 2b4ffa0684f0..e220132f582f 100644
--- a/sympy/integrals/tests/test_heurisch.py
+++ b/sympy/integrals/tests/test_heurisch.py
@@ -2,7 +2,7 @@
 from sympy.core.add import Add
 from sympy.core.function import (Derivative, Function, diff)
 from sympy.core.numbers import (I, Rational, pi)
-from sympy.core.relational import Ne
+from sympy.core.relational import Eq, Ne
 from sympy.core.symbol import (Symbol, symbols)
 from sympy.functions.elementary.exponential import (LambertW, exp, log)
 from sympy.functions.elementary.hyperbolic import (asinh, cosh, sinh, tanh)
@@ -12,6 +12,8 @@
 from sympy.functions.special.bessel import (besselj, besselk, bessely, jn)
 from sympy.functions.special.error_functions import erf
 from sympy.integrals.integrals import Integral
+from sympy.logic.boolalg import And
+from sympy.matrices import Matrix
 from sympy.simplify.ratsimp import ratsimp
 from sympy.simplify.simplify import simplify
 from sympy.integrals.heurisch import components, heurisch, heurisch_wrapper
@@ -365,3 +367,33 @@ def f(x):
     Uz = integrate(f(z), z)
     Ut = integrate(f(t), t)
     assert Ut == Uz.subs(z, t)
+
+
+def test_heurisch_complex_erf_issue_26338():
+    r = symbols('r', real=True)
+    a = exp(-r**2/(2*(2 - I)**2))
+    assert heurisch(a, r, hints=[]) is None  # None, not a wrong soln
+    a = sqrt(pi)*erf((1 + I)/2)/2
+    assert integrate(exp(-I*r**2/2), (r, 0, 1)) == a - I*a
+
+    a = exp(-x**2/(2*(2 - I)**2))
+    assert heurisch(a, x, hints=[]) is None  # None, not a wrong soln
+    a = sqrt(pi)*erf((1 + I)/2)/2
+    assert integrate(exp(-I*x**2/2), (x, 0, 1)) == a - I*a
+
+
+def test_issue_15498():
+    Z0 = Function('Z0')
+    k01, k10, t, s= symbols('k01 k10 t s', real=True, positive=True)
+    m = Matrix([[exp(-k10*t)]])
+    _83 = Rational(83, 100)  # 0.83 works, too
+    [a, b, c, d, e, f, g] = [100, 0.5, _83, 50, 0.6, 2, 120]
+    AIF_btf = a*(d*e*(1 - exp(-(t - b)/e)) + f*g*(1 - exp(-(t - b)/g)))
+    AIF_atf = a*(d*e*exp(-(t - b)/e)*(exp((c - b)/e) - 1
+        ) + f*g*exp(-(t - b)/g)*(exp((c - b)/g) - 1))
+    AIF_sym = Piecewise((0, t < b), (AIF_btf, And(b <= t, t < c)), (AIF_atf, c <= t))
+    aif_eq = Eq(Z0(t), AIF_sym)
+    f_vec = Matrix([[k01*Z0(t)]])
+    integrand = m*m.subs(t, s)**-1*f_vec.subs(aif_eq.lhs, aif_eq.rhs).subs(t, s)
+    solution = integrate(integrand[0], (s, 0, t))
+    assert solution is not None  # does not hang and takes less than 10 s
diff --git a/sympy/integrals/tests/test_integrals.py b/sympy/integrals/tests/test_integrals.py
index 6fa978f3bfaa..860580896807 100644
--- a/sympy/integrals/tests/test_integrals.py
+++ b/sympy/integrals/tests/test_integrals.py
@@ -1148,7 +1148,7 @@ def test_issue_3940():
     assert integrate(exp(-x**2 + I*c*x), x) == \
         -sqrt(pi)*exp(-c**2/4)*erf(I*c/2 - x)/2
     assert integrate(exp(a*x**2 + b*x + c), x) == \
-        sqrt(pi)*exp(c)*exp(-b**2/(4*a))*erfi(sqrt(a)*x + b/(2*sqrt(a)))/(2*sqrt(a))
+        sqrt(pi)*exp(c - b**2/(4*a))*erfi((2*a*x + b)/(2*sqrt(a)))/(2*sqrt(a))
 
     from sympy.core.function import expand_mul
     from sympy.abc import k