exception sympy.polys.polyerrors.PolificationFailed: Cannot construct polynomials using integrate in sympy 1.10.1

[nasser1]

I am getting this error. I do not know if this is known issue already or not.

>python
Python 3.10.5 (main, Jun  6 2022, 18:49:26) [GCC 12.1.0] on linux
Type "help", "copyright", "credits" or "license" for more information.
>>>from sympy import *
>>>F, G, H, a, b, c, d, e, f, g, h, i, k, m, n, p, r, s, t, x=symbols('F G H a b c d e f g h i k m n p r s t x')
>>>integrand=((-x*exp(exp(5))*ln(x)+(x**2+13*x+42)*exp(exp(5)))*exp(ln(x)/(6+x))+((-2*x**4-24*x**3-73*x**2-12*x-36)*exp(x**2)+(-x**2-12*x-36)*exp(2))*exp(exp(5)))/(x**2+12*x+36)
>>>integrate(integrand,x)

Gives

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/usr/lib/python3.10/site-packages/sympy/integrals/integrals.py", line 1566, in integrate
    return integral.doit(**doit_flags)
  File "/usr/lib/python3.10/site-packages/sympy/integrals/integrals.py", line 613, in doit
    antideriv = self._eval_integral(
  File "/usr/lib/python3.10/site-packages/sympy/integrals/integrals.py", line 956, in _eval_integral
    result, i = risch_integrate(f, x, separate_integral=True,
  File "/usr/lib/python3.10/site-packages/sympy/integrals/risch.py", line 1829, in risch_integrate
    ans, i, b = integrate_hyperexponential(fa, fd, DE, conds=conds)
  File "/usr/lib/python3.10/site-packages/sympy/integrals/risch.py", line 1578, in integrate_hyperexponential
    (integrate((p - i).subs(DE.t, 1).subs(s), DE.x), True)
  File "/usr/lib/python3.10/site-packages/sympy/integrals/integrals.py", line 1566, in integrate
    return integral.doit(**doit_flags)
  File "/usr/lib/python3.10/site-packages/sympy/integrals/integrals.py", line 613, in doit
    antideriv = self._eval_integral(
  File "/usr/lib/python3.10/site-packages/sympy/integrals/integrals.py", line 956, in _eval_integral
    result, i = risch_integrate(f, x, separate_integral=True,
  File "/usr/lib/python3.10/site-packages/sympy/integrals/risch.py", line 1831, in risch_integrate
    ans, i, b = integrate_primitive(fa, fd, DE)
  File "/usr/lib/python3.10/site-packages/sympy/integrals/risch.py", line 1460, in integrate_primitive
    q, i, b = integrate_primitive_polynomial(p, DE)
  File "/usr/lib/python3.10/site-packages/sympy/integrals/risch.py", line 1409, in integrate_primitive_polynomial
    rv = limited_integrate(aa, ad, [(Dta, Dtb)], DE)
  File "/usr/lib/python3.10/site-packages/sympy/integrals/prde.py", line 837, in limited_integrate
    h, A = param_rischDE(Fa, Fd, G, DE)
  File "/usr/lib/python3.10/site-packages/sympy/integrals/prde.py", line 766, in param_rischDE
    M = Matrix([wl[:m] + wl[-v:] for wl in W])  # excise dj's.
  File "/usr/lib/python3.10/site-packages/sympy/polys/polymatrix.py", line 83, in __new__
    return cls.from_list(rows, cols, items, gens, ring)
  File "/usr/lib/python3.10/site-packages/sympy/polys/polymatrix.py", line 106, in from_list
    items, info = parallel_poly_from_expr(items, gens, field=True)
  File "/usr/lib/python3.10/site-packages/sympy/polys/polytools.py", line 4414, in parallel_poly_from_expr
    return _parallel_poly_from_expr(exprs, opt)
  File "/usr/lib/python3.10/site-packages/sympy/polys/polytools.py", line 4469, in _parallel_poly_from_expr
    raise PolificationFailed(opt, origs, exprs, True)
sympy.polys.polyerrors.PolificationFailed: Cannot construct polynomials from 
>>> 

Versions

>which python
/usr/bin/python
>python --version
Python 3.10.5
>python
Python 3.10.5 (main, Jun  6 2022, 18:49:26) [GCC 12.1.0] on linux
Type "help", "copyright", "credits" or "license" for more information.
>>> import sympy 
>>> sympy.__version__
'1.10.1'
>>> 

[oscarbenjamin]

I tried this with latest master and while it took some time (maybe 10 minutes) an answer was returned:

In [1]: >>> from sympy import *
   ...: >>> F, G, H, a, b, c, d, e, f, g, h, i, k, m, n, p, r, s, t, x=symbols('F G H a b c d e f g h
   ...:  i k m n p r s t x')
   ...: >>> integrand=((-x*exp(exp(5))*ln(x)+(x**2+13*x+42)*exp(exp(5)))*exp(ln(x)/(6+x))+((-2*x**4-2
   ...: 4*x**3-73*x**2-12*x-36)*exp(x**2)+(-x**2-12*x-36)*exp(2))*exp(exp(5)))/(x**2+12*x+36)
   ...: >>> integrate(integrand,x)
Out[1]: 
                     log(x)                   
     ⎛ 2⎞  ⎛ 5⎞      ──────  ⎛ 5⎞         ⎛ 5⎞
     ⎝x ⎠  ⎝ℯ ⎠      x + 6   ⎝ℯ ⎠      2  ⎝ℯ ⎠
- x⋅ℯ    ⋅ℯ     + x⋅ℯ      ⋅ℯ     - x⋅ℯ ⋅ℯ  

Can you try with the master branch from git? Alternatively you could try the 1.11rc1 pre-release which can be installed with:

pip install --pre --upgrade sympy

[oscarbenjamin]

This was fixed in ec8a19d from #23717

[oscarbenjamin]

I do wonder if it should be so slow though. Maybe it's possible to improve that.

[nasser1]

"I do wonder if it should be so slow though. Maybe it's possible to improve that."

10 minutes seems a little long, I agree. It will be nice if performance is improved, yes. But this is not as important as getting correct results and avoiding exceptions.

Fyi, Mathematica 13.1 does this in 0.1 seconds

integrand = ((-x*Exp[Exp[5]]*Log[x] + (x^2 + 13*x + 42)*Exp[Exp[5]])*
     Exp[Log[x]/(6 + x)] + ((-2*x^4 - 24*x^3 - 73*x^2 - 12*x - 36)*
        Exp[x^2] + (-x^2 - 12*x - 36)*Exp[2])*Exp[Exp[5]])/(x^2 + 
    12*x + 36)
Timing[Integrate[integrand, x]]

Gives

{0.15625,-E^E^5 x (E^2+E^x^2-x^(1/(6+x)))}

[oscarbenjamin]

Gives

{0.15625,-E^E^5 x (E^2+E^x^2-x^(1/(6+x)))}

This looks the same apart from x^(1/(6+x)) rather than e^(log(x)/(x + 6)). Those look generically equivalent although I might have missed something.

The log form is invalid for x=0 although it works in the limit as x -> 0+. Otherwise are they both equivalent for all complex numbers?

[oscarbenjamin]

The profiler shows what takes the time. Looks like the main thing is risch/heurisch which are polynomial heavy. Converting from polynomials (from_expr) comes out high. Manual integrate is also slow here.

In [3]: %prun -s cumulative integrate(integrand, x)
 
        764581812 function calls (736746730 primitive calls) in 578.001 seconds

   Ordered by: cumulative time

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
   4856/1    0.570    0.000  615.103  615.103 {built-in method builtins.exec}
      4/1    0.000    0.000  615.103  615.103 integrals.py:1399(integrate)
     11/1    0.001    0.000  615.100  615.100 integrals.py:382(doit)
     11/1    0.001    0.000  615.042  615.042 integrals.py:816(_eval_integral)
      2/1    0.000    0.000  614.869  614.869 risch.py:1706(risch_integrate)
        2    0.000    0.000  613.975  306.988 risch.py:1525(integrate_hyperexponential)
        7    0.000    0.000  500.435   71.491 heurisch.py:109(heurisch_wrapper)
     32/7    1.030    0.032  500.434   71.491 heurisch.py:295(heurisch)
       62    0.046    0.001  493.709    7.963 heurisch.py:635(_integrate)
94349/251   19.221    0.000  371.897    1.482 rings.py:393(from_expr)
94349/251    0.283    0.000  371.895    1.482 rings.py:370(_rebuild_expr)
457511/251    3.160    0.000  371.894    1.482 rings.py:373(_rebuild)
68287/68283    0.826    0.000  318.365    0.005 {built-in method _functools.reduce}
3472664/341224   31.588    0.000  305.135    0.001 rings.py:1071(__mul__)
 10191589    9.682    0.000  199.093    0.000 fields.py:490(__mul__)
  1886408    3.422    0.000  194.602    0.000 fields.py:300(new)
  1928031    3.929    0.000  191.883    0.000 rings.py:2200(cancel)
  1899366    2.816    0.000  175.192    0.000 rings.py:2121(cofactors)
  1898377   14.877    0.000  171.722    0.000 rings.py:2150(_gcd_monom)
5896508/4215105   24.160    0.000  121.278    0.000 cache.py:67(wrapper)
  5701960   29.490    0.000  104.183    0.000 gaussiandomains.py:185(__divmod__)
        1    0.000    0.000  100.805  100.805 manualintegrate.py:1864(manualintegrate)
    642/1    0.031    0.000  100.767  100.767 manualintegrate.py:1404(integral_steps)
   1594/1    0.007    0.000  100.733  100.733 core.py:80(do_one_rl)
   4204/5    0.010    0.000  100.733   20.147 core.py:61(null_safe_rl)
   196/41    0.001    0.000   98.426    2.401 operations.py:453(doit)
   196/41    0.001    0.000   98.421    2.401 operations.py:455(<listcomp>)
        1    0.000    0.000   98.418   98.418 integrals.py:1118(<listcomp>)
    191/1    0.008    0.000   93.428   93.428 manualintegrate.py:347(_alternatives)
    635/4    0.002    0.000   93.376   23.344 core.py:90(switch_rl)
    764/8    0.002    0.000   84.530   10.566 core.py:28(conditioned_rl)
    764/8    0.006    0.000   84.286   10.536 manualintegrate.py:307(_rewriter)
     45/2    0.000    0.000   84.010   42.005 manualintegrate.py:559(add_rule)
     45/2    0.000    0.000   84.001   42.000 manualintegrate.py:561(<listcomp>)
   207/26    0.002    0.000   83.387    3.207 manualintegrate.py:566(mul_rule)
1358419/704035   11.880    0.000   78.564    0.000 rings.py:916(__add__)
  3804654    2.762    0.000   73.458    0.000 ring.py:26(quo)
95995/17769    0.490    0.000   72.997    0.004 basic.py:1702(match)
79760/4174    2.188    0.000   71.422    0.017 operations.py:192(_matches_commutative)
   191/55    0.004    0.000   71.314    1.297 manualintegrate.py:1297(substitution_rule)
  3794390    2.828    0.000   70.696    0.000 gaussiandomains.py:152(__floordiv__)
60041/13268    0.173    0.000   67.693    0.005 mul.py:1018(matches)
19719/5870    0.042    0.000   66.643    0.011 add.py:521(matches)
  1897700    3.958    0.000   65.733    0.000 gaussiandomains.py:459(gcd)
  289/141    0.012    0.000   62.652    0.444 manualintegrate.py:806(quadratic_denom_rule)
  2727649    3.368    0.000   61.222    0.000 fields.py:395(__add__)
7794/5978    0.049    0.000   57.528    0.010 power.py:1579(matches)
    65/30    0.001    0.000   52.945    1.765 manualintegrate.py:670(parts_rule)
11866039/11383887    6.452    0.000   51.696    0.000 domain.py:403(convert)
148481/126706    2.116    0.000   51.184    0.000 operations.py:52(__new__)
     16/1    0.000    0.000   50.045   50.045 manualintegrate.py:722(make_second_step)