Repository: sympy/sympy
Pull request: 18778 fuzzy comparisons
<!-- 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
Fixes #18778 


#### Brief description of what is fixed or changed
Add to relational a function called cmp(lhs, rhs, op). Where op can be one of
In calculus utils,  and under AccumulationBounds, move the logic into _eval_is_lt, _eval_is_gt.
Did the same thing for Indexed in the tensor module.
With this new pull request, any class that inherits from Expr should implement _eval_is_X, if you want to override the default comparison. Returning None if nothing can be evaluated.

Moved the logic from Expr cmp into cmp(lhs, rhs, op).
Moved the logic from Equlaity into cmp - which calls _is_eq
The logic under _InEquality also uses the cmp logic.

 Examples:
    >>> from sympy.core.relational import cmp
    >>> from sympy import sympify
    >>> from sympy import Symbol
    >>> x = Symbol('x')
    >>> two = sympify(2)
    >>> cmp(x,two, "<") is None

    >>> from sympy.core.relational import cmp
    >>> from sympy import sympify
    >>> from sympy import Symbol
    >>> x = Symbol('x')
    >>> two = sympify(2)
    >>> cmp(x,two, "=") is None

    >>> from sympy.core.relational import cmp
    >>> from sympy import sympify
    >>> from sympy import Symbol
    >>> three = sympify(3)
    >>> two = sympify(2)
    >>> cmp(three,two, "<") is False

    >>> from sympy.core.relational import cmp
    >>> from sympy import sympify
    >>> from sympy import Symbol
    >>> three = sympify(3)
    >>> two = sympify(2)
    >>> cmp(three,two, ">") is True

#### Other comments
_eval_Eq doesn't match with _eval_is_gt however,there are already a great number of classes

#### Release Notes
<!-- BEGIN RELEASE NOTES -->
NO ENTRY
<!-- END RELEASE NOTES -->

Problem statement:
Add fuzzy logic comparison functions
There are many bugs in sympy that come from using inequalities in a bool context e.g.:
```julia
In [5]: def f(x): 
   ...:     if x < 0: 
   ...:         return 'negative' 
   ...:     else: 
   ...:         return 'nonnegative' 
   ...:                                                                                                                                        

In [6]: f(1)                                                                                                                                   
Out[6]: 'nonnegative'

In [7]: f(-1)                                                                                                                                  
Out[7]: 'negative'

In [8]: f(Symbol('y'))                                                                                                                         
--------------------------------------------------------------------------
TypeError: cannot determine truth value of Relational
```
The solution to this is to use `x.is_negative` which gives a fuzzy bool (True, False, None) and does not raise. That works for comparisons with zero but not for other comparisons like `x<1`:
```
if x < 1: # raises for symbol x

if (x - 1).is_negative: # works but creates a new expression needlessly

if (x < 1) is S.true: # works but not obvious

if (x < 1) == True: # works but not obvious
```
I suggest adding new functions that can replace the use of `Eq`, `Ne`, `Lt` etc and return fuzzy bools. The new functions can have names analogous to those of the old assumptions that also use fuzzy logic e.g.:
```python
is_equal(x, y)
is_unequal(x, y)
is_less_than(x, y)
is_less_than_or_equal(x, y)
...
```
These can be used in place of e.g. `Eq(x, y) is S.true` and return True, False, None. The implementation of `Eq` should then become
```python
class Eq(Relational):
    def __new__(cls, lhs, rhs):
        value = is_equal(lhs, rhs)
        if value is True:
            return S.true
        elif value is False:
            return S.false
        else:
            # Return unevaluated
            return Relational.__new__(cls, lhs, rhs)
```
This would cleanly separate the evaluation layer from the querying layer so that evaluation is built on top of querying and querying does not in itself involve creating new instances and evaluating them. Internal code should be changed to use these functions rather than creating Relational instances.

Some of the names can be long e.g. `is_less_than_or_equal` so perhaps an alternative is to use LaTeX-style short names:
```
is_eq(x, y)
is_ne(x, y)
is_lt(x, y)
is_le(x, y)
```
Another possibility is something like
```
cmp(x, '<=', y)
```
which might be more readable for someone new to the codebase.

Making this happen is as simple as taking the code out of the various `Relational.__new__` methods and putting them into separate functions.

Hints from issue discussion:
> works but creates a new expression needlessly

Wouldn't that happen internally anyway?
I only have concerns that this can make things like `is_eq` big because it can accept any expressions.
So managing such things would be an another concern.

I was also undecided about how things like `x > 2` should be handled in binary evaluations or should be coerced into `(x-2).is_positive`.
As far as I can see, sympy doesn't have logics for simple inequalities like triangle inequality, so the things should be directed to make such implementations easier.
The `is_eq` function would just be the body of `Equality.__new__` so it won't be any bigger than that would be anyway. Currently `Equality.__new__` defers to `_eval_Eq` where appropriate.

Currently `Expr.__lt__` and `Lt.__new__` will usually end up doing `(lhs-rhs).is_negative` but that's partly just because there is no `is_lt` function to call. Asking whether `Add(a, -b)` is negative is not inherently any easier than asking whether `a < b`. We just only have an API for one of these and not the other. If we had `is_lt` we could move some of the logic from `Add._eval_is_negative` there and use `is_lt` as part of `Add._eval_is_negative`.

For inequalities the current situation is this:
- `LessThan(a, b)` calls `a.__lt__(b)`  (not `a<b` because of #7951)
- `Expr.__lt__(a, b)` usually calls `(a - b).is_extended_negative`
- If `is_extended_negative` is None then `Expr.__lt__` returns `LessThan(a, b, evaluate=False)`

I would propose to change that to:
- `Expr.__lt__` returns `LessThan(a, b)`.
- `LessThan(a, b)` calls `is_lt(a, b)` and evaluates or not from there.
- The `is_lt` function has some logic of its own (taken from current `Expr.__lt__`) and calls `_eval_is_lt` if present.
- At least for now `is_lt` will mostly defer to `(a - b).is_extended_negative`.
- Various `__lt__` methods around sympy should become `_eval_is_lt` instead.

That way all logic always go through the same function `is_lt`.
Ideally `ask(x < y)` would work. 
> Ideally `ask(x < y)` would work.

Yes, but that doesn't help us here unless we plan to stop evaluating `x < y` automatically.
I think it isn't a problem if it gives True or False, since ask(True) and ask(False) do what you'd expect. However, this is a problem

```py
>>> I > 3
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "./sympy/core/expr.py", line 403, in __gt__
    return self._cmp(other, ">", StrictGreaterThan)
  File "./sympy/core/expr.py", line 348, in _cmp
    raise TypeError("Invalid comparison of non-real %s" % me)
TypeError: Invalid comparison of non-real I
```

I think we should strongly consider making that return False instead, for reasons I've outlined in other issues. 
how can I work with this issue?

> [regarding `I > 3`]...I think we should strongly consider making that return False instead, for reasons I've outlined in other issues.

I thought it would be better to leave it unevaluated, but I am now in favor with returning False since a pure imaginary evaluates as False in terms of sign:
```python
>>> ok = Dummy(negative=False, positive=False, zero=False)  # not inconsistent
>>> I.is_negative
False
>>> I.is_positive
False
>>> I.is_zero
False
```
Do we want the same thing for `nan`?
> Some of the names can be long e.g. `is_less_than_or_equal` so perhaps an alternative is to use LaTeX-style short names:
> 
> ```
> is_eq(x, y)
> is_ne(x, y)
> is_lt(x, y)
> is_le(x, y)
> ```
> 
> Another possibility is something like
> 
> ```
> cmp(x, '<=', y)
> ```
> 
> which might be more readable for someone new to the codebase.

What about `x.lt(y)` or just `lt(x,y)` which is very similar to the current Lt(x, y)? The difference being that the lowercase version always returns something while the uppercase can remain unevaluated in keeping with things like `sum vs Sum`.
The `is_*` naming scheme can be extended to the `fuzzy_*` names, too: `is_and(*x) -> fuzzy_and(x)` (note different signature).

Commit message: feat(core): Add fuzzy comparison functions

Remove unneeded imports and blank lines

Patch to review:
diff --git a/sympy/calculus/util.py b/sympy/calculus/util.py
index 50a9be1c8282..71ecb0050754 100644
--- a/sympy/calculus/util.py
+++ b/sympy/calculus/util.py
@@ -13,7 +13,7 @@
 from sympy.sets.fancysets import ImageSet
 from sympy.solvers.inequalities import solve_univariate_inequality
 from sympy.utilities import filldedent
-
+from sympy.multipledispatch import dispatch
 
 def continuous_domain(f, symbol, domain):
     """
@@ -1372,153 +1372,7 @@ def __abs__(self):
         else:
             return self
 
-    def __lt__(self, other):
-        """
-        Returns True if range of values attained by `self` AccumulationBounds
-        object is less than the range of values attained by `other`, where
-        other may be any value of type AccumulationBounds object or extended
-        real number value, False if `other` satisfies the same property, else
-        an unevaluated Relational
 
-        Examples
-        ========
-
-        >>> from sympy import AccumBounds, oo
-        >>> AccumBounds(1, 3) < AccumBounds(4, oo)
-        True
-        >>> AccumBounds(1, 4) < AccumBounds(3, 4)
-        AccumBounds(1, 4) < AccumBounds(3, 4)
-        >>> AccumBounds(1, oo) < -1
-        False
-
-        """
-        other = _sympify(other)
-        if isinstance(other, AccumBounds):
-            if self.max < other.min:
-                return True
-            if self.min >= other.max:
-                return False
-        elif not other.is_extended_real:
-            raise TypeError(
-                "Invalid comparison of %s %s" %
-                (type(other), other))
-        elif other.is_comparable:
-            if self.max < other:
-                return True
-            if self.min >= other:
-                return False
-        return super().__lt__(other)
-
-    def __le__(self, other):
-        """
-        Returns True if range of values attained by `self` AccumulationBounds
-        object is less than or equal to the range of values attained by
-        `other`, where other may be any value of type AccumulationBounds
-        object or extended real number value, False if `other`
-        satisfies the same property, else an unevaluated Relational.
-
-        Examples
-        ========
-
-        >>> from sympy import AccumBounds, oo
-        >>> AccumBounds(1, 3) <= AccumBounds(4, oo)
-        True
-        >>> AccumBounds(1, 4) <= AccumBounds(3, 4)
-        AccumBounds(1, 4) <= AccumBounds(3, 4)
-        >>> AccumBounds(1, 3) <= 0
-        False
-
-        """
-        other = _sympify(other)
-        if isinstance(other, AccumBounds):
-            if self.max <= other.min:
-                return True
-            if self.min > other.max:
-                return False
-        elif not other.is_extended_real:
-            raise TypeError(
-                "Invalid comparison of %s %s" %
-                (type(other), other))
-        elif other.is_comparable:
-            if self.max <= other:
-                return True
-            if self.min > other:
-                return False
-        return super().__le__(other)
-
-    def __gt__(self, other):
-        """
-        Returns True if range of values attained by `self` AccumulationBounds
-        object is greater than the range of values attained by `other`,
-        where other may be any value of type AccumulationBounds object or
-        extended real number value, False if `other` satisfies
-        the same property, else an unevaluated Relational.
-
-        Examples
-        ========
-
-        >>> from sympy import AccumBounds, oo
-        >>> AccumBounds(1, 3) > AccumBounds(4, oo)
-        False
-        >>> AccumBounds(1, 4) > AccumBounds(3, 4)
-        AccumBounds(1, 4) > AccumBounds(3, 4)
-        >>> AccumBounds(1, oo) > -1
-        True
-
-        """
-        other = _sympify(other)
-        if isinstance(other, AccumBounds):
-            if self.min > other.max:
-                return True
-            if self.max <= other.min:
-                return False
-        elif not other.is_extended_real:
-            raise TypeError(
-                "Invalid comparison of %s %s" %
-                (type(other), other))
-        elif other.is_comparable:
-            if self.min > other:
-                return True
-            if self.max <= other:
-                return False
-        return super().__gt__(other)
-
-    def __ge__(self, other):
-        """
-        Returns True if range of values attained by `self` AccumulationBounds
-        object is less that the range of values attained by `other`, where
-        other may be any value of type AccumulationBounds object or extended
-        real number value, False if `other` satisfies the same
-        property, else an unevaluated Relational.
-
-        Examples
-        ========
-
-        >>> from sympy import AccumBounds, oo
-        >>> AccumBounds(1, 3) >= AccumBounds(4, oo)
-        False
-        >>> AccumBounds(1, 4) >= AccumBounds(3, 4)
-        AccumBounds(1, 4) >= AccumBounds(3, 4)
-        >>> AccumBounds(1, oo) >= 1
-        True
-
-        """
-        other = _sympify(other)
-        if isinstance(other, AccumBounds):
-            if self.min >= other.max:
-                return True
-            if self.max < other.min:
-                return False
-        elif not other.is_extended_real:
-            raise TypeError(
-                "Invalid comparison of %s %s" %
-                (type(other), other))
-        elif other.is_comparable:
-            if self.min >= other:
-                return True
-            if self.max < other:
-                return False
-        return super().__ge__(other)
 
     def __contains__(self, other):
         """
@@ -1627,5 +1481,103 @@ def union(self, other):
             return AccumBounds(other.min, Max(self.max, other.max))
 
 
+@dispatch(AccumulationBounds, AccumulationBounds) # type: ignore # noqa:F811
+def _eval_is_le(lhs, rhs): # noqa:F811
+    if lhs.max <= rhs.min:
+        return True
+    if lhs.min > rhs.max:
+        return False
+
+@dispatch(AccumulationBounds, Basic) # type: ignore # noqa:F811
+def _eval_is_le(lhs, rhs): # noqa: F811
+
+    """
+    Returns True if range of values attained by `self` AccumulationBounds
+    object is greater than the range of values attained by `other`,
+    where other may be any value of type AccumulationBounds object or
+    extended real number value, False if `other` satisfies
+    the same property, else an unevaluated Relational.
+
+    Examples
+    ========
+
+    >>> from sympy import AccumBounds, oo
+    >>> AccumBounds(1, 3) > AccumBounds(4, oo)
+    False
+    >>> AccumBounds(1, 4) > AccumBounds(3, 4)
+    AccumBounds(1, 4) > AccumBounds(3, 4)
+    >>> AccumBounds(1, oo) > -1
+    True
+
+    """
+    if not rhs.is_extended_real:
+            raise TypeError(
+                "Invalid comparison of %s %s" %
+                (type(rhs), rhs))
+    elif rhs.is_comparable:
+        if lhs.max <= rhs:
+            return True
+        if lhs.min > rhs:
+            return False
+
+@dispatch(AccumulationBounds, AccumulationBounds)
+def _eval_is_ge(lhs, rhs): # noqa:F811
+    if lhs.min >= rhs.max:
+        return True
+    if lhs.max < rhs.min:
+        return False
+
+@dispatch(AccumulationBounds, Basic)
+def _eval_is_ge(lhs, rhs): # noqa: F811
+    """
+    Returns True if range of values attained by `lhs` AccumulationBounds
+    object is less that the range of values attained by `rhs`, where
+    other may be any value of type AccumulationBounds object or extended
+    real number value, False if `rhs` satisfies the same
+    property, else an unevaluated Relational.
+
+    Examples
+    ========
+
+    >>> from sympy import AccumBounds, oo
+    >>> AccumBounds(1, 3) >= AccumBounds(4, oo)
+    False
+    >>> AccumBounds(1, 4) >= AccumBounds(3, 4)
+    AccumBounds(1, 4) >= AccumBounds(3, 4)
+    >>> AccumBounds(1, oo) >= 1
+    True
+    """
+
+    if not rhs.is_extended_real:
+        raise TypeError(
+            "Invalid comparison of %s %s" %
+            (type(rhs), rhs))
+    elif rhs.is_comparable:
+        if lhs.min >= rhs:
+            return True
+        if lhs.max < rhs:
+            return False
+
+
+@dispatch(Basic, AccumulationBounds)
+def _eval_is_ge(lhs, rhs): # noqa:F811
+    if not lhs.is_extended_real:
+        raise TypeError(
+            "Invalid comparison of %s %s" %
+            (type(lhs), lhs))
+    elif lhs.is_comparable:
+        if rhs.max <= lhs:
+            return True
+        if rhs.min > lhs:
+            return False
+
+
+@dispatch(AccumulationBounds, AccumulationBounds)
+def _eval_is_ge(lhs, rhs): # noqa:F811
+    if lhs.min >= rhs.max:
+        return True
+    if lhs.max < rhs.min:
+        return False
+
 # setting an alias for AccumulationBounds
 AccumBounds = AccumulationBounds
diff --git a/sympy/combinatorics/permutations.py b/sympy/combinatorics/permutations.py
index da19cbf50601..874c093219a7 100644
--- a/sympy/combinatorics/permutations.py
+++ b/sympy/combinatorics/permutations.py
@@ -14,7 +14,7 @@
 from sympy.utilities.iterables import (flatten, has_variety, minlex,
     has_dups, runs)
 from mpmath.libmp.libintmath import ifac
-
+from sympy.multipledispatch import dispatch
 
 def _af_rmul(a, b):
     """
@@ -950,15 +950,7 @@ def __new__(cls, *args, **kwargs):
 
         return cls._af_new(aform)
 
-    def _eval_Eq(self, other):
-        other = _sympify(other)
-        if not isinstance(other, Permutation):
-            return None
-
-        if self._size != other._size:
-            return None
 
-        return as_Boolean(self._array_form == other._array_form)
 
     @classmethod
     def _af_new(cls, perm):
@@ -3020,3 +3012,10 @@ def __new__(cls, perm, x, evaluate=None):
 
         obj = super().__new__(cls, perm, x)
         return obj
+
+
+@dispatch(Permutation, Permutation)
+def _eval_Eq(lhs, rhs):
+    if lhs._size != rhs._size:
+        return None
+    return as_Boolean(lhs._array_form == rhs._array_form)
diff --git a/sympy/core/containers.py b/sympy/core/containers.py
index 5db2885d0a88..8acf6b4a19f9 100644
--- a/sympy/core/containers.py
+++ b/sympy/core/containers.py
@@ -13,7 +13,7 @@
 from sympy.core.compatibility import as_int, MutableSet
 from sympy.core.sympify import sympify, converter
 from sympy.utilities.iterables import iterable
-
+from sympy.multipledispatch import dispatch
 
 
 class Tuple(Basic):
@@ -141,25 +141,27 @@ def index(self, value, start=None, stop=None):
         else:
             return self.args.index(value, start, stop)
 
-    def _eval_Eq(self, other):
-        from sympy.core.function import AppliedUndef
-        from sympy.core.logic import fuzzy_and, fuzzy_bool
-        from sympy.core.relational import Eq
 
-        if other.is_Symbol or isinstance(other, AppliedUndef):
-            return None
+converter[tuple] = lambda tup: Tuple(*tup)
+
 
-        if not isinstance(other, Tuple) or len(self) != len(other):
-            return S.false
+@dispatch(Tuple, Basic)
+def _eval_Eq(lhs, rhs):
+    from sympy.core.function import AppliedUndef
+    from sympy.core.logic import fuzzy_and, fuzzy_bool
+    from sympy.core.relational import Eq
 
-        r = fuzzy_and(fuzzy_bool(Eq(s, o)) for s, o in zip(self, other))
-        if r is True:
-            return S.true
-        elif r is False:
-            return S.false
+    if rhs.is_Symbol or isinstance(rhs, AppliedUndef):
+        return None
 
+    if not isinstance(rhs, Tuple) or len(lhs) != len(rhs):
+        return S.false
 
-converter[tuple] = lambda tup: Tuple(*tup)
+    r = fuzzy_and(fuzzy_bool(Eq(s, o)) for s, o in zip(lhs, rhs))
+    if r is True:
+        return S.true
+    elif r is False:
+        return S.false
 
 
 def tuple_wrapper(method):
diff --git a/sympy/core/expr.py b/sympy/core/expr.py
index 2b16305db3d4..d13715c31445 100644
--- a/sympy/core/expr.py
+++ b/sympy/core/expr.py
@@ -1,3 +1,4 @@
+from __future__ import print_function, division
 from typing import Tuple as tTuple
 
 from .sympify import sympify, _sympify, SympifyError
@@ -23,6 +24,10 @@ class Expr(Basic, EvalfMixin):
     used only for argument storage and expression manipulation, i.e.
     pattern matching, substitutions, etc).
 
+    If you want to override the comparisons of expressions:
+    Should use _eval_is_X (where X is gt, ge, lt, or Eq).
+    _eval_is_X should return None if you want to return the object unevaluated
+
     See Also
     ========
 
@@ -339,78 +344,30 @@ def __complex__(self):
         re, im = result.as_real_imag()
         return complex(float(re), float(im))
 
-    def _cmp(self, other, op, cls):
-        assert op in ("<", ">", "<=", ">=")
-        try:
-            other = _sympify(other)
-        except SympifyError:
-            return NotImplemented
-
-        if not isinstance(other, Expr):
-            return NotImplemented
-
-        for me in (self, other):
-            if me.is_extended_real is False:
-                raise TypeError("Invalid comparison of non-real %s" % me)
-            if me is S.NaN:
-                raise TypeError("Invalid NaN comparison")
-
-        n2 = _n2(self, other)
-        if n2 is not None:
-            # use float comparison for infinity.
-            # otherwise get stuck in infinite recursion
-            if n2 in (S.Infinity, S.NegativeInfinity):
-                n2 = float(n2)
-            if op == "<":
-                return _sympify(n2 < 0)
-            elif op == ">":
-                return _sympify(n2 > 0)
-            elif op == "<=":
-                return _sympify(n2 <= 0)
-            else: # >=
-                return _sympify(n2 >= 0)
-
-        if self.is_extended_real and other.is_extended_real:
-            if op in ("<=", ">") \
-                and ((self.is_infinite and self.is_extended_negative) \
-                     or (other.is_infinite and other.is_extended_positive)):
-                return S.true if op == "<=" else S.false
-            if op in ("<", ">=") \
-                and ((self.is_infinite and self.is_extended_positive) \
-                     or (other.is_infinite and other.is_extended_negative)):
-                return S.true if op == ">=" else S.false
-            diff = self - other
-            if diff is not S.NaN:
-                if op == "<":
-                    test = diff.is_extended_negative
-                elif op == ">":
-                    test = diff.is_extended_positive
-                elif op == "<=":
-                    test = diff.is_extended_nonpositive
-                else: # >=
-                    test = diff.is_extended_nonnegative
-
-                if test is not None:
-                    return S.true if test == True else S.false
-
-        # return unevaluated comparison object
-        return cls(self, other, evaluate=False)
-
     def __ge__(self, other):
-        from sympy import GreaterThan
-        return self._cmp(other, ">=", GreaterThan)
+        from .relational import GreaterThan
+        return self._cmp(other, GreaterThan)
 
     def __le__(self, other):
-        from sympy import LessThan
-        return self._cmp(other, "<=", LessThan)
+        from .relational import LessThan
+        return self._cmp(other, LessThan)
 
     def __gt__(self, other):
-        from sympy import StrictGreaterThan
-        return self._cmp(other, ">", StrictGreaterThan)
+        from .relational import StrictGreaterThan
+        return self._cmp(other, StrictGreaterThan)
 
     def __lt__(self, other):
-        from sympy import StrictLessThan
-        return self._cmp(other, "<", StrictLessThan)
+        from .relational import StrictLessThan
+        return self._cmp(other, StrictLessThan)
+
+    def _cmp(self, other, cls):
+        try:
+            other = _sympify(other)
+        except SympifyError:
+            return NotImplemented
+        if not isinstance(other, Expr):
+            return NotImplemented
+        return cls(self, other)
 
     def __trunc__(self):
         if not self.is_number:
@@ -3955,17 +3912,6 @@ def doit(self, **kwargs):
             return self.args[0]
 
 
-def _n2(a, b):
-    """Return (a - b).evalf(2) if a and b are comparable, else None.
-    This should only be used when a and b are already sympified.
-    """
-    # /!\ it is very important (see issue 8245) not to
-    # use a re-evaluated number in the calculation of dif
-    if a.is_comparable and b.is_comparable:
-        dif = (a - b).evalf(2)
-        if dif.is_comparable:
-            return dif
-
 
 def unchanged(func, *args):
     """Return True if `func` applied to the `args` is unchanged.
diff --git a/sympy/core/numbers.py b/sympy/core/numbers.py
index 27b671e148c9..bf1a7c1a8109 100644
--- a/sympy/core/numbers.py
+++ b/sympy/core/numbers.py
@@ -16,7 +16,7 @@
 from sympy.core.compatibility import (as_int, HAS_GMPY, SYMPY_INTS,
     int_info, gmpy)
 from sympy.core.cache import lru_cache
-
+from sympy.multipledispatch import dispatch
 import mpmath
 import mpmath.libmp as mlib
 from mpmath.libmp import bitcount
@@ -276,6 +276,8 @@ def _igcd2_python(a, b):
 
 
 # Use Lehmer's algorithm only for very large numbers.
+# The limit could be different on Python 2.7 and 3.x.
+# If so, then this could be defined in compatibility.py.
 BIGBITS = 5000
 def igcd_lehmer(a, b):
     """Computes greatest common divisor of two integers.
@@ -3264,10 +3266,6 @@ def __eq__(self, other):
     def __ne__(self, other):
         return other is not S.NaN
 
-    def _eval_Eq(self, other):
-        # NaN is not mathematically equal to anything, even NaN
-        return S.false
-
     # Expr will _sympify and raise TypeError
     __gt__ = Expr.__gt__
     __ge__ = Expr.__ge__
@@ -3276,6 +3274,9 @@ def _eval_Eq(self, other):
 
 nan = S.NaN
 
+@dispatch(NaN, Expr)
+def _eval_Eq(a, b): # noqa:F811
+    return S.false
 
 class ComplexInfinity(AtomicExpr, metaclass=Singleton):
     r"""Complex infinity.
@@ -3413,7 +3414,6 @@ def __long__(self):
     def __hash__(self):
         return super().__hash__()
 
-
 class Exp1(NumberSymbol, metaclass=Singleton):
     r"""The `e` constant.
 
@@ -3889,6 +3889,9 @@ def _mpc_(self):
 
 I = S.ImaginaryUnit
 
+@dispatch(Tuple, Number)
+def _eval_Eq(self, other): # noqa: F811
+    return S.false
 
 def sympify_fractions(f):
     return Rational(f.numerator, f.denominator, 1)
diff --git a/sympy/core/relational.py b/sympy/core/relational.py
index c65577b561b9..38819b77d9fe 100644
--- a/sympy/core/relational.py
+++ b/sympy/core/relational.py
@@ -1,15 +1,13 @@
 from typing import Dict, Type, Union
 
 from sympy.utilities.exceptions import SymPyDeprecationWarning
-from .add import _unevaluated_Add, Add
 from .basic import S, Atom
 from .compatibility import ordered
 from .basic import Basic
-from .expr import Expr
 from .evalf import EvalfMixin
-from .sympify import _sympify
+from .sympify import _sympify, SympifyError
 from .parameters import global_parameters
-
+from sympy.core.logic import fuzzy_bool, fuzzy_xor, fuzzy_and, fuzzy_not
 from sympy.logic.boolalg import Boolean, BooleanAtom
 
 __all__ = (
@@ -18,14 +16,21 @@
     'StrictGreaterThan', 'GreaterThan',
 )
 
+from .expr import Expr
+from sympy.multipledispatch import dispatch
+
+
 
 def _nontrivBool(side):
     return isinstance(side, Boolean) and \
-        not isinstance(side, (BooleanAtom, Atom))
+        not isinstance(side, Atom)
+
 
 
 # Note, see issue 4986.  Ideally, we wouldn't want to subclass both Boolean
 # and Expr.
+# from .. import Expr
+
 
 def _canonical(cond):
     # return a condition in which all relationals are canonical
@@ -291,6 +296,7 @@ def equals(self, other, failing_expression=False):
                 return left
 
     def _eval_simplify(self, **kwargs):
+        from .add import Add
         r = self
         r = r.func(*[i.simplify(**kwargs) for i in r.args])
         if r.is_Relational:
@@ -318,9 +324,9 @@ def _eval_simplify(self, **kwargs):
                         if m.is_negative:
                             # Dividing with a negative number, so change order of arguments
                             # canonical will put the symbol back on the lhs later
-                            r = r.func(-b/m, x)
+                            r = r.func(-b / m, x)
                         else:
-                            r = r.func(x, -b/m)
+                            r = r.func(x, -b / m)
                     else:
                         r = r.func(b, S.zero)
                 except ValueError:
@@ -332,8 +338,8 @@ def _eval_simplify(self, **kwargs):
                         constant = c[-1]
                         c[-1] = 0
                         scale = gcd(c)
-                        c = [ctmp/scale for ctmp in c]
-                        r = r.func(Poly.from_list(c, x).as_expr(), -constant/scale)
+                        c = [ctmp / scale for ctmp in c]
+                        r = r.func(Poly.from_list(c, x).as_expr(), -constant / scale)
                     except PolynomialError:
                         pass
             elif len(free) >= 2:
@@ -346,18 +352,18 @@ def _eval_simplify(self, **kwargs):
                     constant = m[-1]
                     del m[-1]
                     scale = gcd(m)
-                    m = [mtmp/scale for mtmp in m]
+                    m = [mtmp / scale for mtmp in m]
                     nzm = list(filter(lambda f: f[0] != 0, list(zip(m, free))))
                     if scale.is_zero is False:
                         if constant != 0:
                             # lhs: expression, rhs: constant
-                            newexpr = Add(*[i*j for i, j in nzm])
-                            r = r.func(newexpr, -constant/scale)
+                            newexpr = Add(*[i * j for i, j in nzm])
+                            r = r.func(newexpr, -constant / scale)
                         else:
                             # keep first term on lhs
-                            lhsterm = nzm[0][0]*nzm[0][1]
+                            lhsterm = nzm[0][0] * nzm[0][1]
                             del nzm[0]
-                            newexpr = Add(*[i*j for i, j in nzm])
+                            newexpr = Add(*[i * j for i, j in nzm])
                             r = r.func(lhsterm, -newexpr)
 
                     else:
@@ -367,7 +373,7 @@ def _eval_simplify(self, **kwargs):
         # Did we get a simplified result?
         r = r.canonical
         measure = kwargs['measure']
-        if measure(r) < kwargs['ratio']*measure(self):
+        if measure(r) < kwargs['ratio'] * measure(self):
             return r
         else:
             return self
@@ -475,12 +481,6 @@ class Equality(Relational):
     is_Equality = True
 
     def __new__(cls, lhs, rhs=None, **options):
-        from sympy.core.add import Add
-        from sympy.core.logic import fuzzy_bool, fuzzy_xor, fuzzy_and, fuzzy_not
-        from sympy.core.expr import _n2
-        from sympy.functions.elementary.complexes import arg
-        from sympy.simplify.simplify import clear_coefficients
-        from sympy.utilities.iterables import sift
 
         if rhs is None:
             SymPyDeprecationWarning(
@@ -490,122 +490,25 @@ def __new__(cls, lhs, rhs=None, **options):
                 deprecated_since_version="1.5"
             ).warn()
             rhs = 0
-
+        evaluate = options.pop('evaluate', global_parameters.evaluate)
         lhs = _sympify(lhs)
         rhs = _sympify(rhs)
-
-        evaluate = options.pop('evaluate', global_parameters.evaluate)
-
         if evaluate:
-            if isinstance(lhs, Boolean) != isinstance(lhs, Boolean):
-                # e.g. 0/1 not recognized as Boolean in SymPy
-                return S.false
-
-            # If one expression has an _eval_Eq, return its results.
-            if hasattr(lhs, '_eval_Eq'):
-                r = lhs._eval_Eq(rhs)
-                if r is not None:
-                    return r
-            if hasattr(rhs, '_eval_Eq'):
-                r = rhs._eval_Eq(lhs)
-                if r is not None:
-                    return r
-            # If expressions have the same structure, they must be equal.
-            if lhs == rhs:
-                return S.true  # e.g. True == True
-            elif all(isinstance(i, BooleanAtom) for i in (rhs, lhs)):
-                return S.false  # True != False
-            elif not (lhs.is_Symbol or rhs.is_Symbol) and (
-                    isinstance(lhs, Boolean) !=
-                    isinstance(rhs, Boolean)):
-                return S.false  # only Booleans can equal Booleans
-
-            if lhs.is_infinite or rhs.is_infinite:
-                if fuzzy_xor([lhs.is_infinite, rhs.is_infinite]):
-                    return S.false
-                if fuzzy_xor([lhs.is_extended_real, rhs.is_extended_real]):
-                    return S.false
-                if fuzzy_and([lhs.is_extended_real, rhs.is_extended_real]):
-                    r = fuzzy_xor([lhs.is_extended_positive, fuzzy_not(rhs.is_extended_positive)])
-                    return S(r)
-
-                # Try to split real/imaginary parts and equate them
-                I = S.ImaginaryUnit
-
-                def split_real_imag(expr):
-                    real_imag = lambda t: (
-                            'real' if t.is_extended_real else
-                            'imag' if (I*t).is_extended_real else None)
-                    return sift(Add.make_args(expr), real_imag)
-
-                lhs_ri = split_real_imag(lhs)
-                if not lhs_ri[None]:
-                    rhs_ri = split_real_imag(rhs)
-                    if not rhs_ri[None]:
-                        eq_real = Eq(Add(*lhs_ri['real']), Add(*rhs_ri['real']))
-                        eq_imag = Eq(I*Add(*lhs_ri['imag']), I*Add(*rhs_ri['imag']))
-                        res = fuzzy_and(map(fuzzy_bool, [eq_real, eq_imag]))
-                        if res is not None:
-                            return S(res)
-
-                # Compare e.g. zoo with 1+I*oo by comparing args
-                arglhs = arg(lhs)
-                argrhs = arg(rhs)
-                # Guard against Eq(nan, nan) -> False
-                if not (arglhs == S.NaN and argrhs == S.NaN):
-                    res = fuzzy_bool(Eq(arglhs, argrhs))
-                    if res is not None:
-                        return S(res)
-
-                return Relational.__new__(cls, lhs, rhs, **options)
-
-            if all(isinstance(i, Expr) for i in (lhs, rhs)):
-                # see if the difference evaluates
-                dif = lhs - rhs
-                z = dif.is_zero
-                if z is not None:
-                    if z is False and dif.is_commutative:  # issue 10728
-                        return S.false
-                    if z:
-                        return S.true
-                # evaluate numerically if possible
-                n2 = _n2(lhs, rhs)
-                if n2 is not None:
-                    return _sympify(n2 == 0)
-                # see if the ratio evaluates
-                n, d = dif.as_numer_denom()
-                rv = None
-                if n.is_zero:
-                    rv = d.is_nonzero
-                elif n.is_finite:
-                    if d.is_infinite:
-                        rv = S.true
-                    elif n.is_zero is False:
-                        rv = d.is_infinite
-                        if rv is None:
-                            # if the condition that makes the denominator
-                            # infinite does not make the original expression
-                            # True then False can be returned
-                            l, r = clear_coefficients(d, S.Infinity)
-                            args = [_.subs(l, r) for _ in (lhs, rhs)]
-                            if args != [lhs, rhs]:
-                                rv = fuzzy_bool(Eq(*args))
-                                if rv is True:
-                                    rv = None
-                elif any(a.is_infinite for a in Add.make_args(n)):
-                    # (inf or nan)/x != 0
-                    rv = S.false
-                if rv is not None:
-                    return _sympify(rv)
+            val = is_eq(lhs, rhs)
+            if val is None:
+                return cls(lhs, rhs, evaluate=False)
+            else:
+                return _sympify(val)
 
-        return Relational.__new__(cls, lhs, rhs, **options)
+        return Relational.__new__(cls, lhs, rhs)
 
     @classmethod
     def _eval_relation(cls, lhs, rhs):
         return _sympify(lhs == rhs)
 
     def _eval_rewrite_as_Add(self, *args, **kwargs):
-        """return Eq(L, R) as L - R. To control the evaluation of
+        """
+        return Eq(L, R) as L - R. To control the evaluation of
         the result set pass `evaluate=True` to give L - R;
         if `evaluate=None` then terms in L and R will not cancel
         but they will be listed in canonical order; otherwise
@@ -624,6 +527,7 @@ def _eval_rewrite_as_Add(self, *args, **kwargs):
         >>> eq.rewrite(Add, evaluate=False).args
         (b, x, b, -x)
         """
+        from .add import _unevaluated_Add, Add
         L, R = args
         evaluate = kwargs.get('evaluate', True)
         if evaluate:
@@ -646,6 +550,7 @@ def binary_symbols(self):
         return set()
 
     def _eval_simplify(self, **kwargs):
+        from .add import Add
         from sympy.solvers.solveset import linear_coeffs
         # standard simplify
         e = super()._eval_simplify(**kwargs)
@@ -658,11 +563,11 @@ def _eval_simplify(self, **kwargs):
                 m, b = linear_coeffs(
                     e.rewrite(Add, evaluate=False), x)
                 if m.is_zero is False:
-                    enew = e.func(x, -b/m)
+                    enew = e.func(x, -b / m)
                 else:
-                    enew = e.func(m*x, -b)
+                    enew = e.func(m * x, -b)
                 measure = kwargs['measure']
-                if measure(enew) <= kwargs['ratio']*measure(e):
+                if measure(enew) <= kwargs['ratio'] * measure(e):
                     e = enew
             except ValueError:
                 pass
@@ -727,17 +632,13 @@ class Unequality(Relational):
     def __new__(cls, lhs, rhs, **options):
         lhs = _sympify(lhs)
         rhs = _sympify(rhs)
-
         evaluate = options.pop('evaluate', global_parameters.evaluate)
-
         if evaluate:
-            if isinstance(lhs, Boolean) != isinstance(lhs, Boolean):
-                # e.g. 0/1 not recognized as Boolean in SymPy
-                return S.true
-
-            is_equal = Equality(lhs, rhs)
-            if isinstance(is_equal, BooleanAtom):
-                return is_equal.negated
+            val = is_neq(lhs, rhs)
+            if val is None:
+                return cls(lhs, rhs, evaluate=False)
+            else:
+                return _sympify(val)
 
         return Relational.__new__(cls, lhs, rhs, **options)
 
@@ -776,12 +677,20 @@ class _Inequality(Relational):
     __slots__ = ()
 
     def __new__(cls, lhs, rhs, **options):
-        lhs = _sympify(lhs)
-        rhs = _sympify(rhs)
 
-        evaluate = options.pop('evaluate', global_parameters.evaluate)
+        try:
+            lhs = _sympify(lhs)
+            rhs = _sympify(rhs)
+        except SympifyError:
+            return NotImplemented
 
+        evaluate = options.pop('evaluate', global_parameters.evaluate)
         if evaluate:
+            for me in (lhs, rhs):
+                if me.is_extended_real is False:
+                    raise TypeError("Invalid comparison of non-real %s" % me)
+                if me is S.NaN:
+                    raise TypeError("Invalid NaN comparison")
             # First we invoke the appropriate inequality method of `lhs`
             # (e.g., `lhs.__lt__`).  That method will try to reduce to
             # boolean or raise an exception.  It may keep calling
@@ -790,16 +699,17 @@ def __new__(cls, lhs, rhs, **options):
             # nor a subclass was able to reduce to boolean or raise an
             # exception).  In that case, it must call us with
             # `evaluate=False` to prevent infinite recursion.
-            r = cls._eval_relation(lhs, rhs)
-            if r is not None:
-                return r
-            # Note: not sure r could be None, perhaps we never take this
-            # path?  In principle, could use this to shortcut out if a
-            # class realizes the inequality cannot be evaluated further.
+            return cls._eval_relation(lhs, rhs, **options)
 
         # make a "non-evaluated" Expr for the inequality
         return Relational.__new__(cls, lhs, rhs, **options)
 
+    @classmethod
+    def _eval_relation(cls, lhs, rhs, **options):
+        return cls._eval_fuzzy_relation(lhs, rhs)
+
+
+
 class _Greater(_Inequality):
     """Not intended for general use
 
@@ -1066,9 +976,12 @@ class GreaterThan(_Greater):
     rel_op = '>='
 
     @classmethod
-    def _eval_relation(cls, lhs, rhs):
-        # We don't use the op symbol here: workaround issue #7951
-        return _sympify(lhs.__ge__(rhs))
+    def _eval_fuzzy_relation(cls, lhs, rhs):
+        val = is_ge(lhs, rhs)
+        if val is None:
+            return cls(lhs, rhs, evaluate=False)
+        else:
+            return _sympify(val)
 
 
 Ge = GreaterThan
@@ -1081,9 +994,12 @@ class LessThan(_Less):
     rel_op = '<='
 
     @classmethod
-    def _eval_relation(cls, lhs, rhs):
-        # We don't use the op symbol here: workaround issue #7951
-        return _sympify(lhs.__le__(rhs))
+    def _eval_fuzzy_relation(cls, lhs, rhs):
+        val = is_le(lhs, rhs)
+        if val is None:
+            return cls(lhs, rhs, evaluate=False)
+        else:
+            return _sympify(val)
 
 
 Le = LessThan
@@ -1096,9 +1012,12 @@ class StrictGreaterThan(_Greater):
     rel_op = '>'
 
     @classmethod
-    def _eval_relation(cls, lhs, rhs):
-        # We don't use the op symbol here: workaround issue #7951
-        return _sympify(lhs.__gt__(rhs))
+    def _eval_fuzzy_relation(cls, lhs, rhs):
+        val = is_gt(lhs, rhs)
+        if val is None:
+            return cls(lhs, rhs, evaluate=False)
+        else:
+            return _sympify(val)
 
 
 Gt = StrictGreaterThan
@@ -1111,13 +1030,16 @@ class StrictLessThan(_Less):
     rel_op = '<'
 
     @classmethod
-    def _eval_relation(cls, lhs, rhs):
-        # We don't use the op symbol here: workaround issue #7951
-        return _sympify(lhs.__lt__(rhs))
+    def _eval_fuzzy_relation(cls, lhs, rhs):
+        val = is_lt(lhs, rhs)
+        if val is None:
+            return cls(lhs, rhs, evaluate=False)
+        else:
+            return _sympify(val)
 
 
-Lt = StrictLessThan
 
+Lt = StrictLessThan
 
 # A class-specific (not object-specific) data item used for a minor speedup.
 # It is defined here, rather than directly in the class, because the classes
@@ -1138,3 +1060,169 @@ def _eval_relation(cls, lhs, rhs):
     '<': StrictLessThan,
     'lt': StrictLessThan,
 }
+
+
+def _n2(a, b):
+    """Return (a - b).evalf(2) if a and b are comparable, else None.
+    This should only be used when a and b are already sympified.
+    """
+    # /!\ it is very important (see issue 8245) not to
+    # use a re-evaluated number in the calculation of dif
+    if a.is_comparable and b.is_comparable:
+        dif = (a - b).evalf(2)
+        if dif.is_comparable:
+            return dif
+
+
+@dispatch(Basic, Basic)
+def _eval_is_le(lhs, rhs):
+    return None
+
+
+@dispatch(Basic, Basic)
+def _eval_is_ge(lhs, rhs):
+    return None
+
+
+@dispatch(Basic, Basic)
+def _eval_Eq(lhs, rhs):
+    return None
+
+
+
+def is_lt(lhs, rhs):
+    return fuzzy_not(is_ge(lhs, rhs))
+
+
+def is_gt(lhs, rhs):
+    return fuzzy_not(is_le(lhs, rhs))
+
+
+def is_le(lhs, rhs):
+    return is_ge(rhs, lhs)
+
+
+def is_ge(lhs, rhs):
+    retval = _eval_is_ge(lhs, rhs)
+    if retval is not None:
+        return retval
+    else:
+        n2 = _n2(lhs, rhs)
+        if n2 is not None:
+            # use float comparison for infinity.
+            # otherwise get stuck in infinite recursion
+            if n2 in (S.Infinity, S.NegativeInfinity):
+                n2 = float(n2)
+            return _sympify(n2 >= 0)
+        if lhs.is_extended_real and rhs.is_extended_real:
+            if (lhs.is_infinite and lhs.is_extended_positive) \
+                     or (rhs.is_infinite and rhs.is_extended_negative):
+                return S.true
+            diff = lhs - rhs
+            if diff is not S.NaN:
+                rv = diff.is_extended_nonnegative
+                if rv is not None:
+                   return rv
+
+
+def is_neq(lhs, rhs):
+    from .logic import fuzzy_not
+    return fuzzy_not(is_eq(lhs, rhs))
+
+
+def is_eq(lhs, rhs):
+    from sympy.core.add import Add
+    from sympy.functions.elementary.complexes import arg
+    from sympy.simplify.simplify import clear_coefficients
+    from sympy.utilities.iterables import sift
+
+    retval = _eval_Eq(lhs, rhs)
+    if retval is not None:
+        return retval
+
+    retval = _eval_Eq(rhs, lhs)
+    if retval is not None:
+        return retval
+    else:
+        # retval is still None, so go through the equality logic
+        # If expressions have the same structure, they must be equal.
+        if lhs == rhs:
+            return S.true  # e.g. True == True
+        elif all(isinstance(i, BooleanAtom) for i in (rhs, lhs)):
+            return S.false # True != False
+        elif not (lhs.is_Symbol or rhs.is_Symbol) and (
+            isinstance(lhs, Boolean) !=
+            isinstance(rhs, Boolean)):
+            return S.false  # only Booleans can equal Booleans
+
+        if lhs.is_infinite or rhs.is_infinite:
+            if fuzzy_xor([lhs.is_infinite, rhs.is_infinite]):
+                return S.false
+            if fuzzy_xor([lhs.is_extended_real, rhs.is_extended_real]):
+                return S.false
+            if fuzzy_and([lhs.is_extended_real, rhs.is_extended_real]):
+                return fuzzy_xor([lhs.is_extended_positive, fuzzy_not(rhs.is_extended_positive)])
+
+            # Try to split real/imaginary parts and equate them
+            I = S.ImaginaryUnit
+
+            def split_real_imag(expr):
+                real_imag = lambda t: (
+                    'real' if t.is_extended_real else
+                    'imag' if (I * t).is_extended_real else None)
+                return sift(Add.make_args(expr), real_imag)
+
+            lhs_ri = split_real_imag(lhs)
+            if not lhs_ri[None]:
+                rhs_ri = split_real_imag(rhs)
+                if not rhs_ri[None]:
+                    eq_real = Eq(Add(*lhs_ri['real']), Add(*rhs_ri['real']))
+                    eq_imag = Eq(I * Add(*lhs_ri['imag']), I * Add(*rhs_ri['imag']))
+                    return fuzzy_and(map(fuzzy_bool, [eq_real, eq_imag]))
+
+            # Compare e.g. zoo with 1+I*oo by comparing args
+            arglhs = arg(lhs)
+            argrhs = arg(rhs)
+            # Guard against Eq(nan, nan) -> False
+            if not (arglhs == S.NaN and argrhs == S.NaN):
+                return fuzzy_bool(Eq(arglhs, argrhs))
+
+        if all(isinstance(i, Expr) for i in (lhs, rhs)):
+            # see if the difference evaluates
+            dif = lhs - rhs
+            z = dif.is_zero
+            if z is not None:
+                if z is False and dif.is_commutative:  # issue 10728
+                    return S.false
+                if z:
+                    return S.true
+
+            n2 = _n2(lhs, rhs)
+            if n2 is not None:
+                return _sympify(n2 == 0)
+
+            # see if the ratio evaluates
+            n, d = dif.as_numer_denom()
+            rv = None
+            if n.is_zero:
+                rv = d.is_nonzero
+            elif n.is_finite:
+                if d.is_infinite:
+                    rv = True
+                elif n.is_zero is False:
+                    rv = d.is_infinite
+                    if rv is None:
+                        # if the condition that makes the denominator
+                        # infinite does not make the original expression
+                        # True then False can be returned
+                        l, r = clear_coefficients(d, S.Infinity)
+                        args = [_.subs(l, r) for _ in (lhs, rhs)]
+                        if args != [lhs, rhs]:
+                            rv = fuzzy_bool(Eq(*args))
+                            if rv is True:
+                                rv = None
+            elif any(a.is_infinite for a in Add.make_args(n)):
+                # (inf or nan)/x != 0
+                rv = S.false
+            if rv is not None:
+                return rv
diff --git a/sympy/core/tests/test_relational.py b/sympy/core/tests/test_relational.py
index 915af7fe4b8f..4ee5ab4810af 100644
--- a/sympy/core/tests/test_relational.py
+++ b/sympy/core/tests/test_relational.py
@@ -690,7 +690,7 @@ def test_issue_8245():
     r = Rational(str(a.n(29)))
     assert rel_check(a, r)
 
-    assert Eq(log(cos(2)**2 + sin(2)**2), 0) == True
+    assert Eq(log(cos(2)**2 + sin(2)**2), 0) is S.true
 
 
 def test_issue_8449():
diff --git a/sympy/functions/elementary/integers.py b/sympy/functions/elementary/integers.py
index 960121c8aa73..ac91311c1715 100644
--- a/sympy/functions/elementary/integers.py
+++ b/sympy/functions/elementary/integers.py
@@ -1,3 +1,8 @@
+
+from __future__ import print_function, division
+
+from sympy import Basic
+
 from sympy.core import Add, S
 from sympy.core.evalf import get_integer_part, PrecisionExhausted
 from sympy.core.function import Function
@@ -6,7 +11,7 @@
 from sympy.core.relational import Gt, Lt, Ge, Le, Relational
 from sympy.core.symbol import Symbol
 from sympy.core.sympify import _sympify
-
+from sympy.multipledispatch import dispatch
 
 ###############################################################################
 ######################### FLOOR and CEILING FUNCTIONS #########################
@@ -154,11 +159,7 @@ def _eval_rewrite_as_ceiling(self, arg, **kwargs):
     def _eval_rewrite_as_frac(self, arg, **kwargs):
         return arg - frac(arg)
 
-    def _eval_Eq(self, other):
-        if isinstance(self, floor):
-            if (self.rewrite(ceiling) == other) or \
-                    (self.rewrite(frac) == other):
-                return S.true
+
 
     def __le__(self, other):
         other = S(other)
@@ -216,6 +217,12 @@ def __lt__(self, other):
 
         return Lt(self, other, evaluate=False)
 
+@dispatch(floor, Basic)
+def _eval_Eq(lhs, rhs): # noqa:F811
+    if (lhs.rewrite(ceiling) == rhs) or \
+        (lhs.rewrite(frac) == rhs):
+        return S.true
+
 class ceiling(RoundFunction):
     """
     Ceiling is a univariate function which returns the smallest integer
@@ -289,11 +296,7 @@ def _eval_is_positive(self):
     def _eval_is_nonpositive(self):
         return self.args[0].is_nonpositive
 
-    def _eval_Eq(self, other):
-        if isinstance(self, ceiling):
-            if (self.rewrite(floor) == other) or \
-                    (self.rewrite(frac) == other):
-                return S.true
+
 
     def __lt__(self, other):
         other = S(other)
@@ -351,6 +354,12 @@ def __le__(self, other):
 
         return Le(self, other, evaluate=False)
 
+@dispatch(ceiling, Basic)
+def _eval_Eq(lhs, rhs): # noqa:F811
+    if (lhs.rewrite(floor) == rhs) or \
+        (lhs.rewrite(frac) == rhs):
+        return S.true
+
 class frac(Function):
     r"""Represents the fractional part of x
 
@@ -442,18 +451,7 @@ def _eval_rewrite_as_floor(self, arg, **kwargs):
     def _eval_rewrite_as_ceiling(self, arg, **kwargs):
         return arg + ceiling(-arg)
 
-    def _eval_Eq(self, other):
-        if isinstance(self, frac):
-            if (self.rewrite(floor) == other) or \
-                    (self.rewrite(ceiling) == other):
-                return S.true
-            # Check if other < 0
-            if other.is_extended_negative:
-                return S.false
-            # Check if other >= 1
-            res = self._value_one_or_more(other)
-            if res is not None:
-                return S.false
+
 
     def _eval_is_finite(self):
         return True
@@ -529,3 +527,16 @@ def _value_one_or_more(self, other):
                     return S.true
             if other.is_integer and other.is_positive:
                 return S.true
+
+@dispatch(frac, Basic)
+def _eval_Eq(lhs, rhs): # noqa:F811
+    if (lhs.rewrite(floor) == rhs) or \
+        (lhs.rewrite(ceiling) == rhs):
+        return S.true
+    # Check if other < 0
+    if rhs.is_extended_negative:
+        return S.false
+    # Check if other >= 1
+    res = lhs._value_one_or_more(rhs)
+    if res is not None:
+        return S.false
diff --git a/sympy/matrices/expressions/matexpr.py b/sympy/matrices/expressions/matexpr.py
index 31f1226cce69..f1c1f021f499 100644
--- a/sympy/matrices/expressions/matexpr.py
+++ b/sympy/matrices/expressions/matexpr.py
@@ -14,6 +14,7 @@
 from sympy.simplify import simplify
 from sympy.utilities.misc import filldedent
 from sympy.assumptions.ask import ask, Q
+from sympy.multipledispatch import dispatch
 
 
 def _sympifyit(arg, retval=None):
@@ -593,15 +594,16 @@ def applyfunc(self, func):
         from .applyfunc import ElementwiseApplyFunction
         return ElementwiseApplyFunction(func, self)
 
-    def _eval_Eq(self, other):
-        if not isinstance(other, MatrixExpr):
-            return False
-        if self.shape != other.shape:
-            return False
-        if (self - other).is_ZeroMatrix:
-            return True
-        return Eq(self, other, evaluate=False)
+@dispatch(MatrixExpr, Basic)
+def _eval_Eq(lhs, rhs): # noqa:F811
+    return False
 
+@dispatch(MatrixExpr, MatrixExpr)
+def _eval_Eq(self, other): # noqa:F811
+    if self.shape != other.shape:
+        return False
+    if (self - other).is_ZeroMatrix:
+        return True
 
 def get_postprocessor(cls):
     def _postprocessor(expr):
diff --git a/sympy/matrices/immutable.py b/sympy/matrices/immutable.py
index b559180b3ac2..175a34a0d660 100644
--- a/sympy/matrices/immutable.py
+++ b/sympy/matrices/immutable.py
@@ -7,20 +7,26 @@
 from sympy.matrices.expressions import MatrixExpr
 from sympy.matrices.matrices import MatrixBase
 from sympy.matrices.sparse import MutableSparseMatrix, SparseMatrix
+from sympy.multipledispatch import dispatch
 
 
 def sympify_matrix(arg):
     return arg.as_immutable()
+
+
 sympify_converter[MatrixBase] = sympify_matrix
 
 
 def sympify_mpmath_matrix(arg):
     mat = [_sympify(x) for x in arg]
     return ImmutableDenseMatrix(arg.rows, arg.cols, mat)
+
+
 sympify_converter[_matrix] = sympify_mpmath_matrix
 
 
-class ImmutableDenseMatrix(DenseMatrix, MatrixExpr):
+
+class ImmutableDenseMatrix(DenseMatrix, MatrixExpr):  # type: ignore
     """Create an immutable version of a matrix.
 
     Examples
@@ -77,27 +83,6 @@ def _entry(self, i, j, **kwargs):
     def __setitem__(self, *args):
         raise TypeError("Cannot set values of {}".format(self.__class__))
 
-    def _eval_Eq(self, other):
-        """Helper method for Equality with matrices.
-
-        Relational automatically converts matrices to ImmutableDenseMatrix
-        instances, so this method only applies here.  Returns True if the
-        matrices are definitively the same, False if they are definitively
-        different, and None if undetermined (e.g. if they contain Symbols).
-        Returning None triggers default handling of Equalities.
-
-        """
-        if not hasattr(other, 'shape') or self.shape != other.shape:
-            return S.false
-        if isinstance(other, MatrixExpr) and not isinstance(
-                other, ImmutableDenseMatrix):
-            return None
-        diff = (self - other).is_zero_matrix
-        if diff is True:
-            return S.true
-        elif diff is False:
-            return S.false
-
     def _eval_extract(self, rowsList, colsList):
         # self._mat is a Tuple.  It is slightly faster to index a
         # tuple over a Tuple, so grab the internal tuple directly
@@ -125,6 +110,7 @@ def as_immutable(self):
     def is_diagonalizable(self, reals_only=False, **kwargs):
         return super().is_diagonalizable(
             reals_only=reals_only, **kwargs)
+
     is_diagonalizable.__doc__ = DenseMatrix.is_diagonalizable.__doc__
     is_diagonalizable = cacheit(is_diagonalizable)
 
@@ -179,7 +165,7 @@ def __setitem__(self, *args):
     def _entry(self, i, j, **kwargs):
         return SparseMatrix.__getitem__(self, (i, j))
 
-    _eval_Eq = ImmutableDenseMatrix._eval_Eq
+    # _eval_Eq = ImmutableDenseMatrix._eval_Eq
 
     @property
     def cols(self):
@@ -199,5 +185,27 @@ def as_immutable(self):
     def is_diagonalizable(self, reals_only=False, **kwargs):
         return super().is_diagonalizable(
             reals_only=reals_only, **kwargs)
+
     is_diagonalizable.__doc__ = SparseMatrix.is_diagonalizable.__doc__
     is_diagonalizable = cacheit(is_diagonalizable)
+
+
+@dispatch(ImmutableDenseMatrix, MatrixExpr)
+def _eval_Eq(lhs, rhs): # noqa:F811
+    return None
+
+
+@dispatch(ImmutableDenseMatrix, ImmutableDenseMatrix)
+def _eval_Eq(lhs, rhs): # noqa:F811
+    """Helper method for Equality with matrices.sympy.
+
+    Relational automatically converts matrices to ImmutableDenseMatrix
+    instances, so this method only applies here.  Returns True if the
+    matrices are definitively the same, False if they are definitively
+    different, and None if undetermined (e.g. if they contain Symbols).
+    Returning None triggers default handling of Equalities.
+
+    """
+    if lhs.shape != rhs.shape:
+        return S.false
+    return (lhs - rhs).is_zero_matrix
diff --git a/sympy/polys/domains/finitefield.py b/sympy/polys/domains/finitefield.py
index 0b1cceee49a1..d60bbd990b24 100644
--- a/sympy/polys/domains/finitefield.py
+++ b/sympy/polys/domains/finitefield.py
@@ -3,7 +3,7 @@
 from __future__ import print_function, division
 
 from sympy.polys.domains.field import Field
-from sympy.polys.domains.groundtypes import SymPyInteger
+
 from sympy.polys.domains.modularinteger import ModularIntegerFactory
 from sympy.polys.domains.simpledomain import SimpleDomain
 from sympy.polys.polyerrors import CoercionFailed
@@ -57,6 +57,7 @@ def get_field(self):
         return self
 
     def to_sympy(self, a):
+        from sympy.polys.domains.groundtypes import SymPyInteger
         """Convert ``a`` to a SymPy object. """
         return SymPyInteger(int(a))
 
diff --git a/sympy/polys/rootoftools.py b/sympy/polys/rootoftools.py
index 8d8a199ff826..e0ad45663ac4 100644
--- a/sympy/polys/rootoftools.py
+++ b/sympy/polys/rootoftools.py
@@ -2,6 +2,7 @@
 
 from __future__ import print_function, division
 
+from sympy import Basic
 from sympy.core import (S, Expr, Integer, Float, I, oo, Add, Lambda,
     symbols, sympify, Rational, Dummy)
 from sympy.core.cache import cacheit
@@ -25,7 +26,7 @@
 
 from mpmath import mpf, mpc, findroot, workprec
 from mpmath.libmp.libmpf import dps_to_prec, prec_to_dps
-
+from sympy.multipledispatch import dispatch
 from itertools import chain
 
 
@@ -970,42 +971,49 @@ def eval_rational(self, dx=None, dy=None, n=15):
         self._set_interval(interval)
         return real + I*imag
 
-    def _eval_Eq(self, other):
-        # CRootOf represents a Root, so if other is that root, it should set
-        # the expression to zero *and* it should be in the interval of the
-        # CRootOf instance. It must also be a number that agrees with the
-        # is_real value of the CRootOf instance.
-        if type(self) == type(other):
-            return sympify(self == other)
-        if not other.is_number:
-            return None
-        if not other.is_finite:
-            return S.false
-        z = self.expr.subs(self.expr.free_symbols.pop(), other).is_zero
-        if z is False:  # all roots will make z True but we don't know
-                        # whether this is the right root if z is True
-            return S.false
-        o = other.is_real, other.is_imaginary
-        s = self.is_real, self.is_imaginary
-        assert None not in s  # this is part of initial refinement
-        if o != s and None not in o:
-            return S.false
-        re, im = other.as_real_imag()
-        if self.is_real:
-            if im:
-                return S.false
-            i = self._get_interval()
-            a, b = [Rational(str(_)) for _ in (i.a, i.b)]
-            return sympify(a <= other and other <= b)
-        i = self._get_interval()
-        r1, r2, i1, i2 = [Rational(str(j)) for j in (
-            i.ax, i.bx, i.ay, i.by)]
-        return sympify((
-            r1 <= re and re <= r2) and (
-            i1 <= im and im <= i2))
 
 CRootOf = ComplexRootOf
 
+
+@dispatch(ComplexRootOf, ComplexRootOf)
+def _eval_Eq(lhs, rhs): # noqa:F811
+    return sympify(lhs == rhs)
+
+
+@dispatch(ComplexRootOf, Basic)
+def _eval_Eq(lhs, rhs): # noqa:F811
+    # CRootOf represents a Root, so if rhs is that root, it should set
+    # the expression to zero *and* it should be in the interval of the
+    # CRootOf instance. It must also be a number that agrees with the
+    # is_real value of the CRootOf instance.
+    if not rhs.is_number:
+        return None
+    if not rhs.is_finite:
+        return S.false
+    z = lhs.expr.subs(lhs.expr.free_symbols.pop(), rhs).is_zero
+    if z is False:  # all roots will make z True but we don't know
+        # whether this is the right root if z is True
+        return S.false
+    o = rhs.is_real, rhs.is_imaginary
+    s = lhs.is_real, lhs.is_imaginary
+    assert None not in s  # this is part of initial refinement
+    if o != s and None not in o:
+        return S.false
+    re, im = rhs.as_real_imag()
+    if lhs.is_real:
+        if im:
+            return S.false
+        i = lhs._get_interval()
+        a, b = [Rational(str(_)) for _ in (i.a, i.b)]
+        return sympify(a <= rhs and rhs <= b)
+    i = lhs._get_interval()
+    r1, r2, i1, i2 = [Rational(str(j)) for j in (
+        i.ax, i.bx, i.ay, i.by)]
+    return sympify((
+                       r1 <= re and re <= r2) and (
+                       i1 <= im and im <= i2))
+
+
 @public
 class RootSum(Expr):
     """Represents a sum of all roots of a univariate polynomial. """
diff --git a/sympy/sets/sets.py b/sympy/sets/sets.py
index 64d52737d9f5..49abd73462f2 100644
--- a/sympy/sets/sets.py
+++ b/sympy/sets/sets.py
@@ -27,7 +27,7 @@
 from sympy.utilities.exceptions import SymPyDeprecationWarning
 from sympy.utilities.iterables import iproduct, sift, roundrobin
 from sympy.utilities.misc import func_name, filldedent
-
+from sympy.multipledispatch import dispatch
 from mpmath import mpi, mpf
 
 
@@ -757,15 +757,7 @@ def _flatten(sets):
                     yield s
         return ProductSet(*_flatten(self.sets))
 
-    def _eval_Eq(self, other):
-        if not other.is_ProductSet:
-            return
 
-        if len(self.sets) != len(other.sets):
-            return false
-
-        eqs = (Eq(x, y) for x, y in zip(self.sets, other.sets))
-        return tfn[fuzzy_and(map(fuzzy_bool, eqs))]
 
     def _contains(self, element):
         """
@@ -920,7 +912,7 @@ def __new__(cls, start, end, left_open=False, right_open=False):
             raise ValueError("Non-real intervals are not supported")
 
         # evaluate if possible
-        if (end < start) == True:
+        if (end >= start) == False:
             return S.EmptySet
         elif (end - start).is_negative:
             return S.EmptySet
@@ -1124,11 +1116,6 @@ def _eval_Eq(self, other):
                 return None
             return false
 
-        return And(Eq(self.left, other.left),
-                   Eq(self.right, other.right),
-                   self.left_open == other.left_open,
-                   self.right_open == other.right_open)
-
 
 class Union(Set, LatticeOp, EvalfMixin):
     """
@@ -1787,24 +1774,6 @@ def __new__(cls, *args, **kwargs):
         obj._args_set = _args_set
         return obj
 
-    def _eval_Eq(self, other):
-        if not isinstance(other, FiniteSet):
-            # XXX: If Interval(x, x, evaluate=False) worked then the line
-            # below would mean that
-            #     FiniteSet(x) & Interval(x, x, evaluate=False) -> false
-            if isinstance(other, Interval):
-                return false
-            elif isinstance(other, Set):
-                return None
-            return false
-
-        def all_in_both():
-            s_set = set(self.args)
-            o_set = set(other.args)
-            yield fuzzy_and(self._contains(e) for e in o_set - s_set)
-            yield fuzzy_and(other._contains(e) for e in s_set - o_set)
-
-        return tfn[fuzzy_and(all_in_both())]
 
     def __iter__(self):
         return iter(self.args)
@@ -2516,3 +2485,46 @@ def set_pow(x, y):
 def set_function(f, x):
     from sympy.sets.handlers.functions import _set_function
     return _set_function(f, x)
+
+
+@dispatch(Interval, FiniteSet)
+def _eval_Eq(lhs, rhs): # noqa: F811
+    return false
+
+
+@dispatch(FiniteSet, Interval)
+def _eval_Eq(lhs, rhs): # noqa: F811
+    return false
+
+
+@dispatch(Interval, Interval)
+def _eval_Eq(lhs, rhs): # noqa: F811
+    return And(Eq(lhs.left, rhs.left),
+               Eq(lhs.right, rhs.right),
+               lhs.left_open == rhs.left_open,
+               lhs.right_open == rhs.right_open)
+
+
+@dispatch(FiniteSet, Interval)
+def _eval_Eq(lhs, rhs): # noqa: F811
+    return False
+
+
+@dispatch(FiniteSet, FiniteSet)
+def _eval_Eq(lhs, rhs): # noqa: F811
+    def all_in_both():
+        s_set = set(lhs.args)
+        o_set = set(rhs.args)
+        yield fuzzy_and(lhs._contains(e) for e in o_set - s_set)
+        yield fuzzy_and(rhs._contains(e) for e in s_set - o_set)
+
+    return tfn[fuzzy_and(all_in_both())]
+
+
+@dispatch(ProductSet, ProductSet)
+def _eval_Eq(lhs, rhs): # noqa: F811
+    if len(lhs.sets) != len(rhs.sets):
+        return false
+
+    eqs = (Eq(x, y) for x, y in zip(lhs.sets, rhs.sets))
+    return tfn[fuzzy_and(map(fuzzy_bool, eqs))]
diff --git a/sympy/tensor/indexed.py b/sympy/tensor/indexed.py
index b8b4da9461e7..30deb7a6d94d 100644
--- a/sympy/tensor/indexed.py
+++ b/sympy/tensor/indexed.py
@@ -106,6 +106,7 @@
 
 from __future__ import print_function, division
 
+from sympy import Number
 from sympy.core.assumptions import StdFactKB
 from sympy.core import Expr, Tuple, sympify, S
 from sympy.core.symbol import _filter_assumptions, Symbol
@@ -114,6 +115,7 @@
 from sympy.core.logic import fuzzy_bool, fuzzy_not
 from sympy.core.sympify import _sympify
 from sympy.functions.special.tensor_functions import KroneckerDelta
+from sympy.multipledispatch import dispatch
 
 
 class IndexException(Exception):
@@ -751,58 +753,41 @@ def name(self):
     def free_symbols(self):
         return {self}
 
-    def __le__(self, other):
-        if isinstance(other, Idx):
-            other_upper = other if other.upper is None else other.upper
-            other_lower = other if other.lower is None else other.lower
-        else:
-            other_upper = other
-            other_lower = other
-
-        if self.upper is not None and (self.upper <= other_lower) == True:
-            return True
-        if self.lower is not None and (self.lower > other_upper) == True:
-            return False
-        return super(Idx, self).__le__(other)
-
-    def __ge__(self, other):
-        if isinstance(other, Idx):
-            other_upper = other if other.upper is None else other.upper
-            other_lower = other if other.lower is None else other.lower
-        else:
-            other_upper = other
-            other_lower = other
-
-        if self.lower is not None and (self.lower >= other_upper) == True:
-            return True
-        if self.upper is not None and (self.upper < other_lower) == True:
-            return False
-        return super(Idx, self).__ge__(other)
-
-    def __lt__(self, other):
-        if isinstance(other, Idx):
-            other_upper = other if other.upper is None else other.upper
-            other_lower = other if other.lower is None else other.lower
-        else:
-            other_upper = other
-            other_lower = other
-
-        if self.upper is not None and (self.upper < other_lower) == True:
-            return True
-        if self.lower is not None and (self.lower >= other_upper) == True:
-            return False
-        return super(Idx, self).__lt__(other)
-
-    def __gt__(self, other):
-        if isinstance(other, Idx):
-            other_upper = other if other.upper is None else other.upper
-            other_lower = other if other.lower is None else other.lower
-        else:
-            other_upper = other
-            other_lower = other
-
-        if self.lower is not None and (self.lower > other_upper) == True:
-            return True
-        if self.upper is not None and (self.upper <= other_lower) == True:
-            return False
-        return super(Idx, self).__gt__(other)
+
+@dispatch(Idx, Idx)
+def _eval_is_ge(lhs, rhs): # noqa:F811
+
+    other_upper = rhs if rhs.upper is None else rhs.upper
+    other_lower = rhs if rhs.lower is None else rhs.lower
+
+    if lhs.lower is not None and (lhs.lower >= other_upper) == True:
+        return True
+    if lhs.upper is not None and (lhs.upper < other_lower) == True:
+        return False
+    return None
+
+
+@dispatch(Idx, Number)
+def _eval_is_ge(lhs, rhs): # noqa:F811
+
+    other_upper = rhs
+    other_lower = rhs
+
+    if lhs.lower is not None and (lhs.lower >= other_upper) == True:
+        return True
+    if lhs.upper is not None and (lhs.upper < other_lower) == True:
+        return False
+    return None
+
+
+@dispatch(Number, Idx)
+def _eval_is_ge(lhs, rhs): # noqa:F811
+
+    other_upper = lhs
+    other_lower = lhs
+
+    if rhs.upper is not None and (rhs.upper <= other_lower) == True:
+        return True
+    if rhs.lower is not None and (rhs.lower > other_upper) == True:
+        return False
+    return None


Code review comments:
diff_hunk: @@ -141,25 +141,27 @@ def index(self, value, start=None, stop=None):
         else:
             return self.args.index(value, start, stop)
 
-    def _eval_Eq(self, other):
-        from sympy.core.function import AppliedUndef
-        from sympy.core.logic import fuzzy_and, fuzzy_bool
-        from sympy.core.relational import Eq
 
-        if other.is_Symbol or isinstance(other, AppliedUndef):
-            return None
+converter[tuple] = lambda tup: Tuple(*tup)
+
 
-        if not isinstance(other, Tuple) or len(self) != len(other):
-            return S.false
+@dispatch(Tuple, Basic)
+def _eval_Eq(lhs, rhs):
+    from sympy.core.function import AppliedUndef
+    from sympy.core.logic import fuzzy_and, fuzzy_bool
+    from sympy.core.relational import Eq
 
-        r = fuzzy_and(fuzzy_bool(Eq(s, o)) for s, o in zip(self, other))
-        if r is True:
-            return S.true
-        elif r is False:
-            return S.false
+    if rhs.is_Symbol or isinstance(rhs, AppliedUndef):
+        return None
path: sympy/core/containers.py
text: @user1:
Maybe it would make sense to do the type checking in dispatch by having one handler for `Symbol` and one for `AppliedUndef` and then another for `Tuple`.

@author:
Is every object with is_Symbol guaranteed to be of type Symbol? I think there is one case, I can find out more details for you, but if I do find that there is class that does not subclass Symbol but has is_Symbol set to true should I edit to subclass Symbol

@user1:
I'm not sure. What is the example?

@author:
nevermind, it did actually work.

diff_hunk: @@ -141,25 +141,27 @@ def index(self, value, start=None, stop=None):
         else:
             return self.args.index(value, start, stop)
 
-    def _eval_Eq(self, other):
-        from sympy.core.function import AppliedUndef
-        from sympy.core.logic import fuzzy_and, fuzzy_bool
-        from sympy.core.relational import Eq
 
-        if other.is_Symbol or isinstance(other, AppliedUndef):
-            return None
+converter[tuple] = lambda tup: Tuple(*tup)
+
 
-        if not isinstance(other, Tuple) or len(self) != len(other):
-            return S.false
+@dispatch(Tuple, Basic)
+def _eval_Eq(lhs, rhs):
+    from sympy.core.function import AppliedUndef
+    from sympy.core.logic import fuzzy_and, fuzzy_bool
+    from sympy.core.relational import Eq
 
-        r = fuzzy_and(fuzzy_bool(Eq(s, o)) for s, o in zip(self, other))
-        if r is True:
-            return S.true
-        elif r is False:
-            return S.false
+    if rhs.is_Symbol or isinstance(rhs, AppliedUndef):
+        return None
 
+    if not isinstance(rhs, Tuple) or len(lhs) != len(rhs):
+        return S.false
 
-converter[tuple] = lambda tup: Tuple(*tup)
+    r = fuzzy_and(fuzzy_bool(Eq(s, o)) for s, o in zip(lhs, rhs))
path: sympy/core/containers.py
text: @user1:
Should this use `is_eq` rather than `Eq`?

@author:
Yes.

diff_hunk: @@ -1,3 +1,4 @@
+from __future__ import print_function, division
path: sympy/core/expr.py
text: @user1:
We don't need future imports for py3 things as python 2 isn't supported any more.

diff_hunk: @@ -3020,3 +3012,10 @@ def __new__(cls, perm, x, evaluate=None):
 
         obj = super().__new__(cls, perm, x)
         return obj
+
+
+@dispatch(Permutation, Permutation)
+def _eval_Eq(lhs, rhs):
+    if lhs._size != rhs._size:
+        return None
+    return as_Boolean(lhs._array_form == rhs._array_form)
path: sympy/combinatorics/permutations.py
text: @user1:
I wonder if `_eval_Eq` should return a fuzzy bool rather than a `Boolean`...

@author:
that makes sense to me

diff_hunk: @@ -276,6 +276,8 @@ def _igcd2_python(a, b):
 
 
 # Use Lehmer's algorithm only for very large numbers.
+# The limit could be different on Python 2.7 and 3.x.
+# If so, then this could be defined in compatibility.py.
path: sympy/core/numbers.py
text: @user1:
I think this is a hangover from a merge. These comments were probably removed because Python 2 is no longer supported.

diff_hunk: @@ -339,78 +344,30 @@ def __complex__(self):
         re, im = result.as_real_imag()
         return complex(float(re), float(im))
 
-    def _cmp(self, other, op, cls):
-        assert op in ("<", ">", "<=", ">=")
-        try:
-            other = _sympify(other)
-        except SympifyError:
-            return NotImplemented
-
-        if not isinstance(other, Expr):
-            return NotImplemented
-
-        for me in (self, other):
-            if me.is_extended_real is False:
-                raise TypeError("Invalid comparison of non-real %s" % me)
-            if me is S.NaN:
-                raise TypeError("Invalid NaN comparison")
-
-        n2 = _n2(self, other)
-        if n2 is not None:
-            # use float comparison for infinity.
-            # otherwise get stuck in infinite recursion
-            if n2 in (S.Infinity, S.NegativeInfinity):
-                n2 = float(n2)
-            if op == "<":
-                return _sympify(n2 < 0)
-            elif op == ">":
-                return _sympify(n2 > 0)
-            elif op == "<=":
-                return _sympify(n2 <= 0)
-            else: # >=
-                return _sympify(n2 >= 0)
-
-        if self.is_extended_real and other.is_extended_real:
-            if op in ("<=", ">") \
-                and ((self.is_infinite and self.is_extended_negative) \
-                     or (other.is_infinite and other.is_extended_positive)):
-                return S.true if op == "<=" else S.false
-            if op in ("<", ">=") \
-                and ((self.is_infinite and self.is_extended_positive) \
-                     or (other.is_infinite and other.is_extended_negative)):
-                return S.true if op == ">=" else S.false
-            diff = self - other
-            if diff is not S.NaN:
-                if op == "<":
-                    test = diff.is_extended_negative
-                elif op == ">":
-                    test = diff.is_extended_positive
-                elif op == "<=":
-                    test = diff.is_extended_nonpositive
-                else: # >=
-                    test = diff.is_extended_nonnegative
-
-                if test is not None:
-                    return S.true if test == True else S.false
-
-        # return unevaluated comparison object
-        return cls(self, other, evaluate=False)
-
     def __ge__(self, other):
-        from sympy import GreaterThan
-        return self._cmp(other, ">=", GreaterThan)
+        from .relational import GreaterThan
+        return self._cmp(other, GreaterThan)
path: sympy/core/expr.py
text: @user1:
If the purpose of the `_cmp` method is just to sympy the argument can we use `sympify_return` instead?

That's what's used for e.g. `Expr.__add__`.

diff_hunk: @@ -1138,3 +1060,169 @@ def _eval_relation(cls, lhs, rhs):
     '<': StrictLessThan,
     'lt': StrictLessThan,
 }
+
+
+def _n2(a, b):
+    """Return (a - b).evalf(2) if a and b are comparable, else None.
+    This should only be used when a and b are already sympified.
+    """
+    # /!\ it is very important (see issue 8245) not to
+    # use a re-evaluated number in the calculation of dif
+    if a.is_comparable and b.is_comparable:
+        dif = (a - b).evalf(2)
+        if dif.is_comparable:
+            return dif
+
+
+@dispatch(Basic, Basic)
+def _eval_is_le(lhs, rhs):
+    return None
+
+
+@dispatch(Basic, Basic)
+def _eval_is_ge(lhs, rhs):
+    return None
+
+
+@dispatch(Basic, Basic)
+def _eval_Eq(lhs, rhs):
+    return None
+
+
+
+def is_lt(lhs, rhs):
+    return fuzzy_not(is_ge(lhs, rhs))
+
+
+def is_gt(lhs, rhs):
+    return fuzzy_not(is_le(lhs, rhs))
+
+
+def is_le(lhs, rhs):
+    return is_ge(rhs, lhs)
+
+
+def is_ge(lhs, rhs):
+    retval = _eval_is_ge(lhs, rhs)
+    if retval is not None:
+        return retval
+    else:
+        n2 = _n2(lhs, rhs)
+        if n2 is not None:
+            # use float comparison for infinity.
+            # otherwise get stuck in infinite recursion
+            if n2 in (S.Infinity, S.NegativeInfinity):
+                n2 = float(n2)
+            return _sympify(n2 >= 0)
+        if lhs.is_extended_real and rhs.is_extended_real:
+            if (lhs.is_infinite and lhs.is_extended_positive) \
+                     or (rhs.is_infinite and rhs.is_extended_negative):
+                return S.true
+            diff = lhs - rhs
+            if diff is not S.NaN:
+                rv = diff.is_extended_nonnegative
+                if rv is not None:
+                   return rv
+
+
+def is_neq(lhs, rhs):
+    from .logic import fuzzy_not
+    return fuzzy_not(is_eq(lhs, rhs))
path: sympy/core/relational.py
text: @user1:
I think `fuzzy_not` is already imported.

diff_hunk: @@ -1066,9 +976,12 @@ class GreaterThan(_Greater):
     rel_op = '>='
 
     @classmethod
-    def _eval_relation(cls, lhs, rhs):
-        # We don't use the op symbol here: workaround issue #7951
-        return _sympify(lhs.__ge__(rhs))
+    def _eval_fuzzy_relation(cls, lhs, rhs):
+        val = is_ge(lhs, rhs)
+        if val is None:
path: sympy/core/relational.py
text: @user1:
There are a few of these `_eval_fuzzy_relation` methods and they differ only in which function they call (here `is_ge` for evaluating). Perhaps it would make more sense to move the other logic to the superclass `_eval_relation` method keeping only the call to `is_ge` here since that's the only bit that it unique to this class...



@author:
Yes, I agree that would be quite a bite better

diff_hunk: @@ -1138,3 +1060,169 @@ def _eval_relation(cls, lhs, rhs):
     '<': StrictLessThan,
     'lt': StrictLessThan,
 }
+
+
+def _n2(a, b):
+    """Return (a - b).evalf(2) if a and b are comparable, else None.
+    This should only be used when a and b are already sympified.
+    """
+    # /!\ it is very important (see issue 8245) not to
+    # use a re-evaluated number in the calculation of dif
+    if a.is_comparable and b.is_comparable:
+        dif = (a - b).evalf(2)
+        if dif.is_comparable:
+            return dif
+
+
+@dispatch(Basic, Basic)
+def _eval_is_le(lhs, rhs):
+    return None
+
+
+@dispatch(Basic, Basic)
+def _eval_is_ge(lhs, rhs):
+    return None
+
+
+@dispatch(Basic, Basic)
+def _eval_Eq(lhs, rhs):
+    return None
+
+
+
+def is_lt(lhs, rhs):
+    return fuzzy_not(is_ge(lhs, rhs))
+
+
+def is_gt(lhs, rhs):
+    return fuzzy_not(is_le(lhs, rhs))
+
+
+def is_le(lhs, rhs):
+    return is_ge(rhs, lhs)
+
+
+def is_ge(lhs, rhs):
+    retval = _eval_is_ge(lhs, rhs)
+    if retval is not None:
+        return retval
+    else:
+        n2 = _n2(lhs, rhs)
+        if n2 is not None:
+            # use float comparison for infinity.
+            # otherwise get stuck in infinite recursion
+            if n2 in (S.Infinity, S.NegativeInfinity):
+                n2 = float(n2)
+            return _sympify(n2 >= 0)
+        if lhs.is_extended_real and rhs.is_extended_real:
+            if (lhs.is_infinite and lhs.is_extended_positive) \
+                     or (rhs.is_infinite and rhs.is_extended_negative):
+                return S.true
+            diff = lhs - rhs
+            if diff is not S.NaN:
+                rv = diff.is_extended_nonnegative
+                if rv is not None:
+                   return rv
+
+
+def is_neq(lhs, rhs):
+    from .logic import fuzzy_not
+    return fuzzy_not(is_eq(lhs, rhs))
+
+
+def is_eq(lhs, rhs):
+    from sympy.core.add import Add
+    from sympy.functions.elementary.complexes import arg
+    from sympy.simplify.simplify import clear_coefficients
+    from sympy.utilities.iterables import sift
+
+    retval = _eval_Eq(lhs, rhs)
+    if retval is not None:
+        return retval
path: sympy/core/relational.py
text: @user1:
I think this method will still need to check e.g. `lhs._eval_Eq` for backward compatibility. Users and downstream libraries have probably defined `_eval_Eq` methods and we need to try to ensure that their code would still work with these changes.

diff_hunk: @@ -216,6 +217,12 @@ def __lt__(self, other):
 
         return Lt(self, other, evaluate=False)
 
+@dispatch(floor, Basic)
+def _eval_Eq(lhs, rhs): # noqa:F811
path: sympy/functions/elementary/integers.py
text: @user1:
I don't think a `floor` could ever compare equal to a non-`Expr` so this should probably by `(floor, Expr)`.

diff_hunk: @@ -1627,5 +1481,103 @@ def union(self, other):
             return AccumBounds(other.min, Max(self.max, other.max))
 
 
+@dispatch(AccumulationBounds, AccumulationBounds) # type: ignore # noqa:F811
+def _eval_is_le(lhs, rhs): # noqa:F811
+    if lhs.max <= rhs.min:
+        return True
+    if lhs.min > rhs.max:
+        return False
+
+@dispatch(AccumulationBounds, Basic) # type: ignore # noqa:F811
+def _eval_is_le(lhs, rhs): # noqa: F811
+
+    """
+    Returns True if range of values attained by `self` AccumulationBounds
+    object is greater than the range of values attained by `other`,
+    where other may be any value of type AccumulationBounds object or
+    extended real number value, False if `other` satisfies
+    the same property, else an unevaluated Relational.
+
+    Examples
+    ========
+
+    >>> from sympy import AccumBounds, oo
+    >>> AccumBounds(1, 3) > AccumBounds(4, oo)
+    False
+    >>> AccumBounds(1, 4) > AccumBounds(3, 4)
+    AccumBounds(1, 4) > AccumBounds(3, 4)
+    >>> AccumBounds(1, oo) > -1
+    True
+
+    """
+    if not rhs.is_extended_real:
+            raise TypeError(
+                "Invalid comparison of %s %s" %
+                (type(rhs), rhs))
+    elif rhs.is_comparable:
+        if lhs.max <= rhs:
+            return True
+        if lhs.min > rhs:
+            return False
path: sympy/calculus/util.py
text: @user1:
What exactly does the use of inequalities here do? Would it make more sense to use `is_le`?

@author:
Actually yeah it would make more sense, to use is_le on these functions

diff_hunk: @@ -2516,3 +2485,46 @@ def set_pow(x, y):
 def set_function(f, x):
     from sympy.sets.handlers.functions import _set_function
     return _set_function(f, x)
+
+
+@dispatch(Interval, FiniteSet)
+def _eval_Eq(lhs, rhs): # noqa: F811
+    return false
path: sympy/sets/sets.py
text: @user1:
Rather than have these things here I think they should go in a module in `sympy/sets/handlers/`

diff_hunk: @@ -199,5 +185,27 @@ def as_immutable(self):
     def is_diagonalizable(self, reals_only=False, **kwargs):
         return super().is_diagonalizable(
             reals_only=reals_only, **kwargs)
+
     is_diagonalizable.__doc__ = SparseMatrix.is_diagonalizable.__doc__
     is_diagonalizable = cacheit(is_diagonalizable)
+
+
+@dispatch(ImmutableDenseMatrix, MatrixExpr)
+def _eval_Eq(lhs, rhs): # noqa:F811
+    return None
path: sympy/matrices/immutable.py
text: @user1:
Is this handler needed? Can it not fall back to a general `Basic, Basic` handler that returns None?