(Integration) regressions?

[mpeychev]

Hello, I am trying to integrate the following expressions (f) which work under SymPy v1.4 but fail on both the master version and 1.7.1.

Behavior in v1.4 is:

>>> from sympy import *
>>> x = Symbol('x', real=True, nonzero=True)

>>> f1 = x*(-x**4/asin(5)**4 - x*sinh(x + log(asin(5))) + 5)
>>> integrate(f1, x)
-x**6/(6*asin(5)**4) - x**2*cosh(x + log(asin(5))) + 5*x**2/2 + 2*x*sinh(x + log(asin(5))) - 2*cosh(x + log(asin(5)))

>>> f2 = (x + cosh(cos(4)))/(x*(x + 1/(12*x)))
>>> integrate(f2, x)
log(x**2 + 1/12)/2 + 2*sqrt(3)*cosh(cos(4))*atan(2*sqrt(3)*x)

For f1: v1.7.1 and master raise the following error: RecursionError: maximum recursion depth exceeded
For f2: v1.7.1 raises CoercionFailed error (issue possibly related to #21024), while master fails with RecursionError: maximum recursion depth exceeded

[oscarbenjamin]

For f1 I've bisected this to fc9a7dd from #17770

[oscarbenjamin]

The problem is this leads to infinite recursion:

In [2]: e = x*(-x**4/asin(5)**4 - x*sinh(x + log(asin(5))) + 5)                                                                   

In [3]: e.is_zero                                                                                                                 

That gets stuck in

sinh(x + log(asin(5))).is_extended_positive 

It seems that this freezes:

In [3]: e = -I*log((re(asin(5)) + I*im(asin(5)))/sqrt(re(asin(5))**2 + im(asin(5))**2))                                           

In [4]: e % pi                                                                                                                    
^C

That freezes because this does:

In [7]: e2 = -I*log((re(asin(5)) + I*im(asin(5)))/sqrt(re(asin(5))**2 + im(asin(5))**2))/pi                                       

In [8]: e2.round(2)                                                                                                               

There seems to be an infinite recursion on these lines:

sympy/sympy/core/expr.py

Lines 3770 to 3772 in 4c8422b

	if not x.is_extended_real:
	r, i = x.as_real_imag()
	return r.round(n) + S.ImaginaryUnit*i.round(n)

The problem is self==r. This expression seems to be real but fails the is_extended_real test.

The reason the is_extended_real test fails is because it gives None. The bug here is that we should not use of not x.is_y with three-valued-logic as in the old assumptions because not None gives True. The test should be

if x.is_extended_real is False:

That gives this diff:

diff --git a/sympy/core/expr.py b/sympy/core/expr.py
index 38c423a70f..1c7ef70b6a 100644
--- a/sympy/core/expr.py
+++ b/sympy/core/expr.py
@@ -3767,7 +3767,7 @@ def round(self, n=None):
                     'Expected a number but got %s:' % func_name(x))
         elif x in (S.NaN, S.Infinity, S.NegativeInfinity, S.ComplexInfinity):
             return x
-        if not x.is_extended_real:
+        if x.is_extended_real is False:
             r, i = x.as_real_imag()
             return r.round(n) + S.ImaginaryUnit*i.round(n)
         if not x:

With that diff we get:

In [1]: f1 = x*(-x**4/asin(5)**4 - x*sinh(x + log(asin(5))) + 5)

In [2]: integrate(f1, x)
Out[2]: 
       6                                      2                                                        
      x         2                          5⋅x                                                         
- ────────── - x ⋅cosh(x + log(asin(5))) + ──── + 2⋅x⋅sinh(x + log(asin(5))) - 2⋅cosh(x + log(asin(5)))
        4                                   2                                                          
  6⋅asin (5) 

I'm going to mark this as easy to fix although the issue with f2 might be something else...

[oscarbenjamin]

To fix this issue the diff above should be applied and also tests should be added for all of the issues that I highlighted (as well as the integral for f1).

For f2 more consideration is needed...

[nisargsc]

Hey, @oscarbenjamin, I have applied the diff you mentioned. Where should I add tests? I have added in sympy/core/test_expr.py but while running ./bin/test core it gives NameError: name 'asin' is not defined
This is my first contribution please help.

[oscarbenjamin]

You'll need to import the asin function:

from sympy.functions.elementary.trigonometric import asin

I would add the test in sympy/core/tests/test_arit.py.

[nisargsc]

Thank you

[prajwalranjan]

Hi, I'd like to try to fix this issue

[mpeychev]

Dear Prajwal,

To sum up, as far as I understood, the current status is:

• The problem with f1is fixed by RecursionError in integrate fixed #21053.
• The problem with f2 is still open.

[oscarbenjamin]

A first step would be to use git bisect to find when in between 1.7.1 and the current sympy master branch the CoercionFailed changed into a RecursionError.

[oscarbenjamin]

The RecursionError for f1 comes from this:

In [7]: Mod(4*x, I*pi)
---------------------------------------------------------------------------
RecursionError 

That gives a RecursionError in sympy 1.7.1 as well as master so the change between 1.7.1 and current master must be due to something else but the bug is there in both.

The infinite recursion comes from here:

sympy/sympy/core/mod.py

Line 221 in b564d9b

return G*cls(p, q, evaluate=(p, q) != (pwas, qwas))

The guard evaluate=(p, q) != (pwas, qwas) is intended to prevent infinite recursion but it doesn't work because in this example p = -pwas and q = -qwas:

ipdb> p p
4*x
ipdb> p q
I*pi
ipdb> p pwas
-4*x
ipdb> p qwas
-I*pi

[akshanshbhatt]

Integration of f2 was failing on master because of the changes made in the sympy.solvers.solvers.solve_linear_system. When compared with ver 1.4, sympy.solvers.solvers.solve_linear_system used to have the flag simplify, which was removed in PR #18814. This flag was important because when this keyword argument was passed True, the elements of the system matrix were simplified before being solved by being converted into rref. For this specific expression (f2), the solution set returned by solve_linear_system was fairly complex due to sympy.solvers.solvers.solve (refer to the blob below) unnecessarily converting cosh(cos(4)) into 6*exp(exp(4*I)/2 + exp(-4*I)/2) + 6*exp(-exp(-4*I)/2 - exp(4*I)/2)). Although this was happening in ver 1.4 too, but it was later simplified back to cosh(cos(4)) by the simplify flag pointed above. The recursion error was due to this fairly complex expression passed on for the later stages of computation.

sympy/sympy/solvers/solvers.py

Lines 913 to 917 in 7053279

	# rewrite hyperbolics in terms of exp
	f[i] = f[i].replace(lambda w: isinstance(w, HyperbolicFunction),
	lambda w: w.rewrite(exp))

I tested the integral by introducing this flag and necessary code it works fine again.

[abbyjng]

I'm interested in working on this issue. Is the problem with f2 still available?

Edit: Nevermind, my apologies, I didn't see the most recent comment when I posted this.

[oscarbenjamin]

I think that the bug in Mod(4*x, I*pi) should be fixed. It isn't clear that the simplify flag would be needed if that was fixed.

[akshanshbhatt]

I saw this behavior of the Mod function, but RecursionError is particularly for this case. I tried a different integral on master this time, and the result was still not comprehensible.

>>> from sympy import *
>>> x = Symbol('x', real=True, nonzero=True)
>>> f2 = (x + cosh(sin(4)))/(x*(x + 1/(12*x)))
>>> print(integrate(f2, x))
(1/2 - sqrt(3)*sqrt(-exp(I*exp(-4*I))*exp(I*exp(4*I)))*(exp(I*exp(4*I)) + exp(I*exp(-4*I)))*exp(-I*exp(-4*I))*exp(-I*exp(4*I))/2)*log(x + (-exp(7*I*exp(-4*I)/2)*exp(I*exp(4*I)/2) + 2*(1/2 - sqrt(3)*sqrt(-exp(I*exp(-4*I))*exp(I*exp(4*I)))*(exp(I*exp(4*I)) + exp(I*exp(-4*I)))*exp(-I*exp(-4*I))*exp(-I*exp(4*I))/2)*exp(7*I*exp(-4*I)/2)*exp(I*exp(4*I)/2) - 3*exp(5*I*exp(-4*I)/2)*exp(3*I*exp(4*I)/2) - 3*exp(3*I*exp(-4*I)/2)*exp(5*I*exp(4*I)/2) + 6*(1/2 - sqrt(3)*sqrt(-exp(I*exp(-4*I))*exp(I*exp(4*I)))*(exp(I*exp(4*I)) + exp(I*exp(-4*I)))*exp(-I*exp(-4*I))*exp(-I*exp(4*I))/2)*exp(5*I*exp(-4*I)/2)*exp(3*I*exp(4*I)/2) - exp(I*exp(-4*I)/2)*exp(7*I*exp(4*I)/2) + 6*(1/2 - sqrt(3)*sqrt(-exp(I*exp(-4*I))*exp(I*exp(4*I)))*(exp(I*exp(4*I)) + exp(I*exp(-4*I)))*exp(-I*exp(-4*I))*exp(-I*exp(4*I))/2)*exp(3*I*exp(-4*I)/2)*exp(5*I*exp(4*I)/2) + 2*(1/2 - sqrt(3)*sqrt(-exp(I*exp(-4*I))*exp(I*exp(4*I)))*(exp(I*exp(4*I)) + exp(I*exp(-4*I)))*exp(-I*exp(-4*I))*exp(-I*exp(4*I))/2)*exp(I*exp(-4*I)/2)*exp(7*I*exp(4*I)/2))/(6*exp(4*I*exp(4*I)) + 24*exp(I*exp(-4*I))*exp(3*I*exp(4*I)) + 36*exp(2*I*exp(-4*I))*exp(2*I*exp(4*I)) + 24*exp(3*I*exp(-4*I))*exp(I*exp(4*I)) + 6*exp(4*I*exp(-4*I)))) + (1/2 + sqrt(3)*sqrt(-exp(I*exp(-4*I))*exp(I*exp(4*I)))*(exp(I*exp(4*I)) + exp(I*exp(-4*I)))*exp(-I*exp(-4*I))*exp(-I*exp(4*I))/2)*log(x + (2*(1/2 + sqrt(3)*sqrt(-exp(I*exp(-4*I))*exp(I*exp(4*I)))*(exp(I*exp(4*I)) + exp(I*exp(-4*I)))*exp(-I*exp(-4*I))*exp(-I*exp(4*I))/2)*exp(I*exp(-4*I)/2)*exp(7*I*exp(4*I)/2) + 6*(1/2 + sqrt(3)*sqrt(-exp(I*exp(-4*I))*exp(I*exp(4*I)))*(exp(I*exp(4*I)) + exp(I*exp(-4*I)))*exp(-I*exp(-4*I))*exp(-I*exp(4*I))/2)*exp(3*I*exp(-4*I)/2)*exp(5*I*exp(4*I)/2) + 6*(1/2 + sqrt(3)*sqrt(-exp(I*exp(-4*I))*exp(I*exp(4*I)))*(exp(I*exp(4*I)) + exp(I*exp(-4*I)))*exp(-I*exp(-4*I))*exp(-I*exp(4*I))/2)*exp(5*I*exp(-4*I)/2)*exp(3*I*exp(4*I)/2) + 2*(1/2 + sqrt(3)*sqrt(-exp(I*exp(-4*I))*exp(I*exp(4*I)))*(exp(I*exp(4*I)) + exp(I*exp(-4*I)))*exp(-I*exp(-4*I))*exp(-I*exp(4*I))/2)*exp(7*I*exp(-4*I)/2)*exp(I*exp(4*I)/2) - exp(7*I*exp(-4*I)/2)*exp(I*exp(4*I)/2) - 3*exp(5*I*exp(-4*I)/2)*exp(3*I*exp(4*I)/2) - 3*exp(3*I*exp(-4*I)/2)*exp(5*I*exp(4*I)/2) - exp(I*exp(-4*I)/2)*exp(7*I*exp(4*I)/2))/(6*exp(4*I*exp(4*I)) + 24*exp(I*exp(-4*I))*exp(3*I*exp(4*I)) + 36*exp(2*I*exp(-4*I))*exp(2*I*exp(4*I)) + 24*exp(3*I*exp(-4*I))*exp(I*exp(4*I)) + 6*exp(4*I*exp(-4*I))))

The actual result should have been,

log(x**2 + 1/12)/2 + 2*sqrt(3)*cosh(sin(4))*atan(2*sqrt(3)*x).

So, in my opinion, the problem with Mod might be specific to this problem. The snippet in my previous comment (solvers.py blob) is causing this issue. The author who made that commit might be thinking about expanding those hyperbolic functions to exp that contained the symbols for which the system was to be solved. Expanding constants like cosh(cos(4)) not only convolute the system further but also takes additional computation time.

[oscarbenjamin]

It would be better to avoid/limit the rewrite than to simplify later on. The rewrite could be changed to something like:

In [5]: f = cosh(x) + cosh(4)

In [6]: f.replace(lambda e: isinstance(e, cosh) and x in e.free_symbols, lambda e: e.rewrite(cosh, exp))
Out[6]: 
 x              -x
ℯ              ℯ  
── + cosh(4) + ───
2               2