Inefficiency in the Integrator with a Rational Expression

[andreubernadserra]

Hello everyone, I'm facing long run times when integrating a rational expression. The code to reproduce the issue is the following:

from sympy import symbols, Function, integrate, diff

p,T,nu = symbols("p T nu", positive=True)
R,b    = symbols("R b"   , positive=True)
a      = Function("a"    , positive=True)(T)

srk_p = R*T/(nu-b) - a/(nu**2 + b*nu)
integrand = p - T*diff(srk_p, T)

integrate(integrand, nu)

The traceback after 5min is the following:

Traceback (most recent call last):
File "", line 1, in
File "sympy\integrals\integrals.py", line 1566, in integrate
return integral.doit(**doit_flags)
File "sympy\integrals\integrals.py", line 613, in doit
antideriv = self._eval_integral(
File "sympy\integrals\integrals.py", line 956, in _eval_integral
result, i = risch_integrate(f, x, separate_integral=True,
File "sympy\integrals\risch.py", line 1835, in risch_integrate
ans = integrate(fa.as_expr()/fd.as_expr(), DE.x, risch=False)
File "sympy\integrals\integrals.py", line 1566, in integrate
return integral.doit(**doit_flags)
File "sympy\integrals\integrals.py", line 613, in doit
antideriv = self._eval_integral(
File "sympy\integrals\integrals.py", line 1041, in _eval_integral
parts.append(coeff * ratint(g, x))
File "sympy\integrals\rationaltools.py", line 112, in ratint
R = log_to_real(h, q, x, t)
File "sympy\integrals\rationaltools.py", line 375, in log_to_real
if len(R_u) != R.count_roots():
File "sympy\polys\polytools.py", line 3547, in count_roots
count = f.rep.count_real_roots(inf=inf, sup=sup)
File "sympy\polys\polyclasses.py", line 863, in count_real_roots
return dup_count_real_roots(f.rep, f.dom, inf=inf, sup=sup)
File "sympy\polys\rootisolation.py", line 788, in dup_count_real_roots
sturm = dup_sturm(f, K)
File "sympy\polys\rootisolation.py", line 65, in dup_sturm
s = dup_rem(sturm[-2], sturm[-1], K)
File "sympy\polys\densearith.py", line 1557, in dup_rem
return dup_div(f, g, K)[1]
File "sympy\polys\densearith.py", line 1534, in dup_div
return dup_ff_div(f, g, K)
File "sympy\polys\densearith.py", line 1444, in dup_ff_div
r = dup_sub(r, h, K)
File "sympy\polys\densearith.py", line 613, in dup_sub
return dup_strip([ a - b for a, b in zip(f, g) ])
File "sympy\polys\densearith.py", line 613, in
return dup_strip([ a - b for a, b in zip(f, g) ])
File "sympy\polys\domains\expressiondomain.py", line 89, in sub
return f.simplify(f.ex - g.ex)
File "sympy\polys\domains\expressiondomain.py", line 50, in simplify
return f.class(ex.cancel().expand(**eflags))
File "sympy\core\expr.py", line 3733, in cancel
return cancel(self, *gens, **args)
File "sympy\polys\polytools.py", line 6801, in cancel
c, (P, Q) = 1, F.cancel(G)
File "sympy\polys\rings.py", line 2223, in cancel
_, p, q = f.cofactors(g)
File "sympy\polys\rings.py", line 2139, in cofactors
h, cff, cfg = f._gcd(g)
File "sympy\polys\rings.py", line 2172, in _gcd
return f._gcd_ZZ(g)
File "sympy\polys\rings.py", line 2177, in _gcd_ZZ
return heugcd(f, g)
File "sympy\polys\heuristicgcd.py", line 80, in heugcd
h, cff, cfg = heugcd(ff, gg)
File "sympy\polys\heuristicgcd.py", line 80, in heugcd
h, cff, cfg = heugcd(ff, gg)
File "sympy\polys\heuristicgcd.py", line 80, in heugcd
h, cff, cfg = heugcd(ff, gg)
File "sympy\polys\heuristicgcd.py", line 61, in heugcd
gcd, f, g = f.extract_ground(g)
File "sympy\polys\rings.py", line 2040, in extract_ground
gcd = f.ring.domain.gcd(fc, gc)
File "sympy\polys\domains\integerring.py", line 212, in gcd
return gcd(a, b)
KeyboardInterrupt

Thanks a lot!

This is posted here as suggested in the Google Groups.

[oscarbenjamin]

I'm not really sure what count_roots is supposed to do here. The polynomial in question is this:

Poly(64*_u**9 + 192*R*T*_u**8 + (224*R**2*T**2 - 96*T**2*Derivative(a(T), T)**2/b**2)*_u**7 + (128*R**3
*T**3 - 256*R*T**3*Derivative(a(T), T)**2/b**2)*_u**6 + (36*R**4*T**4 - 264*R**2*T**4*Derivative(a(T), 
T)**2/b**2 + 36*T**4*Derivative(a(T), T)**4/b**4)*_u**5 + (4*R**5*T**5 - 136*R**3*T**5*Derivative(a(T),
 T)**2/b**2 + 68*R*T**5*Derivative(a(T), T)**4/b**4)*_u**4 + (-36*R**4*T**6*Derivative(a(T), T)**2/b**2
 + 40*R**2*T**6*Derivative(a(T), T)**4/b**4 - 4*T**6*Derivative(a(T), T)**6/b**6)*_u**3 + (-4*R**5*T**7
*Derivative(a(T), T)**2/b**2 + 8*R**3*T**7*Derivative(a(T), T)**4/b**4 - 4*R*T**7*Derivative(a(T), T)**
6/b**6)*_u**2, _u, domain='EX')

I don't see how you could meaningfully count the roots.

[oscarbenjamin]

I just ran this again with latest master and got a correct result:

In [2]: from sympy import symbols, Function, integrate, diff
   ...: 
   ...: p,T,nu = symbols("p T nu", positive=True)
   ...: R,b    = symbols("R b"   , positive=True)
   ...: a      = Function("a"    , positive=True)(T)
   ...: 
   ...: srk_p = R*T/(nu-b) - a/(nu**2 + b*nu)
   ...: integrand = p - T*diff(srk_p, T)
   ...: 
   ...: integrate(integrand, nu)
Out[2]: 
             ⎛                                                                  ⎛                  
             ⎜                                                                  ⎜   ⎛              
             ⎜                                   2                  2           ⎜   ⎜              
             ⎜ 3  2    2      2      2 ⎛d       ⎞       3 ⎛d       ⎞            ⎜ 2 ⎜              
ν⋅p + RootSum⎜t ⋅b  + t ⋅R⋅T⋅b  - t⋅T ⋅⎜──(a(T))⎟  - R⋅T ⋅⎜──(a(T))⎟ , t ↦ t⋅log⎜t ⋅⎜- ────────────
             ⎜                         ⎝dT      ⎠         ⎝dT      ⎠            ⎜   ⎜              
             ⎜                                                                  ⎜   ⎜     2  2  2 d
             ⎜                                                                  ⎜   ⎜  2⋅R ⋅T ⋅b ⋅─
             ⎝                                                                  ⎝   ⎝             d

                                                                                                   
                                            3 d                    ⎞                               
        4                                3⋅b ⋅──(a(T))             ⎟          2                 2  
     R⋅b                                      dT                   ⎟       t⋅b                 R ⋅b
────────────────────────── - ──────────────────────────────────────⎟ - ──────────── + ─────────────
                         3                                        3⎟       d                       
             2 ⎛d       ⎞       2  2  2 d             2 ⎛d       ⎞ ⎟   2⋅T⋅──(a(T))      2  2     ⎛
─(a(T)) - 2⋅T ⋅⎜──(a(T))⎟    2⋅R ⋅T ⋅b ⋅──(a(T)) - 2⋅T ⋅⎜──(a(T))⎟ ⎟       dT         2⋅R ⋅b  - 2⋅⎜
T              ⎝dT      ⎠               dT              ⎝dT      ⎠ ⎠                              ⎝

                                                         2        ⎞⎞
                     2 d                       ⎛d       ⎞         ⎟⎟
3                 R⋅b ⋅──(a(T))            2⋅b⋅⎜──(a(T))⎟         ⎟⎟
                       dT                      ⎝dT      ⎠         ⎟⎟
────────── + ─────────────────────── + ─────────────────────── + ν⎟⎟
         2                         2                         2    ⎟⎟
d       ⎞       2  2     ⎛d       ⎞       2  2     ⎛d       ⎞     ⎟⎟
──(a(T))⎟    2⋅R ⋅b  - 2⋅⎜──(a(T))⎟    2⋅R ⋅b  - 2⋅⎜──(a(T))⎟     ⎟⎟
dT      ⎠                ⎝dT      ⎠                ⎝dT      ⎠     ⎠⎠

Simplifying that gives:

In [9]: f.doit().simplify()
Out[9]: 
                            ⎛ T        -T⎞ d       
                         log⎝ν ⋅(b + ν)  ⎠⋅──(a(T))
                                           dT      
-R⋅T⋅log(-b + ν) + ν⋅p + ──────────────────────────
                                     b  

Both simplified and unsimplified results are correct:

In [10]: ratsimp(f.diff(nu) - integrand)
Out[10]: 0

In [14]: ratsimp(f.doit().simplify().diff(nu) - integrand)
Out[14]: 0

It is still slow though (it took 2 minutes). We should bisect to find what fixed this.

Also this should be kept open as a performance issue.

[oscarbenjamin]

We should bisect to find what fixed this.

I can't reproduce the original problem with older sympy versions (I tried 1.6 to 1.12).

[asmeurer]

The derivative constant term is really throwing it off. If you replace it with a dummy it gives an answer much faster:

>>> var('da')
da
>>> integrate(integrand.subs(diff(a, T), da), nu)
T*(-R*log(nu + (-2*R**2*T**2*b**3*da + 2*T**2*b*da**3)/(2*R**2*T**2*b**2*da - 2*T**2*da**3)) + da*log(nu)/b - da*log(nu + (2*R**2*T**2*b**3*da - 2*T**2*b*da**3)/(2*R**2*T**2*b**2*da - 2*T**2*da**3))/b) + nu*p

Also notice from the traceback that it's going through risch_integrate. It should go directly to ratint.

[oscarbenjamin]

Maybe the problem is with apart:

In [6]: integrand
Out[6]: 
    ⎛         d       ⎞    
    ⎜         ──(a(T))⎟    
    ⎜  R      dT      ⎟    
- T⋅⎜────── - ────────⎟ + p
    ⎜-b + ν          2⎟    
    ⎝         b⋅ν + ν ⎠    

In [7]: integrand.apart(nu)
Out[7]: 
  ⎛           2     d            d       ⎞    
T⋅⎜R⋅b⋅ν + R⋅ν  + b⋅──(a(T)) - ν⋅──(a(T))⎟    
  ⎝                 dT           dT      ⎠    
────────────────────────────────────────── + p
            ν⋅(b - ν)⋅(b + ν)                 

In [8]: integrand.subs(T, 1).apart(nu)
Out[8]: 
            ⎛d       ⎞│      ⎛d       ⎞│   
            ⎜──(a(T))⎟│      ⎜──(a(T))⎟│   
  R         ⎝dT      ⎠│T=1   ⎝dT      ⎠│T=1
───── + p - ────────────── + ──────────────
b - ν         b⋅(b + ν)           b⋅ν 

Although integrating that expression is slow as well for some reason. Profile:

In [10]: %prun -s cumulative integrand.subs(T, 1).apart(nu).integrate(nu)

         8999807 function calls (7727813 primitive calls) in 18.752 seconds

   Ordered by: cumulative time

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
    148/1    0.018    0.000   18.791   18.791 {built-in method builtins.exec}
        1    0.000    0.000   18.790   18.790 <string>:1(<module>)
        1    0.000    0.000   18.647   18.647 expr.py:3723(integrate)
      2/1    0.000    0.000   18.647   18.647 integrals.py:1399(integrate)
      2/1    0.000    0.000   18.646   18.646 integrals.py:382(doit)
      2/1    0.000    0.000   18.631   18.631 integrals.py:816(_eval_integral)
        1    0.000    0.000   18.629   18.629 risch.py:1706(risch_integrate)
        1    0.000    0.000   18.580   18.580 rationaltools.py:15(ratint)
       33    0.001    0.000   11.893    0.360 densearith.py:1410(dup_ff_div)
        1    0.000    0.000   11.414   11.414 rationaltools.py:327(log_to_real)
      245    0.000    0.000   11.282    0.046 densearith.py:1515(dup_div)
      897    0.002    0.000   11.042    0.012 densearith.py:140(dup_mul_term)
      897    0.001    0.000   11.040    0.012 densearith.py:157(<listcomp>)
     1032    0.003    0.000   10.665    0.010 rings.py:2219(cancel)
  977/972    0.003    0.000   10.656    0.011 rings.py:2140(cofactors)
      135    0.000    0.000   10.600    0.079 rings.py:2185(_gcd)
      134    0.001    0.000   10.598    0.079 rings.py:2195(_gcd_ZZ)
  507/134    0.026    0.000   10.597    0.079 heuristicgcd.py:7(heugcd)
        1    0.000    0.000   10.570   10.570 polytools.py:3501(count_roots)
        1    0.000    0.000   10.570   10.570 polyclasses.py:861(count_real_roots)
        1    0.000    0.000   10.570   10.570 rootisolation.py:779(dup_count_real_roots)
      675    0.002    0.000   10.535    0.016 fields.py:300(new)
        1    0.000    0.000   10.388   10.388 rootisolation.py:32(dup_sturm)
        7    0.000    0.000    9.858    1.408 densearith.py:1539(dup_rem)
      456    0.002    0.000    9.816    0.022 fields.py:490(__mul__)
      367    0.010    0.000    6.674    0.018 polytools.py:6804(cancel)
      360    0.002    0.000    6.643    0.018 expressiondomain.py:49(simplify)
      360    0.003    0.000    6.633    0.018 expr.py:3788(cancel)
        1    0.000    0.000    6.370    6.370 rationaltools.py:187(ratint_logpart)
     1182    5.736    0.005    5.785    0.005 rings.py:2310(evaluate)
      275    0.001    0.000    4.549    0.017 exprtools.py:1156(factor_terms)