sympy hangs on evaluating expression

[zachetienne]

From at least sympy 0.7.6 through the latest checkout (Nov 27, 2017 1.1.2-dev), the below minimal-ish example causes sympy to hang indefinitely.

from sympy import sqrt, symbols, sin, cos, log, exp

r,r_m  = symbols('r r_m', positive=True)
th = symbols('th', real=True)
a,M = symbols('a M', real=True)
kappa,gamma = symbols('kappa gamma', real=True)

l  = sqrt(M/r_m**3) * (r_m**4 + r_m**2*a**2 - 2*M*r_m*a**2 - a*sqrt(M*r_m)*(r_m**2-a**2))
l /= r_m**2 - 3*M*r_m + 2*a*sqrt(M*r_m)

Delta = r**2 - 2*M*r + a**2
Sigma = r**2 + a**2*cos(th)**2
A = (r**2 + a**2)**2 - Delta*a**2*sin(th)**2

tmp3 = sqrt(1 + l**2*Sigma**2*Delta/(sin(th))**2)
ln_h  = log( (1 + tmp3) / (Sigma*Delta/A) ) - tmp3
hm1 = exp(ln_h) - 1

print("sympy will hang here, computing rho0:")
rho0 = ( hm1*(gamma-1)/(kappa*gamma) )**(1/(gamma - 1))

[gxyd]

Thanks for reporting the issue. I verified, this still hangs on the development version. A fix is welcome.

[gxyd]

Here is the keyboard interrupt:

^C---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
<ipython-input-5-14f1d7d831a6> in <module>()
     18 
     19 print("sympy will hang here, computing rho0:")
---> 20 rho0 = ( hm1*(gamma-1)/(kappa*gamma) )**(1/(gamma - 1))

~/dev/sympy/sympy/core/expr.py in __pow__(self, other, mod)
    152     def __pow__(self, other, mod=None):
    153         if mod is None:
--> 154             return self._pow(other)
    155         try:
    156             _self, other, mod = as_int(self), as_int(other), as_int(mod)

~/dev/sympy/sympy/core/decorators.py in __sympifyit_wrapper(a, b)
     89                 if not hasattr(b, '_op_priority'):
     90                     b = sympify(b, strict=True)
---> 91                 return func(a, b)
     92             except SympifyError:
     93                 return retval

~/dev/sympy/sympy/core/decorators.py in binary_op_wrapper(self, other)
    130                     else:
    131                         return f(self)
--> 132             return func(self, other)
    133         return binary_op_wrapper
    134     return priority_decorator

~/dev/sympy/sympy/core/expr.py in _pow(self, other)
    148     @call_highest_priority('__rpow__')
    149     def _pow(self, other):
--> 150         return Pow(self, other)
    151 
    152     def __pow__(self, other, mod=None):

~/dev/sympy/sympy/core/cache.py in wrapper(*args, **kwargs)
     91             def wrapper(*args, **kwargs):
     92                 try:
---> 93                     retval = cfunc(*args, **kwargs)
     94                 except TypeError:
     95                     retval = func(*args, **kwargs)

/usr/lib/python3.4/functools.py in wrapper(*args, **kwds)
    470                         hits += 1
    471                         return result
--> 472                 result = user_function(*args, **kwds)
    473                 with lock:
    474                     if key in cache:

~/dev/sympy/sympy/core/power.py in __new__(cls, b, e, evaluate)
    220                         return S.Exp1**(c*numer(ex))
    221                     elif den.is_Add:
--> 222                         s = sign(im(b))
    223                         if s.is_Number and s and den == \
    224                                 log(-factor_terms(b, sign=False)) + s*S.ImaginaryUnit*S.Pi:

~/dev/sympy/sympy/core/cache.py in wrapper(*args, **kwargs)
     91             def wrapper(*args, **kwargs):
     92                 try:
---> 93                     retval = cfunc(*args, **kwargs)
     94                 except TypeError:
     95                     retval = func(*args, **kwargs)

/usr/lib/python3.4/functools.py in wrapper(*args, **kwds)
    470                         hits += 1
    471                         return result
--> 472                 result = user_function(*args, **kwds)
    473                 with lock:
    474                     if key in cache:

~/dev/sympy/sympy/core/function.py in __new__(cls, *args, **options)
    437 
    438         evaluate = options.get('evaluate', global_evaluate[0])
--> 439         result = super(Function, cls).__new__(cls, *args, **options)
    440         if evaluate and isinstance(result, cls) and result.args:
    441             pr2 = min(cls._should_evalf(a) for a in result.args)

~/dev/sympy/sympy/core/cache.py in wrapper(*args, **kwargs)
     91             def wrapper(*args, **kwargs):
     92                 try:
---> 93                     retval = cfunc(*args, **kwargs)
     94                 except TypeError:
     95                     retval = func(*args, **kwargs)

/usr/lib/python3.4/functools.py in wrapper(*args, **kwds)
    470                         hits += 1
    471                         return result
--> 472                 result = user_function(*args, **kwds)
    473                 with lock:
    474                     if key in cache:

~/dev/sympy/sympy/core/function.py in __new__(cls, *args, **options)
    248 
    249         if evaluate:
--> 250             evaluated = cls.eval(*args)
    251             if evaluated is not None:
    252                 return evaluated

~/dev/sympy/sympy/functions/elementary/complexes.py in eval(cls, arg)
    174                     # impossible, don't try to do im(arg) again
    175                     # (because this is what we are trying to do now).
--> 176                     real_imag = term.as_real_imag(ignore=arg)
    177                     if real_imag:
    178                         excluded.append(real_imag[1])

~/dev/sympy/sympy/core/mul.py in as_real_imag(self, deep, **hints)
    781         addterms = S.One
    782         for a in self.args:
--> 783             r, i = a.as_real_imag()
    784             if i.is_zero:
    785                 coeffr.append(r)

~/dev/sympy/sympy/core/add.py in as_real_imag(self, deep, **hints)
    821         re_part, im_part = [], []
    822         for term in sargs:
--> 823             re, im = term.as_real_imag(deep=deep)
    824             re_part.append(re)
    825             im_part.append(im)

~/dev/sympy/sympy/core/mul.py in as_real_imag(self, deep, **hints)
    781         addterms = S.One
    782         for a in self.args:
--> 783             r, i = a.as_real_imag()
    784             if i.is_zero:
    785                 coeffr.append(r)

~/dev/sympy/sympy/functions/elementary/exponential.py in as_real_imag(self, deep, **hints)
    354         re, im = self.args[0].as_real_imag()
    355         if deep:
--> 356             re = re.expand(deep, **hints)
    357             im = im.expand(deep, **hints)
    358         cos, sin = sympy.cos(im), sympy.sin(im)

~/dev/sympy/sympy/core/cache.py in wrapper(*args, **kwargs)
     91             def wrapper(*args, **kwargs):
     92                 try:
---> 93                     retval = cfunc(*args, **kwargs)
     94                 except TypeError:
     95                     retval = func(*args, **kwargs)

/usr/lib/python3.4/functools.py in wrapper(*args, **kwds)
    470                         hits += 1
    471                         return result
--> 472                 result = user_function(*args, **kwds)
    473                 with lock:
    474                     if key in cache:

~/dev/sympy/sympy/core/expr.py in expand(self, deep, modulus, power_base, power_exp, mul, log, multinomial, basic, **hints)
   3111             if hints.get('mul', False):
   3112                 expr, _ = Expr._expand_hint(
-> 3113                     expr, '_eval_expand_mul', deep=deep, **hints)
   3114             if hints.get('log', False):
   3115                 expr, _ = Expr._expand_hint(

~/dev/sympy/sympy/core/expr.py in _expand_hint(expr, hint, deep, **hints)
   3035             sargs = []
   3036             for arg in expr.args:
-> 3037                 arg, arghit = Expr._expand_hint(arg, hint, **hints)
   3038                 hit |= arghit
   3039                 sargs.append(arg)

~/dev/sympy/sympy/core/expr.py in _expand_hint(expr, hint, deep, **hints)
   3035             sargs = []
   3036             for arg in expr.args:
-> 3037                 arg, arghit = Expr._expand_hint(arg, hint, **hints)
   3038                 hit |= arghit
   3039                 sargs.append(arg)

~/dev/sympy/sympy/core/expr.py in _expand_hint(expr, hint, deep, **hints)
   3035             sargs = []
   3036             for arg in expr.args:
-> 3037                 arg, arghit = Expr._expand_hint(arg, hint, **hints)
   3038                 hit |= arghit
   3039                 sargs.append(arg)

~/dev/sympy/sympy/core/expr.py in _expand_hint(expr, hint, deep, **hints)
   3043 
   3044         if hasattr(expr, hint):
-> 3045             newexpr = getattr(expr, hint)(**hints)
   3046             if newexpr != expr:
   3047                 return (newexpr, True)

~/dev/sympy/sympy/core/mul.py in _eval_expand_mul(self, **hints)
    855         if d.is_Mul:
    856             n, d = [i._eval_expand_mul(**hints) if i.is_Mul else i
--> 857                 for i in (n, d)]
    858             expr = n/d
    859             if not expr.is_Mul:

~/dev/sympy/sympy/core/mul.py in <listcomp>(.0)
    855         if d.is_Mul:
    856             n, d = [i._eval_expand_mul(**hints) if i.is_Mul else i
--> 857                 for i in (n, d)]
    858             expr = n/d
    859             if not expr.is_Mul:

~/dev/sympy/sympy/core/mul.py in _eval_expand_mul(self, **hints)
    877             if sums:
    878                 deep = hints.get("deep", False)
--> 879                 terms = self.func._expandsums(sums)
    880                 args = []
    881                 for term in terms:

~/dev/sympy/sympy/core/mul.py in _expandsums(sums)
    842         right = Mul._expandsums(sums[L//2:])
    843 
--> 844         terms = [Mul(a, b) for a in left for b in right]
    845         added = Add(*terms)
    846         return Add.make_args(added)  # it may have collapsed down to one term

~/dev/sympy/sympy/core/mul.py in <listcomp>(.0)
    842         right = Mul._expandsums(sums[L//2:])
    843 
--> 844         terms = [Mul(a, b) for a in left for b in right]
    845         added = Add(*terms)
    846         return Add.make_args(added)  # it may have collapsed down to one term

~/dev/sympy/sympy/core/cache.py in wrapper(*args, **kwargs)
     91             def wrapper(*args, **kwargs):
     92                 try:
---> 93                     retval = cfunc(*args, **kwargs)
     94                 except TypeError:
     95                     retval = func(*args, **kwargs)

/usr/lib/python3.4/functools.py in wrapper(*args, **kwds)
    470                         hits += 1
    471                         return result
--> 472                 result = user_function(*args, **kwds)
    473                 with lock:
    474                     if key in cache:

~/dev/sympy/sympy/core/operations.py in __new__(cls, *args, **options)
     43         is_commutative = not nc_part
     44         obj = cls._from_args(c_part + nc_part, is_commutative)
---> 45         obj = cls._exec_constructor_postprocessors(obj)
     46 
     47         if order_symbols is not None:

~/dev/sympy/sympy/core/basic.py in _exec_constructor_postprocessors(cls, obj)
   1672                         if cls in Basic._constructor_postprocessor_mapping
   1673                     )
-> 1674                     for k, v in chain.from_iterable(postprocessor_mappings):
   1675                         postprocessors[k].extend([j for j in v if j not in postprocessors[k]])
   1676             except TypeError:

~/dev/sympy/sympy/core/basic.py in <genexpr>(.0)
   1670                         Basic._constructor_postprocessor_mapping[cls].items()
   1671                         for cls in type(i).mro()
-> 1672                         if cls in Basic._constructor_postprocessor_mapping
   1673                     )
   1674                     for k, v in chain.from_iterable(postprocessor_mappings):

KeyboardInterrupt: 

[smichr]

A more minimal working example:

>>> ok = hm1**(1/gamma)
>>> ok = x**(1/(gamma-1))
>>> bad = hm1**(1/(gamma - 1))

[oscarbenjamin]

This is due to

from sympy import *

a, r, M, r_m, th = symbols('a, r, M, r_m, th')

expr = (-a**2*(-2*M*r + a**2 + r**2)*sin(th)**2 + (a**2 + r**2)**2)*(sqrt(M*(a**2*cos(th)**2 + r**2)**2*(-2*M*r + a**2 + r**2)*(-sqrt(M)*a*sqrt(r_m)*(-a**2 + r_m**2) - 2*M*a**2*r_m + a**2*r_m**2 + r_m**4)**2/(r_m**3*(2*sqrt(M)*a*sqrt(r_m) - 3*M*r_m + r_m**2)**2*sin(th)**2) + 1) + 1)*exp(-sqrt(M*(a**2*cos(th)**2 + r**2)**2*(-2*M*r + a**2 + r**2)*(-sqrt(M)*a*sqrt(r_m)*(-a**2 + r_m**2) - 2*M*a**2*r_m + a**2*r_m**2 + r_m**4)**2/(r_m**3*(2*sqrt(M)*a*sqrt(r_m) - 3*M*r_m + r_m**2)**2*sin(th)**2) + 1))/((a**2*cos(th)**2 + r**2)*(-2*M*r + a**2 + r**2))

print(expr.as_real_imag()) # hangs