Repository: sympy/sympy

Issue:
Matching failure for simple linear expression
This match works:
```julia
In [20]: x, y, z = symbols('x y z')                                                                                                            

In [21]: a = Wild('a')                                                                                                                         

In [22]: b = Wild('b')                                                                                                                         

In [23]: (x + (y-1)*z).match(a*x + b*z)                                                                                                        
Out[23]: {a: 1, b: y - 1}
```
However removing the `a` Wild from the pattern means it longer matches:
```julia
In [24]: (x + (y-1)*z).match(x + b*z) 
```
`dummy_eq` failing for equal Sum
`dummy_eq` should treat a Dummy variable in `self` as a Wild (actually or idiomatically) to see if two expressions are equivalent or not. The following test fails:
```python
from sympy import Sum, Dummy
from sympy.abc import x, a c
assert not Sum(x, (x, a, Dummy())).dummy_eq(Sum(x, (x, a, c)))
```
If `Dummy -> c` then the two are equivalent (and that is the purpose and intended output of the method)), so the current test is wrong. This is corrected in #19261.

And, of course, if that is fixed then the following is wrong (and fixed by #19261, too) since `x` and `y` are structurally bound symbols whose identity does not matter mathematically:
```python
from sympy import Sum, Dummy
from sympy.abc import x, a c
assert not Sum(x, (x, a, Dummy())).dummy_eq(Sum(y, (y, a, c)))
```


Discussion:
sympy-bot: :white_check_mark:

Hi, I am the [SymPy bot](https://github.com/sympy/sympy-bot) (v158). 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).



Your release notes are in good order.

Here is what the release notes will look like:
* core
  * `Sum` is now included in `count_ops` results ([#19261](https://github.com/sympy/sympy/pull/19261) by [@smichr](https://github.com/smichr))

  * `match` has been optimized to improve term-matching in Add ([#19261](https://github.com/sympy/sympy/pull/19261) by [@smichr](https://github.com/smichr))

  * `dummy_eq` can treat more than one symbol as a wild-type ([#19261](https://github.com/sympy/sympy/pull/19261) by [@smichr](https://github.com/smichr))

This will be added to https://github.com/sympy/sympy/wiki/Release-Notes-for-1.7.

Note: This comment will be updated with the latest check if you edit the pull request. You need to reload the page to see it. <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. -->


    #### Brief description of what is fixed or changed

    * count_ops

    `Sum`, unlike `Integral` did not show up in `count_ops` results. After reviewing the code I tidied it up a bit and included `Sum`.

    before:
    ```python
    >>> Sum(x, (x, 1, 2)).count_ops()
    0
    ```
    now:
    ```python
    >>> Sum(x, (x, 1, 2)).count_ops(visual=1)
    SUM
    ```

    * matching

    The use of `ordered` while matching sometimes orders terms of Add in a way that leads to sub-optimal matching:

    before:
    ```python
    >>> X, Y = map(Wild, "XY")
    >>> x, y = symbols('x y')
    >>> (5*y - x).match(5*X - Y)
    {Y_: -5*y, X_: -x/5}
    ```
    now:
    ```python
    >>> (5*y - x).match(5*X - Y)
    {X: y, Y: x}
    ```

    * dummy_eq

    `dummy_eq` was very limited in its ability to check expressions with Dummy symbols. In anticipation of #18884 which generates expressions with with more than one Dummy in the result, `dummy_eq` has been enhanced to allow equality checking when there is more than one Dummy present.

    ```python
    >>> (Dummy() + 2*Dummy()).dummy_eq(x + 2*y)
    True
    ```

    Specifying `symbol` -- now permitting more than one symbol -- is only needed when wanting to indicate which symbols in `other` are to be treated as having a structurally equivalent symbol in `self`:
    ```python
    >>> (b+2*a).dummy_eq(x + 2*y)  # no Dummies on the left (default behavior)
    False
    >>> (b+2*a).dummy_eq(x + 2*y, x, y)  # see if x and y can match to a and b
    True
    ```
    #### Other comments


    #### Release Notes

    <!-- Write the release notes for this release below. 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
        * `Sum` is now included in `count_ops` results
        * `match` has been optimized to improve term-matching in Add
        * `dummy_eq` can treat more than one symbol as a wild-type
    <!-- END RELEASE NOTES -->

</details><p>

codecov[bot]: # [Codecov](https://codecov.io/gh/sympy/sympy/pull/19261?src=pr&el=h1) Report
> Merging [#19261](https://codecov.io/gh/sympy/sympy/pull/19261?src=pr&el=desc) into [master](https://codecov.io/gh/sympy/sympy/commit/b561fee8a4810924f9481a5d78a83436b544a672&el=desc) will **increase** coverage by `0.027%`.
> The diff coverage is `100.000%`.

```diff
@@              Coverage Diff              @@
##            master    #19261       +/-   ##
=============================================
+ Coverage   75.649%   75.676%   +0.027%     
=============================================
  Files          651       651               
  Lines       169487    169499       +12     
  Branches     40021     40029        +8     
=============================================
+ Hits        128216    128271       +55     
+ Misses       35654     35614       -40     
+ Partials      5617      5614        -3     
```

codecov[bot]: # [Codecov](https://codecov.io/gh/sympy/sympy/pull/19261?src=pr&el=h1) Report
> Merging [#19261](https://codecov.io/gh/sympy/sympy/pull/19261?src=pr&el=desc) into [master](https://codecov.io/gh/sympy/sympy/commit/b561fee8a4810924f9481a5d78a83436b544a672&el=desc) will **increase** coverage by `0.027%`.
> The diff coverage is `100.000%`.

```diff
@@              Coverage Diff              @@
##            master    #19261       +/-   ##
=============================================
+ Coverage   75.649%   75.676%   +0.027%     
=============================================
  Files          651       651               
  Lines       169487    169499       +12     
  Branches     40021     40029        +8     
=============================================
+ Hits        128216    128271       +55     
+ Misses       35654     35614       -40     
+ Partials      5617      5614        -3     
```

oscarbenjamin: I still don't understand why `dummy_eq` needs to be changed at all...
smichr: > why `dummy_eq` needs to be changed at all.

* It can't test more than 1 dummy symbol -- tests in your domain solver generate 2 Dummies (`trigintegrate`). How will those be tested?
oscarbenjamin: > * It can't test more than 1 dummy symbol

Okay I see that now.

Maybe that part of `dummy_eq` is just a bad idea anyway though. Really if we want the Dummies to compare equal is that not because they are in fact bound?
oscarbenjamin: Are there any *good* examples using `dummy_eq` with a non-bound `Dummy`?
oscarbenjamin: I tried this diff, ran the tests, and saw no particularly compelling failures:
```diff
diff --git a/sympy/core/basic.py b/sympy/core/basic.py
index bbb285a72f..8dc8aa9fde 100644
--- a/sympy/core/basic.py
+++ b/sympy/core/basic.py
@@ -388,28 +388,7 @@ def dummy_eq(self, other, symbol=None):
         False
 
         """
-        s = self.as_dummy()
-        o = _sympify(other)
-        o = o.as_dummy()
-
-        dummy_symbols = [i for i in s.free_symbols if i.is_Dummy]
-
-        if len(dummy_symbols) == 1:
-            dummy = dummy_symbols.pop()
-        else:
-            return s == o
-
-        if symbol is None:
-            symbols = o.free_symbols
-
-            if len(symbols) == 1:
-                symbol = symbols.pop()
-            else:
-                return s == o
-
-        tmp = dummy.__class__()
-
-        return s.subs(dummy, tmp) == o.subs(symbol, tmp)
+        return self.as_dummy() == _sympify(other).as_dummy()
 
     # Note, we always use the default ordering (lex) in __str__ and __repr__,
     # regardless of the global setting.  See issue 5487.
```
smichr: > Are there any _good_ examples using `dummy_eq` with a non-bound `Dummy`?

Searching through codebase:
* comparing Lambdas - obsolete since Lambda(x,x+1)==Lambda(y,y+1) now
* Polys that need to create a generator, e.g. `gegenbauer_poly(1, S.Half).dummy_eq(x)`

All the rest appear (looking quickly) to be bound-symbol related.

What about the 2 constants of integration that were being returned by your branch? Should we always be looking for a way to send back expressions with uniquely named symbols or should just use Dummy symbols?
oscarbenjamin: > * comparing Lambdas - obsolete since `Lambda(x,x+1)==Lambda(y,y+1)` now

That should be changed. Either `Lambda` should canonicalise the symbols or those should compare unequal.

> * Polys that need to create a generator, e.g. `gegenbauer_poly(1, S.Half).dummy_eq(x)`

Maybe that API should be fixed to require a generator or return a `PurePoly`.

> What about the 2 constants of integration that were being returned by your branch?

I'm not sure which branch you are referring to.
smichr: > I'm not sure which branch you are referring to.

#15396, e.g.:
```python
>>> from sympy.integrals.trigonometry import *
>>> trigintegrate(cos(x)*cot(x), x)
-(_b1 - _b2)*log(cos(x) + 1)/2 + (_b1 + _b2)*log(cos(x) + (_b2 - (_b1 + _b2))/_b1)/2 + cos(x) + 1
```


> That should be changed.

#17498

Lots of discussion at #15240 and your conclusions stated well [here](https://github.com/sympy/sympy/issues/5539#issuecomment-524504574).
smichr: This is a bit of a mystery: how can a doctest pass and the identical test in the test suite fail?
```
Traceback (most recent call last):
  File "/home/travis/virtualenv/python3.8.0/lib/python3.8/site-packages/sympy-1.7.dev0-py3.8.egg/sympy/testing/runtests.py", line 1335, in _timeout
    function()
  File "/home/travis/virtualenv/python3.8.0/lib/python3.8/site-packages/sympy-1.7.dev0-py3.8.egg/sympy/physics/quantum/tests/test_cg.py", line 150, in test_cg_simp_sum
    assert cg_simp(x * Sum(CG(a, alpha, b, 0, a, alpha), (alpha, -a, a)
AssertionError
```
vs this passing doctest 
```
    """
    >>> from sympy import symbols, Sum, KroneckerDelta
    >>> from sympy.physics.quantum.cg import CG, cg_simp
    >>> x, a, b, alpha = symbols('x a b alpha')
    >>> k = cg_simp(x * Sum(CG(a, alpha, b, 0, a, alpha), (alpha, -a, a)))
    >>> x*(2*a + 1)*KroneckerDelta(0, b)
    x*(2*a + 1)*KroneckerDelta(0, b)
    >>> _ == k
    True
```
Anyone have any ideas?
oscarbenjamin: > Anyone have any ideas? I reverted the `count_ops` change and the problem still remains, so it has to do with the ordering of args during matching.

No idea. Ideally the matching should be done so that the ordering of the args doesn't affect the result...
smichr: It is system/version dependent because it passes locally under 3.7 with x86 Windows.
smichr: Wow -- would not have expected this and haven't tracked down how it is happening, but an extra symbol is finding its way into the match result.
```
pat, expr = Sum(CG(a_, _0, b_, 0, a_, _0), (_0, -a_, a_)), Sum(CG(a, _0, b, 0, a, _0), (_0, -a, a))
pat.match(expr) -> {V_: x, a_: a, b_: b}
```
Where is that `V_` coming from?? It doesn't happen under 3.7 (Windows) but is happening under 3.8 during travis tests. I am now sanitizing the match dictionary to correct for this. I imagine that it gets in when the dictionary is treated like a defaultdict and a query like `d[V_]` is made:
```
>>> d
defaultdict(<class 'int'>, {})
>>> d[x]
0
>>> d
defaultdict(<class 'int'>, {x: 0})
```
But how it ends up giving an `x` for the value...must be a special defaultdict somewhere...
smichr: @oscarbenjamin , every issue has been resolved here. The injection of the "V_" into the match results -- a system/Python-dependent incident -- has been worked around. I don't know the root cause of that, however, as I mention [above](https://github.com/sympy/sympy/pull/19261#issuecomment-630589553).
smichr: @iammosespaulr , if you could take a look through these many small changes that would be helpful.
smichr: @iammosespaulr , each commit is self contained if that helps you work through this.
oscarbenjamin: All looks good except not sure about the last commit. If I go back one commit from the tip of this branch how would I potentially reproduce the problem?
smichr: > how would I potentially reproduce the problem?

push the commit to travis - the first assertion in `test_cg_simp_sum` will fail.  The helper for the simplifcation routine checks the length of the match dictionary and, because a spurious symbol formerly came back in the dictionary, the assumption that only the requested symbols would be present is violated. It failed every time on travis after the `302bdfd` commit ("match ignores bound symbols unless requested"); it never fails locally under Windows with Python 3.7.
smichr: @oscarbenjamin , I moved the test that will fail to `basic` and put an assertion in the match function. It should fail on my `debug` branch travis [run](https://travis-ci.org/github/smichr/sympy/builds/688918208?utm_source=github_status&utm_medium=notification).
oscarbenjamin: Looks like lots of things are failing on the debug branch
smichr: > debug branch

https://travis-ci.org/github/smichr/sympy/builds/688918208
iammosespaulr: > @iammosespaulr , each commit is self contained if that helps you work through this.

It's kinda a big pr, I'm starting from the issues fixed 😄, also checking out your `debug` branch.
smichr: > Looks like lots of things are failing on the debug branch

The assertion that there be no extraneous symbols in the match is what is causing others to fail. But the code using those results wasn't influenced by the presence of extra symbols.

VERY CURIOUS that the same test, moved to basic, did not fail but failed in the cg context. I somewhat think this is a Python issue and not a SymPy issue (directly).
smichr: > also checking out your `debug` branch.

My humble recommendation is to not worry about that branch. I'm not sure that it's really a SymPy issue. Going through the commits as they appear here in order -- each is pretty small -- would be a good way to go. But your choice. :-)
oscarbenjamin: > VERY CURIOUS that the same test, moved to basic, did not fail but failed in the cg context.

I think it's because it depends on one test having been run before another.
smichr: > I think it's because it depends on one test having been run before another.

Is there some issue that we are abusing by providing a default argument to the match routine:

`def matches(...repl_dict={}...)` Should that be None and then internally provide a dictionary if it is None?
oscarbenjamin: > Is there some issue that we are abusing by providing a default argument to the match routine:
> 
> `def matches(...repl_dict={}...)` Should that be None and then internally provide a dictionary if it is None?

Oh, yes. That could be the cause of this. It should be fixed either way. Don't use mutable defaults!

The dict is directly returned if `self==expr`:
https://github.com/sympy/sympy/blob/0180df08428d4c9ccff124fd2984720e8409f46f/sympy/core/basic.py#L1556-L1589
That allows the caller to modify it which will affect subsequent calls to matches.
smichr: > Don't use mutable defaults!

this should be a TODO -- there are 70 cases of `={}` in the codebase and 71 of `=[]`.
oscarbenjamin: Here's a demo:
```julia
In [1]: exp(x).matches(exp(y))                                                                                                    

In [2]: d = exp(x).matches(exp(x))                                                                                                

In [3]: d[x] = y                                                                                                                  

In [4]: exp(x).matches(exp(y))                                                                                                    
Out[4]: {x: y}
```
smichr: So should None be used or should we just make sure that each routine that uses it starts with `repl_dict = repl_dict.copy()`?
smichr: @oscarbenjamin , assuming you don't want this in 1.6, I will look to merge as soon as 1.6 is released.
oscarbenjamin: 1.6 is already on a separate branch so this can go into master.

I think we should always use None to avoid mutable defaults so that there is no confusion. Then the rule is simple: no mutable defaults.
smichr: @iammosespaulr, Oscar has already looked this over and the mystery issue has been solved. I'll wait to see if you have any comments and, if not, will commit in about 24 hours.
oscarbenjamin: Yes, this looks good. I think it's better to use None as the default value rather than calling copy but that can be fixed in a more general cleanup for #19372 
iammosespaulr: @smichr saw that a list was used and then converted back to dict, it was in `core/basic` so I thought one less step makes it faster right ?
iammosespaulr: ( also you've made extensive use of `dumeq`, does it enhance `dummy_eq`? Would it be useful to have it's functionality inside `dummy_eq` so that you don't have to declare a seperate func ? ) just a thought
smichr: dumeq just handles the iterable case. Nonbasic list and tuple can't use `dummy_eq` method.
oscarbenjamin: > dumeq just handles the iterable case. Nonbasic list and tuple can't use `dummy_eq` method.

This is one reason why I think it would be better if we used functions more rather than methods e.g. we could have:
```
>>> subs([x+y, x-y], x, 2)
[2+y, 2-y]
```
Same for evalf etc

Patch:
diff --git a/sympy/core/basic.py b/sympy/core/basic.py
index bbb285a72f..53beba2396 100644
--- a/sympy/core/basic.py
+++ b/sympy/core/basic.py
@@ -409,7 +409,7 @@ def dummy_eq(self, other, symbol=None):
 
         tmp = dummy.__class__()
 
-        return s.subs(dummy, tmp) == o.subs(symbol, tmp)
+        return s.xreplace({dummy: tmp}) == o.xreplace({symbol: tmp})
 
     # Note, we always use the default ordering (lex) in __str__ and __repr__,
     # regardless of the global setting.  See issue 5487.
@@ -546,7 +546,9 @@ def as_dummy(self):
         """Return the expression with any objects having structurally
         bound symbols replaced with unique, canonical symbols within
         the object in which they appear and having only the default
-        assumption for commutativity being True.
+        assumption for commutativity being True. When applied to a
+        symbol a new symbol having only the same commutativity will be
+        returned.
 
         Examples
         ========
@@ -558,41 +560,40 @@ def as_dummy(self):
         Integral(_0, (_0, x))
         >>> _.variables[0].is_real is None
         True
+        >>> r.as_dummy()
+        _r
 
         Notes
         =====
 
-        Any object that has structural dummy variables should have
-        a property, `bound_symbols` that returns a list of structural
-        dummy symbols of the object itself.
+        Any object that has structurally bound variables should have
+        a property, `bound_symbols` that returns those symbols
+        appearing in the object.
 
-        Lambda and Subs have bound symbols, but because of how they
-        are cached, they already compare the same regardless of their
-        bound symbols:
-
-        >>> from sympy import Lambda
-        >>> Lambda(x, x + 1) == Lambda(y, y + 1)
-        True
         """
+        from sympy.core.symbol import Dummy, Symbol
         def can(x):
-            d = {i: i.as_dummy() for i in x.bound_symbols}
             # mask free that shadow bound
+            free = x.free_symbols
+            bound = set(x.bound_symbols)
+            d = {i: Dummy() for i in bound & free}
             x = x.subs(d)
-            c = x.canonical_variables
-            # replace bound
-            x = x.xreplace(c)
-            # undo masking
-            x = x.xreplace({v: k for k, v in d.items()})
-            return x
+            # replace bound with canonical names
+            x = x.xreplace(x.canonical_variables)
+            # return after undoing masking
+            return x.xreplace({v: k for k, v in d.items()})
+        if not self.has(Symbol):
+            return self
         return self.replace(
             lambda x: hasattr(x, 'bound_symbols'),
-            lambda x: can(x))
+            lambda x: can(x),
+            simultaneous=False)
 
     @property
     def canonical_variables(self):
         """Return a dictionary mapping any variable defined in
         ``self.bound_symbols`` to Symbols that do not clash
-        with any existing symbol in the expression.
+        with any free symbols in the expression.
 
         Examples
         ========
@@ -602,22 +603,21 @@ def canonical_variables(self):
         >>> Lambda(x, 2*x).canonical_variables
         {x: _0}
         """
-        from sympy.core.symbol import Symbol
         from sympy.utilities.iterables import numbered_symbols
         if not hasattr(self, 'bound_symbols'):
             return {}
         dums = numbered_symbols('_')
         reps = {}
-        v = self.bound_symbols
-        # this free will include bound symbols that are not part of
-        # self's bound symbols
-        free = {i.name for i in self.atoms(Symbol) - set(v)}
-        for v in v:
+        # watch out for free symbol that are not in bound symbols;
+        # those that are in bound symbols are about to get changed
+        bound = self.bound_symbols
+        names = {i.name for i in self.free_symbols - set(bound)}
+        for b in bound:
             d = next(dums)
-            if v.is_Symbol:
-                while v.name == d.name or d.name in free:
+            if b.is_Symbol:
+                while d.name in names:
                     d = next(dums)
-            reps[v] = d
+            reps[b] = d
         return reps
 
     def rcall(self, *args):
@@ -900,9 +900,10 @@ def subs(self, *args, **kwargs):
         sympy.core.evalf.EvalfMixin.evalf: calculates the given formula to a desired level of precision
 
         """
+        from sympy.core.compatibility import _nodes, default_sort_key
         from sympy.core.containers import Dict
-        from sympy.utilities.iterables import sift
-        from sympy import Dummy, Symbol
+        from sympy.core.symbol import Dummy, Symbol
+        from sympy.utilities.misc import filldedent
 
         unordered = False
         if len(args) == 1:
@@ -913,7 +914,6 @@ def subs(self, *args, **kwargs):
                 unordered = True
                 sequence = sequence.items()
             elif not iterable(sequence):
-                from sympy.utilities.misc import filldedent
                 raise ValueError(filldedent("""
                    When a single argument is passed to subs
                    it should be a dictionary of old: new pairs or an iterable
@@ -941,10 +941,17 @@ def subs(self, *args, **kwargs):
 
         if unordered:
             sequence = dict(sequence)
-            atoms, nonatoms = sift(list(sequence),
-                lambda x: x.is_Atom, binary=True)
-            sequence = [(k, sequence[k]) for k in
-                list(reversed(list(ordered(nonatoms)))) + list(ordered(atoms))]
+            # order so more complex items are first and items
+            # of identical complexity are ordered so
+            # f(x) < f(y) < x < y
+            # \___ 2 __/    \_1_/  <- number of nodes
+            #
+            # For more complex ordering use an unordered sequence.
+            k = list(ordered(sequence, default=False, keys=(
+                lambda x: -_nodes(x),
+                lambda x: default_sort_key(x),
+                )))
+            sequence = [(k, sequence[k]) for k in k]
 
         if kwargs.pop('simultaneous', False):  # XXX should this be the default for dict subs?
             reps = {}
@@ -1391,11 +1398,11 @@ def replace(self, query, value, map=False, simultaneous=True, exact=None):
 
         >>> e = x**(1 + y)
         >>> (x**(1 + y)).replace(x**(1 + a), lambda a: x**-a, exact=False)
-        1
+        x
         >>> (x**(1 + y)).replace(x**(1 + a), lambda a: x**-a, exact=True)
         x**(-x - y + 1)
         >>> (x**y).replace(x**(1 + a), lambda a: x**-a, exact=False)
-        1
+        x
         >>> (x**y).replace(x**(1 + a), lambda a: x**-a, exact=True)
         x**(1 - y)
 
@@ -1417,8 +1424,8 @@ def replace(self, query, value, map=False, simultaneous=True, exact=None):
                   using matching rules
 
         """
-        from sympy.core.symbol import Dummy, Wild
-        from sympy.simplify.simplify import bottom_up
+        from sympy.core.symbol import Wild
+
 
         try:
             query = _sympify(query)
@@ -1480,55 +1487,40 @@ def replace(self, query, value, map=False, simultaneous=True, exact=None):
                 "first argument to replace() must be a "
                 "type, an expression or a callable")
 
+        def walk(rv, F):
+            """Apply ``F`` to args and then to result.
+            """
+            args = getattr(rv, 'args', None)
+            if args is not None:
+                if args:
+                    newargs = tuple([walk(a, F) for a in args])
+                    if args != newargs:
+                        rv = rv.func(*newargs)
+                        if simultaneous:
+                            # if rv is something that was already
+                            # matched (that was changed) then skip
+                            # applying F again
+                            for i, e in enumerate(args):
+                                if rv == e and e != newargs[i]:
+                                    return rv
+                rv = F(rv)
+            return rv
+
+
         mapping = {}  # changes that took place
-        mask = []  # the dummies that were used as change placeholders
 
         def rec_replace(expr):
             result = _query(expr)
             if result or result == {}:
-                new = _value(expr, result)
-                if new is not None and new != expr:
-                    mapping[expr] = new
-                    if simultaneous:
-                        # don't let this change during rebuilding;
-                        # XXX this may fail if the object being replaced
-                        # cannot be represented as a Dummy in the expression
-                        # tree, e.g. an ExprConditionPair in Piecewise
-                        # cannot be represented with a Dummy
-                        com = getattr(new, 'is_commutative', True)
-                        if com is None:
-                            com = True
-                        d = Dummy('rec_replace', commutative=com)
-                        mask.append((d, new))
-                        expr = d
-                    else:
-                        expr = new
+                v = _value(expr, result)
+                if v is not None and v != expr:
+                    if map:
+                        mapping[expr] = v
+                    expr = v
             return expr
 
-        rv = bottom_up(self, rec_replace, atoms=True)
-
-        # restore original expressions for Dummy symbols
-        if simultaneous:
-            mask = list(reversed(mask))
-            for o, n in mask:
-                r = {o: n}
-                # if a sub-expression could not be replaced with
-                # a Dummy then this will fail; either filter
-                # against such sub-expressions or figure out a
-                # way to carry out simultaneous replacement
-                # in this situation.
-                rv = rv.xreplace(r)  # if this fails, see above
-
-        if not map:
-            return rv
-        else:
-            if simultaneous:
-                # restore subexpressions in mapping
-                for o, n in mask:
-                    r = {o: n}
-                    mapping = {k.xreplace(r): v.xreplace(r)
-                        for k, v in mapping.items()}
-            return rv, mapping
+        rv = walk(self, rec_replace)
+        return (rv, mapping) if map else rv
 
     def find(self, query, group=False):
         """Find all subexpressions matching a query. """
@@ -1569,6 +1561,7 @@ def matches(self, expr, repl_dict={}, old=False):
         >>> Basic(a + x, x).matches(Basic(a + b + c, b + c))
         {x_: b + c}
         """
+        repl_dict = repl_dict.copy()
         expr = sympify(expr)
         if not isinstance(expr, self.__class__):
             return None
@@ -1602,7 +1595,7 @@ def match(self, pattern, old=False):
         Examples
         ========
 
-        >>> from sympy import Wild
+        >>> from sympy import Wild, Sum
         >>> from sympy.abc import x, y
         >>> p = Wild("p")
         >>> q = Wild("q")
@@ -1618,6 +1611,16 @@ def match(self, pattern, old=False):
         >>> (p*q**r).xreplace(e.match(p*q**r))
         4*x**2
 
+        Structurally bound symbols are ignored during matching:
+
+        >>> Sum(x, (x, 1, 2)).match(Sum(y, (y, 1, p)))
+        {p_: 2}
+
+        But they can be identified if desired:
+
+        >>> Sum(x, (x, 1, 2)).match(Sum(q, (q, 1, p)))
+        {p_: 2, q_: x}
+
         The ``old`` flag will give the old-style pattern matching where
         expressions and patterns are essentially solved to give the
         match. Both of the following give None unless ``old=True``:
@@ -1628,8 +1631,37 @@ def match(self, pattern, old=False):
         {p_: 2/x**2}
 
         """
+        from sympy.core.symbol import Wild
+        from sympy.core.function import WildFunction
+        from sympy.utilities.misc import filldedent
+
         pattern = sympify(pattern)
-        return pattern.matches(self, old=old)
+        # match non-bound symbols
+        canonical = lambda x: x if x.is_Symbol else x.as_dummy()
+        m = canonical(pattern).matches(canonical(self), old=old)
+        if m is None:
+            return m
+        wild = pattern.atoms(Wild, WildFunction)
+        # sanity check
+        if set(m) - wild:
+            raise ValueError(filldedent('''
+            Some `matches` routine did not use a copy of repl_dict
+            and injected unexpected symbols. Report this as an
+            error at https://github.com/sympy/sympy/issues'''))
+        # now see if bound symbols were requested
+        bwild = wild - set(m)
+        if not bwild:
+            return m
+        # replace free-Wild symbols in pattern with match result
+        # so they will match but not be in the next match
+        wpat = pattern.xreplace(m)
+        # identify remaining bound wild
+        w = wpat.matches(self, old=old)
+        # add them to m
+        if w:
+            m.update(w)
+        # done
+        return m
 
     def count_ops(self, visual=None):
         """wrapper for count_ops that returns the operation count."""
@@ -1839,7 +1871,7 @@ class Atom(Basic):
 
     def matches(self, expr, repl_dict={}, old=False):
         if self == expr:
-            return repl_dict
+            return repl_dict.copy()
 
     def xreplace(self, rule, hack2=False):
         return rule.get(self, self)
diff --git a/sympy/core/compatibility.py b/sympy/core/compatibility.py
index 18836ab7ce..34cf359cd2 100644
--- a/sympy/core/compatibility.py
+++ b/sympy/core/compatibility.py
@@ -553,8 +553,11 @@ def _nodes(e):
     for which the sum of nodes is returned).
     """
     from .basic import Basic
+    from .function import Derivative
 
     if isinstance(e, Basic):
+        if isinstance(e, Derivative):
+            return _nodes(e.expr) + len(e.variables)
         return e.count(Basic)
     elif iterable(e):
         return 1 + sum(_nodes(ei) for ei in e)
diff --git a/sympy/core/function.py b/sympy/core/function.py
index d3c85215ee..d769951163 100644
--- a/sympy/core/function.py
+++ b/sympy/core/function.py
@@ -1918,6 +1918,9 @@ class Lambda(Expr):
     'lambda x: expr'. A function of several variables is written as
     Lambda((x, y, ...), expr).
 
+    Examples
+    ========
+
     A simple example:
 
     >>> from sympy import Lambda
@@ -2057,24 +2060,6 @@ def rmatch(pars, args):
 
         return symargmap
 
-    def __eq__(self, other):
-        if not isinstance(other, Lambda):
-            return False
-        if self.nargs != other.nargs:
-            return False
-
-        try:
-            d = self._match_signature(other.signature, self.signature)
-        except BadArgumentsError:
-            return False
-        return self.args == other.xreplace(d).args
-
-    def __hash__(self):
-        return super().__hash__()
-
-    def _hashable_content(self):
-        return (self.expr.xreplace(self.canonical_variables),)
-
     @property
     def is_identity(self):
         """Return ``True`` if this ``Lambda`` is an identity function. """
@@ -3088,7 +3073,7 @@ def count_ops(expr, visual=False):
     2*ADD + SIN
 
     """
-    from sympy import Integral, Symbol
+    from sympy import Integral, Sum, Symbol
     from sympy.core.relational import Relational
     from sympy.simplify.radsimp import fraction
     from sympy.logic.boolalg import BooleanFunction
@@ -3156,18 +3141,22 @@ def count_ops(expr, visual=False):
                 ops.append(DIV)
                 args.append(a.base)  # won't be -Mul but could be Add
                 continue
-            if (a.is_Mul or
-                a.is_Pow or
-                a.is_Function or
-                isinstance(a, Derivative) or
-                    isinstance(a, Integral)):
-
+            if a.is_Mul or isinstance(a, LatticeOp):
                 o = Symbol(a.func.__name__.upper())
                 # count the args
-                if (a.is_Mul or isinstance(a, LatticeOp)):
-                    ops.append(o*(len(a.args) - 1))
-                else:
-                    ops.append(o)
+                ops.append(o*(len(a.args) - 1))
+            elif a.args and (
+                    a.is_Pow or
+                    a.is_Function or
+                    isinstance(a, Derivative) or
+                    isinstance(a, Integral) or
+                    isinstance(a, Sum)):
+                # if it's not in the list above we don't
+                # consider a.func something to count, e.g.
+                # Tuple, MatrixSymbol, etc...
+                o = Symbol(a.func.__name__.upper())
+                ops.append(o)
+
             if not a.is_Symbol:
                 args.extend(a.args)
 
diff --git a/sympy/core/mul.py b/sympy/core/mul.py
index 506b0f0204..7b2f28cc8d 100644
--- a/sympy/core/mul.py
+++ b/sympy/core/mul.py
@@ -952,6 +952,7 @@ def _matches_simple(self, expr, repl_dict):
 
     def matches(self, expr, repl_dict={}, old=False):
         expr = sympify(expr)
+        repl_dict = repl_dict.copy()
         if self.is_commutative and expr.is_commutative:
             return self._matches_commutative(expr, repl_dict, old)
         elif self.is_commutative is not expr.is_commutative:
@@ -1000,6 +1001,7 @@ def _matches_noncomm(nodes, targets, repl_dict={}):
         expression, while `targets` is a list of arguments in the
         multiplication expression being matched against.
         """
+        repl_dict = repl_dict.copy()
         # List of possible future states to be considered
         agenda = []
         # The current matching state, storing index in nodes and targets
diff --git a/sympy/core/operations.py b/sympy/core/operations.py
index bcf9c8ad23..da5337c1a1 100644
--- a/sympy/core/operations.py
+++ b/sympy/core/operations.py
@@ -182,6 +182,7 @@ def _matches_commutative(self, expr, repl_dict={}, old=False):
         # make sure expr is Expr if pattern is Expr
         from .expr import Add, Expr
         from sympy import Mul
+        repl_dict = repl_dict.copy()
         if isinstance(self, Expr) and not isinstance(expr, Expr):
             return None
 
@@ -201,6 +202,13 @@ def _matches_commutative(self, expr, repl_dict={}, old=False):
             binary=True)
         if not exact_part:
             wild_part = list(ordered(wild_part))
+            if self.is_Add:
+                # in addition to normal ordered keys, impose
+                # sorting on Muls with leading Number to put
+                # them in order
+                wild_part = sorted(wild_part, key=lambda x:
+                    x.args[0] if x.is_Mul and x.args[0].is_Number else
+                    0)
         else:
             exact = self._new_rawargs(*exact_part)
             free = expr.free_symbols
@@ -221,7 +229,15 @@ def _matches_commutative(self, expr, repl_dict={}, old=False):
         saw = set()
         while expr not in saw:
             saw.add(expr)
-            expr_list = (self.identity,) + tuple(ordered(self.make_args(expr)))
+            args = tuple(ordered(self.make_args(expr)))
+            if self.is_Add and expr.is_Add:
+                # in addition to normal ordered keys, impose
+                # sorting on Muls with leading Number to put
+                # them in order
+                args = tuple(sorted(args, key=lambda x:
+                    x.args[0] if x.is_Mul and x.args[0].is_Number else
+                    0))
+            expr_list = (self.identity,) + args
             for last_op in reversed(expr_list):
                 for w in reversed(wild_part):
                     d1 = w.matches(last_op, repl_dict)
diff --git a/sympy/core/power.py b/sympy/core/power.py
index 1397e5eeae..239b6221e8 100644
--- a/sympy/core/power.py
+++ b/sympy/core/power.py
@@ -1433,11 +1433,11 @@ def as_numer_denom(self):
 
     def matches(self, expr, repl_dict={}, old=False):
         expr = _sympify(expr)
+        repl_dict = repl_dict.copy()
 
         # special case, pattern = 1 and expr.exp can match to 0
         if expr is S.One:
-            d = repl_dict.copy()
-            d = self.exp.matches(S.Zero, d)
+            d = self.exp.matches(S.Zero, repl_dict)
             if d is not None:
                 return d
 
diff --git a/sympy/core/symbol.py b/sympy/core/symbol.py
index 258bacd440..1839feccc0 100644
--- a/sympy/core/symbol.py
+++ b/sympy/core/symbol.py
@@ -273,7 +273,9 @@ def sort_key(self, order=None):
         return self.class_key(), (1, (self.name,)), S.One.sort_key(), S.One
 
     def as_dummy(self):
-        return Dummy(self.name)
+        # only put commutativity in explicitly if it is False
+        return Dummy(self.name) if self.is_commutative is not False \
+            else Dummy(self.name, commutative=self.is_commutative)
 
     def as_real_imag(self, deep=True, **hints):
         from sympy import im, re
diff --git a/sympy/functions/special/tensor_functions.py b/sympy/functions/special/tensor_functions.py
index 92a6f2aecf..3b3f86065e 100644
--- a/sympy/functions/special/tensor_functions.py
+++ b/sympy/functions/special/tensor_functions.py
@@ -189,7 +189,7 @@ def eval(cls, i, j, delta_range=None):
         # to make KroneckerDelta canonical
         # following lines will check if inputs are in order
         # if not, will return KroneckerDelta with correct order
-        if i is not min(i, j, key=default_sort_key):
+        if i != min(i, j, key=default_sort_key):
             if delta_range:
                 return cls(j, i, delta_range)
             else:
diff --git a/sympy/matrices/expressions/applyfunc.py b/sympy/matrices/expressions/applyfunc.py
index e6b44f8038..36870d4410 100644
--- a/sympy/matrices/expressions/applyfunc.py
+++ b/sympy/matrices/expressions/applyfunc.py
@@ -49,7 +49,7 @@ def __new__(cls, function, expr):
         if not expr.is_Matrix:
             raise ValueError("{} must be a matrix instance.".format(expr))
 
-        if not isinstance(function, FunctionClass):
+        if not isinstance(function, (FunctionClass, Lambda)):
             d = Dummy('d')
             function = Lambda(d, function(d))
 
diff --git a/sympy/physics/quantum/cg.py b/sympy/physics/quantum/cg.py
index 6aa000c4ba..d89ba5a2e6 100644
--- a/sympy/physics/quantum/cg.py
+++ b/sympy/physics/quantum/cg.py
@@ -681,21 +681,21 @@ def _check_varsh_sum_871_2(e):
 
 
 def _check_varsh_sum_872_4(e):
+    alpha = symbols('alpha')
+    beta = symbols('beta')
     a = Wild('a')
-    alpha = Wild('alpha')
     b = Wild('b')
-    beta = Wild('beta')
     c = Wild('c')
     cp = Wild('cp')
     gamma = Wild('gamma')
     gammap = Wild('gammap')
-    match1 = e.match(Sum(CG(a, alpha, b, beta, c, gamma)*CG(
-        a, alpha, b, beta, cp, gammap), (alpha, -a, a), (beta, -b, b)))
-    if match1 is not None and len(match1) == 8:
+    cg1 = CG(a, alpha, b, beta, c, gamma)
+    cg2 = CG(a, alpha, b, beta, cp, gammap)
+    match1 = e.match(Sum(cg1*cg2, (alpha, -a, a), (beta, -b, b)))
+    if match1 is not None and len(match1) == 6:
         return (KroneckerDelta(c, cp)*KroneckerDelta(gamma, gammap)).subs(match1)
-    match2 = e.match(Sum(
-        CG(a, alpha, b, beta, c, gamma)**2, (alpha, -a, a), (beta, -b, b)))
-    if match2 is not None and len(match2) == 6:
+    match2 = e.match(Sum(cg1**2, (alpha, -a, a), (beta, -b, b)))
+    if match2 is not None and len(match2) == 4:
         return 1
     return e
 
diff --git a/sympy/sets/sets.py b/sympy/sets/sets.py
index 004e47d2ad..1e3c9d292f 100644
--- a/sympy/sets/sets.py
+++ b/sympy/sets/sets.py
@@ -1770,7 +1770,18 @@ def __new__(cls, *args, **kwargs):
         else:
             args = list(map(sympify, args))
 
-        _args_set = set(args)
+        # keep the form of the first canonical arg
+        dargs = {}
+        for i in reversed(list(ordered(args))):
+            if i.is_Symbol:
+                dargs[i] = i
+            else:
+                try:
+                    dargs[i.as_dummy()] = i
+                except TypeError:
+                    # e.g. i = class without args like `Interval`
+                    dargs[i] = i
+        _args_set = set(dargs.values())
         args = list(ordered(_args_set, Set._infimum_key))
         obj = Basic.__new__(cls, *args)
         obj._args_set = _args_set


Test patch:
diff --git a/sympy/codegen/tests/test_array_utils.py b/sympy/codegen/tests/test_array_utils.py
index d3e1bef324..7e0275b8f8 100644
--- a/sympy/codegen/tests/test_array_utils.py
+++ b/sympy/codegen/tests/test_array_utils.py
@@ -465,13 +465,13 @@ def test_recognize_diagonalized_vectors():
     assert recognize_matrix_expression(cg) == A*DiagMatrix(a)*C*DiagMatrix(a)*B
 
     cg = CodegenArrayContraction(CodegenArrayTensorProduct(a, I1, b, I1, (a.T*b).applyfunc(cos)), (1, 2, 8), (5, 6, 9))
-    assert cg.split_multiple_contractions() == CodegenArrayContraction(CodegenArrayTensorProduct(a, I1, b, I1, (a.T*b).applyfunc(cos)), (1, 2), (3, 8), (5, 6), (7, 9))
-    assert recognize_matrix_expression(cg) == MatMul(a, I1, (a.T*b).applyfunc(cos), Transpose(I1), b.T)
+    assert cg.split_multiple_contractions().dummy_eq(CodegenArrayContraction(CodegenArrayTensorProduct(a, I1, b, I1, (a.T*b).applyfunc(cos)), (1, 2), (3, 8), (5, 6), (7, 9)))
+    assert recognize_matrix_expression(cg).dummy_eq(MatMul(a, I1, (a.T*b).applyfunc(cos), Transpose(I1), b.T))
 
     cg = CodegenArrayContraction(CodegenArrayTensorProduct(A.T, a, b, b.T, (A*X*b).applyfunc(cos)), (1, 2, 8), (5, 6, 9))
-    assert cg.split_multiple_contractions() == CodegenArrayContraction(
+    assert cg.split_multiple_contractions().dummy_eq(CodegenArrayContraction(
         CodegenArrayTensorProduct(A.T, DiagMatrix(a), b, b.T, (A*X*b).applyfunc(cos)),
-                               (1, 2), (3, 8), (5, 6, 9))
+                               (1, 2), (3, 8), (5, 6, 9)))
     # assert recognize_matrix_expression(cg)
 
     # Check no overlap of lines:
diff --git a/sympy/core/tests/test_basic.py b/sympy/core/tests/test_basic.py
index 4c4dc187a3..f686f0f65d 100644
--- a/sympy/core/tests/test_basic.py
+++ b/sympy/core/tests/test_basic.py
@@ -7,7 +7,7 @@
 from sympy.core.basic import (Basic, Atom, preorder_traversal, as_Basic,
     _atomic, _aresame)
 from sympy.core.singleton import S
-from sympy.core.symbol import symbols, Symbol
+from sympy.core.symbol import symbols, Symbol, Dummy
 from sympy.core.sympify import SympifyError
 from sympy.core.function import Function, Lambda
 from sympy.core.compatibility import default_sort_key
@@ -281,12 +281,27 @@ def test_as_dummy():
     u, v, x, y, z, _0, _1 = symbols('u v x y z _0 _1')
     assert Lambda(x, x + 1).as_dummy() == Lambda(_0, _0 + 1)
     assert Lambda(x, x + _0).as_dummy() == Lambda(_1, _0 + _1)
-    assert (1 + Sum(x, (x, 1, x))).as_dummy() == 1 + Sum(_0, (_0, 1, x))
+    eq = (1 + Sum(x, (x, 1, x)))
+    ans = 1 + Sum(_0, (_0, 1, x))
+    once = eq.as_dummy()
+    assert once == ans
+    twice = once.as_dummy()
+    assert twice == ans
+    assert Integral(x + _0, (x, x + 1), (_0, 1, 2)
+        ).as_dummy() == Integral(_0 + _1, (_0, x + 1), (_1, 1, 2))
+    for T in (Symbol, Dummy):
+        d = T('x', real=True)
+        D = d.as_dummy()
+        assert D != d and D.func == Dummy and D.is_real is None
+    assert Dummy().as_dummy().is_commutative
+    assert Dummy(commutative=False).as_dummy().is_commutative is False
 
 
 def test_canonical_variables():
     x, i0, i1 = symbols('x _:2')
     assert Integral(x, (x, x + 1)).canonical_variables == {x: i0}
+    assert Integral(x, (x, x + 1), (i0, 1, 2)).canonical_variables == {
+        x: i0, i0: i1}
     assert Integral(x, (x, x + i0)).canonical_variables == {x: i1}
 
 
diff --git a/sympy/core/tests/test_count_ops.py b/sympy/core/tests/test_count_ops.py
index 6e377848f6..80fab1caf2 100644
--- a/sympy/core/tests/test_count_ops.py
+++ b/sympy/core/tests/test_count_ops.py
@@ -1,6 +1,6 @@
 from sympy import symbols, sin, exp, cos, Derivative, Integral, Basic, \
     count_ops, S, And, I, pi, Eq, Or, Not, Xor, Nand, Nor, Implies, \
-    Equivalent, MatrixSymbol, Symbol, ITE, Rel, Rational
+    Equivalent, MatrixSymbol, Symbol, ITE, Rel, Rational, Sum
 from sympy.core.containers import Tuple
 
 x, y, z = symbols('x,y,z')
@@ -30,8 +30,8 @@ def count(val):
 
 
 def test_count_ops_visual():
-    ADD, MUL, POW, SIN, COS, EXP, AND, D, G = symbols(
-        'Add Mul Pow sin cos exp And Derivative Integral'.upper())
+    ADD, MUL, POW, SIN, COS, EXP, AND, D, G, M = symbols(
+        'Add Mul Pow sin cos exp And Derivative Integral Sum'.upper())
     DIV, SUB, NEG = symbols('DIV SUB NEG')
     LT, LE, GT, GE, EQ, NE = symbols('LT LE GT GE EQ NE')
     NOT, OR, AND, XOR, IMPLIES, EQUIVALENT, _ITE, BASIC, TUPLE = symbols(
@@ -88,6 +88,7 @@ def count(val):
 
     assert count(Derivative(x, x)) == D
     assert count(Integral(x, x) + 2*x/(1 + x)) == G + DIV + MUL + 2*ADD
+    assert count(Sum(x, (x, 1, x + 1)) + 2*x/(1 + x)) == M + DIV + MUL + 3*ADD
     assert count(Basic()) is S.Zero
 
     assert count({x + 1: sin(x)}) == ADD + SIN
diff --git a/sympy/core/tests/test_expr.py b/sympy/core/tests/test_expr.py
index b0701d95d1..b9924f4d9e 100644
--- a/sympy/core/tests/test_expr.py
+++ b/sympy/core/tests/test_expr.py
@@ -909,7 +909,8 @@ def test_replace():
     # if not simultaneous then y*sin(x) -> y*sin(x)/y = sin(x) -> sin(x)/y
     assert (y*sin(x)).replace(sin, lambda expr: sin(expr)/y,
         simultaneous=False) == sin(x)/y
-    assert (x**2 + O(x**3)).replace(Pow, lambda b, e: b**e/e) == O(1, x)
+    assert (x**2 + O(x**3)).replace(Pow, lambda b, e: b**e/e
+        ) == x**2/2 + O(x**3)
     assert (x**2 + O(x**3)).replace(Pow, lambda b, e: b**e/e,
         simultaneous=False) == x**2/2 + O(x**3)
     assert (x*(x*y + 3)).replace(lambda x: x.is_Mul, lambda x: 2 + x) == \
diff --git a/sympy/core/tests/test_function.py b/sympy/core/tests/test_function.py
index c45ff01744..97979fa80f 100644
--- a/sympy/core/tests/test_function.py
+++ b/sympy/core/tests/test_function.py
@@ -218,9 +218,8 @@ def test_Lambda():
 
     assert unchanged(Lambda, (x,), x**2)
     assert Lambda(x, x**2) == Lambda((x,), x**2)
-    assert Lambda(x, x**2) == Lambda(y, y**2)
-    assert Lambda(x, x**2) != Lambda(y, y**2 + 1)
-    assert Lambda((x, y), x**y) == Lambda((y, x), y**x)
+    assert Lambda(x, x**2) != Lambda(x, x**2 + 1)
+    assert Lambda((x, y), x**y) != Lambda((y, x), y**x)
     assert Lambda((x, y), x**y) != Lambda((x, y), y**x)
 
     assert Lambda((x, y), x**y)(x, y) == x**y
@@ -239,7 +238,9 @@ def test_Lambda():
     p = x, y, z, t
     assert Lambda(p, t*(x + y + z))(*p) == t * (x + y + z)
 
-    assert Lambda(x, 2*x) + Lambda(y, 2*y) == 2*Lambda(x, 2*x)
+    eq = Lambda(x, 2*x) + Lambda(y, 2*y)
+    assert eq != 2*Lambda(x, 2*x)
+    assert eq.as_dummy() == 2*Lambda(x, 2*x).as_dummy()
     assert Lambda(x, 2*x) not in [ Lambda(x, x) ]
     raises(BadSignatureError, lambda: Lambda(1, x))
     assert Lambda(x, 1)(1) is S.One
@@ -300,12 +301,19 @@ def test_Lambda_arguments():
 
 
 def test_Lambda_equality():
-    assert Lambda(x, 2*x) == Lambda(y, 2*y)
-    # although variables are casts as Dummies, the expressions
-    # should still compare equal
     assert Lambda((x, y), 2*x) == Lambda((x, y), 2*x)
+    # these, of course, should never be equal
     assert Lambda(x, 2*x) != Lambda((x, y), 2*x)
     assert Lambda(x, 2*x) != 2*x
+    # But it is tempting to want expressions that differ only
+    # in bound symbols to compare the same.  But this is not what
+    # Python's `==` is intended to do; two objects that compare
+    # as equal means that they are indistibguishable and cache to the
+    # same value.  We wouldn't want to expression that are
+    # mathematically the same but written in different variables to be
+    # interchanged else what is the point of allowing for different
+    # variable names?
+    assert Lambda(x, 2*x) != Lambda(y, 2*y)
 
 
 def test_Subs():
diff --git a/sympy/core/tests/test_match.py b/sympy/core/tests/test_match.py
index aefad7b83b..255ce674d5 100644
--- a/sympy/core/tests/test_match.py
+++ b/sympy/core/tests/test_match.py
@@ -1,6 +1,6 @@
 from sympy import (abc, Add, cos, collect, Derivative, diff, exp, Float, Function,
     I, Integer, log, Mul, oo, Poly, Rational, S, sin, sqrt, Symbol, symbols,
-    Wild, pi, meijerg
+    Wild, pi, meijerg, Sum
 )
 
 from sympy.testing.pytest import XFAIL
@@ -508,7 +508,7 @@ def test_issue_3883():
     a, b, c = symbols('a b c', cls=Wild, exclude=(gamma,))
 
     assert f.match(a * log(gamma) + b * gamma + c) == \
-        {a: Rational(-1, 2), b: -(x - mu)**2/2, c: log(2*pi)/2}
+        {a: Rational(-1, 2), b: -(mu - x)**2/2, c: log(2*pi)/2}
     assert f.expand().collect(gamma).match(a * log(gamma) + b * gamma + c) == \
         {a: Rational(-1, 2), b: (-(x - mu)**2/2).expand(), c: (log(2*pi)/2).expand()}
     g1 = Wild('g1', exclude=[gamma])
@@ -706,3 +706,19 @@ def test_match_issue_17397():
         - 4*Derivative(f(x), (x, 2)) - 2*Derivative(f(x), x)/x + 4*Derivative(f(x), (x, 2))/x
     r = collect(eq, [f(x).diff(x, 2), f(x).diff(x), f(x)]).match(deq)
     assert r == {a3: x - 4 + 4/x, b3: 1 - 2/x, c3: x - 4}
+
+
+def test_match_terms():
+    X, Y = map(Wild, "XY")
+    x, y, z = symbols('x y z')
+    assert (5*y - x).match(5*X - Y) == {X: y, Y: x}
+    # 15907
+    assert (x + (y - 1)*z).match(x + X*z) == {X: y - 1}
+
+
+def test_match_bound():
+    V, W = map(Wild, "VW")
+    x, y = symbols('x y')
+    assert Sum(x, (x, 1, 2)).match(Sum(y, (y, 1, W))) == {W: 2}
+    assert Sum(x, (x, 1, 2)).match(Sum(V, (V, 1, W))) == {W: 2, V:x}
+    assert Sum(x, (x, 1, 2)).match(Sum(V, (V, 1, 2))) == {V:x}
diff --git a/sympy/core/tests/test_subs.py b/sympy/core/tests/test_subs.py
index 53d428269a..d2f4eeb329 100644
--- a/sympy/core/tests/test_subs.py
+++ b/sympy/core/tests/test_subs.py
@@ -850,3 +850,8 @@ def test_issue_17823():
     expr = q1.diff().diff()**2*q1 + q1.diff()*q2.diff()
     reps={q1: a, q1.diff(): a*x*y, q1.diff().diff(): z}
     assert expr.subs(reps) == a*x*y*Derivative(q2, t) + a*z**2
+
+
+def test_issue_19326():
+    x, y = [i(t) for i in map(Function, 'xy')]
+    assert (x*y).subs({x: 1 + x, y: x}) == (1 + x)*x
diff --git a/sympy/integrals/tests/test_rationaltools.py b/sympy/integrals/tests/test_rationaltools.py
index fb13ca381b..ff2b0529bb 100644
--- a/sympy/integrals/tests/test_rationaltools.py
+++ b/sympy/integrals/tests/test_rationaltools.py
@@ -1,4 +1,5 @@
-from sympy import S, symbols, I, atan, log, Poly, sqrt, simplify, integrate, Rational
+from sympy import (S, symbols, I, atan, log, Poly, sqrt, simplify,
+    integrate, Rational, Dummy)
 
 from sympy.integrals.rationaltools import ratint, ratint_logpart, log_to_atan
 
@@ -99,6 +100,10 @@ def test_ratint():
     assert ratint(1/(x**2 + 1), x, symbol=x) == ans
     assert ratint(1/(x**2 + 1), x, symbol='x') == ans
     assert ratint(1/(x**2 + 1), x, symbol=a) == ans
+    # this asserts that as_dummy must return a unique symbol
+    # even if the symbol is already a Dummy
+    d = Dummy()
+    assert ratint(1/(d**2 + 1), d, symbol=d) == atan(d)
 
 
 def test_ratint_logpart():
diff --git a/sympy/matrices/expressions/tests/test_applyfunc.py b/sympy/matrices/expressions/tests/test_applyfunc.py
index c1e923d2f9..0d62274977 100644
--- a/sympy/matrices/expressions/tests/test_applyfunc.py
+++ b/sympy/matrices/expressions/tests/test_applyfunc.py
@@ -36,8 +36,8 @@ def test_applyfunc_matrix():
 
     expr = ElementwiseApplyFunction(exp, X*Y)
     assert expr.expr == X*Y
-    assert expr.function == Lambda(x, exp(x))
-    assert expr == (X*Y).applyfunc(exp)
+    assert expr.function.dummy_eq(Lambda(x, exp(x)))
+    assert expr.dummy_eq((X*Y).applyfunc(exp))
     assert expr.func(*expr.args) == expr
 
     assert isinstance(X*expr, MatMul)
@@ -55,7 +55,7 @@ def test_applyfunc_matrix():
     M = Matrix([[x, y], [z, t]])
     expr = ElementwiseApplyFunction(sin, M)
     assert isinstance(expr, ElementwiseApplyFunction)
-    assert expr.function == Lambda(x, sin(x))
+    assert expr.function.dummy_eq(Lambda(x, sin(x)))
     assert expr.expr == M
     assert expr.doit() == M.applyfunc(sin)
     assert expr.doit() == Matrix([[sin(x), sin(y)], [sin(z), sin(t)]])
diff --git a/sympy/matrices/expressions/tests/test_derivatives.py b/sympy/matrices/expressions/tests/test_derivatives.py
index 1f0f357bc7..c0089f0821 100644
--- a/sympy/matrices/expressions/tests/test_derivatives.py
+++ b/sympy/matrices/expressions/tests/test_derivatives.py
@@ -56,7 +56,8 @@ def test_matrix_derivative_by_scalar():
     assert (x*x.T).diff(i) == ZeroMatrix(k, k)
     assert (x + y).diff(i) == ZeroMatrix(k, 1)
     assert hadamard_power(x, 2).diff(i) == ZeroMatrix(k, 1)
-    assert hadamard_power(x, i).diff(i) == HadamardProduct(x.applyfunc(log), HadamardPower(x, i))
+    assert hadamard_power(x, i).diff(i).dummy_eq(
+        HadamardProduct(x.applyfunc(log), HadamardPower(x, i)))
     assert hadamard_product(x, y).diff(i) == ZeroMatrix(k, 1)
     assert hadamard_product(i*OneMatrix(k, 1), x, y).diff(i) == hadamard_product(x, y)
     assert (i*x).diff(i) == x
@@ -360,41 +361,46 @@ def test_derivatives_elementwise_applyfunc():
     from sympy.matrices.expressions.diagonal import DiagMatrix
 
     expr = x.applyfunc(tan)
-    assert expr.diff(x) == DiagMatrix(x.applyfunc(lambda x: tan(x)**2 + 1))
+    assert expr.diff(x).dummy_eq(
+        DiagMatrix(x.applyfunc(lambda x: tan(x)**2 + 1)))
     assert expr[i, 0].diff(x[m, 0]).doit() == (tan(x[i, 0])**2 + 1)*KDelta(i, m)
     _check_derivative_with_explicit_matrix(expr, x, expr.diff(x))
 
     expr = (i**2*x).applyfunc(sin)
-    assert expr.diff(i) == HadamardProduct((2*i)*x, (i**2*x).applyfunc(cos))
+    assert expr.diff(i).dummy_eq(
+        HadamardProduct((2*i)*x, (i**2*x).applyfunc(cos)))
     assert expr[i, 0].diff(i).doit() == 2*i*x[i, 0]*cos(i**2*x[i, 0])
     _check_derivative_with_explicit_matrix(expr, i, expr.diff(i))
 
     expr = (log(i)*A*B).applyfunc(sin)
-    assert expr.diff(i) == HadamardProduct(A*B/i, (log(i)*A*B).applyfunc(cos))
+    assert expr.diff(i).dummy_eq(
+        HadamardProduct(A*B/i, (log(i)*A*B).applyfunc(cos)))
     _check_derivative_with_explicit_matrix(expr, i, expr.diff(i))
 
     expr = A*x.applyfunc(exp)
-    assert expr.diff(x) == DiagMatrix(x.applyfunc(exp))*A.T
+    assert expr.diff(x).dummy_eq(DiagMatrix(x.applyfunc(exp))*A.T)
     _check_derivative_with_explicit_matrix(expr, x, expr.diff(x))
 
     expr = x.T*A*x + k*y.applyfunc(sin).T*x
-    assert expr.diff(x) == A.T*x + A*x + k*y.applyfunc(sin)
+    assert expr.diff(x).dummy_eq(A.T*x + A*x + k*y.applyfunc(sin))
     _check_derivative_with_explicit_matrix(expr, x, expr.diff(x))
 
     expr = x.applyfunc(sin).T*y
-    assert expr.diff(x) == DiagMatrix(x.applyfunc(cos))*y
+    assert expr.diff(x).dummy_eq(DiagMatrix(x.applyfunc(cos))*y)
     _check_derivative_with_explicit_matrix(expr, x, expr.diff(x))
 
     expr = (a.T * X * b).applyfunc(sin)
-    assert expr.diff(X) == a*(a.T*X*b).applyfunc(cos)*b.T
+    assert expr.diff(X).dummy_eq(a*(a.T*X*b).applyfunc(cos)*b.T)
     _check_derivative_with_explicit_matrix(expr, X, expr.diff(X))
 
     expr = a.T * X.applyfunc(sin) * b
-    assert expr.diff(X) == DiagMatrix(a)*X.applyfunc(cos)*DiagMatrix(b)
+    assert expr.diff(X).dummy_eq(
+        DiagMatrix(a)*X.applyfunc(cos)*DiagMatrix(b))
     _check_derivative_with_explicit_matrix(expr, X, expr.diff(X))
 
     expr = a.T * (A*X*B).applyfunc(sin) * b
-    assert expr.diff(X) == A.T*DiagMatrix(a)*(A*X*B).applyfunc(cos)*DiagMatrix(b)*B.T
+    assert expr.diff(X).dummy_eq(
+        A.T*DiagMatrix(a)*(A*X*B).applyfunc(cos)*DiagMatrix(b)*B.T)
     _check_derivative_with_explicit_matrix(expr, X, expr.diff(X))
 
     expr = a.T * (A*X*b).applyfunc(sin) * b.T
@@ -403,7 +409,8 @@ def test_derivatives_elementwise_applyfunc():
     #_check_derivative_with_explicit_matrix(expr, X, expr.diff(X))
 
     expr = a.T*A*X.applyfunc(sin)*B*b
-    assert expr.diff(X) == DiagMatrix(A.T*a)*X.applyfunc(cos)*DiagMatrix(B*b)
+    assert expr.diff(X).dummy_eq(
+        DiagMatrix(A.T*a)*X.applyfunc(cos)*DiagMatrix(B*b))
 
     expr = a.T * (A*X.applyfunc(sin)*B).applyfunc(log) * b
     # TODO: wrong
diff --git a/sympy/matrices/expressions/tests/test_funcmatrix.py b/sympy/matrices/expressions/tests/test_funcmatrix.py
index 7ec810b01a..4479b33dfc 100644
--- a/sympy/matrices/expressions/tests/test_funcmatrix.py
+++ b/sympy/matrices/expressions/tests/test_funcmatrix.py
@@ -28,7 +28,7 @@ def test_funcmatrix_creation():
 
     m = FunctionMatrix(2, 2, KroneckerDelta)
     assert m.as_explicit() == Identity(2).as_explicit()
-    assert m.args[2] == Lambda((i, j), KroneckerDelta(i, j))
+    assert m.args[2].dummy_eq(Lambda((i, j), KroneckerDelta(i, j)))
 
     n = symbols('n')
     assert FunctionMatrix(n, n, Lambda((i, j), 0))
diff --git a/sympy/polys/tests/test_partfrac.py b/sympy/polys/tests/test_partfrac.py
index aa697ec448..bcf2411199 100644
--- a/sympy/polys/tests/test_partfrac.py
+++ b/sympy/polys/tests/test_partfrac.py
@@ -90,24 +90,24 @@ def test_apart_full():
     f = 1/(x**2 + 1)
 
     assert apart(f, full=False) == f
-    assert apart(f, full=True) == \
-        -RootSum(x**2 + 1, Lambda(a, a/(x - a)), auto=False)/2
+    assert apart(f, full=True).dummy_eq(
+        -RootSum(x**2 + 1, Lambda(a, a/(x - a)), auto=False)/2)
 
     f = 1/(x**3 + x + 1)
 
     assert apart(f, full=False) == f
-    assert apart(f, full=True) == \
+    assert apart(f, full=True).dummy_eq(
         RootSum(x**3 + x + 1,
-        Lambda(a, (a**2*Rational(6, 31) - a*Rational(9, 31) + Rational(4, 31))/(x - a)), auto=False)
+        Lambda(a, (a**2*Rational(6, 31) - a*Rational(9, 31) + Rational(4, 31))/(x - a)), auto=False))
 
     f = 1/(x**5 + 1)
 
     assert apart(f, full=False) == \
         (Rational(-1, 5))*((x**3 - 2*x**2 + 3*x - 4)/(x**4 - x**3 + x**2 -
          x + 1)) + (Rational(1, 5))/(x + 1)
-    assert apart(f, full=True) == \
+    assert apart(f, full=True).dummy_eq(
         -RootSum(x**4 - x**3 + x**2 - x + 1,
-        Lambda(a, a/(x - a)), auto=False)/5 + (Rational(1, 5))/(x + 1)
+        Lambda(a, a/(x - a)), auto=False)/5 + (Rational(1, 5))/(x + 1))
 
 
 def test_apart_undetermined_coeffs():
@@ -126,27 +126,33 @@ def test_apart_undetermined_coeffs():
 
 def test_apart_list():
     from sympy.utilities.iterables import numbered_symbols
+    def dummy_eq(i, j):
+        if type(i) in (list, tuple):
+            return all(dummy_eq(i, j) for i, j in zip(i, j))
+        return i == j or i.dummy_eq(j)
 
     w0, w1, w2 = Symbol("w0"), Symbol("w1"), Symbol("w2")
     _a = Dummy("a")
 
     f = (-2*x - 2*x**2) / (3*x**2 - 6*x)
-    assert apart_list(f, x, dummies=numbered_symbols("w")) == (-1,
-        Poly(Rational(2, 3), x, domain='QQ'),
+    got = apart_list(f, x, dummies=numbered_symbols("w"))
+    ans = (-1, Poly(Rational(2, 3), x, domain='QQ'),
         [(Poly(w0 - 2, w0, domain='ZZ'), Lambda(_a, 2), Lambda(_a, -_a + x), 1)])
+    assert dummy_eq(got, ans)
 
-    assert apart_list(2/(x**2-2), x, dummies=numbered_symbols("w")) == (1,
-                                      Poly(0, x, domain='ZZ'),
-                                      [(Poly(w0**2 - 2, w0, domain='ZZ'),
-                                        Lambda(_a, _a/2),
-                                        Lambda(_a, -_a + x), 1)])
+    got = apart_list(2/(x**2-2), x, dummies=numbered_symbols("w"))
+    ans = (1, Poly(0, x, domain='ZZ'), [(Poly(w0**2 - 2, w0, domain='ZZ'),
+        Lambda(_a, _a/2),
+        Lambda(_a, -_a + x), 1)])
+    assert dummy_eq(got, ans)
 
     f = 36 / (x**5 - 2*x**4 - 2*x**3 + 4*x**2 + x - 2)
-    assert apart_list(f, x, dummies=numbered_symbols("w")) == (1,
-                             Poly(0, x, domain='ZZ'),
-                             [(Poly(w0 - 2, w0, domain='ZZ'), Lambda(_a, 4), Lambda(_a, -_a + x), 1),
-                              (Poly(w1**2 - 1, w1, domain='ZZ'), Lambda(_a, -3*_a - 6), Lambda(_a, -_a + x), 2),
-                              (Poly(w2 + 1, w2, domain='ZZ'), Lambda(_a, -4), Lambda(_a, -_a + x), 1)])
+    got = apart_list(f, x, dummies=numbered_symbols("w"))
+    ans = (1, Poly(0, x, domain='ZZ'),
+        [(Poly(w0 - 2, w0, domain='ZZ'), Lambda(_a, 4), Lambda(_a, -_a + x), 1),
+        (Poly(w1**2 - 1, w1, domain='ZZ'), Lambda(_a, -3*_a - 6), Lambda(_a, -_a + x), 2),
+        (Poly(w2 + 1, w2, domain='ZZ'), Lambda(_a, -4), Lambda(_a, -_a + x), 1)])
+    assert dummy_eq(got, ans)
 
 
 def test_assemble_partfrac_list():
diff --git a/sympy/printing/pretty/tests/test_pretty.py b/sympy/printing/pretty/tests/test_pretty.py
index 09f751c759..1cfd0476fe 100644
--- a/sympy/printing/pretty/tests/test_pretty.py
+++ b/sympy/printing/pretty/tests/test_pretty.py
@@ -3248,13 +3248,13 @@ def test_MatrixExpressions():
     expr = (n*X).applyfunc(lamda)
     ascii_str = """\
 /     1\\      \n\
-|d -> -|.(n*X)\n\
-\\     d/      \
+|x -> -|.(n*X)\n\
+\\     x/      \
 """
     ucode_str = u("""\
 ⎛    1⎞      \n\
-⎜d ↦ ─⎟˳(n⋅X)\n\
-⎝    d⎠      \
+⎜x ↦ ─⎟˳(n⋅X)\n\
+⎝    x⎠      \
 """)
     assert pretty(expr) == ascii_str
     assert upretty(expr) == ucode_str
diff --git a/sympy/printing/tests/test_latex.py b/sympy/printing/tests/test_latex.py
index 5d7c7927c1..4900a85de2 100644
--- a/sympy/printing/tests/test_latex.py
+++ b/sympy/printing/tests/test_latex.py
@@ -1866,7 +1866,7 @@ def test_ElementwiseApplyFunction():
     expr = (X.T*X).applyfunc(sin)
     assert latex(expr) == r"{\left( d \mapsto \sin{\left(d \right)} \right)}_{\circ}\left({X^{T} X}\right)"
     expr = X.applyfunc(Lambda(x, 1/x))
-    assert latex(expr) == r'{\left( d \mapsto \frac{1}{d} \right)}_{\circ}\left({X}\right)'
+    assert latex(expr) == r'{\left( x \mapsto \frac{1}{x} \right)}_{\circ}\left({X}\right)'
 
 
 def test_ZeroMatrix():
diff --git a/sympy/printing/tests/test_repr.py b/sympy/printing/tests/test_repr.py
index d3298da71d..1173fee44d 100644
--- a/sympy/printing/tests/test_repr.py
+++ b/sympy/printing/tests/test_repr.py
@@ -66,6 +66,7 @@ def test_Function():
     sT(sin(x), "sin(Symbol('x'))")
     sT(sin, "sin")
 
+
 def test_Geometry():
     sT(Point(0, 0), "Point2D(Integer(0), Integer(0))")
     sT(Ellipse(Point(0, 0), 5, 1),
@@ -206,12 +207,14 @@ def test_Mul():
     assert srepr(3*x**3*y, order='old') == "Mul(Integer(3), Symbol('y'), Pow(Symbol('x'), Integer(3)))"
     assert srepr(sympify('(x+4)*2*x*7', evaluate=False), order='none') == "Mul(Add(Symbol('x'), Integer(4)), Integer(2), Symbol('x'), Integer(7))"
 
+
 def test_AlgebraicNumber():
     a = AlgebraicNumber(sqrt(2))
     sT(a, "AlgebraicNumber(Pow(Integer(2), Rational(1, 2)), [Integer(1), Integer(0)])")
     a = AlgebraicNumber(root(-2, 3))
     sT(a, "AlgebraicNumber(Pow(Integer(-2), Rational(1, 3)), [Integer(1), Integer(0)])")
 
+
 def test_PolyRing():
     assert srepr(ring("x", ZZ, lex)[0]) == "PolyRing((Symbol('x'),), ZZ, lex)"
     assert srepr(ring("x,y", QQ, grlex)[0]) == "PolyRing((Symbol('x'), Symbol('y')), QQ, grlex)"
@@ -233,6 +236,7 @@ def test_FracElement():
     F, x, y = field("x,y", ZZ)
     assert srepr((3*x**2*y + 1)/(x - y**2)) == "FracElement(FracField((Symbol('x'), Symbol('y')), ZZ, lex), [((2, 1), 3), ((0, 0), 1)], [((1, 0), 1), ((0, 2), -1)])"
 
+
 def test_FractionField():
     assert srepr(QQ.frac_field(x)) == \
         "FractionField(FracField((Symbol('x'),), QQ, lex))"
@@ -308,6 +312,7 @@ def test_Permutation():
     import_stmt = "from sympy.combinatorics import Permutation"
     sT(Permutation(1, 2), "Permutation(1, 2)", import_stmt)
 
+
 def test_diffgeom():
     from sympy.diffgeom import Manifold, Patch, CoordSystem, BaseScalarField
     m = Manifold('M', 2)
diff --git a/sympy/printing/tests/test_str.py b/sympy/printing/tests/test_str.py
index 082b7844d4..287c2145f3 100644
--- a/sympy/printing/tests/test_str.py
+++ b/sympy/printing/tests/test_str.py
@@ -867,7 +867,7 @@ def test_MatrixExpressions():
 
     lamda = Lambda(x, 1/x)
     expr = (n*X).applyfunc(lamda)
-    assert str(expr) == 'Lambda(_d, 1/_d).(n*X)'
+    assert str(expr) == 'Lambda(x, 1/x).(n*X)'
 
 
 def test_Subs_printing():
diff --git a/sympy/sets/tests/test_fancysets.py b/sympy/sets/tests/test_fancysets.py
index f89c571742..a01f3101fa 100644
--- a/sympy/sets/tests/test_fancysets.py
+++ b/sympy/sets/tests/test_fancysets.py
@@ -562,8 +562,8 @@ def test_infinitely_indexed_set_1():
             imageset(Lambda(n, 2*n + 1), S.Integers)) is S.EmptySet
 
     assert imageset(Lambda(m, 2*m), S.Integers).intersect(
-                imageset(Lambda(n, 3*n), S.Integers)) == \
-            ImageSet(Lambda(t, 6*t), S.Integers)
+                imageset(Lambda(n, 3*n), S.Integers)).dummy_eq(
+            ImageSet(Lambda(t, 6*t), S.Integers))
 
     assert imageset(x, x/2 + Rational(1, 3), S.Integers).intersect(S.Integers) is S.EmptySet
     assert imageset(x, x/2 + S.Half, S.Integers).intersect(S.Integers) is S.Integers
@@ -649,37 +649,38 @@ def test_imageset_intersect_diophantine():
             ImageSet(Lambda(n, -(n - 3)**2), S.Integers)) == FiniteSet(0)
     # Single parametric solution for diophantine solution:
     assert ImageSet(Lambda(n, n**2 + 5), S.Integers).intersect(
-            ImageSet(Lambda(m, 2*m), S.Integers)) == ImageSet(
-            Lambda(n, 4*n**2 + 4*n + 6), S.Integers)
+            ImageSet(Lambda(m, 2*m), S.Integers)).dummy_eq(ImageSet(
+            Lambda(n, 4*n**2 + 4*n + 6), S.Integers))
     # 4 non-parametric solution couples for dioph. equation:
     assert ImageSet(Lambda(n, n**2 - 9), S.Integers).intersect(
             ImageSet(Lambda(m, -m**2), S.Integers)) == FiniteSet(-9, 0)
     # Double parametric solution for diophantine solution:
     assert ImageSet(Lambda(m, m**2 + 40), S.Integers).intersect(
-            ImageSet(Lambda(n, 41*n), S.Integers)) == Intersection(
+            ImageSet(Lambda(n, 41*n), S.Integers)).dummy_eq(Intersection(
             ImageSet(Lambda(m, m**2 + 40), S.Integers),
-            ImageSet(Lambda(n, 41*n), S.Integers))
+            ImageSet(Lambda(n, 41*n), S.Integers)))
     # Check that diophantine returns *all* (8) solutions (permute=True)
     assert ImageSet(Lambda(n, n**4 - 2**4), S.Integers).intersect(
             ImageSet(Lambda(m, -m**4 + 3**4), S.Integers)) == FiniteSet(0, 65)
     assert ImageSet(Lambda(n, pi/12 + n*5*pi/12), S.Integers).intersect(
-            ImageSet(Lambda(n, 7*pi/12 + n*11*pi/12), S.Integers)) == ImageSet(
-            Lambda(n, 55*pi*n/12 + 17*pi/4), S.Integers)
+            ImageSet(Lambda(n, 7*pi/12 + n*11*pi/12), S.Integers)).dummy_eq(ImageSet(
+            Lambda(n, 55*pi*n/12 + 17*pi/4), S.Integers))
     # TypeError raised by diophantine (#18081)
-    assert ImageSet(Lambda(n, n*log(2)), S.Integers).intersection(S.Integers) \
-            == Intersection(ImageSet(Lambda(n, n*log(2)), S.Integers), S.Integers)
+    assert ImageSet(Lambda(n, n*log(2)), S.Integers).intersection(
+        S.Integers).dummy_eq(Intersection(ImageSet(
+        Lambda(n, n*log(2)), S.Integers), S.Integers))
     # NotImplementedError raised by diophantine (no solver for cubic_thue)
     assert ImageSet(Lambda(n, n**3 + 1), S.Integers).intersect(
-            ImageSet(Lambda(n, n**3), S.Integers)) == Intersection(
+            ImageSet(Lambda(n, n**3), S.Integers)).dummy_eq(Intersection(
             ImageSet(Lambda(n, n**3 + 1), S.Integers),
-            ImageSet(Lambda(n, n**3), S.Integers))
+            ImageSet(Lambda(n, n**3), S.Integers)))
 
 
 def test_infinitely_indexed_set_3():
     from sympy.abc import n, m, t
     assert imageset(Lambda(m, 2*pi*m), S.Integers).intersect(
-            imageset(Lambda(n, 3*pi*n), S.Integers)) == \
-        ImageSet(Lambda(t, 6*pi*t), S.Integers)
+            imageset(Lambda(n, 3*pi*n), S.Integers)).dummy_eq(
+        ImageSet(Lambda(t, 6*pi*t), S.Integers))
     assert imageset(Lambda(n, 2*n + 1), S.Integers) == \
         imageset(Lambda(n, 2*n - 1), S.Integers)
     assert imageset(Lambda(n, 3*n + 2), S.Integers) == \
diff --git a/sympy/sets/tests/test_setexpr.py b/sympy/sets/tests/test_setexpr.py
index b079429e38..e6490e7757 100644
--- a/sympy/sets/tests/test_setexpr.py
+++ b/sympy/sets/tests/test_setexpr.py
@@ -297,13 +297,17 @@ def test_SetExpr_Interval_pow():
 
 def test_SetExpr_Integers():
     assert SetExpr(S.Integers) + 1 == SetExpr(S.Integers)
-    assert SetExpr(S.Integers) + I == SetExpr(ImageSet(Lambda(_d, _d + I), S.Integers))
+    assert (SetExpr(S.Integers) + I).dummy_eq(
+        SetExpr(ImageSet(Lambda(_d, _d + I), S.Integers)))
     assert SetExpr(S.Integers)*(-1) == SetExpr(S.Integers)
-    assert SetExpr(S.Integers)*2 == SetExpr(ImageSet(Lambda(_d, 2*_d), S.Integers))
-    assert SetExpr(S.Integers)*I == SetExpr(ImageSet(Lambda(_d, I*_d), S.Integers))
+    assert (SetExpr(S.Integers)*2).dummy_eq(
+        SetExpr(ImageSet(Lambda(_d, 2*_d), S.Integers)))
+    assert (SetExpr(S.Integers)*I).dummy_eq(
+        SetExpr(ImageSet(Lambda(_d, I*_d), S.Integers)))
     # issue #18050:
-    assert SetExpr(S.Integers)._eval_func(Lambda(x, I*x + 1)) == SetExpr(
-            ImageSet(Lambda(_d, I*_d + 1), S.Integers))
+    assert SetExpr(S.Integers)._eval_func(Lambda(x, I*x + 1)).dummy_eq(
+        SetExpr(ImageSet(Lambda(_d, I*_d + 1), S.Integers)))
     # needs improvement:
-    assert SetExpr(S.Integers)*I + 1 == SetExpr(
-            ImageSet(Lambda(x, x + 1), ImageSet(Lambda(_d, _d*I), S.Integers)))
+    assert (SetExpr(S.Integers)*I + 1).dummy_eq(
+        SetExpr(ImageSet(Lambda(x, x + 1),
+        ImageSet(Lambda(_d, _d*I), S.Integers))))
diff --git a/sympy/sets/tests/test_sets.py b/sympy/sets/tests/test_sets.py
index aa4a373577..61c058c05e 100644
--- a/sympy/sets/tests/test_sets.py
+++ b/sympy/sets/tests/test_sets.py
@@ -47,8 +47,10 @@ def f(x):
     assert (str(imageset(lambda x: x + clash, Interval(-2, 1)).lamda.expr)
         in ('_x + x', 'x + _x'))
     x1, x2 = symbols("x1, x2")
-    assert imageset(lambda x, y: Add(x, y), Interval(1, 2), Interval(2, 3)) == \
-        ImageSet(Lambda((x1, x2), x1+x2), Interval(1, 2), Interval(2, 3))
+    assert imageset(lambda x, y:
+        Add(x, y), Interval(1, 2), Interval(2, 3)).dummy_eq(
+        ImageSet(Lambda((x1, x2), x1 + x2),
+        Interval(1, 2), Interval(2, 3)))
 
 
 def test_is_empty():
diff --git a/sympy/solvers/tests/test_solveset.py b/sympy/solvers/tests/test_solveset.py
index 133634fd96..82f4c58696 100644
--- a/sympy/solvers/tests/test_solveset.py
+++ b/sympy/solvers/tests/test_solveset.py
@@ -30,8 +30,7 @@
 from sympy.tensor.indexed import Indexed
 from sympy.utilities.iterables import numbered_symbols
 
-from sympy.testing.pytest import (XFAIL, raises, skip, slow, SKIP,
-    nocache_fail)
+from sympy.testing.pytest import (XFAIL, raises, skip, slow, SKIP)
 from sympy.testing.randtest import verify_numerically as tn
 from sympy.physics.units import cm
 
@@ -47,6 +46,12 @@
     t, w, x, y, z)
 
 
+def dumeq(i, j):
+    if type(i) in (list, tuple):
+        return all(dumeq(i, j) for i, j in zip(i, j))
+    return i == j or i.dummy_eq(j)
+
+
 def test_invert_real():
     x = Symbol('x', real=True)
 
@@ -92,48 +97,48 @@ def ireal(x, s=S.Reals):
     base_values =  FiniteSet(y - 1, -y - 1)
     assert invert_real(Abs(x**31 + x + 1), y, x) == (lhs, base_values)
 
-    assert invert_real(sin(x), y, x) == \
-        (x, imageset(Lambda(n, n*pi + (-1)**n*asin(y)), S.Integers))
+    assert dumeq(invert_real(sin(x), y, x),
+        (x, imageset(Lambda(n, n*pi + (-1)**n*asin(y)), S.Integers)))
 
-    assert invert_real(sin(exp(x)), y, x) == \
-        (x, imageset(Lambda(n, log((-1)**n*asin(y) + n*pi)), S.Integers))
+    assert dumeq(invert_real(sin(exp(x)), y, x),
+        (x, imageset(Lambda(n, log((-1)**n*asin(y) + n*pi)), S.Integers)))
 
-    assert invert_real(csc(x), y, x) == \
-        (x, imageset(Lambda(n, n*pi + (-1)**n*acsc(y)), S.Integers))
+    assert dumeq(invert_real(csc(x), y, x),
+        (x, imageset(Lambda(n, n*pi + (-1)**n*acsc(y)), S.Integers)))
 
-    assert invert_real(csc(exp(x)), y, x) == \
-        (x, imageset(Lambda(n, log((-1)**n*acsc(y) + n*pi)), S.Integers))
+    assert dumeq(invert_real(csc(exp(x)), y, x),
+        (x, imageset(Lambda(n, log((-1)**n*acsc(y) + n*pi)), S.Integers)))
 
-    assert invert_real(cos(x), y, x) == \
+    assert dumeq(invert_real(cos(x), y, x),
         (x, Union(imageset(Lambda(n, 2*n*pi + acos(y)), S.Integers), \
-                imageset(Lambda(n, 2*n*pi - acos(y)), S.Integers)))
+                imageset(Lambda(n, 2*n*pi - acos(y)), S.Integers))))
 
-    assert invert_real(cos(exp(x)), y, x) == \
+    assert dumeq(invert_real(cos(exp(x)), y, x),
         (x, Union(imageset(Lambda(n, log(2*n*pi + acos(y))), S.Integers), \
-                imageset(Lambda(n, log(2*n*pi - acos(y))), S.Integers)))
+                imageset(Lambda(n, log(2*n*pi - acos(y))), S.Integers))))
 
-    assert invert_real(sec(x), y, x) == \
+    assert dumeq(invert_real(sec(x), y, x),
         (x, Union(imageset(Lambda(n, 2*n*pi + asec(y)), S.Integers), \
-                imageset(Lambda(n, 2*n*pi - asec(y)), S.Integers)))
+                imageset(Lambda(n, 2*n*pi - asec(y)), S.Integers))))
 
-    assert invert_real(sec(exp(x)), y, x) == \
+    assert dumeq(invert_real(sec(exp(x)), y, x),
         (x, Union(imageset(Lambda(n, log(2*n*pi + asec(y))), S.Integers), \
-                imageset(Lambda(n, log(2*n*pi - asec(y))), S.Integers)))
+                imageset(Lambda(n, log(2*n*pi - asec(y))), S.Integers))))
 
-    assert invert_real(tan(x), y, x) == \
-        (x, imageset(Lambda(n, n*pi + atan(y)), S.Integers))
+    assert dumeq(invert_real(tan(x), y, x),
+        (x, imageset(Lambda(n, n*pi + atan(y)), S.Integers)))
 
-    assert invert_real(tan(exp(x)), y, x) == \
-        (x, imageset(Lambda(n, log(n*pi + atan(y))), S.Integers))
+    assert dumeq(invert_real(tan(exp(x)), y, x),
+        (x, imageset(Lambda(n, log(n*pi + atan(y))), S.Integers)))
 
-    assert invert_real(cot(x), y, x) == \
-        (x, imageset(Lambda(n, n*pi + acot(y)), S.Integers))
+    assert dumeq(invert_real(cot(x), y, x),
+        (x, imageset(Lambda(n, n*pi + acot(y)), S.Integers)))
 
-    assert invert_real(cot(exp(x)), y, x) == \
-        (x, imageset(Lambda(n, log(n*pi + acot(y))), S.Integers))
+    assert dumeq(invert_real(cot(exp(x)), y, x),
+        (x, imageset(Lambda(n, log(n*pi + acot(y))), S.Integers)))
 
-    assert invert_real(tan(tan(x)), y, x) == \
-        (tan(x), imageset(Lambda(n, n*pi + atan(y)), S.Integers))
+    assert dumeq(invert_real(tan(tan(x)), y, x),
+        (tan(x), imageset(Lambda(n, n*pi + atan(y)), S.Integers)))
 
     x = Symbol('x', positive=True)
     assert invert_real(x**pi, y, x) == (x, FiniteSet(y**(1/pi)))
@@ -143,8 +148,8 @@ def test_invert_complex():
     assert invert_complex(x + 3, y, x) == (x, FiniteSet(y - 3))
     assert invert_complex(x*3, y, x) == (x, FiniteSet(y / 3))
 
-    assert invert_complex(exp(x), y, x) == \
-        (x, imageset(Lambda(n, I*(2*pi*n + arg(y)) + log(Abs(y))), S.Integers))
+    assert dumeq(invert_complex(exp(x), y, x),
+        (x, imageset(Lambda(n, I*(2*pi*n + arg(y)) + log(Abs(y))), S.Integers)))
 
     assert invert_complex(log(x), y, x) == (x, FiniteSet(exp(y)))
 
@@ -622,16 +627,16 @@ def test_solve_abs():
     eqab = eq.subs(reps)
     for si in sol.subs(reps):
         assert not eqab.subs(x, si)
-    assert solveset(Eq(sin(Abs(x)), 1), x, domain=S.Reals) == Union(
+    assert dumeq(solveset(Eq(sin(Abs(x)), 1), x, domain=S.Reals), Union(
         Intersection(Interval(0, oo),
             ImageSet(Lambda(n, (-1)**n*pi/2 + n*pi), S.Integers)),
         Intersection(Interval(-oo, 0),
-            ImageSet(Lambda(n, n*pi - (-1)**(-n)*pi/2), S.Integers)))
+            ImageSet(Lambda(n, n*pi - (-1)**(-n)*pi/2), S.Integers))))
 
 
 def test_issue_9824():
-    assert solveset(sin(x)**2 - 2*sin(x) + 1, x) == ImageSet(Lambda(n, 2*n*pi + pi/2), S.Integers)
-    assert solveset(cos(x)**2 - 2*cos(x) + 1, x) == ImageSet(Lambda(n, 2*n*pi), S.Integers)
+    assert dumeq(solveset(sin(x)**2 - 2*sin(x) + 1, x), ImageSet(Lambda(n, 2*n*pi + pi/2), S.Integers))
+    assert dumeq(solveset(cos(x)**2 - 2*cos(x) + 1, x), ImageSet(Lambda(n, 2*n*pi), S.Integers))
 
 
 def test_issue_9565():
@@ -735,13 +740,13 @@ def test_solve_quintics():
 
 def test_solveset_complex_exp():
     from sympy.abc import x, n
-    assert solveset_complex(exp(x) - 1, x) == \
-        imageset(Lambda(n, I*2*n*pi), S.Integers)
-    assert solveset_complex(exp(x) - I, x) == \
-        imageset(Lambda(n, I*(2*n*pi + pi/2)), S.Integers)
+    assert dumeq(solveset_complex(exp(x) - 1, x),
+        imageset(Lambda(n, I*2*n*pi), S.Integers))
+    assert dumeq(solveset_complex(exp(x) - I, x),
+        imageset(Lambda(n, I*(2*n*pi + pi/2)), S.Integers))
     assert solveset_complex(1/exp(x), x) == S.EmptySet
-    assert solveset_complex(sinh(x).rewrite(exp), x) == \
-        imageset(Lambda(n, n*pi*I), S.Integers)
+    assert dumeq(solveset_complex(sinh(x).rewrite(exp), x),
+        imageset(Lambda(n, n*pi*I), S.Integers))
 
 
 def test_solveset_real_exp():
@@ -774,62 +779,60 @@ def test_solve_complex_sqrt():
 
 def test_solveset_complex_tan():
     s = solveset_complex(tan(x).rewrite(exp), x)
-    assert s == imageset(Lambda(n, pi*n), S.Integers) - \
-        imageset(Lambda(n, pi*n + pi/2), S.Integers)
+    assert dumeq(s, imageset(Lambda(n, pi*n), S.Integers) - \
+        imageset(Lambda(n, pi*n + pi/2), S.Integers))
 
 
-@nocache_fail
 def test_solve_trig():
     from sympy.abc import n
-    assert solveset_real(sin(x), x) == \
+    assert dumeq(solveset_real(sin(x), x),
         Union(imageset(Lambda(n, 2*pi*n), S.Integers),
-              imageset(Lambda(n, 2*pi*n + pi), S.Integers))
+              imageset(Lambda(n, 2*pi*n + pi), S.Integers)))
 
-    assert solveset_real(sin(x) - 1, x) == \
-        imageset(Lambda(n, 2*pi*n + pi/2), S.Integers)
+    assert dumeq(solveset_real(sin(x) - 1, x),
+        imageset(Lambda(n, 2*pi*n + pi/2), S.Integers))
 
-    assert solveset_real(cos(x), x) == \
+    assert dumeq(solveset_real(cos(x), x),
         Union(imageset(Lambda(n, 2*pi*n + pi/2), S.Integers),
-              imageset(Lambda(n, 2*pi*n + pi*Rational(3, 2)), S.Integers))
+              imageset(Lambda(n, 2*pi*n + pi*Rational(3, 2)), S.Integers)))
 
-    assert solveset_real(sin(x) + cos(x), x) == \
+    assert dumeq(solveset_real(sin(x) + cos(x), x),
         Union(imageset(Lambda(n, 2*n*pi + pi*Rational(3, 4)), S.Integers),
-              imageset(Lambda(n, 2*n*pi + pi*Rational(7, 4)), S.Integers))
+              imageset(Lambda(n, 2*n*pi + pi*Rational(7, 4)), S.Integers)))
 
     assert solveset_real(sin(x)**2 + cos(x)**2, x) == S.EmptySet
 
-    # This fails when running with the cache off:
-    assert solveset_complex(cos(x) - S.Half, x) == \
+    assert dumeq(solveset_complex(cos(x) - S.Half, x),
         Union(imageset(Lambda(n, 2*n*pi + pi*Rational(5, 3)), S.Integers),
-              imageset(Lambda(n, 2*n*pi + pi/3), S.Integers))
+              imageset(Lambda(n, 2*n*pi + pi/3), S.Integers)))
 
-    assert solveset(sin(y + a) - sin(y), a, domain=S.Reals) == \
+    assert dumeq(solveset(sin(y + a) - sin(y), a, domain=S.Reals),
         Union(ImageSet(Lambda(n, 2*n*pi), S.Integers),
         Intersection(ImageSet(Lambda(n, -I*(I*(
         2*n*pi + arg(-exp(-2*I*y))) +
-        2*im(y))), S.Integers), S.Reals))
+        2*im(y))), S.Integers), S.Reals)))
 
-    assert solveset_real(sin(2*x)*cos(x) + cos(2*x)*sin(x)-1, x) == \
-                            ImageSet(Lambda(n, n*pi*Rational(2, 3) + pi/6), S.Integers)
+    assert dumeq(solveset_real(sin(2*x)*cos(x) + cos(2*x)*sin(x)-1, x),
+        ImageSet(Lambda(n, n*pi*Rational(2, 3) + pi/6), S.Integers))
 
     # Tests for _solve_trig2() function
-    assert solveset_real(2*cos(x)*cos(2*x) - 1, x) == \
+    assert dumeq(solveset_real(2*cos(x)*cos(2*x) - 1, x),
           Union(ImageSet(Lambda(n, 2*n*pi + 2*atan(sqrt(-2*2**Rational(1, 3)*(67 +
                   9*sqrt(57))**Rational(2, 3) + 8*2**Rational(2, 3) + 11*(67 +
                   9*sqrt(57))**Rational(1, 3))/(3*(67 + 9*sqrt(57))**Rational(1, 6)))), S.Integers),
                   ImageSet(Lambda(n, 2*n*pi - 2*atan(sqrt(-2*2**Rational(1, 3)*(67 +
                   9*sqrt(57))**Rational(2, 3) + 8*2**Rational(2, 3) + 11*(67 +
                   9*sqrt(57))**Rational(1, 3))/(3*(67 + 9*sqrt(57))**Rational(1, 6))) +
-                  2*pi), S.Integers))
+                  2*pi), S.Integers)))
 
-    assert solveset_real(2*tan(x)*sin(x) + 1, x) == Union(
+    assert dumeq(solveset_real(2*tan(x)*sin(x) + 1, x), Union(
         ImageSet(Lambda(n, 2*n*pi + atan(sqrt(2)*sqrt(-1 +sqrt(17))/
             (1 - sqrt(17))) + pi), S.Integers),
         ImageSet(Lambda(n, 2*n*pi - atan(sqrt(2)*sqrt(-1 + sqrt(17))/
-            (1 - sqrt(17))) + pi), S.Integers))
+            (1 - sqrt(17))) + pi), S.Integers)))
 
-    assert solveset_real(cos(2*x)*cos(4*x) - 1, x) == \
-                            ImageSet(Lambda(n, n*pi), S.Integers)
+    assert dumeq(solveset_real(cos(2*x)*cos(4*x) - 1, x),
+                            ImageSet(Lambda(n, n*pi), S.Integers))
 
 
 def test_solve_hyperbolic():
@@ -852,19 +855,19 @@ def test_solve_hyperbolic():
         log(sqrt(4 + sqrt(17))))
     assert solveset_real(sinh(x) + tanh(x) - 1, x) == FiniteSet(
         log(sqrt(2)/2 + sqrt(-S(1)/2 + sqrt(2))))
-    assert solveset_complex(sinh(x) - I/2, x) == Union(
+    assert dumeq(solveset_complex(sinh(x) - I/2, x), Union(
         ImageSet(Lambda(n, I*(2*n*pi + 5*pi/6)), S.Integers),
-        ImageSet(Lambda(n, I*(2*n*pi + pi/6)), S.Integers))
-    assert solveset_complex(sinh(x) + sech(x), x) == Union(
+        ImageSet(Lambda(n, I*(2*n*pi + pi/6)), S.Integers)))
+    assert dumeq(solveset_complex(sinh(x) + sech(x), x), Union(
         ImageSet(Lambda(n, 2*n*I*pi + log(sqrt(-2 + sqrt(5)))), S.Integers),
         ImageSet(Lambda(n, I*(2*n*pi + pi/2) + log(sqrt(2 + sqrt(5)))), S.Integers),
         ImageSet(Lambda(n, I*(2*n*pi + pi) + log(sqrt(-2 + sqrt(5)))), S.Integers),
-        ImageSet(Lambda(n, I*(2*n*pi - pi/2) + log(sqrt(2 + sqrt(5)))), S.Integers))
+        ImageSet(Lambda(n, I*(2*n*pi - pi/2) + log(sqrt(2 + sqrt(5)))), S.Integers)))
     # issues #9606 / #9531:
     assert solveset(sinh(x), x, S.Reals) == FiniteSet(0)
-    assert solveset(sinh(x), x, S.Complexes) == Union(
+    assert dumeq(solveset(sinh(x), x, S.Complexes), Union(
         ImageSet(Lambda(n, I*(2*n*pi + pi)), S.Integers),
-        ImageSet(Lambda(n, 2*n*I*pi), S.Integers))
+        ImageSet(Lambda(n, 2*n*I*pi), S.Integers)))
 
 
 def test_solve_invalid_sol():
@@ -875,14 +878,14 @@ def test_solve_invalid_sol():
 @XFAIL
 def test_solve_trig_simplified():
     from sympy.abc import n
-    assert solveset_real(sin(x), x) == \
-        imageset(Lambda(n, n*pi), S.Integers)
+    assert dumeq(solveset_real(sin(x), x),
+        imageset(Lambda(n, n*pi), S.Integers))
 
-    assert solveset_real(cos(x), x) == \
-        imageset(Lambda(n, n*pi + pi/2), S.Integers)
+    assert dumeq(solveset_real(cos(x), x),
+        imageset(Lambda(n, n*pi + pi/2), S.Integers))
 
-    assert solveset_real(cos(x) + sin(x), x) == \
-        imageset(Lambda(n, n*pi - pi/4), S.Integers)
+    assert dumeq(solveset_real(cos(x) + sin(x), x),
+        imageset(Lambda(n, n*pi - pi/4), S.Integers))
 
 
 @XFAIL
@@ -1005,9 +1008,9 @@ def test_solveset():
     assert solveset(x - 1 >= 0, x, S.Reals) == Interval(1, oo)
     assert solveset(exp(x) - 1 >= 0, x, S.Reals) == Interval(0, oo)
 
-    assert solveset(exp(x) - 1, x) == imageset(Lambda(n, 2*I*pi*n), S.Integers)
-    assert solveset(Eq(exp(x), 1), x) == imageset(Lambda(n, 2*I*pi*n),
-                                                  S.Integers)
+    assert dumeq(solveset(exp(x) - 1, x), imageset(Lambda(n, 2*I*pi*n), S.Integers))
+    assert dumeq(solveset(Eq(exp(x), 1), x), imageset(Lambda(n, 2*I*pi*n),
+                                                  S.Integers))
     # issue 13825
     assert solveset(x**2 + f(0) + 1, x) == {-sqrt(-f(0) - 1), sqrt(-f(0) - 1)}
 
@@ -1026,10 +1029,10 @@ def test__solveset_multi():
     assert _solveset_multi([x+y, x+1], [y, x], [Reals, Reals]) == FiniteSet((1, -1))
     assert _solveset_multi([x+y, x-y-1], [x, y], [Reals, Reals]) == FiniteSet((S(1)/2, -S(1)/2))
     assert _solveset_multi([x-1, y-2], [x, y], [Reals, Reals]) == FiniteSet((1, 2))
-    #assert _solveset_multi([x+y], [x, y], [Reals, Reals]) == ImageSet(Lambda(x, (x, -x)), Reals)
-    assert _solveset_multi([x+y], [x, y], [Reals, Reals]) == Union(
+    # assert dumeq(_solveset_multi([x+y], [x, y], [Reals, Reals]), ImageSet(Lambda(x, (x, -x)), Reals))
+    assert dumeq(_solveset_multi([x+y], [x, y], [Reals, Reals]), Union(
             ImageSet(Lambda(((x,),), (x, -x)), ProductSet(Reals)),
-            ImageSet(Lambda(((y,),), (-y, y)), ProductSet(Reals)))
+            ImageSet(Lambda(((y,),), (-y, y)), ProductSet(Reals))))
     assert _solveset_multi([x+y, x+y+1], [x, y], [Reals, Reals]) == S.EmptySet
     assert _solveset_multi([x+y, x-y, x-1], [x, y], [Reals, Reals]) == S.EmptySet
     assert _solveset_multi([x+y, x-y, x-1], [y, x], [Reals, Reals]) == S.EmptySet
@@ -1045,21 +1048,21 @@ def test__solveset_multi():
     assert _solveset_multi([x**2-1, y], [x, y], [Reals, Reals]) == FiniteSet((1, 0), (-1, 0))
     #assert _solveset_multi([x**2-y**2], [x, y], [Reals, Reals]) == Union(
     #        ImageSet(Lambda(x, (x, -x)), Reals), ImageSet(Lambda(x, (x, x)), Reals))
-    assert _solveset_multi([x**2-y**2], [x, y], [Reals, Reals]) == Union(
+    assert dumeq(_solveset_multi([x**2-y**2], [x, y], [Reals, Reals]), Union(
             ImageSet(Lambda(((x,),), (x, -Abs(x))), ProductSet(Reals)),
             ImageSet(Lambda(((x,),), (x, Abs(x))), ProductSet(Reals)),
             ImageSet(Lambda(((y,),), (-Abs(y), y)), ProductSet(Reals)),
-            ImageSet(Lambda(((y,),), (Abs(y), y)), ProductSet(Reals)))
+            ImageSet(Lambda(((y,),), (Abs(y), y)), ProductSet(Reals))))
     assert _solveset_multi([r*cos(theta)-1, r*sin(theta)], [theta, r],
             [Interval(0, pi), Interval(-1, 1)]) == FiniteSet((0, 1), (pi, -1))
     assert _solveset_multi([r*cos(theta)-1, r*sin(theta)], [r, theta],
             [Interval(0, 1), Interval(0, pi)]) == FiniteSet((1, 0))
     #assert _solveset_multi([r*cos(theta)-r, r*sin(theta)], [r, theta],
     #        [Interval(0, 1), Interval(0, pi)]) == ?
-    assert _solveset_multi([r*cos(theta)-r, r*sin(theta)], [r, theta],
-            [Interval(0, 1), Interval(0, pi)]) == Union(
+    assert dumeq(_solveset_multi([r*cos(theta)-r, r*sin(theta)], [r, theta],
+            [Interval(0, 1), Interval(0, pi)]), Union(
             ImageSet(Lambda(((r,),), (r, 0)), ImageSet(Lambda(r, (r,)), Interval(0, 1))),
-            ImageSet(Lambda(((theta,),), (0, theta)), ImageSet(Lambda(theta, (theta,)), Interval(0, pi))))
+            ImageSet(Lambda(((theta,),), (0, theta)), ImageSet(Lambda(theta, (theta,)), Interval(0, pi)))))
 
 
 def test_conditionset():
@@ -1069,8 +1072,8 @@ def test_conditionset():
     assert solveset(Eq(x**2 + x*sin(x), 1), x, domain=S.Reals
         ) == ConditionSet(x, Eq(x**2 + x*sin(x) - 1, 0), S.Reals)
 
-    assert solveset(Eq(-I*(exp(I*x) - exp(-I*x))/2, 1), x
-        ) == imageset(Lambda(n, 2*n*pi + pi/2), S.Integers)
+    assert dumeq(solveset(Eq(-I*(exp(I*x) - exp(-I*x))/2, 1), x
+        ), imageset(Lambda(n, 2*n*pi + pi/2), S.Integers))
 
     assert solveset(x + sin(x) > 1, x, domain=S.Reals
         ) == ConditionSet(x, x + sin(x) > 1, S.Reals)
@@ -1296,9 +1299,9 @@ def test_solve_decomposition():
     s5 = ImageSet(Lambda(n, 2*n*pi - 1 + pi), S.Integers)
 
     assert solve_decomposition(f1, x, S.Reals) == FiniteSet(0, log(2), log(3))
-    assert solve_decomposition(f2, x, S.Reals) == s3
-    assert solve_decomposition(f3, x, S.Reals) == Union(s1, s2, s3)
-    assert solve_decomposition(f4, x, S.Reals) == Union(s4, s5)
+    assert dumeq(solve_decomposition(f2, x, S.Reals), s3)
+    assert dumeq(solve_decomposition(f3, x, S.Reals), Union(s1, s2, s3))
+    assert dumeq(solve_decomposition(f4, x, S.Reals), Union(s4, s5))
     assert solve_decomposition(f5, x, S.Reals) == FiniteSet(-2)
     assert solve_decomposition(f6, x, S.Reals) == S.EmptySet
     assert solve_decomposition(f7, x, S.Reals) == S.EmptySet
@@ -1313,7 +1316,7 @@ def test_nonlinsolve_basic():
     assert nonlinsolve([x],[x, y]) == FiniteSet((0, y))
     assert nonlinsolve(system, [y]) == FiniteSet((x + 5,))
     soln = (ImageSet(Lambda(n, 2*n*pi + pi/2), S.Integers),)
-    assert nonlinsolve([sin(x) - 1], [x]) == FiniteSet(tuple(soln))
+    assert dumeq(nonlinsolve([sin(x) - 1], [x]), FiniteSet(tuple(soln)))
     assert nonlinsolve([x**2 - 1], [x]) == FiniteSet((-1,), (1,))
 
     soln = FiniteSet((y, y))
@@ -1349,7 +1352,7 @@ def test_trig_system():
     assert nonlinsolve([sin(x) - 1, cos(x) -1 ], x) == S.EmptySet
     soln1 = (ImageSet(Lambda(n, 2*n*pi + pi/2), S.Integers),)
     soln = FiniteSet(soln1)
-    assert nonlinsolve([sin(x) - 1, cos(x)], x) == soln
+    assert dumeq(nonlinsolve([sin(x) - 1, cos(x)], x), soln)
 
 
 @XFAIL
@@ -1364,7 +1367,7 @@ def test_trig_system_fail():
         ImageSet(Lambda(n, n*pi+ pi/2), S.Integers))
     soln_2 = FiniteSet(soln_2)
     soln = soln_1 + soln_2
-    assert nonlinsolve(sys, [x, y]) == soln
+    assert dumeq(nonlinsolve(sys, [x, y]), soln)
 
     # Add more cases from here
     # http://www.vitutor.com/geometry/trigonometry/equations_systems.html#uno
@@ -1373,7 +1376,7 @@ def test_trig_system_fail():
         ImageSet(Lambda(n, 2*n*pi + pi*Rational(2, 3)), S.Integers))
     soln_y = Union(ImageSet(Lambda(n, 2*n*pi + pi/6), S.Integers),
         ImageSet(Lambda(n, 2*n*pi + pi*Rational(5, 6)), S.Integers))
-    assert nonlinsolve(sys, [x, y]) == FiniteSet((soln_x, soln_y))
+    assert dumeq(nonlinsolve(sys, [x, y]), FiniteSet((soln_x, soln_y)))
 
 
 def test_nonlinsolve_positive_dimensional():
@@ -1440,21 +1443,21 @@ def test_nonlinsolve_using_substitution():
 
 def test_nonlinsolve_complex():
     n = Dummy('n')
-    assert nonlinsolve([exp(x) - sin(y), 1/y - 3], [x, y]) == {
-        (ImageSet(Lambda(n, 2*n*I*pi + log(sin(Rational(1, 3)))), S.Integers), Rational(1, 3))}
+    assert dumeq(nonlinsolve([exp(x) - sin(y), 1/y - 3], [x, y]), {
+        (ImageSet(Lambda(n, 2*n*I*pi + log(sin(Rational(1, 3)))), S.Integers), Rational(1, 3))})
 
     system = [exp(x) - sin(y), 1/exp(y) - 3]
-    assert nonlinsolve(system, [x, y]) == {
+    assert dumeq(nonlinsolve(system, [x, y]), {
         (ImageSet(Lambda(n, I*(2*n*pi + pi)
                          + log(sin(log(3)))), S.Integers), -log(3)),
         (ImageSet(Lambda(n, I*(2*n*pi + arg(sin(2*n*I*pi - log(3))))
                          + log(Abs(sin(2*n*I*pi - log(3))))), S.Integers),
-        ImageSet(Lambda(n, 2*n*I*pi - log(3)), S.Integers))}
+        ImageSet(Lambda(n, 2*n*I*pi - log(3)), S.Integers))})
 
     system = [exp(x) - sin(y), y**2 - 4]
-    assert nonlinsolve(system, [x, y]) == {
+    assert dumeq(nonlinsolve(system, [x, y]), {
         (ImageSet(Lambda(n, I*(2*n*pi + pi) + log(sin(2))), S.Integers), -2),
-        (ImageSet(Lambda(n, 2*n*I*pi + log(sin(2))), S.Integers), 2)}
+        (ImageSet(Lambda(n, 2*n*I*pi + log(sin(2))), S.Integers), 2)})
 
 
 @XFAIL
@@ -1491,7 +1494,7 @@ def test_issue_5132_1():
                                             (s_complex_y, s_complex_z_2)
                                         )
     soln = soln_real + soln_complex
-    assert nonlinsolve(eqs, [y, z]) == soln
+    assert dumeq(nonlinsolve(eqs, [y, z]), soln)
 
 
 def test_issue_5132_2():
@@ -1504,7 +1507,7 @@ def test_issue_5132_2():
     # not sure about the complex soln. But it looks correct.
     soln_complex = (img, z)
     soln = FiniteSet(soln_real, soln_complex)
-    assert nonlinsolve(eqs, [x, z]) == soln
+    assert dumeq(nonlinsolve(eqs, [x, z]), soln)
 
     system = [r - x**2 - y**2, tan(t) - y/x]
     s_x = sqrt(r/(tan(t)**2 + 1))
@@ -1647,7 +1650,7 @@ def test_issue_5132_substitution():
         (s_complex_y, s_complex_z_1),
         (s_complex_y, s_complex_z_2))
     soln = soln_real + soln_complex
-    assert substitution(eqs, [y, z]) == soln
+    assert dumeq(substitution(eqs, [y, z]), soln)
 
 
 def test_raises_substitution():
@@ -1872,9 +1875,9 @@ def test_issue_14300():
 
     assert solveset(f, x, S.Reals) == \
         Intersection(S.Reals, a1)
-    assert solveset(f, x) == \
+    assert dumeq(solveset(f, x),
         ImageSet(Lambda(n, -I*(2*n*pi + arg(-y + 1))/18000000 -
-            log(Abs(y - 1))/18000000), S.Integers)
+            log(Abs(y - 1))/18000000), S.Integers))
 
 
 def test_issue_14454():
@@ -1953,14 +1956,14 @@ def test_exponential_complex():
     from sympy import Dummy
     n = Dummy('n')
 
-    assert solveset_complex(2**x + 4**x, x) == imageset(
-        Lambda(n, I*(2*n*pi + pi)/log(2)), S.Integers)
+    assert dumeq(solveset_complex(2**x + 4**x, x),imageset(
+        Lambda(n, I*(2*n*pi + pi)/log(2)), S.Integers))
     assert solveset_complex(x**z*y**z - 2, z) == FiniteSet(
         log(2)/(log(x) + log(y)))
-    assert solveset_complex(4**(x/2) - 2**(x/3), x) == imageset(
-        Lambda(n, 3*n*I*pi/log(2)), S.Integers)
-    assert solveset(2**x + 32, x) == imageset(
-        Lambda(n, (I*(2*n*pi + pi) + 5*log(2))/log(2)), S.Integers)
+    assert dumeq(solveset_complex(4**(x/2) - 2**(x/3), x), imageset(
+        Lambda(n, 3*n*I*pi/log(2)), S.Integers))
+    assert dumeq(solveset(2**x + 32, x), imageset(
+        Lambda(n, (I*(2*n*pi + pi) + 5*log(2))/log(2)), S.Integers))
 
     eq = (2**exp(y**2/x) + 2)/(x**2 + 15)
     a = sqrt(x)*sqrt(-log(log(2)) + log(log(2) + 2*n*I*pi))
@@ -1968,14 +1971,14 @@ def test_exponential_complex():
 
     union1 = imageset(Lambda(n, I*(2*n*pi - pi*Rational(2, 3))/log(2)), S.Integers)
     union2 = imageset(Lambda(n, I*(2*n*pi + pi*Rational(2, 3))/log(2)), S.Integers)
-    assert solveset(2**x + 4**x + 8**x, x) == Union(union1, union2)
+    assert dumeq(solveset(2**x + 4**x + 8**x, x), Union(union1, union2))
 
     eq = 4**(x + 1) + 4**(x + 2) + 4**(x - 1) - 3**(x + 2) - 3**(x + 3)
     res = solveset(eq, x)
     num = 2*n*I*pi - 4*log(2) + 2*log(3)
     den = -2*log(2) + log(3)
     ans = imageset(Lambda(n, num/den), S.Integers)
-    assert res == ans
+    assert dumeq(res, ans)
 
 
 def test_expo_conditionset():
@@ -2149,25 +2152,25 @@ def test_invert_modular():
     assert invert_modular(Mod(exp(x), 7), S(5), n, x) == (Mod(exp(x), 7), 5)
     assert invert_modular(Mod(log(x), 7), S(5), n, x) == (Mod(log(x), 7), 5)
     # a is symbol
-    assert invert_modular(Mod(x, 7), S(5), n, x) == \
-            (x, ImageSet(Lambda(n, 7*n + 5), S.Integers))
+    assert dumeq(invert_modular(Mod(x, 7), S(5), n, x),
+            (x, ImageSet(Lambda(n, 7*n + 5), S.Integers)))
     # a.is_Add
-    assert invert_modular(Mod(x + 8, 7), S(5), n, x) == \
-            (x, ImageSet(Lambda(n, 7*n + 4), S.Integers))
+    assert dumeq(invert_modular(Mod(x + 8, 7), S(5), n, x),
+            (x, ImageSet(Lambda(n, 7*n + 4), S.Integers)))
     assert invert_modular(Mod(x**2 + x, 7), S(5), n, x) == \
             (Mod(x**2 + x, 7), 5)
     # a.is_Mul
-    assert invert_modular(Mod(3*x, 7), S(5), n, x) == \
-            (x, ImageSet(Lambda(n, 7*n + 4), S.Integers))
+    assert dumeq(invert_modular(Mod(3*x, 7), S(5), n, x),
+            (x, ImageSet(Lambda(n, 7*n + 4), S.Integers)))
     assert invert_modular(Mod((x + 1)*(x + 2), 7), S(5), n, x) == \
             (Mod((x + 1)*(x + 2), 7), 5)
     # a.is_Pow
     assert invert_modular(Mod(x**4, 7), S(5), n, x) == \
             (x, EmptySet())
-    assert invert_modular(Mod(3**x, 4), S(3), n, x) == \
-            (x, ImageSet(Lambda(n, 2*n + 1), S.Naturals0))
-    assert invert_modular(Mod(2**(x**2 + x + 1), 7), S(2), n, x) == \
-            (x**2 + x + 1, ImageSet(Lambda(n, 3*n + 1), S.Naturals0))
+    assert dumeq(invert_modular(Mod(3**x, 4), S(3), n, x),
+            (x, ImageSet(Lambda(n, 2*n + 1), S.Naturals0)))
+    assert dumeq(invert_modular(Mod(2**(x**2 + x + 1), 7), S(2), n, x),
+            (x**2 + x + 1, ImageSet(Lambda(n, 3*n + 1), S.Naturals0)))
     assert invert_modular(Mod(sin(x)**4, 7), S(5), n, x) == (x, EmptySet())
 
 
@@ -2188,17 +2191,18 @@ def test_solve_modular():
     assert solveset(7 - Mod(x, 5), x, S.Integers) == EmptySet()
     assert solveset(5 - Mod(x, 5), x, S.Integers) == EmptySet()
     # Negative m
-    assert solveset(2 + Mod(x, -3), x, S.Integers) == \
-            ImageSet(Lambda(n, -3*n - 2), S.Integers)
+    assert dumeq(solveset(2 + Mod(x, -3), x, S.Integers),
+            ImageSet(Lambda(n, -3*n - 2), S.Integers))
     assert solveset(4 + Mod(x, -3), x, S.Integers) == EmptySet()
     # linear expression in Mod
-    assert solveset(3 - Mod(x, 5), x, S.Integers) == ImageSet(Lambda(n, 5*n + 3), S.Integers)
-    assert solveset(3 - Mod(5*x - 8, 7), x, S.Integers) == \
-                ImageSet(Lambda(n, 7*n + 5), S.Integers)
-    assert solveset(3 - Mod(5*x, 7), x, S.Integers) == \
-                ImageSet(Lambda(n, 7*n + 2), S.Integers)
+    assert dumeq(solveset(3 - Mod(x, 5), x, S.Integers),
+        ImageSet(Lambda(n, 5*n + 3), S.Integers))
+    assert dumeq(solveset(3 - Mod(5*x - 8, 7), x, S.Integers),
+                ImageSet(Lambda(n, 7*n + 5), S.Integers))
+    assert dumeq(solveset(3 - Mod(5*x, 7), x, S.Integers),
+                ImageSet(Lambda(n, 7*n + 2), S.Integers))
     # higher degree expression in Mod
-    assert solveset(Mod(x**2, 160) - 9, x, S.Integers) == \
+    assert dumeq(solveset(Mod(x**2, 160) - 9, x, S.Integers),
             Union(ImageSet(Lambda(n, 160*n + 3), S.Integers),
             ImageSet(Lambda(n, 160*n + 13), S.Integers),
             ImageSet(Lambda(n, 160*n + 67), S.Integers),
@@ -2206,35 +2210,35 @@ def test_solve_modular():
             ImageSet(Lambda(n, 160*n + 83), S.Integers),
             ImageSet(Lambda(n, 160*n + 93), S.Integers),
             ImageSet(Lambda(n, 160*n + 147), S.Integers),
-            ImageSet(Lambda(n, 160*n + 157), S.Integers))
+            ImageSet(Lambda(n, 160*n + 157), S.Integers)))
     assert solveset(3 - Mod(x**4, 7), x, S.Integers) == EmptySet()
-    assert solveset(Mod(x**4, 17) - 13, x, S.Integers) == \
+    assert dumeq(solveset(Mod(x**4, 17) - 13, x, S.Integers),
             Union(ImageSet(Lambda(n, 17*n + 3), S.Integers),
             ImageSet(Lambda(n, 17*n + 5), S.Integers),
             ImageSet(Lambda(n, 17*n + 12), S.Integers),
-            ImageSet(Lambda(n, 17*n + 14), S.Integers))
+            ImageSet(Lambda(n, 17*n + 14), S.Integers)))
     # a.is_Pow tests
-    assert solveset(Mod(7**x, 41) - 15, x, S.Integers) == \
-            ImageSet(Lambda(n, 40*n + 3), S.Naturals0)
-    assert solveset(Mod(12**x, 21) - 18, x, S.Integers) == \
-            ImageSet(Lambda(n, 6*n + 2), S.Naturals0)
-    assert solveset(Mod(3**x, 4) - 3, x, S.Integers) == \
-            ImageSet(Lambda(n, 2*n + 1), S.Naturals0)
-    assert solveset(Mod(2**x, 7) - 2 , x, S.Integers) == \
-            ImageSet(Lambda(n, 3*n + 1), S.Naturals0)
-    assert solveset(Mod(3**(3**x), 4) - 3, x, S.Integers) == \
+    assert dumeq(solveset(Mod(7**x, 41) - 15, x, S.Integers),
+            ImageSet(Lambda(n, 40*n + 3), S.Naturals0))
+    assert dumeq(solveset(Mod(12**x, 21) - 18, x, S.Integers),
+            ImageSet(Lambda(n, 6*n + 2), S.Naturals0))
+    assert dumeq(solveset(Mod(3**x, 4) - 3, x, S.Integers),
+            ImageSet(Lambda(n, 2*n + 1), S.Naturals0))
+    assert dumeq(solveset(Mod(2**x, 7) - 2 , x, S.Integers),
+            ImageSet(Lambda(n, 3*n + 1), S.Naturals0))
+    assert dumeq(solveset(Mod(3**(3**x), 4) - 3, x, S.Integers),
             Intersection(ImageSet(Lambda(n, Intersection({log(2*n + 1)/log(3)},
-            S.Integers)), S.Naturals0), S.Integers)
+            S.Integers)), S.Naturals0), S.Integers))
     # Implemented for m without primitive root
     assert solveset(Mod(x**3, 7) - 2, x, S.Integers) == EmptySet()
-    assert solveset(Mod(x**3, 8) - 1, x, S.Integers) == \
-            ImageSet(Lambda(n, 8*n + 1), S.Integers)
-    assert solveset(Mod(x**4, 9) - 4, x, S.Integers) == \
+    assert dumeq(solveset(Mod(x**3, 8) - 1, x, S.Integers),
+            ImageSet(Lambda(n, 8*n + 1), S.Integers))
+    assert dumeq(solveset(Mod(x**4, 9) - 4, x, S.Integers),
             Union(ImageSet(Lambda(n, 9*n + 4), S.Integers),
-            ImageSet(Lambda(n, 9*n + 5), S.Integers))
+            ImageSet(Lambda(n, 9*n + 5), S.Integers)))
     # domain intersection
-    assert solveset(3 - Mod(5*x - 8, 7), x, S.Naturals0) == \
-            Intersection(ImageSet(Lambda(n, 7*n + 5), S.Integers), S.Naturals0)
+    assert dumeq(solveset(3 - Mod(5*x - 8, 7), x, S.Naturals0),
+            Intersection(ImageSet(Lambda(n, 7*n + 5), S.Integers), S.Naturals0))
     # Complex args
     assert solveset(Mod(x, 3) - I, x, S.Integers) == \
             EmptySet()
@@ -2244,11 +2248,11 @@ def test_solve_modular():
             ConditionSet(x, Eq(Mod(x + I, 3) - 2, 0), S.Integers)
 
     # issue 17373 (https://github.com/sympy/sympy/issues/17373)
-    assert solveset(Mod(x**4, 14) - 11, x, S.Integers) == \
+    assert dumeq(solveset(Mod(x**4, 14) - 11, x, S.Integers),
             Union(ImageSet(Lambda(n, 14*n + 3), S.Integers),
-            ImageSet(Lambda(n, 14*n + 11), S.Integers))
-    assert solveset(Mod(x**31, 74) - 43, x, S.Integers) == \
-            ImageSet(Lambda(n, 74*n + 31), S.Integers)
+            ImageSet(Lambda(n, 14*n + 11), S.Integers)))
+    assert dumeq(solveset(Mod(x**31, 74) - 43, x, S.Integers),
+            ImageSet(Lambda(n, 74*n + 31), S.Integers))
 
     # issue 13178
     n = symbols('n', integer=True)
@@ -2256,18 +2260,18 @@ def test_solve_modular():
     b = 1898888478
     m = 2**31 - 1
     c = 20170816
-    assert solveset(c - Mod(a**n*b, m), n, S.Integers) == \
-            ImageSet(Lambda(n, 2147483646*n + 100), S.Naturals0)
-    assert solveset(c - Mod(a**n*b, m), n, S.Naturals0) == \
+    assert dumeq(solveset(c - Mod(a**n*b, m), n, S.Integers),
+            ImageSet(Lambda(n, 2147483646*n + 100), S.Naturals0))
+    assert dumeq(solveset(c - Mod(a**n*b, m), n, S.Naturals0),
             Intersection(ImageSet(Lambda(n, 2147483646*n + 100), S.Naturals0),
-            S.Naturals0)
-    assert solveset(c - Mod(a**(2*n)*b, m), n, S.Integers) == \
+            S.Naturals0))
+    assert dumeq(solveset(c - Mod(a**(2*n)*b, m), n, S.Integers),
             Intersection(ImageSet(Lambda(n, 1073741823*n + 50), S.Naturals0),
-            S.Integers)
+            S.Integers))
     assert solveset(c - Mod(a**(2*n + 7)*b, m), n, S.Integers) == EmptySet()
-    assert solveset(c - Mod(a**(n - 4)*b, m), n, S.Integers) == \
+    assert dumeq(solveset(c - Mod(a**(n - 4)*b, m), n, S.Integers),
             Intersection(ImageSet(Lambda(n, 2147483646*n + 104), S.Naturals0),
-            S.Integers)
+            S.Integers))
 
 # end of modular tests
 
diff --git a/sympy/stats/tests/test_rv.py b/sympy/stats/tests/test_rv.py
index 6d91e3097d..8891fb3ea3 100644
--- a/sympy/stats/tests/test_rv.py
+++ b/sympy/stats/tests/test_rv.py
@@ -44,6 +44,7 @@ def test_random_symbols():
     assert set(random_symbols(2*X + Y.symbol)) == set((X,))
     assert set(random_symbols(2)) == set()
 
+
 def test_characteristic_function():
     #  Imports I from sympy
     from sympy import I
@@ -56,9 +57,10 @@ def test_characteristic_function():
     R = Lambda(t, exp(2 * exp(t*I) - 2))
 
 
-    assert characteristic_function(X) == P
-    assert characteristic_function(Y) == Q
-    assert characteristic_function(Z) == R
+    assert characteristic_function(X).dummy_eq(P)
+    assert characteristic_function(Y).dummy_eq(Q)
+    assert characteristic_function(Z).dummy_eq(R)
+
 
 def test_moment_generating_function():
 
@@ -71,9 +73,9 @@ def test_moment_generating_function():
     R = Lambda(t, exp(2 * exp(t) - 2))
 
 
-    assert moment_generating_function(X) == P
-    assert moment_generating_function(Y) == Q
-    assert moment_generating_function(Z) == R
+    assert moment_generating_function(X).dummy_eq(P)
+    assert moment_generating_function(Y).dummy_eq(Q)
+    assert moment_generating_function(Z).dummy_eq(R)
 
 def test_sample_iter():
 
diff --git a/sympy/utilities/tests/test_wester.py b/sympy/utilities/tests/test_wester.py
index 28cc5f2224..1e6d676803 100644
--- a/sympy/utilities/tests/test_wester.py
+++ b/sympy/utilities/tests/test_wester.py
@@ -28,7 +28,7 @@
 from sympy.functions.special.error_functions import Ci, Si, erf
 from sympy.functions.special.zeta_functions import zeta
 from sympy.testing.pytest import (XFAIL, slow, SKIP, skip, ON_TRAVIS,
-    raises, nocache_fail)
+    raises)
 from sympy.utilities.iterables import partitions
 from mpmath import mpi, mpc
 from sympy.matrices import Matrix, GramSchmidt, eye
@@ -934,14 +934,14 @@ def test_M14():
     assert solveset_real(tan(x) - 1, x) == ImageSet(Lambda(n, n*pi + pi/4), S.Integers)
 
 
-@nocache_fail
 def test_M15():
     n = Dummy('n')
-    # This test fails when running with the cache off:
-    assert solveset(sin(x) - S.Half) in (Union(ImageSet(Lambda(n, 2*n*pi + pi/6), S.Integers),
-                                           ImageSet(Lambda(n, 2*n*pi + pi*R(5, 6)), S.Integers)),
-                                           Union(ImageSet(Lambda(n, 2*n*pi + pi*R(5, 6)), S.Integers),
-                                           ImageSet(Lambda(n, 2*n*pi + pi/6), S.Integers)))
+    got = solveset(sin(x) - S.Half)
+    assert any(got.dummy_eq(i) for i in (
+        Union(ImageSet(Lambda(n, 2*n*pi + pi/6), S.Integers),
+        ImageSet(Lambda(n, 2*n*pi + pi*R(5, 6)), S.Integers)),
+        Union(ImageSet(Lambda(n, 2*n*pi + pi*R(5, 6)), S.Integers),
+        ImageSet(Lambda(n, 2*n*pi + pi/6), S.Integers))))
 
 
 @XFAIL