Repository: sympy/sympy

Issue:
(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.


Discussion:
sympy-bot: 

Hi, I am the [SymPy bot](https://github.com/sympy/sympy-bot) (v161). I'm here to help you write a release notes entry. Please read the [guide on how to write release notes](https://github.com/sympy/sympy/wiki/Writing-Release-Notes).

:x: There was an issue with the release notes. **Please do not close this pull request;** instead edit the description after reading the [guide on how to write release notes](https://github.com/sympy/sympy/wiki/Writing-Release-Notes).

* No release notes were found. Please edit the PR description and write the release notes under `<!-- BEGIN RELEASE NOTES -->`.

<details><summary>Click here to see the pull request description that was parsed.</summary>

    <!-- 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


    #### Other comments


    #### 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 -->

    <!-- END RELEASE NOTES -->


</details><p>

github-actions[bot]: Benchmark results from GitHub Actions

Lower numbers are good, higher numbers are bad. A ratio less than 1
means a speed up and greater than 1 means a slowdown. Green lines
beginning with `+` are slowdowns (the PR is slower then master or
master is slower than the previous release). Red lines beginning
with `-` are speedups.

Significantly changed benchmark results (PR vs master)
```diff

```
Significantly changed benchmark results (master vs previous release)
```diff

```
Full benchmark results can be found as artifacts in GitHub Actions
(click on checks at the top of the PR).
eagleoflqj: > I don't want to add any type hints until mypy is running in CI

Then I’d wait #22212 to be merged first
oscarbenjamin: Can you split this PR into two PRs where one fixes the cited issues and the other adds type hints?

Patch:
diff --git a/sympy/codegen/ast.py b/sympy/codegen/ast.py
index a02714b174..80133758c7 100644
--- a/sympy/codegen/ast.py
+++ b/sympy/codegen/ast.py
@@ -27,36 +27,38 @@
        |
        |
        |--->Token
-                |--->Attribute
-                |--->For
-                |--->String
-                |       |--->QuotedString
-                |       |--->Comment
-                |--->Type
-                |       |--->IntBaseType
-                |       |              |--->_SizedIntType
-                |       |                               |--->SignedIntType
-                |       |                               |--->UnsignedIntType
-                |       |--->FloatBaseType
-                |                        |--->FloatType
-                |                        |--->ComplexBaseType
-                |                                           |--->ComplexType
-                |--->Node
-                |       |--->Variable
-                |       |           |---> Pointer
-                |       |--->FunctionPrototype
-                |                            |--->FunctionDefinition
-                |--->Element
-                |--->Declaration
-                |--->While
-                |--->Scope
-                |--->Stream
-                |--->Print
-                |--->FunctionCall
-                |--->BreakToken
-                |--->ContinueToken
-                |--->NoneToken
-                |--->Return
+       |        |--->Attribute
+       |        |--->For
+       |        |--->String
+       |        |       |--->QuotedString
+       |        |       |--->Comment
+       |        |--->Type
+       |        |       |--->IntBaseType
+       |        |       |              |--->_SizedIntType
+       |        |       |                               |--->SignedIntType
+       |        |       |                               |--->UnsignedIntType
+       |        |       |--->FloatBaseType
+       |        |                        |--->FloatType
+       |        |                        |--->ComplexBaseType
+       |        |                                           |--->ComplexType
+       |        |--->Node
+       |        |       |--->Variable
+       |        |       |           |---> Pointer
+       |        |       |--->FunctionPrototype
+       |        |                            |--->FunctionDefinition
+       |        |--->Element
+       |        |--->Declaration
+       |        |--->While
+       |        |--->Scope
+       |        |--->Stream
+       |        |--->Print
+       |        |--->FunctionCall
+       |        |--->BreakToken
+       |        |--->ContinueToken
+       |        |--->NoneToken
+       |        |--->Return
+       |
+       |--->Statement
 
 
 Predefined types
diff --git a/sympy/core/evalf.py b/sympy/core/evalf.py
index c9313200a4..c39e52f463 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 c9c9e8150c..fb010a5331 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/facts.py b/sympy/core/facts.py
index c333a97965..a3e1bb538b 100644
--- a/sympy/core/facts.py
+++ b/sympy/core/facts.py
@@ -7,7 +7,7 @@
 
 For rule-based inference engines, the classical work is RETE algorithm [1],
 [2] Although we are not implementing it in full (or even significantly)
-it's still worth a read to understand the underlying ideas.
+it's still still worth a read to understand the underlying ideas.
 
 In short, every rule in a system of rules is one of two forms:
 
diff --git a/sympy/functions/combinatorial/numbers.py b/sympy/functions/combinatorial/numbers.py
index d8e47a8ef6..a3f686271c 100644
--- a/sympy/functions/combinatorial/numbers.py
+++ b/sympy/functions/combinatorial/numbers.py
@@ -15,12 +15,11 @@
 from sympy.core.expr import Expr
 from sympy.core.function import Function, expand_mul
 from sympy.core.logic import fuzzy_not
-from sympy.core.mul import Mul, prod
-from sympy.core.numbers import E, pi, oo, Rational, Integer
-from sympy.core.relational import is_le, is_gt
+from sympy.core.mul import prod
+from sympy.core.numbers import E, pi, Rational, Integer
+from sympy.core.relational import LessThan, StrictGreaterThan
 from sympy.external.gmpy import SYMPY_INTS
-from sympy.functions.combinatorial.factorials import (binomial,
-    factorial, subfactorial)
+from sympy.functions.combinatorial.factorials import binomial, factorial
 from sympy.functions.elementary.exponential import log
 from sympy.functions.elementary.integers import floor
 from sympy.functions.elementary.miscellaneous import sqrt, cbrt
@@ -831,15 +830,16 @@ def eval(cls, n, m=None):
         if m.is_zero:
             return n
 
-        if n is S.Infinity:
+        if n is S.Infinity and m.is_Number:
+            # TODO: Fix for symbolic values of m
             if m.is_negative:
                 return S.NaN
-            elif is_le(m, S.One):
+            elif LessThan(m, S.One):
                 return S.Infinity
-            elif is_gt(m, S.One):
+            elif StrictGreaterThan(m, S.One):
                 return zeta(m)
             else:
-                return
+                return cls
 
         if n == 0:
             return S.Zero
@@ -2169,123 +2169,3 @@ def eval(cls, n):
         if n < 0:
             raise ValueError('The provided number must be a positive integer')
         return Integer(cls._motzkin(n - 1))
-
-
-def nD(i=None, brute=None, *, n=None, m=None):
-    """return the number of derangements for: ``n`` unique items, ``i``
-    items (as a sequence or multiset), or multiplicities, ``m`` given
-    as a sequence or multiset.
-
-    Examples
-    ========
-
-    >>> from sympy.utilities.iterables import generate_derangements as enum
-    >>> from sympy.functions.combinatorial.numbers import nD
-
-    A derangement ``d`` of sequence ``s`` has all ``d[i] != s[i]``:
-
-    >>> set([''.join(i) for i in enum('abc')])
-    {'bca', 'cab'}
-    >>> nD('abc')
-    2
-
-    Input as iterable or dictionary (multiset form) is accepted:
-
-    >>> assert nD([1, 2, 2, 3, 3, 3]) == nD({1: 1, 2: 2, 3: 3})
-
-    By default, a brute-force enumeration and count of multiset permutations
-    is only done if there are fewer than 9 elements. There may be cases when
-    there is high multiplicty with few unique elements that will benefit
-    from a brute-force enumeration, too. For this reason, the `brute`
-    keyword (default None) is provided. When False, the brute-force
-    enumeration will never be used. When True, it will always be used.
-
-    >>> nD('1111222233', brute=True)
-    44
-
-    For convenience, one may specify ``n`` distinct items using the
-    ``n`` keyword:
-
-    >>> assert nD(n=3) == nD('abc') == 2
-
-    Since the number of derangments depends on the multiplicity of the
-    elements and not the elements themselves, it may be more convenient
-    to give a list or multiset of multiplicities using keyword ``m``:
-
-    >>> assert nD('abc') == nD(m=(1,1,1)) == nD(m={1:3}) == 2
-
-    """
-    from sympy.utilities.iterables import multiset, multiset_derangements, iterable
-    from sympy.integrals.integrals import integrate
-    from sympy.functions.elementary.exponential import exp
-    from sympy.functions.special.polynomials import laguerre
-    from sympy.abc import x
-    def ok(x):
-        if not isinstance(x, SYMPY_INTS):
-            raise TypeError('expecting integer values')
-        if x < 0:
-            raise ValueError('value must not be negative')
-        return True
-
-    if (i, n, m).count(None) != 2:
-        raise ValueError('enter only 1 of i, n, or m')
-    if i is not None:
-        if isinstance(i, SYMPY_INTS):
-            raise TypeError('items must be a list or dictionary')
-        if not i:
-            return S.Zero
-        if type(i) is not dict:
-            s = list(i)
-            ms = multiset(s)
-        elif type(i) is dict:
-            all(ok(_) for _ in i.values())
-            ms = {k: v for k, v in i.items() if v}
-            s = None
-        if not ms:
-            return S.Zero
-        N = sum(ms.values())
-        counts = multiset(ms.values())
-        nkey = len(ms)
-    elif n is not None:
-        ok(n)
-        if not n:
-            return S.Zero
-        return subfactorial(n)
-    elif m is not None:
-        if isinstance(m, dict):
-            all(ok(i) and ok(j) for i, j in m.items())
-            counts = {k: v for k, v in m.items() if k*v}
-        elif iterable(m) or isinstance(m, str):
-            m = list(m)
-            all(ok(i) for i in m)
-            counts = multiset([i for i in m if i])
-        else:
-            raise TypeError('expecting iterable')
-        if not counts:
-            return S.Zero
-        N = sum(k*v for k, v in counts.items())
-        nkey = sum(counts.values())
-        s = None
-    big = int(max(counts))
-    if big == 1:  # no repetition
-        return subfactorial(nkey)
-    nval = len(counts)
-    if big*2 > N:
-        return S.Zero
-    if big*2 == N:
-        if nkey == 2 and nval == 1:
-            return S.One  # aaabbb
-        if nkey - 1 == big:  # one element repeated
-            return factorial(big)  # e.g. abc part of abcddd
-    if N < 9 and brute is None or brute:
-        # for all possibilities, this was found to be faster
-        if s is None:
-            s = []
-            i = 0
-            for m, v in counts.items():
-                for j in range(v):
-                    s.extend([i]*m)
-                    i += 1
-        return Integer(sum(1 for i in multiset_derangements(s)))
-    return Integer(abs(integrate(exp(-x)*Mul(*[
-        laguerre(i, x)**m for i, m in counts.items()]), (x, 0, oo))))
diff --git a/sympy/printing/codeprinter.py b/sympy/printing/codeprinter.py
index 8c2f803ca3..dff63a56b5 100644
--- a/sympy/printing/codeprinter.py
+++ b/sympy/printing/codeprinter.py
@@ -391,6 +391,10 @@ def _print_FunctionCall(self, expr):
     def _print_Variable(self, expr):
         return self._print(expr.symbol)
 
+    def _print_Statement(self, expr):
+        arg, = expr.args
+        return self._get_statement(self._print(arg))
+
     def _print_Symbol(self, expr):
 
         name = super()._print_Symbol(expr)
diff --git a/sympy/tensor/array/expressions/array_expressions.py b/sympy/tensor/array/expressions/array_expressions.py
index 51ae5c1af4..57e3ce84b5 100644
--- a/sympy/tensor/array/expressions/array_expressions.py
+++ b/sympy/tensor/array/expressions/array_expressions.py
@@ -899,11 +899,6 @@ def _normalize(self):
         if isinstance(expr, ArrayAdd):
             return self._ArrayContraction_denest_ArrayAdd(expr, *contraction_indices)
 
-        # Check single index contractions on 1-dimensional axes:
-        contraction_indices = [i for i in contraction_indices if len(i) > 1 or get_shape(expr)[i[0]] != 1]
-        if len(contraction_indices) == 0:
-            return expr
-
         return self.func(expr, *contraction_indices, normalize=False)
 
     def __mul__(self, other):
diff --git a/sympy/tensor/array/expressions/conv_array_to_matrix.py b/sympy/tensor/array/expressions/conv_array_to_matrix.py
index c4fde81c92..930e15bee1 100644
--- a/sympy/tensor/array/expressions/conv_array_to_matrix.py
+++ b/sympy/tensor/array/expressions/conv_array_to_matrix.py
@@ -157,11 +157,7 @@ def _(expr: ArrayContraction):
         return _array2matrix(expr)
     subexpr = expr.expr
     contraction_indices: tTuple[tTuple[int]] = expr.contraction_indices
-    if contraction_indices == ((0,), (1,)) or (
-        contraction_indices == ((0,),) and subexpr.shape[1] == 1
-    ) or (
-        contraction_indices == ((1,),) and subexpr.shape[0] == 1
-    ):
+    if contraction_indices == ((0,), (1,)):
         shape = subexpr.shape
         subexpr = _array2matrix(subexpr)
         if isinstance(subexpr, MatrixExpr):
diff --git a/sympy/tensor/array/expressions/conv_matrix_to_array.py b/sympy/tensor/array/expressions/conv_matrix_to_array.py
index aca690c30a..a49d12b62b 100644
--- a/sympy/tensor/array/expressions/conv_matrix_to_array.py
+++ b/sympy/tensor/array/expressions/conv_matrix_to_array.py
@@ -75,7 +75,6 @@ def convert_matrix_to_array(expr: Basic) -> Basic:
         if isinstance(exp, Integer) and exp > 0:
             return convert_matrix_to_array(HadamardProduct.fromiter(base for i in range(exp)))
         else:
-            d = Dummy("d")
-            return ArrayElementwiseApplyFunc(Lambda(d, d**exp), base)
+            raise NotImplementedError("conversion of Hadamard symbolic power is currently not supported")
     else:
         return expr
diff --git a/sympy/utilities/iterables.py b/sympy/utilities/iterables.py
index 111d3f31bd..4f557cf02a 100644
--- a/sympy/utilities/iterables.py
+++ b/sympy/utilities/iterables.py
@@ -2168,238 +2168,13 @@ def generate_involutions(n):
             yield p
 
 
-def multiset_derangements(s):
-    """Generate derangements of the elements of s *in place*.
-
-    Examples
-    ========
-
-    >>> from sympy.utilities.iterables import multiset_derangements, uniq
-
-    Because the derangements of multisets (not sets) are generated
-    in place, copies of the return value must be made if a collection
-    of derangements is desired or else all values will be the same:
-
-    >>> list(uniq([i for i in multiset_derangements('1233')]))
-    [['3', '3', '2', '1']]
-    >>> [i.copy() for i in multiset_derangements('1233')]
-    [['3', '3', '1', '2'], ['3', '3', '2', '1']]
-    """
-    ms = multiset(s)
-    mx = max(ms.values())
-    n = len(s)
-    # special cases
-
-    # 0) impossible case
-    if mx*2 > n:
-        return
-
-    # 1) singletons
-    if len(ms) == n:
-        for p in generate_derangements(s):
-            yield p
-        return
-
-    for M in ms:
-        if ms[M] == mx:
-            break
-    inonM = [i for i in range(n) if s[i] != M]
-    iM = [i for i in range(n) if s[i] == M]
-    rv = [None]*n
-
-    # 2) half are the same
-    if 2*mx == n:
-        for i in inonM:
-            rv[i] = M
-        for p in multiset_permutations([s[i] for i in inonM]):
-            for i, pi in zip(iM, p):
-                rv[i] = pi
-            yield rv
-        return
-
-    # 3) single repeat covers all but 1 of the non-repeats
-    if n - 2*mx == 1 and len(ms.values()) - 1 == n - mx:
-        for i in range(len(inonM)):
-            i1 = inonM[i]
-            ifill = inonM[:i] + inonM[i+1:]
-            for j in ifill:
-                rv[j] = M
-            rv[i1] = s[i1]
-            for p in permutations([s[j] for j in ifill]):
-                for j, pi in zip(iM, p):
-                    rv[j] = pi
-                for j in iM:
-                    rv[j], rv[i1] = rv[i1], rv[j]
-                    yield rv
-                    i1 = j
-        return
-
-    def finish_derangements():
-        """Place the last two elements into the partially completed
-        derangement, and yield the results.
-        In non-recursive version, this will be inlined, but a little
-        easier to understand as a function for now.
-        """
-
-        a = take[1][0]  # penultimate element
-        a_ct = take[1][1]
-        b = take[0][0]  # last element to be placed
-        b_ct = take[0][1]
-
-        # split the indexes of the not-already-assigned elemements of rv into
-        # three categories
-        forced_a = []  # positions which must have an a
-        forced_b = []  # positions which must have a b
-        open_free = []  # positions which could take either
-        for i in range(len(s)):
-            if rv[i] is None:
-                if s[i] == a:
-                    forced_b.append(i)
-                elif s[i] == b:
-                    forced_a.append(i)
-                else:
-                    open_free.append(i)
-
-        if len(forced_a) > a_ct or len(forced_b) > b_ct:
-            # No derangement possible
-            return
-        for i in forced_a:
-            rv[i] = a
-        for i in forced_b:
-            rv[i] = b
-        for a_place in subsets(open_free, a_ct - len(forced_a)):
-            for a_pos in a_place:
-                rv[a_pos] = a
-            for i in open_free:
-                if rv[i] is None:  # anything not in the subset is set to b
-                    rv[i] = b
-            yield rv
-            # Clean up/undo the final placements
-            for i in open_free:
-                rv[i] = None
-        # additional cleanup - clear forced_a, forced_b
-        for i in forced_a:
-            rv[i] = None
-        for i in forced_b:
-            rv[i] = None
-
-    def iopen(v):
-        return [i for i in range(n) if rv[i] is None and s[i] != v]
-
-    def do(j):
-        if j == -1:
-            yield rv
-        else:
-            M, mx = take[j]
-            for i in subsets(iopen(M), mx):
-                for ii in i:
-                    rv[ii] = M
-                yield from do(j - 1)
-                for ii in i:
-                    rv[ii] = None
-    take = sorted(ms.items(), key=lambda x:(x[1], x[0]))
-    yield from do(len(take) - 1)
-
-
-def random_derangement(t, choice=None, strict=True):
-    """Return a list of elements in which none are in the same positions
-    as they were originally. If an element fills more than half of the positions
-    then an error will be raised since no derangement is possible. To obtain
-    a derangement of as many items as possible--with some of the most numerous
-    remaining in their original positions--pass `strict=False`. To produce a
-    pseudorandom derangment, pass a pseudorandom selector like `Random(seed).choice`.
-
-    Examples
-    ========
-
-    >>> from sympy.utilities.iterables import random_derangement
-    >>> from random import Random
-    >>> t = 'SymPy: a CAS in pure Python'
-    >>> d = random_derangement(t)
-    >>> all(i != j for i, j in zip(d, t))
-    True
-
-    A predictable result can be obtained by using a pseudorandom
-    generator for the choice:
-
-    >>> c = Random(1).choice
-    >>> d = [''.join(random_derangement(t, c)) for i in range(5)]
-    >>> assert len(set(d)) != 1  # we got different values
-
-    By resetting c, the same sequence can be obtained:
-
-    >>> c = Random(1).choice
-    >>> d2 = [''.join(random_derangement(t, c)) for i in range(5)]
-    >>> assert d == d2
-    """
-    if choice is None:
-        import secrets
-        choice = secrets.choice
-    def shuffle(rv):
-        '''Knuth shuffle'''
-        for i in range(len(rv) - 1, 0, -1):
-            x = choice(rv[:i + 1])
-            j = rv.index(x)
-            rv[i], rv[j] = rv[j], rv[i]
-    def pick(rv, n):
-        '''shuffle rv and return the first n values
-        '''
-        shuffle(rv)
-        return rv[:n]
-    ms = multiset(t)
-    tot = len(t)
-    ms = sorted(ms.items(), key=lambda x: x[1])
-  # if there are not enough spaces for the most
-  # plentiful element to move to then some of them
-  # will have to stay in place
-    M, mx = ms[-1]
-    n = len(t)
-    xs = 2*mx - tot
-    if xs > 0:
-        if strict:
-            raise ValueError('no derangement possible')
-        opts = [i for (i, c) in enumerate(t) if c == ms[-1][0]]
-        pick(opts, xs)
-        stay = sorted(opts[:xs])
-        rv = list(t)
-        for i in reversed(stay):
-            rv.pop(i)
-        rv = random_derangement(rv, choice)
-        for i in stay:
-            rv.insert(i, ms[-1][0])
-        return ''.join(rv) if type(t) is str else rv
-  # the normal derangement calculated from here
-    if n == len(ms):
-      # approx 1/3 will succeed
-        rv = list(t)
-        while True:
-            shuffle(rv)
-            if all(i != j for i,j in zip(rv, t)):
-                break
-    else:
-      # general case
-        rv = [None]*n
-        while True:
-            j = 0
-            while j > -len(ms):  # do most numerous first
-                j -= 1
-                e, c = ms[j]
-                opts = [i for i in range(n) if rv[i] is None and t[i] != e]
-                if len(opts) < c:
-                    for i in range(n):
-                        rv[i] = None
-                    break # try again
-                pick(opts, c)
-                for i in range(c):
-                    rv[opts[i]] = e
-            else:
-                return rv
-    return rv
-
-
 def generate_derangements(perm):
     """
-    Routine to generate unique derangements or sets or multisets.
+    Routine to generate unique derangements.
+
+    TODO: This will be rewritten to use the
+    ECO operator approach once the permutations
+    branch is in master.
 
     Examples
     ========
@@ -2420,21 +2195,9 @@ def generate_derangements(perm):
     sympy.functions.combinatorial.factorials.subfactorial
 
     """
-    if not has_dups(perm):
-        s = perm
-        if len(perm) == 2:
-            yield [s[1],s[0]]
-            return
-        if len(perm) == 3:
-            yield [s[1],s[2],s[0]]
-            yield [s[2],s[0],s[1]]
-            return
-        for p in permutations(s):
-            if not any(i == j for i, j in zip(p, s)):
-                yield list(p)
-    else:
-        for p in multiset_derangements(perm):
-            yield list(p)
+    for p in multiset_permutations(perm):
+        if not any(i == j for i, j in zip(perm, p)):
+            yield p
 
 
 def necklaces(n, k, free=False):


Test patch:
diff --git a/sympy/core/tests/test_evalf.py b/sympy/core/tests/test_evalf.py
index cd099e7d05..f2c429c3ea 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}) == zoo
+    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 6982c827e7..19458f30c7 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/combinatorial/tests/test_comb_numbers.py b/sympy/functions/combinatorial/tests/test_comb_numbers.py
index 607bb0362a..54f7cd5275 100644
--- a/sympy/functions/combinatorial/tests/test_comb_numbers.py
+++ b/sympy/functions/combinatorial/tests/test_comb_numbers.py
@@ -22,7 +22,9 @@
 from sympy.core.numbers import GoldenRatio, Integer
 
 from sympy.testing.pytest import XFAIL, raises, nocache_fail
-from sympy.abc import x
+
+
+x = Symbol('x')
 
 
 def test_carmichael():
@@ -216,13 +218,6 @@ def test_harmonic():
     assert harmonic(oo, 1) is oo
     assert harmonic(oo, 2) == (pi**2)/6
     assert harmonic(oo, 3) == zeta(3)
-    assert harmonic(oo, Dummy(negative=True)) is S.NaN
-    ip = Dummy(integer=True, positive=True)
-    if (1/ip <= 1) is True:  #---------------------------------+
-        assert None, 'delete this if-block and the next line' #|
-    ip = Dummy(even=True, positive=True)  #--------------------+
-    assert harmonic(oo, 1/ip) is oo
-    assert harmonic(oo, 1 + ip) is zeta(1 + ip)
 
     assert harmonic(0, m) == 0
 
@@ -690,55 +685,3 @@ def test_motzkin():
     raises(ValueError, lambda: motzkin.find_motzkin_numbers_in_range(-2,7))
     raises(ValueError, lambda: motzkin.find_motzkin_numbers_in_range(13,7))
     raises(ValueError, lambda: motzkin.find_first_n_motzkins(112.8))
-
-
-def test_nD_derangements():
-    from sympy.utilities.iterables import (partitions, multiset,
-        multiset_derangements, multiset_permutations)
-    from sympy.functions.combinatorial.numbers import nD
-
-    got = []
-    for i in partitions(8, k=4):
-        s = []
-        it = 0
-        for k, v in i.items():
-            for i in range(v):
-                s.extend([it]*k)
-                it += 1
-        ms = multiset(s)
-        c1 = sum(1 for i in multiset_permutations(s) if
-            all(i != j for i, j in zip(i, s)))
-        assert c1 == nD(ms) == nD(ms, 0) == nD(ms, 1)
-        v = [tuple(i) for i in multiset_derangements(s)]
-        c2 = len(v)
-        assert c2 == len(set(v))
-        assert c1 == c2
-        got.append(c1)
-    assert got == [1, 4, 6, 12, 24, 24, 61, 126, 315, 780, 297, 772,
-        2033, 5430, 14833]
-
-    assert nD('1112233456', brute=True) == nD('1112233456') == 16356
-    assert nD('') == nD([]) == nD({}) == 0
-    assert nD({1: 0}) == 0
-    raises(ValueError, lambda: nD({1: -1}))
-    assert nD('112') == 0
-    assert nD(i='112') == 0
-    assert [nD(n=i) for i in range(6)] == [0, 0, 1, 2, 9, 44]
-    assert nD((i for i in range(4))) == nD('0123') == 9
-    assert nD(m=(i for i in range(4))) == 3
-    assert nD(m={0: 1, 1: 1, 2: 1, 3: 1}) == 3
-    assert nD(m=[0, 1, 2, 3]) == 3
-    raises(TypeError, lambda: nD(m=0))
-    raises(TypeError, lambda: nD(-1))
-    assert nD({-1: 1, -2: 1}) == 1
-    assert nD(m={0: 3}) == 0
-    raises(ValueError, lambda: nD(i='123', n=3))
-    raises(ValueError, lambda: nD(i='123', m=(1,2)))
-    raises(ValueError, lambda: nD(n=0, m=(1,2)))
-    raises(ValueError, lambda: nD({1: -1}))
-    raises(ValueError, lambda: nD(m={-1: 1, 2: 1}))
-    raises(ValueError, lambda: nD(m={1: -1, 2: 1}))
-    raises(ValueError, lambda: nD(m=[-1, 2]))
-    raises(TypeError, lambda: nD({1: x}))
-    raises(TypeError, lambda: nD(m={1: x}))
-    raises(TypeError, lambda: nD(m={x: 1}))
diff --git a/sympy/series/tests/test_limits.py b/sympy/series/tests/test_limits.py
index 996e5d1fbc..22201e479d 100644
--- a/sympy/series/tests/test_limits.py
+++ b/sympy/series/tests/test_limits.py
@@ -607,16 +607,13 @@ def test_issue_8481():
     limit(lamda**k * exp(-lamda) / factorial(k), k, oo) == 0
 
 
-def test_issue_8635_18176():
+def test_issue_8635():
     x = Symbol('x', real=True)
     k = Symbol('k', positive=True)
     assert limit(x**n - x**(n - 0), x, oo) == 0
     assert limit(x**n - x**(n - 5), x, oo) == oo
     assert limit(x**n - x**(n - 2.5), x, oo) == oo
     assert limit(x**n - x**(n - k - 1), x, oo) == oo
-    x = Symbol('x', positive=True)
-    assert limit(x**n - x**(n - 1), x, oo) == oo
-    assert limit(x**n - x**(n + 2), x, oo) == -oo
 
 
 def test_issue_8730():
diff --git a/sympy/tensor/array/expressions/tests/test_array_expressions.py b/sympy/tensor/array/expressions/tests/test_array_expressions.py
index 2ba5241509..a5d31e8cd7 100644
--- a/sympy/tensor/array/expressions/tests/test_array_expressions.py
+++ b/sympy/tensor/array/expressions/tests/test_array_expressions.py
@@ -112,12 +112,6 @@ def test_arrayexpr_contraction_construction():
     assert indtup == [[(0, 0), (1, 1)], [(0, 1), (2, 0)]]
     assert cg._contraction_tuples_to_contraction_indices(cg.expr, indtup) == [(0, 3), (1, 4)]
 
-    # Test removal of trivial contraction:
-    assert ArrayContraction(a, (1,)) == a
-    assert ArrayContraction(
-        ArrayTensorProduct(a, b), (0, 2), (1,), (3,)) == ArrayContraction(
-        ArrayTensorProduct(a, b), (0, 2))
-
 
 def test_arrayexpr_array_flatten():
 
diff --git a/sympy/tensor/array/expressions/tests/test_convert_matrix_to_array.py b/sympy/tensor/array/expressions/tests/test_convert_matrix_to_array.py
index 6042a59531..deb0ec50a2 100644
--- a/sympy/tensor/array/expressions/tests/test_convert_matrix_to_array.py
+++ b/sympy/tensor/array/expressions/tests/test_convert_matrix_to_array.py
@@ -1,5 +1,4 @@
-from sympy import Lambda
-from sympy.core.symbol import symbols, Dummy
+from sympy.core.symbol import symbols
 from sympy.matrices.expressions.hadamard import (HadamardPower, HadamardProduct)
 from sympy.matrices.expressions.inverse import Inverse
 from sympy.matrices.expressions.matexpr import MatrixSymbol
@@ -8,7 +7,7 @@
 from sympy.matrices.expressions.trace import Trace
 from sympy.matrices.expressions.transpose import Transpose
 from sympy.tensor.array.expressions.array_expressions import ArrayTensorProduct, ArrayContraction, \
-    PermuteDims, ArrayDiagonal, ArrayElementwiseApplyFunc
+    PermuteDims, ArrayDiagonal
 from sympy.tensor.array.expressions.conv_array_to_matrix import convert_array_to_matrix
 from sympy.tensor.array.expressions.conv_matrix_to_array import convert_matrix_to_array
 
@@ -106,11 +105,6 @@ def test_arrayexpr_convert_matrix_to_array():
     result = ArrayDiagonal(ArrayContraction(ArrayTensorProduct(M, N, M, N), (1, 2), (5, 6)), (0, 2), (1, 3))
     assert convert_matrix_to_array(expr) == result
 
-    expr = HadamardPower(M, n)
-    d0 = Dummy("d0")
-    result = ArrayElementwiseApplyFunc(Lambda(d0, d0**n), M)
-    assert convert_matrix_to_array(expr).dummy_eq(result)
-
     expr = M**2
     assert isinstance(expr, MatPow)
     assert convert_matrix_to_array(expr) == ArrayContraction(ArrayTensorProduct(M, M), (1, 2))
diff --git a/sympy/utilities/tests/test_iterables.py b/sympy/utilities/tests/test_iterables.py
index c2c22e6295..d74d271815 100644
--- a/sympy/utilities/tests/test_iterables.py
+++ b/sympy/utilities/tests/test_iterables.py
@@ -19,7 +19,7 @@
     prefixes, reshape, rotate_left, rotate_right, runs, sift,
     strongly_connected_components, subsets, take, topological_sort, unflatten,
     uniq, variations, ordered_partitions, rotations, is_palindromic, iterable,
-    NotIterable, multiset_derangements)
+    NotIterable)
 from sympy.utilities.enumerative import (
     factoring_visitor, multiset_partitions_taocp )
 
@@ -547,19 +547,6 @@ def test_derangements():
     assert list(generate_derangements([0, 1, 2, 2])) == [
         [2, 2, 0, 1], [2, 2, 1, 0]]
     assert list(generate_derangements('ba')) == [list('ab')]
-    # multiset_derangements
-    D = multiset_derangements
-    assert list(D('abb')) == []
-    assert [''.join(i) for i in D('ab')] == ['ba']
-    assert [''.join(i) for i in D('abc')] == ['bca', 'cab']
-    assert [''.join(i) for i in D('aabb')] == ['bbaa']
-    assert [''.join(i) for i in D('aabbcccc')] == [
-        'ccccaabb', 'ccccabab', 'ccccabba', 'ccccbaab', 'ccccbaba',
-        'ccccbbaa']
-    assert [''.join(i) for i in D('aabbccc')] == [
-        'cccabba', 'cccabab', 'cccaabb', 'ccacbba', 'ccacbab',
-        'ccacabb', 'cbccbaa', 'cbccaba', 'cbccaab', 'bcccbaa',
-        'bcccaba', 'bcccaab']
 
 
 def test_necklaces():