Repository: sympy/sympy

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.

Discussion:
sympy-bot: :white_check_mark:

Hi, I am the [SymPy bot](https://github.com/sympy/sympy-bot). I'm here to help you write a release notes entry. Please read the [guide on how to write release notes](https://github.com/sympy/sympy/wiki/Writing-Release-Notes).



* No release notes entry will be added for this pull request.

<details><summary>Click here to see the pull request description that was parsed.</summary>

    <!-- Your title above should be a short description of what
    was changed. Do not include the issue number in the title. -->

    #### References to other Issues or PRs
    <!-- If this pull request fixes an issue, write "Fixes #NNNN" in that exact
    format, e.g. "Fixes #1234" (see
    https://tinyurl.com/auto-closing for more information). Also, please
    write a comment on that issue linking back to this pull request once it is
    open. -->


    #### Brief description of what is fixed or changed


    #### Other comments


    #### Release Notes

    <!-- Write the release notes for this release below between the BEGIN and END
    statements. The basic format is a bulleted list with the name of the subpackage
    and the release note for this PR. For example:

    * solvers
      * Added a new solver for logarithmic equations.

    * functions
      * Fixed a bug with log of integers. Formerly, `log(-x)` incorrectly gave `-log(x)`.

    * physics.units
      * Corrected a semantical error in the conversion between volt and statvolt which
        reported the volt as being larger than the statvolt.

    or if no release note(s) should be included use:

    NO ENTRY

    See https://github.com/sympy/sympy/wiki/Writing-Release-Notes for more
    information on how to write release notes. The bot will check your release
    notes automatically to see if they are formatted correctly. -->

    <!-- BEGIN RELEASE NOTES -->
    NO ENTRY
    <!-- END RELEASE NOTES -->


</details><p>


Patch:
diff --git a/sympy/integrals/heurisch.py b/sympy/integrals/heurisch.py
index 344edf250a..2a1b61c27d 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 2b4ffa0684..e220132f58 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 6fa978f3bf..8605808968 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