Repository: sympy/sympy
Pull request: Fix handling of zoo in evalf
<!-- 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. -->

Fixes #8242
Fixes #21147
Closes #21183
#### Brief description of what is fixed or changed

Fix by allowing `evalf`s to return `zoo` in addition to mpf value tuple, as there's no equivalent mpc representation in mpmath. Check this return value in each `evalf`s and handle it.
#### Other comments

Some additional logics are added to `equals` and `is_constant` in order to fix newly broken tests. One case is leaving unfixed, which imho is acceptable because it’s a bit complicated and previously based on wrong evalf intermediate result.

#### 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.

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 -->
* core
  * Evaluation of expressions involving zoo is now handled properly in evalf without raising obscure exceptions.
<!-- END RELEASE NOTES -->


Problem statement:
(1/x).evalf(subs={x: 0}) produces 0
Here's how it looks like:

``` python
In [1]: sp.__version__
Out[1]: '0.7.5-git'

In [2]: x = sp.abc.x

In [3]: (1/x).evalf(subs={x: 0})
Out[3]: 0

In [4]: (1/x).subs({x: 0})
Out[4]: zoo

In [5]: (1/x).xreplace({x: 0})
Out[5]: zoo
```

I did investigate it some time ago and, AFAIR, it happens in `Pow` evaluation that short-circuits to 0 whenever the base is 0 while it should return infs if the exponent is negative.


Evalf and evaluate=False issues
`Mul(2, Pow(0,-1, evaluate=False), evaluate=False).evalf()` returns 0 instead of complex infinity.

Hints from issue discussion:
Another interesting and related issue: while

``` python
In [3]: x = Symbol('x', real=True)

In [4]: (1/x).subs({x: 0.})
Out[4]: +inf

In [5]: (1/x).xreplace({x: 0.})
Out[5]: +inf
```

the following

``` python
(1/x).evalf(subs={x: 0.})
```

raises

``` pytb
Traceback (most recent call last):
  File "<ipython-input-6-226388149257>", line 1, in <module>
    (1/x).evalf(subs={x: 0.})
  File "sympy/core/evalf.py", line 1326, in evalf
    result = evalf(self, prec + 4, options)
  File "sympy/core/evalf.py", line 1220, in evalf
    r = rf(x, prec, options)
  File "sympy/core/evalf.py", line 622, in evalf_pow
    return mpf_pow_int(re, p, target_prec), None, target_prec, None
  File "sympy/mpmath/libmp/libmpf.py", line 1044, in mpf_pow_int
    if n == -1: return mpf_div(fone, s, prec, rnd)
  File "sympy/mpmath/libmp/libmpf.py", line 934, in mpf_div
    raise ZeroDivisionError
ZeroDivisionError
```

Valid for sympy-1.0 and sympy-1.1-58-gd250ff0, too, as noted in #12904. There are many related issues:

#10662
#9807
#7867
#7192
#6974


I get these results

```
>>> (x*y).n(subs={y:0})
0
>>> log(x).evalf(subs={x: 0})
zoo
>>> (x-1).evalf(subs={x:1})
0.e-125
>>> (x**I).evalf(subs={x:0})
nan
```

with these mods (which are suggestive as to where the problem lies)

```diff
diff --git a/sympy/core/evalf.py b/sympy/core/evalf.py
index c794656..233bbe5 100644
--- a/sympy/core/evalf.py
+++ b/sympy/core/evalf.py
@@ -821,6 +821,8 @@ def evalf_log(expr, prec, options):
 
     imaginary_term = (mpf_cmp(xre, fzero) < 0)
 
+    if xre == (0,)*4:
+        return S.ComplexInfinity
     re = mpf_log(mpf_abs(xre), prec, rnd)
     size = fastlog(re)
     if prec - size > workprec and re != fzero:
@@ -1284,13 +1286,25 @@ def evalf(x, prec, options):
     try:
         rf = evalf_table[x.func]
         r = rf(x, prec, options)
-    except KeyError:
+        if r is S.ComplexInfinity:
+            return r
+        if x and r[0] == r[1] == None:
+            assert not 'subs' in options
+    except (KeyError, AssertionError):
         try:
             # Fall back to ordinary evalf if possible
             if 'subs' in options:
                 x = x.subs(evalf_subs(prec, options['subs']))
-            xe = x._eval_evalf(prec)
-            re, im = xe.as_real_imag()
+            if x is S.ComplexInfinity:
+                return x
+            if x:
+                if not x.is_number:
+                    return (None,)*4
+                x = evalf(x, prec, options)
+                xe = x._eval_evalf(prec)
+                re, im = xe.as_real_imag()
+            else:
+                return (0, 0, 0, 0), None, None, None
             if re.has(re_) or im.has(im_):
                 raise NotImplementedError
             if re == 0:
@@ -1306,6 +1320,8 @@ def evalf(x, prec, options):
                 im = im._to_mpmath(prec, allow_ints=False)._mpf_
                 imprec = prec
             r = re, im, reprec, imprec
+            if re == im == None:
+                return (0, 0, 0, 0), None, None, None
         except AttributeError:
             raise NotImplementedError
     if options.get("verbose"):
@@ -1392,6 +1408,8 @@ def evalf(self, n=15, subs=None, maxn=100, chop=False, strict=False, quad=None,
             options['quad'] = quad
         try:
             result = evalf(self, prec + 4, options)
+            if result is S.ComplexInfinity:
+                return result
         except NotImplementedError:
             # Fall back to the ordinary evalf
             v = self._eval_evalf(prec)
```

Whoever tackles this beware: there is the `evalf` method and the `evalf` function (which references the `evalf_table` (of special evalf functions) and the `_eval_evalf` method. This is not a beginner's issue...unless you are a special beginner!
I believe the problem is [here](https://github.com/sympy/sympy/blob/b42c1c6645cc1b0d8e8a6c9c9e86c6979ef2c76d/sympy/core/evalf.py#L673). `evalf_pow` assumes that 0 to an integer power is 0 (`None`s represent exact zeros). But this is only true if the exponent is > 0. The question is, what *should* it return. I guess we want it to give zoo, but what value to we return from that function to give that (there is no floating point version of zoo). I guess your patch is in that direction, although I don't like that it makes the evalf functions return just `zoo` instead of a tuple.

We could also make it return `inf` quite easily:

```patch
diff --git a/sympy/core/evalf.py b/sympy/core/evalf.py
index c79465651..9437c9c17 100644
--- a/sympy/core/evalf.py
+++ b/sympy/core/evalf.py
@@ -670,7 +670,12 @@ def evalf_pow(v, prec, options):
                 return None, mpf_neg(z), None, target_prec
         # Zero raised to an integer power
         if not re:
-            return None, None, None, None
+            if exp > 0:
+                return None, None, None, None
+            elif exp < 0:
+                return mpmath_inf._mpf_, None, target_prec, None
+            else:
+                return mpf(1)._mpf_, None, target_prec, None
         # General complex number to arbitrary integer power
         re, im = libmp.mpc_pow_int((re, im), p, prec)
         # Assumes full accuracy in input
```

Maybe we should apply the above as a simpler fix for now (`inf` is not as correct as zoo, but it's definitely less wrong than 0). 
SO [reports this issue](https://stackoverflow.com/questions/45637615), too.
I want to work on this. Ideally mpmath should provide a tuple that represents `zoo`, but it doesn't. And there is no benefit that we assign a tuple as `zoo` then check this tuple everywhere, so I'd like to follow @smichr 's idea that directly return the `zoo` when appropriate.

Commit message: fix issue 8242

Patch to review:
diff --git a/sympy/core/evalf.py b/sympy/core/evalf.py
index c9313200a434..c39e52f4637b 100644
--- a/sympy/core/evalf.py
+++ b/sympy/core/evalf.py
@@ -12,7 +12,7 @@
     make_mpc, make_mpf, mp, mpc, mpf, nsum, quadts, quadosc, workprec)
 from mpmath import inf as mpmath_inf
 from mpmath.libmp import (from_int, from_man_exp, from_rational, fhalf,
-        fnan, fnone, fone, fzero, mpf_abs, mpf_add,
+        fnan, finf, fninf, fnone, fone, fzero, mpf_abs, mpf_add,
         mpf_atan, mpf_atan2, mpf_cmp, mpf_cos, mpf_e, mpf_exp, mpf_log, mpf_lt,
         mpf_mul, mpf_neg, mpf_pi, mpf_pow, mpf_pow_int, mpf_shift, mpf_sin,
         mpf_sqrt, normalize, round_nearest, to_int, to_str)
@@ -206,6 +206,8 @@ def complex_accuracy(result):
     In the worst case (re and im equal), this is wrong by a factor
     sqrt(2), or by log2(sqrt(2)) = 0.5 bit.
     """
+    if result is S.ComplexInfinity:
+        return INF
     re, im, re_acc, im_acc = result
     if not im:
         if not re:
@@ -221,8 +223,10 @@ def complex_accuracy(result):
 
 
 def get_abs(expr, prec, options):
-    re, im, re_acc, im_acc = evalf(expr, prec + 2, options)
-
+    result = evalf(expr, prec + 2, options)
+    if result is S.ComplexInfinity:
+        return finf, None, prec, None
+    re, im, re_acc, im_acc = result
     if not re:
         re, re_acc, im, im_acc = im, im_acc, re, re_acc
     if im:
@@ -246,6 +250,8 @@ def get_complex_part(expr, no, prec, options):
     i = 0
     while 1:
         res = evalf(expr, workprec, options)
+        if res is S.ComplexInfinity:
+            return fnan, None, prec, None
         value, accuracy = res[no::2]
         # XXX is the last one correct? Consider re((1+I)**2).n()
         if (not value) or accuracy >= prec or -value[2] > prec:
@@ -289,6 +295,8 @@ def chop_parts(value, prec):
     """
     Chop off tiny real or complex parts.
     """
+    if value is S.ComplexInfinity:
+        return value
     re, im, re_acc, im_acc = value
     # Method 1: chop based on absolute value
     if re and re not in _infs_nan and (fastlog(re) < -prec + 4):
@@ -323,7 +331,10 @@ def get_integer_part(expr, no, options, return_ints=False):
     from sympy.functions.elementary.complexes import re, im
     # The expression is likely less than 2^30 or so
     assumed_size = 30
-    ire, iim, ire_acc, iim_acc = evalf(expr, assumed_size, options)
+    result = evalf(expr, assumed_size, options)
+    if result is S.ComplexInfinity:
+        raise ValueError("Cannot get integer part of Complex Infinity")
+    ire, iim, ire_acc, iim_acc = result
 
     # We now know the size, so we can calculate how much extra precision
     # (if any) is needed to get within the nearest integer
@@ -532,10 +543,17 @@ def evalf_add(v, prec, options):
         options['maxprec'] = min(oldmaxprec, 2*prec)
 
         terms = [evalf(arg, prec + 10, options) for arg in v.args]
+        n = terms.count(S.ComplexInfinity)
+        if n >= 2:
+            return fnan, None, prec, None
         re, re_acc = add_terms(
-            [a[0::2] for a in terms if a[0]], prec, target_prec)
+            [a[0::2] for a in terms if isinstance(a, tuple) and a[0]], prec, target_prec)
         im, im_acc = add_terms(
-            [a[1::2] for a in terms if a[1]], prec, target_prec)
+            [a[1::2] for a in terms if isinstance(a, tuple) and a[1]], prec, target_prec)
+        if n == 1:
+            if re in (finf, fninf, fnan) or im in (finf, fninf, fnan):
+                return fnan, None, prec, None
+            return S.ComplexInfinity
         acc = complex_accuracy((re, im, re_acc, im_acc))
         if acc >= target_prec:
             if options.get('verbose'):
@@ -568,19 +586,30 @@ def evalf_mul(v, prec, options):
     args = list(v.args)
 
     # see if any argument is NaN or oo and thus warrants a special return
+    has_zero = False
     special = []
     from .numbers import Float
     for arg in args:
-        arg = evalf(arg, prec, options)
-        if arg[0] is None:
+        result = evalf(arg, prec, options)
+        if result is S.ComplexInfinity:
+            special.append(result)
             continue
-        arg = Float._new(arg[0], 1)
-        if arg is S.NaN or arg.is_infinite:
-            special.append(arg)
+        if result[0] is None:
+            if result[1] is None:
+                has_zero = True
+            continue
+        num = Float._new(result[0], 1)
+        if num is S.NaN:
+            return fnan, None, prec, None
+        if num.is_infinite:
+            special.append(num)
     if special:
+        if has_zero:
+            return fnan, None, prec, None
         from .mul import Mul
-        special = Mul(*special)
-        return evalf(special, prec + 4, {})
+        return evalf(Mul(*special), prec + 4, {})
+    if has_zero:
+        return None, None, None, None
 
     # With guard digits, multiplication in the real case does not destroy
     # accuracy. This is also true in the complex case when considering the
@@ -688,7 +717,12 @@ def evalf_pow(v, prec, options):
         # Exponentiation by p magnifies relative error by |p|, so the
         # base must be evaluated with increased precision if p is large
         prec += int(math.log(abs(p), 2))
-        re, im, re_acc, im_acc = evalf(base, prec + 5, options)
+        result = evalf(base, prec + 5, options)
+        if result is S.ComplexInfinity:
+            if p < 0:
+                return None, None, None, None
+            return result
+        re, im, re_acc, im_acc = result
         # Real to integer power
         if re and not im:
             return mpf_pow_int(re, p, target_prec), None, target_prec, None
@@ -706,15 +740,25 @@ def evalf_pow(v, prec, options):
                 return None, mpf_neg(z), None, target_prec
         # Zero raised to an integer power
         if not re:
+            if p < 0:
+                return S.ComplexInfinity
             return None, None, None, None
         # General complex number to arbitrary integer power
         re, im = libmp.mpc_pow_int((re, im), p, prec)
         # Assumes full accuracy in input
         return finalize_complex(re, im, target_prec)
 
+    result = evalf(base, prec + 5, options)
+    if result is S.ComplexInfinity:
+        if exp.is_Rational:
+            if exp < 0:
+                return None, None, None, None
+            return result
+        raise NotImplementedError
+
     # Pure square root
     if exp is S.Half:
-        xre, xim, _, _ = evalf(base, prec + 5, options)
+        xre, xim, _, _ = result
         # General complex square root
         if xim:
             re, im = libmp.mpc_sqrt((xre or fzero, xim), prec)
@@ -730,7 +774,10 @@ def evalf_pow(v, prec, options):
     # We first evaluate the exponent to find its magnitude
     # This determines the working precision that must be used
     prec += 10
-    yre, yim, _, _ = evalf(exp, prec, options)
+    result = evalf(exp, prec, options)
+    if result is S.ComplexInfinity:
+        return fnan, None, prec, None
+    yre, yim, _, _ = result
     # Special cases: x**0
     if not (yre or yim):
         return fone, None, prec, None
@@ -752,6 +799,10 @@ def evalf_pow(v, prec, options):
     xre, xim, _, _ = evalf(base, prec + 5, options)
     # 0**y
     if not (xre or xim):
+        if yim:
+            return fnan, None, prec, None
+        if yre[0] == 1:  # y < 0
+            return S.ComplexInfinity
         return None, None, None, None
 
     # (real ** complex) or (complex ** complex)
@@ -846,7 +897,10 @@ def evalf_log(expr, prec, options):
         return evalf(expr, prec, options)
     arg = expr.args[0]
     workprec = prec + 10
-    xre, xim, xacc, _ = evalf(arg, workprec, options)
+    result = evalf(arg, workprec, options)
+    if result is S.ComplexInfinity:
+        return result
+    xre, xim, xacc, _ = result
 
     # evalf can return NoneTypes if chop=True
     # issue 18516, 19623
@@ -1203,11 +1257,11 @@ def summand(k, _term=[term0]):
 
 def evalf_prod(expr, prec, options):
     if all((l[1] - l[2]).is_Integer for l in expr.limits):
-        re, im, re_acc, im_acc = evalf(expr.doit(), prec=prec, options=options)
+        result = evalf(expr.doit(), prec=prec, options=options)
     else:
         from sympy.concrete.summations import Sum
-        re, im, re_acc, im_acc = evalf(expr.rewrite(Sum), prec=prec, options=options)
-    return re, im, re_acc, im_acc
+        result = evalf(expr.rewrite(Sum), prec=prec, options=options)
+    return result
 
 
 def evalf_sum(expr, prec, options):
@@ -1285,7 +1339,8 @@ def _create_evalf_table():
     from sympy.concrete.summations import Sum
     from .add import Add
     from .mul import Mul
-    from .numbers import Exp1, Float, Half, ImaginaryUnit, Integer, NaN, NegativeOne, One, Pi, Rational, Zero
+    from .numbers import Exp1, Float, Half, ImaginaryUnit, Integer, NaN, NegativeOne, One, Pi, Rational, \
+        Zero, ComplexInfinity
     from .power import Pow
     from .symbol import Dummy, Symbol
     from sympy.functions.elementary.complexes import Abs, im, re
@@ -1307,6 +1362,7 @@ def _create_evalf_table():
         Exp1: lambda x, prec, options: (mpf_e(prec), None, prec, None),
         ImaginaryUnit: lambda x, prec, options: (None, fone, None, prec),
         NegativeOne: lambda x, prec, options: (fnone, None, prec, None),
+        ComplexInfinity: lambda x, prec, options: S.ComplexInfinity,
         NaN: lambda x, prec, options: (fnan, None, prec, None),
 
         exp: lambda x, prec, options: evalf_pow(
@@ -1405,7 +1461,7 @@ def evalf(x, prec, options):
 
     if options.get("verbose"):
         print("### input", x)
-        print("### output", to_str(r[0] or fzero, 50))
+        print("### output", to_str(r[0] or fzero, 50) if isinstance(r, tuple) else r)
         print("### raw", r) # r[0], r[2]
         print()
     chop = options.get('chop', False)
@@ -1536,6 +1592,8 @@ def evalf(self, n=15, subs=None, maxn=100, chop=False, strict=False, quad=None,
             except NotImplementedError:
                 # Probably contains symbols or unknown functions
                 return v
+        if result is S.ComplexInfinity:
+            return result
         re, im, re_acc, im_acc = result
         if re is S.NaN or im is S.NaN:
             return S.NaN
@@ -1571,7 +1629,10 @@ def _to_mpmath(self, prec, allow_ints=True):
         if hasattr(self, '_as_mpf_val'):
             return make_mpf(self._as_mpf_val(prec))
         try:
-            re, im, _, _ = evalf(self, prec, {})
+            result = evalf(self, prec, {})
+            if result is S.ComplexInfinity:
+                raise NotImplementedError
+            re, im, _, _ = result
             if im:
                 if not re:
                     re = fzero
diff --git a/sympy/core/expr.py b/sympy/core/expr.py
index c9c9e8150c6d..fb010a533183 100644
--- a/sympy/core/expr.py
+++ b/sympy/core/expr.py
@@ -651,7 +651,7 @@ def check_denominator_zeros(expression):
         if expr.is_zero:
             return True
 
-        # Don't attempt subsitution or differentiation with non-number symbols
+        # Don't attempt substitution or differentiation with non-number symbols
         wrt_number = {sym for sym in wrt if sym.kind is NumberKind}
 
         # try numerical evaluation to see if we get two different values
@@ -700,9 +700,6 @@ def check_denominator_zeros(expression):
                 if not (pure_complex(deriv, or_real=True)):
                     if flags.get('failing_number', False):
                         return failing_number
-                    elif deriv.free_symbols:
-                        # dead line provided _random returns None in such cases
-                        return None
                 return False
         cd = check_denominator_zeros(self)
         if cd is True:
@@ -752,6 +749,12 @@ def equals(self, other, failing_expression=False):
             # then this can't be a zero
             return False
 
+        factors = diff.as_coeff_mul()[1]
+        if len(factors) > 1:  # avoid infinity recursion
+            fac_zero = [fac.equals(0) for fac in factors]
+            if None not in fac_zero:  # every part can be decided
+                return any(fac_zero)
+
         constant = diff.is_constant(simplify=False, failing_number=True)
 
         if constant is False:
diff --git a/sympy/core/tests/test_evalf.py b/sympy/core/tests/test_evalf.py
index cd099e7d05b8..af8667cf5f98 100644
--- a/sympy/core/tests/test_evalf.py
+++ b/sympy/core/tests/test_evalf.py
@@ -1,3 +1,5 @@
+import math
+
 from sympy.concrete.products import (Product, product)
 from sympy.concrete.summations import Sum
 from sympy.core.add import Add
@@ -5,7 +7,7 @@
 from sympy.core.function import (Function, nfloat)
 from sympy.core.mul import Mul
 from sympy.core import (GoldenRatio)
-from sympy.core.numbers import (E, I, Rational, oo, pi)
+from sympy.core.numbers import (E, I, Rational, oo, zoo, nan, pi)
 from sympy.core.power import Pow
 from sympy.core.relational import Eq
 from sympy.core.singleton import S
@@ -13,7 +15,7 @@
 from sympy.core.sympify import sympify
 from sympy.functions.combinatorial.factorials import factorial
 from sympy.functions.combinatorial.numbers import fibonacci
-from sympy.functions.elementary.complexes import Abs
+from sympy.functions.elementary.complexes import (Abs, re, im)
 from sympy.functions.elementary.exponential import (exp, log)
 from sympy.functions.elementary.hyperbolic import (acosh, cosh)
 from sympy.functions.elementary.integers import (ceiling, floor)
@@ -38,6 +40,7 @@ def NS(e, n=15, **options):
 
 
 def test_evalf_helpers():
+    from mpmath.libmp import finf
     assert complex_accuracy((from_float(2.0), None, 35, None)) == 35
     assert complex_accuracy((from_float(2.0), from_float(10.0), 35, 100)) == 37
     assert complex_accuracy(
@@ -45,6 +48,9 @@ def test_evalf_helpers():
     assert complex_accuracy((from_float(2.0), from_float(10.0), 100, 35)) == 35
     assert complex_accuracy(
         (from_float(2.0), from_float(1000.0), 100, 35)) == 35
+    assert complex_accuracy(finf) == math.inf
+    assert complex_accuracy(zoo) == math.inf
+    raises(ValueError, lambda: get_integer_part(zoo, 1, {}))
 
 
 def test_evalf_basic():
@@ -496,8 +502,7 @@ def test_issue_6632_evalf():
 
 def test_issue_4945():
     from sympy.abc import H
-    from sympy.core.numbers import zoo
-    assert (H/0).evalf(subs={H:1}) == zoo*H
+    assert (H/0).evalf(subs={H:1}) == zoo
 
 
 def test_evalf_integral():
@@ -628,3 +633,28 @@ def test_issue_20291():
 
     sol = Complement(Intersection(FiniteSet(-b/2 - sqrt(b**2-4*pi)/2), Reals), FiniteSet(0))
     assert sol.evalf(subs={b: 1}) == EmptySet
+
+
+def test_evalf_with_zoo():
+    assert (1/x).evalf(subs={x: 0}) == zoo  # issue 8242
+    assert (-1/x).evalf(subs={x: 0}) == zoo  # PR 16150
+    assert (0 ** x).evalf(subs={x: -1}) == zoo  # PR 16150
+    assert (0 ** x).evalf(subs={x: -1 + I}) == nan
+    assert Mul(2, Pow(0, -1, evaluate=False), evaluate=False).evalf() == zoo  # issue 21147
+    assert Mul(x, 1/x, evaluate=False).evalf(subs={x: 0}) == Mul(x, 1/x, evaluate=False).subs(x, 0) == nan
+    assert Mul(1/x, 1/x, evaluate=False).evalf(subs={x: 0})
+    assert Mul(1/x, Abs(1/x), evaluate=False).evalf(subs={x: 0}) == zoo
+    assert Abs(zoo, evaluate=False).evalf() == oo
+    assert re(zoo, evaluate=False).evalf() == nan
+    assert im(zoo, evaluate=False).evalf() == nan
+    assert Add(zoo, zoo, evaluate=False).evalf() == nan
+    assert Add(oo, zoo, evaluate=False).evalf() == nan
+    assert Pow(zoo, -1, evaluate=False).evalf() == 0
+    assert Pow(zoo, Rational(-1, 3), evaluate=False).evalf() == 0
+    assert Pow(zoo, Rational(1, 3), evaluate=False).evalf() == zoo
+    assert Pow(zoo, S.Half, evaluate=False).evalf() == zoo
+    assert Pow(zoo, 2, evaluate=False).evalf() == zoo
+    assert Pow(0, zoo, evaluate=False).evalf() == nan
+    assert log(zoo, evaluate=False).evalf() == zoo
+    assert zoo.evalf(chop=True) == zoo
+    assert x.evalf(subs={x: zoo}) == zoo
diff --git a/sympy/core/tests/test_expr.py b/sympy/core/tests/test_expr.py
index 6982c827e70a..19458f30c751 100644
--- a/sympy/core/tests/test_expr.py
+++ b/sympy/core/tests/test_expr.py
@@ -1797,16 +1797,17 @@ def test_issue_5843():
 
 def test_is_constant():
     from sympy.solvers.solvers import checksol
-    Sum(x, (x, 1, 10)).is_constant() is True
-    Sum(x, (x, 1, n)).is_constant() is False
-    Sum(x, (x, 1, n)).is_constant(y) is True
-    Sum(x, (x, 1, n)).is_constant(n) is False
-    Sum(x, (x, 1, n)).is_constant(x) is True
+    assert Sum(x, (x, 1, 10)).is_constant() is True
+    assert Sum(x, (x, 1, n)).is_constant() is False
+    assert Sum(x, (x, 1, n)).is_constant(y) is True
+    assert Sum(x, (x, 1, n)).is_constant(n) is False
+    assert Sum(x, (x, 1, n)).is_constant(x) is True
     eq = a*cos(x)**2 + a*sin(x)**2 - a
-    eq.is_constant() is True
+    assert eq.is_constant() is True
     assert eq.subs({x: pi, a: 2}) == eq.subs({x: pi, a: 3}) == 0
     assert x.is_constant() is False
     assert x.is_constant(y) is True
+    assert log(x/y).is_constant() is False
 
     assert checksol(x, x, Sum(x, (x, 1, n))) is False
     assert checksol(x, x, Sum(x, (x, 1, n))) is False
@@ -1852,7 +1853,7 @@ def test_equals():
         2*sqrt(2)*x**(S(3)/2)*(1 + 1/(2*x))**(S(5)/2)/(-6 - 3/x)
     ans = sqrt(2*x + 1)*(6*x**2 + x - 1)/15
     diff = i - ans
-    assert diff.equals(0) is False
+    assert diff.equals(0) is None  # should be False, but previously this was False due to wrong intermediate result
     assert diff.subs(x, Rational(-1, 2)/2) == 7*sqrt(2)/120
     # there are regions for x for which the expression is True, for
     # example, when x < -1/2 or x > 0 the expression is zero
@@ -2084,7 +2085,7 @@ def test_issue_6325():
         (a + b*t)**2 + (c + t*z)**2))/sqrt((a + b*t)**2 + (c + t*z)**2)
     e = sqrt((a + b*t)**2 + (c + z*t)**2)
     assert diff(e, t, 2) == ans
-    e.diff(t, 2) == ans
+    assert e.diff(t, 2) == ans
     assert diff(e, t, 2, simplify=False) != ans
 
 
diff --git a/sympy/functions/elementary/exponential.py b/sympy/functions/elementary/exponential.py
index 1b5dd7c27ba5..cc4aa6000730 100644
--- a/sympy/functions/elementary/exponential.py
+++ b/sympy/functions/elementary/exponential.py
@@ -1055,6 +1055,9 @@ def _eval_as_leading_term(self, x, logx=None, cdir=0):
             return self.func(x0) - 2*I*S.Pi
         return self.func(arg)
 
+    def is_constant(self, *wrt, **flags):
+        return self.args[0].is_constant(*wrt, **flags)
+
 
 class LambertW(Function):
     r"""


Code review comments:
diff_hunk: @@ -1055,6 +1055,9 @@ def _eval_as_leading_term(self, x, logx=None, cdir=0):
             return self.func(x0) - 2*I*S.Pi
         return self.func(arg)
 
+    def is_constant(self, *wrt, **flags):
+        return self.args[0].is_constant(*wrt, **flags)
path: sympy/functions/elementary/exponential.py
text: @user1:
Why is this needed for `log`? No other function class has this.

@author:
It was imported when fixing test_issue_22020 in test_expr.py
Seems without it that case can still pass now. I'll remove it to see what CI says.

diff_hunk: @@ -628,3 +633,28 @@ def test_issue_20291():
 
     sol = Complement(Intersection(FiniteSet(-b/2 - sqrt(b**2-4*pi)/2), Reals), FiniteSet(0))
     assert sol.evalf(subs={b: 1}) == EmptySet
+
+
+def test_evalf_with_zoo():
+    assert (1/x).evalf(subs={x: 0}) == zoo  # issue 8242
+    assert (-1/x).evalf(subs={x: 0}) == zoo  # PR 16150
+    assert (0 ** x).evalf(subs={x: -1}) == zoo  # PR 16150
+    assert (0 ** x).evalf(subs={x: -1 + I}) == nan
+    assert Mul(2, Pow(0, -1, evaluate=False), evaluate=False).evalf() == zoo  # issue 21147
+    assert Mul(x, 1/x, evaluate=False).evalf(subs={x: 0}) == Mul(x, 1/x, evaluate=False).subs(x, 0) == nan
+    assert Mul(1/x, 1/x, evaluate=False).evalf(subs={x: 0})
path: sympy/core/tests/test_evalf.py
text: @user1:
There seems to be a missing `== zoo` here