diff --git a/sympy/assumptions/ask.py b/sympy/assumptions/ask.py
index f42c8c37ccea..813ca7c9784a 100644
--- a/sympy/assumptions/ask.py
+++ b/sympy/assumptions/ask.py
@@ -135,7 +135,7 @@ def ask(proposition, assumptions=True, context=global_assumptions):
     # direct resolution method, no logic
     res = key(expr)._eval_ask(assumptions)
     if res is not None:
-        return res
+        return bool(res)
 
     if assumptions == True:
         return
diff --git a/sympy/assumptions/handlers/calculus.py b/sympy/assumptions/handlers/calculus.py
index fa91a4ecd8ef..164444b9ca36 100644
--- a/sympy/assumptions/handlers/calculus.py
+++ b/sympy/assumptions/handlers/calculus.py
@@ -254,11 +254,11 @@ def Pow(expr, assumptions):
             return False
         if base_bounded and exp_bounded:
             return True
-        if (abs(expr.base) <= 1) is True and ask(Q.positive(expr.exp), assumptions):
+        if (abs(expr.base) <= 1) == True and ask(Q.positive(expr.exp), assumptions):
             return True
-        if (abs(expr.base) >= 1) is True and ask(Q.negative(expr.exp), assumptions):
+        if (abs(expr.base) >= 1) == True and ask(Q.negative(expr.exp), assumptions):
             return True
-        if (abs(expr.base) >= 1) is True and exp_bounded is False:
+        if (abs(expr.base) >= 1) == True and exp_bounded is False:
             return False
         return None
 
diff --git a/sympy/concrete/delta.py b/sympy/concrete/delta.py
index 80e33c0be0ff..eb8c37d98566 100644
--- a/sympy/concrete/delta.py
+++ b/sympy/concrete/delta.py
@@ -166,7 +166,7 @@ def deltaproduct(f, limit):
     """
     from sympy.concrete.products import product
 
-    if ((limit[2] - limit[1]) < 0) is True:
+    if ((limit[2] - limit[1]) < 0) == True:
         return S.One
 
     if not f.has(KroneckerDelta):
@@ -284,7 +284,7 @@ def deltasummation(f, limit, no_piecewise=False):
     from sympy.concrete.summations import summation
     from sympy.solvers import solve
 
-    if ((limit[2] - limit[1]) < 0) is True:
+    if ((limit[2] - limit[1]) < 0) == True:
         return S.Zero
 
     if not f.has(KroneckerDelta):
diff --git a/sympy/concrete/summations.py b/sympy/concrete/summations.py
index b12e735247da..570ba23aa1bc 100644
--- a/sympy/concrete/summations.py
+++ b/sympy/concrete/summations.py
@@ -194,7 +194,7 @@ def doit(self, **hints):
         for n, limit in enumerate(self.limits):
             i, a, b = limit
             dif = b - a
-            if dif.is_integer and (dif < 0) is True:
+            if dif.is_integer and (dif < 0) == True:
                 a, b = b + 1, a - 1
                 f = -f
 
@@ -308,7 +308,7 @@ def euler_maclaurin(self, m=0, n=0, eps=0, eval_integral=True):
         if len(self.limits) != 1:
             raise ValueError("More than 1 limit")
         i, a, b = self.limits[0]
-        if (a > b) is True:
+        if (a > b) == True:
             if a - b == 1:
                 return S.Zero,S.Zero
             a, b = b + 1, a - 1
@@ -324,10 +324,9 @@ def euler_maclaurin(self, m=0, n=0, eps=0, eval_integral=True):
                 term = f.subs(i, a)
                 if term:
                     test = abs(term.evalf(3)) < eps
-                    if isinstance(test, bool):
-                        if test is True:
-                            return s, abs(term)
-                    else:
+                    if test == True:
+                        return s, abs(term)
+                    elif not (test == False):
                         # a symbolic Relational class, can't go further
                         return term, S.Zero
                 s += term
@@ -758,7 +757,7 @@ def eval_sum_hyper(f, i_a_b):
                 return None
             (res1, cond1), (res2, cond2) = res1, res2
             cond = And(cond1, cond2)
-            if cond is False:
+            if cond == False:
                 return None
         return Piecewise((res1 - res2, cond), (old_sum, True))
 
@@ -770,7 +769,7 @@ def eval_sum_hyper(f, i_a_b):
         res1, cond1 = res1
         res2, cond2 = res2
         cond = And(cond1, cond2)
-        if cond is False:
+        if cond == False:
             return None
         return Piecewise((res1 + res2, cond), (old_sum, True))
 
diff --git a/sympy/core/containers.py b/sympy/core/containers.py
index 23d9339819bb..0f16ec19241f 100644
--- a/sympy/core/containers.py
+++ b/sympy/core/containers.py
@@ -94,10 +94,10 @@ def _to_mpmath(self, prec):
         return tuple([a._to_mpmath(prec) for a in self.args])
 
     def __lt__(self, other):
-        return self.args < other.args
+        return sympify(self.args < other.args)
 
     def __le__(self, other):
-        return self.args <= other.args
+        return sympify(self.args <= other.args)
 
     # XXX: Basic defines count() as something different, so we can't
     # redefine it here. Originally this lead to cse() test failure.
@@ -247,7 +247,7 @@ def __contains__(self, key):
         return sympify(key) in self._dict
 
     def __lt__(self, other):
-        return self.args < other.args
+        return sympify(self.args < other.args)
 
     @property
     def _sorted_args(self):
diff --git a/sympy/core/expr.py b/sympy/core/expr.py
index 3c2011efe066..96e4b1ebf29d 100644
--- a/sympy/core/expr.py
+++ b/sympy/core/expr.py
@@ -218,7 +218,7 @@ def __ge__(self, other):
             raise TypeError("Invalid comparison of complex %s" % dif)
         if dif.is_nonnegative is not None and \
                 dif.is_nonnegative is not dif.is_negative:
-            return dif.is_nonnegative
+            return sympify(dif.is_nonnegative)
         return C.GreaterThan(self, other)
 
     @_sympifyit('other', False)  # sympy >  other
@@ -228,7 +228,7 @@ def __le__(self, other):
             raise TypeError("Invalid comparison of complex %s" % dif)
         if dif.is_nonpositive is not None and \
                 dif.is_nonpositive is not dif.is_positive:
-            return dif.is_nonpositive
+            return sympify(dif.is_nonpositive)
         return C.LessThan(self, other)
 
     @_sympifyit('other', False)  # sympy >  other
@@ -238,7 +238,7 @@ def __gt__(self, other):
             raise TypeError("Invalid comparison of complex %s" % dif)
         if dif.is_positive is not None and \
                 dif.is_positive is not dif.is_nonpositive:
-            return dif.is_positive
+            return sympify(dif.is_positive)
         return C.StrictGreaterThan(self, other)
 
     @_sympifyit('other', False)  # sympy >  other
@@ -248,7 +248,7 @@ def __lt__(self, other):
             raise TypeError("Invalid comparison of complex %s" % dif)
         if dif.is_negative is not None and \
                 dif.is_negative is not dif.is_nonnegative:
-            return dif.is_negative
+            return sympify(dif.is_negative)
         return C.StrictLessThan(self, other)
 
     @staticmethod
@@ -2019,7 +2019,7 @@ def could_extract_minus_sign(self):
                 return len(negative_args) % 2 == 1
 
             # As a last resort, we choose the one with greater value of .sort_key()
-            return self.sort_key() < negative_self.sort_key()
+            return bool(self.sort_key() < negative_self.sort_key())
 
     def extract_branch_factor(self, allow_half=False):
         """
diff --git a/sympy/core/exprtools.py b/sympy/core/exprtools.py
index a83c12fcadc2..92f5683cfe9b 100644
--- a/sympy/core/exprtools.py
+++ b/sympy/core/exprtools.py
@@ -534,9 +534,9 @@ def gcd(self, other):  # Factors
         for factor, exp in self.factors.items():
             if factor in other.factors:
                 lt = (exp < other.factors[factor])
-                if lt is True:
+                if lt == True:
                     factors[factor] = exp
-                elif lt is False:
+                elif lt == False:
                     factors[factor] = other.factors[factor]
 
         return Factors(factors)
diff --git a/sympy/core/mod.py b/sympy/core/mod.py
index 5b51e4a24d20..738039f8b18c 100644
--- a/sympy/core/mod.py
+++ b/sympy/core/mod.py
@@ -52,7 +52,7 @@ def doit(p, q):
             else:
                 if type(d) is int:
                     rv = p - d*q
-                    if (rv*q < 0) is True:
+                    if (rv*q < 0) == True:
                         rv += q
                     return rv
 
diff --git a/sympy/core/numbers.py b/sympy/core/numbers.py
index 458b78826c29..0cc780da272f 100644
--- a/sympy/core/numbers.py
+++ b/sympy/core/numbers.py
@@ -874,52 +874,56 @@ def __gt__(self, other):
         try:
             other = _sympify(other)
         except SympifyError:
-            return False    # sympy > other
+            return S.false    # sympy > other
         if isinstance(other, NumberSymbol):
             return other.__le__(self)
         if other.is_comparable:
             other = other.evalf()
         if isinstance(other, Number):
-            return bool(mlib.mpf_gt(self._mpf_, other._as_mpf_val(self._prec)))
+            return _sympify(bool(
+                mlib.mpf_gt(self._mpf_, other._as_mpf_val(self._prec))))
         return Expr.__gt__(self, other)
 
     def __ge__(self, other):
         try:
             other = _sympify(other)
         except SympifyError:
-            return False    # sympy > other  -->  ! <=
+            return S.false    # sympy > other  -->  ! <=
         if isinstance(other, NumberSymbol):
             return other.__lt__(self)
         if other.is_comparable:
             other = other.evalf()
         if isinstance(other, Number):
-            return bool(mlib.mpf_ge(self._mpf_, other._as_mpf_val(self._prec)))
+            return _sympify(bool(
+                mlib.mpf_ge(self._mpf_, other._as_mpf_val(self._prec))))
         return Expr.__ge__(self, other)
 
     def __lt__(self, other):
         try:
             other = _sympify(other)
         except SympifyError:
-            return False    # sympy > other
+            return S.false    # sympy > other
         if isinstance(other, NumberSymbol):
             return other.__ge__(self)
         if other.is_real and other.is_number:
             other = other.evalf()
         if isinstance(other, Number):
-            return bool(mlib.mpf_lt(self._mpf_, other._as_mpf_val(self._prec)))
+            return _sympify(bool(
+                mlib.mpf_lt(self._mpf_, other._as_mpf_val(self._prec))))
         return Expr.__lt__(self, other)
 
     def __le__(self, other):
         try:
             other = _sympify(other)
         except SympifyError:
-            return False    # sympy > other  -->  ! <=
+            return S.false    # sympy > other  -->  ! <=
         if isinstance(other, NumberSymbol):
             return other.__gt__(self)
         if other.is_real and other.is_number:
             other = other.evalf()
         if isinstance(other, Number):
-            return bool(mlib.mpf_le(self._mpf_, other._as_mpf_val(self._prec)))
+            return _sympify(bool(
+                mlib.mpf_le(self._mpf_, other._as_mpf_val(self._prec))))
         return Expr.__le__(self, other)
 
     def __hash__(self):
@@ -1320,17 +1324,17 @@ def __gt__(self, other):
         try:
             other = _sympify(other)
         except SympifyError:
-            return False    # sympy > other  --> not <
+            return S.false    # sympy > other  --> not <
         if isinstance(other, NumberSymbol):
             return other.__le__(self)
         if other.is_real and other.is_number and not isinstance(other, Rational):
             other = other.evalf()
         if isinstance(other, Number):
             if isinstance(other, Rational):
-                return bool(self.p*other.q > self.q*other.p)
+                return _sympify(bool(self.p*other.q > self.q*other.p))
             if isinstance(other, Float):
-                return bool(mlib.mpf_gt(
-                    self._as_mpf_val(other._prec), other._mpf_))
+                return _sympify(bool(mlib.mpf_gt(
+                    self._as_mpf_val(other._prec), other._mpf_)))
             if other is S.NaN:
                 return other.__le__(self)
         return Expr.__gt__(self, other)
@@ -1339,17 +1343,17 @@ def __ge__(self, other):
         try:
             other = _sympify(other)
         except SympifyError:
-            return False    # sympy > other  -->  not <=
+            return S.false    # sympy > other  -->  not <=
         if isinstance(other, NumberSymbol):
             return other.__lt__(self)
         if other.is_real and other.is_number and not isinstance(other, Rational):
             other = other.evalf()
         if isinstance(other, Number):
             if isinstance(other, Rational):
-                 return bool(self.p*other.q >= self.q*other.p)
+                 return _sympify(bool(self.p*other.q >= self.q*other.p))
             if isinstance(other, Float):
-                return bool(mlib.mpf_ge(
-                    self._as_mpf_val(other._prec), other._mpf_))
+                return _sympify(bool(mlib.mpf_ge(
+                    self._as_mpf_val(other._prec), other._mpf_)))
             if other is S.NaN:
                 return other.__lt__(self)
         return Expr.__ge__(self, other)
@@ -1358,17 +1362,17 @@ def __lt__(self, other):
         try:
             other = _sympify(other)
         except SympifyError:
-            return False    # sympy > other  --> not <
+            return S.false    # sympy > other  --> not <
         if isinstance(other, NumberSymbol):
             return other.__ge__(self)
         if other.is_real and other.is_number and not isinstance(other, Rational):
             other = other.evalf()
         if isinstance(other, Number):
             if isinstance(other, Rational):
-                return bool(self.p*other.q < self.q*other.p)
+                return _sympify(bool(self.p*other.q < self.q*other.p))
             if isinstance(other, Float):
-                return bool(mlib.mpf_lt(
-                    self._as_mpf_val(other._prec), other._mpf_))
+                return _sympify(bool(mlib.mpf_lt(
+                    self._as_mpf_val(other._prec), other._mpf_)))
             if other is S.NaN:
                 return other.__ge__(self)
         return Expr.__lt__(self, other)
@@ -1377,17 +1381,17 @@ def __le__(self, other):
         try:
             other = _sympify(other)
         except SympifyError:
-            return False    # sympy > other  -->  not <=
+            return S.false    # sympy > other  -->  not <=
         if isinstance(other, NumberSymbol):
             return other.__gt__(self)
         if other.is_real and other.is_number and not isinstance(other, Rational):
             other = other.evalf()
         if isinstance(other, Number):
             if isinstance(other, Rational):
-                return bool(self.p*other.q <= self.q*other.p)
+                return _sympify(bool(self.p*other.q <= self.q*other.p))
             if isinstance(other, Float):
-                return bool(mlib.mpf_le(
-                    self._as_mpf_val(other._prec), other._mpf_))
+                return _sympify(bool(mlib.mpf_le(
+                    self._as_mpf_val(other._prec), other._mpf_)))
             if other is S.NaN:
                 return other.__gt__(self)
         return Expr.__le__(self, other)
@@ -1679,30 +1683,30 @@ def __ne__(self, other):
 
     def __gt__(self, other):
         if isinstance(other, integer_types):
-            return (self.p > other)
+            return _sympify(self.p > other)
         elif isinstance(other, Integer):
-            return (self.p > other.p)
+            return _sympify(self.p > other.p)
         return Rational.__gt__(self, other)
 
     def __lt__(self, other):
         if isinstance(other, integer_types):
-            return (self.p < other)
+            return _sympify(self.p < other)
         elif isinstance(other, Integer):
-            return (self.p < other.p)
+            return _sympify(self.p < other.p)
         return Rational.__lt__(self, other)
 
     def __ge__(self, other):
         if isinstance(other, integer_types):
-            return (self.p >= other)
+            return _sympify(self.p >= other)
         elif isinstance(other, Integer):
-            return (self.p >= other.p)
+            return _sympify(self.p >= other.p)
         return Rational.__ge__(self, other)
 
     def __le__(self, other):
         if isinstance(other, integer_types):
-            return (self.p <= other)
+            return _sympify(self.p <= other)
         elif isinstance(other, Integer):
-            return (self.p <= other.p)
+            return _sympify(self.p <= other.p)
         return Rational.__le__(self, other)
 
     def __hash__(self):
@@ -2240,25 +2244,25 @@ def __ne__(self, other):
     def __lt__(self, other):
         if other.is_number and other.is_real is False:
             raise TypeError("Invalid comparison of %s and %s" % (self, other))
-        return False
+        return S.false
 
     @_sympifyit('other', NotImplemented)
     def __le__(self, other):
         if other.is_number and other.is_real is False:
             raise TypeError("Invalid comparison of %s and %s" % (self, other))
-        return other is S.Infinity
+        return _sympify(other is S.Infinity)
 
     @_sympifyit('other', NotImplemented)
     def __gt__(self, other):
         if other.is_number and other.is_real is False:
             raise TypeError("Invalid comparison of %s and %s" % (self, other))
-        return other is not S.Infinity
+        return _sympify(other is not S.Infinity)
 
     @_sympifyit('other', NotImplemented)
     def __ge__(self, other):
         if other.is_number and other.is_real is False:
             raise TypeError("Invalid comparison of %s and %s" % (self, other))
-        return True
+        return S.true
 
     def __mod__(self, other):
         return S.NaN
@@ -2433,25 +2437,25 @@ def __ne__(self, other):
     def __lt__(self, other):
         if other.is_number and other.is_real is False:
             raise TypeError("Invalid comparison of %s and %s" % (self, other))
-        return other is not S.NegativeInfinity
+        return _sympify(other is not S.NegativeInfinity)
 
     @_sympifyit('other', NotImplemented)
     def __le__(self, other):
         if other.is_number and other.is_real is False:
             raise TypeError("Invalid comparison of %s and %s" % (self, other))
-        return True
+        return S.true
 
     @_sympifyit('other', NotImplemented)
     def __gt__(self, other):
         if other.is_number and other.is_real is False:
             raise TypeError("Invalid comparison of %s and %s" % (self, other))
-        return False
+        return S.false
 
     @_sympifyit('other', NotImplemented)
     def __ge__(self, other):
         if other.is_number and other.is_real is False:
             raise TypeError("Invalid comparison of %s and %s" % (self, other))
-        return other is S.NegativeInfinity
+        return _sympify(other is S.NegativeInfinity)
 
 
 class NaN(with_metaclass(Singleton, Number)):
@@ -2543,16 +2547,16 @@ def __ne__(self, other):
         return other is not S.NaN
 
     def __gt__(self, other):
-        return False
+        return S.false
 
     def __ge__(self, other):
-        return False
+        return S.false
 
     def __lt__(self, other):
-        return False
+        return S.false
 
     def __le__(self, other):
-        return False
+        return S.false
 
 nan = S.NaN
 
@@ -2665,41 +2669,41 @@ def __lt__(self, other):
         try:
             other = _sympify(other)
         except SympifyError:
-            return False    # sympy > other  --> not <
+            return S.false    # sympy > other  --> not <
         if self is other:
-            return False
+            return S.false
         if isinstance(other, Number):
             approx = self.approximation_interval(other.__class__)
             if approx is not None:
                 l, u = approx
                 if other < l:
-                    return False
+                    return S.false
                 if other > u:
-                    return True
-            return self.evalf() < other
+                    return S.true
+            return _sympify(self.evalf() < other)
         if other.is_real and other.is_number:
             other = other.evalf()
-            return self.evalf() < other
+            return _sympify(self.evalf() < other)
         return Expr.__lt__(self, other)
 
     def __le__(self, other):
         try:
             other = _sympify(other)
         except SympifyError:
-            return False    # sympy > other  --> not <=
+            return S.false    # sympy > other  --> not <=
         if self is other:
-            return True
+            return S.true
         if other.is_real and other.is_number:
             other = other.evalf()
         if isinstance(other, Number):
-            return self.evalf() <= other
+            return _sympify(self.evalf() <= other)
         return Expr.__le__(self, other)
 
     def __gt__(self, other):
-        return (-self) < (-other)
+        return _sympify((-self) < (-other))
 
     def __ge__(self, other):
-        return (-self) <= (-other)
+        return _sympify((-self) <= (-other))
 
     def __int__(self):
         # subclass with appropriate return value
diff --git a/sympy/core/power.py b/sympy/core/power.py
index 9560b6f3bfa3..5fb1d78b4f1c 100644
--- a/sympy/core/power.py
+++ b/sympy/core/power.py
@@ -222,7 +222,7 @@ def _eval_power(self, other):
                 e.is_real is False and smallarg is False):
             return -self.func(b, e*other)
         if (other.is_integer or
-            e.is_real and (b_nneg or (abs(e) < 1) is True) or
+            e.is_real and (b_nneg or (abs(e) < 1) == True) or
             e.is_real is False and smallarg is True or
                 b.is_polar):
             return self.func(b, e*other)
@@ -795,9 +795,8 @@ def _eval_is_polynomial(self, syms):
             return False
 
         if self.base.has(*syms):
-            return self.base._eval_is_polynomial(syms) and \
-                self.exp.is_Integer and \
-                (self.exp >= 0) is True
+            return bool(self.base._eval_is_polynomial(syms) and
+                self.exp.is_Integer and (self.exp >= 0))
         else:
             return True
 
diff --git a/sympy/core/sets.py b/sympy/core/sets.py
index 3edc1a67d4c6..065b51ecd300 100644
--- a/sympy/core/sets.py
+++ b/sympy/core/sets.py
@@ -816,7 +816,7 @@ def as_relational(self, symbol):
             right = other < self.end
         else:
             right = other <= self.end
-        if right is True:
+        if right == True:
             if self.left_open:
                 return other > self.start
             else:
diff --git a/sympy/core/tests/test_containers.py b/sympy/core/tests/test_containers.py
index 03d9bdcdd369..2c35f86f791a 100644
--- a/sympy/core/tests/test_containers.py
+++ b/sympy/core/tests/test_containers.py
@@ -63,10 +63,10 @@ def test_Tuple_equality():
 
 
 def test_Tuple_comparision():
-    assert (Tuple(1, 3) >= Tuple(-10, 30)) is True
-    assert (Tuple(1, 3) <= Tuple(-10, 30)) is False
-    assert (Tuple(1, 3) >= Tuple(1, 3)) is True
-    assert (Tuple(1, 3) <= Tuple(1, 3)) is True
+    assert (Tuple(1, 3) >= Tuple(-10, 30)) is S.true
+    assert (Tuple(1, 3) <= Tuple(-10, 30)) is S.false
+    assert (Tuple(1, 3) >= Tuple(1, 3)) is S.true
+    assert (Tuple(1, 3) <= Tuple(1, 3)) is S.true
 
 
 def test_Tuple_tuple_count():
diff --git a/sympy/core/tests/test_expr.py b/sympy/core/tests/test_expr.py
index 84c177bf6d00..1a82de3085b6 100644
--- a/sympy/core/tests/test_expr.py
+++ b/sympy/core/tests/test_expr.py
@@ -161,15 +161,15 @@ def s(a, b):
 
 
 def test_relational():
-    assert (pi < 3) is False
-    assert (pi <= 3) is False
-    assert (pi > 3) is True
-    assert (pi >= 3) is True
-    assert (-pi < 3) is True
-    assert (-pi <= 3) is True
-    assert (-pi > 3) is False
-    assert (-pi >= 3) is False
-    assert (x - 2 < x - 3) is False
+    assert (pi < 3) is S.false
+    assert (pi <= 3) is S.false
+    assert (pi > 3) is S.true
+    assert (pi >= 3) is S.true
+    assert (-pi < 3) is S.true
+    assert (-pi <= 3) is S.true
+    assert (-pi > 3) is S.false
+    assert (-pi >= 3) is S.false
+    assert (x - 2 < x - 3) is S.false
 
 
 def test_relational_assumptions():
@@ -186,26 +186,26 @@ def test_relational_assumptions():
     m2 = Symbol("m2", positive=False, real=True)
     m3 = Symbol("m3", nonpositive=False, real=True)
     m4 = Symbol("m4", negative=False, real=True)
-    assert (m1 < 0) is True
-    assert (m2 <= 0) is True
-    assert (m3 > 0) is True
-    assert (m4 >= 0) is True
+    assert (m1 < 0) is S.true
+    assert (m2 <= 0) is S.true
+    assert (m3 > 0) is S.true
+    assert (m4 >= 0) is S.true
     m1 = Symbol("m1", negative=True)
     m2 = Symbol("m2", nonpositive=True)
     m3 = Symbol("m3", positive=True)
     m4 = Symbol("m4", nonnegative=True)
-    assert (m1 < 0) is True
-    assert (m2 <= 0) is True
-    assert (m3 > 0) is True
-    assert (m4 >= 0) is True
+    assert (m1 < 0) is S.true
+    assert (m2 <= 0) is S.true
+    assert (m3 > 0) is S.true
+    assert (m4 >= 0) is S.true
     m1 = Symbol("m1", negative=False)
     m2 = Symbol("m2", nonpositive=False)
     m3 = Symbol("m3", positive=False)
     m4 = Symbol("m4", nonnegative=False)
-    assert (m1 < 0) is False
-    assert (m2 <= 0) is False
-    assert (m3 > 0) is False
-    assert (m4 >= 0) is False
+    assert (m1 < 0) is S.false
+    assert (m2 <= 0) is S.false
+    assert (m3 > 0) is S.false
+    assert (m4 >= 0) is S.false
 
 
 def test_relational_noncommutative():
diff --git a/sympy/core/tests/test_numbers.py b/sympy/core/tests/test_numbers.py
index e808e5732c04..d8a8bfb87872 100644
--- a/sympy/core/tests/test_numbers.py
+++ b/sympy/core/tests/test_numbers.py
@@ -335,10 +335,10 @@ def test_Rational_cmp():
     assert not (Rational(-1) > 0)
     assert Rational(-1) < 0
 
-    assert (n1 < S.NaN) is False
-    assert (n1 <= S.NaN) is False
-    assert (n1 > S.NaN) is False
-    assert (n1 <= S.NaN) is False
+    assert (n1 < S.NaN) is S.false
+    assert (n1 <= S.NaN) is S.false
+    assert (n1 > S.NaN) is S.false
+    assert (n1 <= S.NaN) is S.false
 
 
 def test_Float():
diff --git a/sympy/core/tests/test_relational.py b/sympy/core/tests/test_relational.py
index ec8f53046875..bc92d56ac0f6 100644
--- a/sympy/core/tests/test_relational.py
+++ b/sympy/core/tests/test_relational.py
@@ -85,32 +85,31 @@ def test_Eq():
 
 def test_rel_Infinity():
     # NOTE: All of these are actually handled by sympy.core.Number, and do
-    # not create Relational objects.  Therefore, they still return True and
-    # False instead of S.true and S.false.
-    assert (oo > oo) is False
-    assert (oo > -oo) is True
-    assert (oo > 1) is True
-    assert (oo < oo) is False
-    assert (oo < -oo) is False
-    assert (oo < 1) is False
-    assert (oo >= oo) is True
-    assert (oo >= -oo) is True
-    assert (oo >= 1) is True
-    assert (oo <= oo) is True
-    assert (oo <= -oo) is False
-    assert (oo <= 1) is False
-    assert (-oo > oo) is False
-    assert (-oo > -oo) is False
-    assert (-oo > 1) is False
-    assert (-oo < oo) is True
-    assert (-oo < -oo) is False
-    assert (-oo < 1) is True
-    assert (-oo >= oo) is False
-    assert (-oo >= -oo) is True
-    assert (-oo >= 1) is False
-    assert (-oo <= oo) is True
-    assert (-oo <= -oo) is True
-    assert (-oo <= 1) is True
+    # not create Relational objects.
+    assert (oo > oo) is S.false
+    assert (oo > -oo) is S.true
+    assert (oo > 1) is S.true
+    assert (oo < oo) is S.false
+    assert (oo < -oo) is S.false
+    assert (oo < 1) is S.false
+    assert (oo >= oo) is S.true
+    assert (oo >= -oo) is S.true
+    assert (oo >= 1) is S.true
+    assert (oo <= oo) is S.true
+    assert (oo <= -oo) is S.false
+    assert (oo <= 1) is S.false
+    assert (-oo > oo) is S.false
+    assert (-oo > -oo) is S.false
+    assert (-oo > 1) is S.false
+    assert (-oo < oo) is S.true
+    assert (-oo < -oo) is S.false
+    assert (-oo < 1) is S.true
+    assert (-oo >= oo) is S.false
+    assert (-oo >= -oo) is S.true
+    assert (-oo >= 1) is S.false
+    assert (-oo <= oo) is S.true
+    assert (-oo <= -oo) is S.true
+    assert (-oo <= 1) is S.true
 
 
 def test_bool():
diff --git a/sympy/core/tests/test_sets.py b/sympy/core/tests/test_sets.py
index 7e3db205dba2..2d0209b40be1 100644
--- a/sympy/core/tests/test_sets.py
+++ b/sympy/core/tests/test_sets.py
@@ -373,7 +373,7 @@ def test_Interval_as_relational():
     assert Interval(-2, oo, left_open=False).as_relational(x) == Ge(x, -2)
     assert Interval(-2, oo, left_open=True).as_relational(x) == Gt(x, -2)
 
-    assert Interval(-oo, oo).as_relational(x) is True
+    assert Interval(-oo, oo).as_relational(x) is S.true
 
 
 def test_Finite_as_relational():
diff --git a/sympy/functions/elementary/complexes.py b/sympy/functions/elementary/complexes.py
index 5fcb29ce158f..e259d72a2dca 100644
--- a/sympy/functions/elementary/complexes.py
+++ b/sympy/functions/elementary/complexes.py
@@ -832,7 +832,7 @@ def mr(expr):
             if arg == 0:
                 return abs(c)*principal_branch(Mul(*m), period)
             return principal_branch(exp_polar(I*arg)*Mul(*m), period)*abs(c)
-        if arg.is_number and ((abs(arg) < period/2) is True or arg == period/2) \
+        if arg.is_number and ((abs(arg) < period/2) == True or arg == period/2) \
                 and m == ():
             return exp_polar(arg*I)*abs(c)
 
diff --git a/sympy/functions/elementary/exponential.py b/sympy/functions/elementary/exponential.py
index b1318c24fdf7..725853c8a721 100644
--- a/sympy/functions/elementary/exponential.py
+++ b/sympy/functions/elementary/exponential.py
@@ -113,9 +113,9 @@ def _eval_power(b, e):
         if be.is_polar:
             return rv
         besmall = abs(be) <= S.Pi
-        if besmall is True:
+        if besmall == True:
             return rv
-        elif besmall is False and e.is_Rational and e.q == 2:
+        elif besmall == False and e.is_Rational and e.q == 2:
             return -rv
 
     def _eval_expand_power_exp(self, **hints):
diff --git a/sympy/functions/elementary/miscellaneous.py b/sympy/functions/elementary/miscellaneous.py
index f1f1738202d1..33a492c40f21 100644
--- a/sympy/functions/elementary/miscellaneous.py
+++ b/sympy/functions/elementary/miscellaneous.py
@@ -374,7 +374,9 @@ def _is_connected(cls, x, y):
         """
         Check if x and y are connected somehow.
         """
-        if (x == y) or isinstance(x > y, bool) or isinstance(x < y, bool):
+        xy = x > y
+        yx = x < y
+        if (x == y) or xy == True or xy == False or yx == True or yx == False:
             return True
         if x.is_Number and y.is_Number:
             return True
@@ -395,15 +397,11 @@ def _is_asneeded(cls, x, y):
             if cls._rel(x, y):
                 return True
         xy = cls._rel(x, y)
-        if isinstance(xy, bool):
-            if xy:
-                return True
-            return False
+        if xy == True or xy == False:
+            return bool(xy)
         yx = cls._rel_inversed(x, y)
-        if isinstance(yx, bool):
-            if yx:
-                return False  # never occurs?
-            return True
+        if yx == True or yx == False:
+            return not bool(yx)
         return False
 
     def _eval_derivative(self, s):
diff --git a/sympy/functions/elementary/piecewise.py b/sympy/functions/elementary/piecewise.py
index 1252ad9cadec..ddf21d91bacd 100644
--- a/sympy/functions/elementary/piecewise.py
+++ b/sympy/functions/elementary/piecewise.py
@@ -209,41 +209,41 @@ def _eval_interval(self, sym, a, b):
                 super(Piecewise, self)._eval_interval(sym, a, b))
 
         mul = 1
-        if (a == b) is True:
+        if (a == b) == True:
             return S.Zero
-        elif (a > b) is True:
+        elif (a > b) == True:
             a, b, mul = b, a, -1
-        elif (a <= b) is not True:
+        elif (a <= b) != True:
             newargs = []
             for e, c in self.args:
                 intervals = self._sort_expr_cond(
                     sym, S.NegativeInfinity, S.Infinity, c)
                 values = []
                 for lower, upper, expr in intervals:
-                    if (a < lower) is True:
+                    if (a < lower) == True:
                         mid = lower
                         rep = b
                         val = e._eval_interval(sym, mid, b)
                         val += self._eval_interval(sym, a, mid)
-                    elif (a > upper) is True:
+                    elif (a > upper) == True:
                         mid = upper
                         rep = b
                         val = e._eval_interval(sym, mid, b)
                         val += self._eval_interval(sym, a, mid)
-                    elif (a >= lower) is True and (a <= upper) is True:
+                    elif (a >= lower) == True and (a <= upper) == True:
                         rep = b
                         val = e._eval_interval(sym, a, b)
-                    elif (b < lower) is True:
+                    elif (b < lower) == True:
                         mid = lower
                         rep = a
                         val = e._eval_interval(sym, a, mid)
                         val += self._eval_interval(sym, mid, b)
-                    elif (b > upper) is True:
+                    elif (b > upper) == True:
                         mid = upper
                         rep = a
                         val = e._eval_interval(sym, a, mid)
                         val += self._eval_interval(sym, mid, b)
-                    elif ((b >= lower) is True) and ((b <= upper) is True):
+                    elif ((b >= lower) == True) and ((b <= upper) == True):
                         rep = a
                         val = e._eval_interval(sym, a, b)
                     else:
@@ -365,8 +365,8 @@ def _sort_expr_cond(self, sym, a, b, targetcond=None):
                         self.__eval_cond(upper == int_expr[n][1]):
                     upper = Min(upper, int_expr[n][0])
                 elif len(int_expr[n][1].free_symbols) and \
-                        (lower >= int_expr[n][0]) is not True and \
-                        (int_expr[n][1] == Min(lower, upper)) is not True:
+                        (lower >= int_expr[n][0]) != True and \
+                        (int_expr[n][1] == Min(lower, upper)) != True:
                     upper = Min(upper, int_expr[n][0])
                 elif self.__eval_cond(upper > int_expr[n][0]) and \
                         self.__eval_cond(upper <= int_expr[n][1]):
@@ -375,7 +375,7 @@ def _sort_expr_cond(self, sym, a, b, targetcond=None):
                         self.__eval_cond(upper < int_expr[n][1]):
                     int_expr[n][0] = Max(upper, int_expr[n][0])
 
-            if self.__eval_cond(lower >= upper) is not True:  # Is it still an interval?
+            if self.__eval_cond(lower >= upper) != True:  # Is it still an interval?
                 int_expr.append([lower, upper, expr])
             if orig_cond == targetcond:
                 return [(lower, upper, None)]
@@ -409,14 +409,14 @@ def _sort_expr_cond(self, sym, a, b, targetcond=None):
         holes = []
         curr_low = a
         for int_a, int_b, expr in int_expr:
-            if (curr_low < int_a) is True:
+            if (curr_low < int_a) == True:
                 holes.append([curr_low, Min(b, int_a), default])
-            elif (curr_low >= int_a) is not True:
+            elif (curr_low >= int_a) != True:
                 holes.append([curr_low, Min(b, int_a), default])
             curr_low = Min(b, int_b)
-        if (curr_low < b) is True:
+        if (curr_low < b) == True:
             holes.append([Min(b, curr_low), b, default])
-        elif (curr_low >= b) is not True:
+        elif (curr_low >= b) != True:
             holes.append([Min(b, curr_low), b, default])
 
         if holes and default is not None:
diff --git a/sympy/functions/elementary/trigonometric.py b/sympy/functions/elementary/trigonometric.py
index 88f3e392fbce..dd34fc7e5f0d 100644
--- a/sympy/functions/elementary/trigonometric.py
+++ b/sympy/functions/elementary/trigonometric.py
@@ -1433,7 +1433,11 @@ def _eval_rewrite_as_log(self, x):
         return -S.ImaginaryUnit*C.log(S.ImaginaryUnit*x + sqrt(1 - x**2))
 
     def _eval_is_real(self):
-        return self.args[0].is_real and (self.args[0] >= -1 and self.args[0] <= 1)
+        r = self.args[0].is_real and (self.args[0] >= -1 and
+                                      self.args[0] <= 1)
+        if r == True or r == False:
+            r = bool(r)
+        return r
 
     def inverse(self, argindex=1):
         """
@@ -1545,7 +1549,11 @@ def _eval_as_leading_term(self, x):
             return self.func(arg)
 
     def _eval_is_real(self):
-        return self.args[0].is_real and (self.args[0] >= -1 and self.args[0] <= 1)
+        r = self.args[0].is_real and (self.args[0] >= -1 and
+                                      self.args[0] <= 1)
+        if r == True or r == False:
+            r = bool(r)
+        return r
 
     def _eval_rewrite_as_log(self, x):
         return S.Pi/2 + S.ImaginaryUnit * C.log(S.ImaginaryUnit * x + sqrt(1 - x**2))
diff --git a/sympy/functions/special/error_functions.py b/sympy/functions/special/error_functions.py
index aeeb2015dd66..5eef476e161e 100644
--- a/sympy/functions/special/error_functions.py
+++ b/sympy/functions/special/error_functions.py
@@ -1210,7 +1210,7 @@ def eval(cls, nu, z):
         if n == 0:
             return
         if nu.is_integer:
-            if (nu > 0) is not True:
+            if (nu > 0) != True:
                 return
             return expint(nu, z) \
                 - 2*pi*I*n*(-1)**(nu - 1)/factorial(nu - 1)*unpolarify(z)**(nu - 1)
diff --git a/sympy/functions/special/gamma_functions.py b/sympy/functions/special/gamma_functions.py
index 3933c1169d04..e666211261c6 100644
--- a/sympy/functions/special/gamma_functions.py
+++ b/sympy/functions/special/gamma_functions.py
@@ -421,11 +421,11 @@ def eval(cls, a, z):
 
         # We extract branching information here. C/f lowergamma.
         nx, n = z.extract_branch_factor()
-        if a.is_integer and (a > 0) is True:
+        if a.is_integer and (a > 0) == True:
             nx = unpolarify(z)
             if z != nx:
                 return uppergamma(a, nx)
-        elif a.is_integer and (a <= 0) is True:
+        elif a.is_integer and (a <= 0) == True:
             if n != 0:
                 return -2*pi*I*n*(-1)**(-a)/factorial(-a) + uppergamma(a, nx)
         elif n != 0:
@@ -569,11 +569,11 @@ def fdiff(self, argindex=2):
             raise ArgumentIndexError(self, argindex)
 
     def _eval_is_positive(self):
-        if self.args[1].is_positive and (self.args[0] > 0) is True:
+        if self.args[1].is_positive and (self.args[0] > 0) == True:
             return self.args[0].is_odd
 
     def _eval_is_negative(self):
-        if self.args[1].is_positive and (self.args[0] > 0) is True:
+        if self.args[1].is_positive and (self.args[0] > 0) == True:
             return self.args[0].is_even
 
     def _eval_is_real(self):
diff --git a/sympy/functions/special/hyper.py b/sympy/functions/special/hyper.py
index e1dc3c29a989..aa58e4a69844 100644
--- a/sympy/functions/special/hyper.py
+++ b/sympy/functions/special/hyper.py
@@ -253,9 +253,9 @@ def radius_of_convergence(self):
         >>> hyper((1, 2), (3, 4), z).radius_of_convergence
         oo
         """
-        if any(a.is_integer and (a <= 0) is True for a in self.ap + self.bq):
-            aints = [a for a in self.ap if a.is_Integer and (a <= 0) is True]
-            bints = [a for a in self.bq if a.is_Integer and (a <= 0) is True]
+        if any(a.is_integer and (a <= 0) == True for a in self.ap + self.bq):
+            aints = [a for a in self.ap if a.is_Integer and (a <= 0) == True]
+            bints = [a for a in self.bq if a.is_Integer and (a <= 0) == True]
             if len(aints) < len(bints):
                 return S(0)
             popped = False
diff --git a/sympy/functions/special/polynomials.py b/sympy/functions/special/polynomials.py
index 6a2a2d8bcdae..e6b0ae6557be 100644
--- a/sympy/functions/special/polynomials.py
+++ b/sympy/functions/special/polynomials.py
@@ -349,7 +349,7 @@ def eval(cls, n, a, x):
         if not n.is_Number:
             # Handle this before the general sign extraction rule
             if x == S.NegativeOne:
-                if (C.re(a) > S.Half) is True:
+                if (C.re(a) > S.Half) == True:
                     return S.ComplexInfinity
                 else:
                     # No sec function available yet
@@ -666,7 +666,7 @@ class chebyshevt_root(Function):
 
     @classmethod
     def eval(cls, n, k):
-        if not ((0 <= k) is (k < n) is True):
+        if not ((0 <= k) and (k < n)):
             raise ValueError("must have 0 <= k < n, "
                 "got k = %s and n = %s" % (k, n))
         return C.cos(S.Pi*(2*k + 1)/(2*n))
@@ -706,7 +706,7 @@ class chebyshevu_root(Function):
 
     @classmethod
     def eval(cls, n, k):
-        if not ((0 <= k) is (k < n) is True):
+        if not ((0 <= k) and (k < n)):
             raise ValueError("must have 0 <= k < n, "
                 "got k = %s and n = %s" % (k, n))
         return C.cos(S.Pi*(k + 1)/(n + 1))
diff --git a/sympy/geometry/line.py b/sympy/geometry/line.py
index acb0c4caa94c..fad9391ec373 100644
--- a/sympy/geometry/line.py
+++ b/sympy/geometry/line.py
@@ -1434,8 +1434,8 @@ def contains(self, o):
                     rv = o.y >= self.source.y
                 else:
                     rv = o.y <= self.source.y
-                if isinstance(rv, bool):
-                    return rv
+                if rv == True or rv == False:
+                    return bool(rv)
                 raise Undecidable(
                     'Cannot determine if %s is in %s' % (o, self))
             else:
@@ -1504,9 +1504,9 @@ def __new__(cls, p1, p2, **kwargs):
         p2 = Point(p2)
         if p1 == p2:
             return Point(p1)
-        if (p1.x > p2.x) is True:
+        if (p1.x > p2.x) == True:
             p1, p2 = p2, p1
-        elif (p1.x == p2.x) is (p1.y > p2.y) is True:
+        elif (p1.x == p2.x) == True and (p1.y > p2.y) == True:
             p1, p2 = p2, p1
         return LinearEntity.__new__(cls, p1, p2, **kwargs)
 
diff --git a/sympy/integrals/meijerint.py b/sympy/integrals/meijerint.py
index 530eca174919..bcfae07be3d0 100644
--- a/sympy/integrals/meijerint.py
+++ b/sympy/integrals/meijerint.py
@@ -878,10 +878,10 @@ def pb(g):
     _, s = _get_coeff_exp(po, x)
     _, b1 = _get_coeff_exp(g1.argument, x)
     _, b2 = _get_coeff_exp(g2.argument, x)
-    if (b1 < 0) is True:
+    if (b1 < 0) == True:
         b1 = -b1
         g1 = _flip_g(g1)
-    if (b2 < 0) is True:
+    if (b2 < 0) == True:
         b2 = -b2
         g2 = _flip_g(g2)
     if not b1.is_Rational or not b2.is_Rational:
@@ -969,7 +969,7 @@ def lambda_s0(c1, c2):
         for a in g1.an:
             for b in g1.bm:
                 diff = a - b
-                if (diff > 0) is True and diff.is_integer:
+                if (diff > 0) == True and diff.is_integer:
                     c1 = False
 
     tmp = []
@@ -1637,7 +1637,7 @@ def _meijerint_indefinite_1(f, x):
 
         def tr(p):
             return [a + rho + 1 for a in p]
-        if any(b.is_integer and (b <= 0) is True for b in tr(g.bm)):
+        if any(b.is_integer and (b <= 0) == True for b in tr(g.bm)):
             r = -meijerg(
                 tr(g.an), tr(g.aother) + [1], tr(g.bm) + [0], tr(g.bother), t)
         else:
@@ -1770,7 +1770,7 @@ def meijerint_definite(f, x, a, b):
     results = []
     if b == oo:
         for split in _find_splitting_points(f, x):
-            if (a - split >= 0) is True:
+            if (a - split >= 0) == True:
                 _debug('Trying x --> x + %s' % split)
                 res = _meijerint_definite_2(f.subs(x, x + split)
                                             *Heaviside(x + split - a), x)
diff --git a/sympy/integrals/rationaltools.py b/sympy/integrals/rationaltools.py
index 1dff1f869552..5468fd7a9758 100644
--- a/sympy/integrals/rationaltools.py
+++ b/sympy/integrals/rationaltools.py
@@ -216,7 +216,7 @@ def ratint_logpart(f, g, x, t=None):
         R_map[r.degree()] = r
 
     def _include_sign(c, sqf):
-        if (c < 0) is True:
+        if (c < 0) == True:
             h, k = sqf[0]
             sqf[0] = h*c, k
 
diff --git a/sympy/integrals/transforms.py b/sympy/integrals/transforms.py
index 4746660d3225..99a6511aaebc 100644
--- a/sympy/integrals/transforms.py
+++ b/sympy/integrals/transforms.py
@@ -474,9 +474,9 @@ def left(c, is_numer):
             return c < b_
         if b_ is None:
             return c <= a_
-        if (c >= b_) is True:
+        if (c >= b_) == True:
             return False
-        if (c <= a_) is True:
+        if (c <= a_) == True:
             return True
         if is_numer:
             return None
@@ -620,8 +620,8 @@ def linear_arg(arg):
         elif isinstance(fact, gamma):
             a, b = linear_arg(fact.args[0])
             if is_numer:
-                if (a > 0 and (left(-b/a, is_numer) is False)) or \
-                   (a < 0 and (left(-b/a, is_numer) is True)):
+                if (a > 0 and (left(-b/a, is_numer) == False)) or \
+                   (a < 0 and (left(-b/a, is_numer) == True)):
                     raise NotImplementedError(
                         'Gammas partially over the strip.')
             ugammas += [(a, b)]
@@ -918,9 +918,9 @@ def bigger(ex1, ex2):
         n = power(ex2)
         if n is None:
             return None
-        if n > 0 and (abs(ex1) <= abs(a)**n) is True:
+        if n > 0 and (abs(ex1) <= abs(a)**n) == True:
             return False
-        if n < 0 and (abs(ex1) >= abs(a)**n) is True:
+        if n < 0 and (abs(ex1) >= abs(a)**n) == True:
             return True
 
     def replie(x, y):
@@ -934,13 +934,14 @@ def replie(x, y):
         return (x < y)
 
     def replue(x, y):
-        if bigger(x, y) in (True, False):
+        b = bigger(x, y)
+        if b == True or b == False:
             return True
         return Unequality(x, y)
 
     def repl(ex, *args):
-        if isinstance(ex, bool):
-            return ex
+        if ex == True or ex == False:
+            return bool(ex)
         return ex.replace(*args)
     expr = repl(expr, StrictLessThan, replie)
     expr = repl(expr, StrictGreaterThan, lambda x, y: replie(y, x))
@@ -1027,7 +1028,7 @@ def process_conds(conds):
     conds = conds2
 
     def cnt(expr):
-        if isinstance(expr, bool):
+        if expr == True or expr == False:
             return 0
         return expr.count_ops()
     conds.sort(key=lambda x: (-x[0], cnt(x[1])))
@@ -1037,8 +1038,8 @@ def cnt(expr):
     a, aux = conds[0]
 
     def sbs(expr):
-        if isinstance(expr, bool):
-            return expr
+        if expr == S.true or expr == S.false:
+            return bool(expr)
         return expr.subs(s, s_)
     if simplify:
         F = _simplifyconds(F, s, a)
diff --git a/sympy/matrices/expressions/blockmatrix.py b/sympy/matrices/expressions/blockmatrix.py
index 474a939dc08c..859b01acc46e 100644
--- a/sympy/matrices/expressions/blockmatrix.py
+++ b/sympy/matrices/expressions/blockmatrix.py
@@ -150,12 +150,12 @@ def transpose(self):
     def _entry(self, i, j):
         # Find row entry
         for row_block, numrows in enumerate(self.rowblocksizes):
-            if (i < numrows) is not False:
+            if (i < numrows) != False:
                 break
             else:
                 i -= numrows
         for col_block, numcols in enumerate(self.colblocksizes):
-            if (j < numcols) is not False:
+            if (j < numcols) != False:
                 break
             else:
                 j -= numcols
diff --git a/sympy/matrices/expressions/matexpr.py b/sympy/matrices/expressions/matexpr.py
index 234b2b52f0b6..5ced97a80277 100644
--- a/sympy/matrices/expressions/matexpr.py
+++ b/sympy/matrices/expressions/matexpr.py
@@ -196,8 +196,8 @@ def valid_index(self, i, j):
         def is_valid(idx):
             return isinstance(idx, (int, Integer, Symbol, Expr))
         return (is_valid(i) and is_valid(j) and
-                (0 <= i) is not False and (i < self.rows) is not False and
-                (0 <= j) is not False and (j < self.cols) is not False)
+                (0 <= i) != False and (i < self.rows) != False and
+                (0 <= j) != False and (j < self.cols) != False)
 
     def __getitem__(self, key):
         if not isinstance(key, tuple) and isinstance(key, slice):
@@ -209,7 +209,7 @@ def __getitem__(self, key):
                 from sympy.matrices.expressions.slice import MatrixSlice
                 return MatrixSlice(self, i, j)
             i, j = sympify(i), sympify(j)
-            if self.valid_index(i, j) is not False:
+            if self.valid_index(i, j) != False:
                 return self._entry(i, j)
             else:
                 raise IndexError("Invalid indices (%s, %s)" % (i, j))
diff --git a/sympy/physics/quantum/state.py b/sympy/physics/quantum/state.py
index 78984f737c68..00b371980bff 100644
--- a/sympy/physics/quantum/state.py
+++ b/sympy/physics/quantum/state.py
@@ -735,7 +735,7 @@ def __call__(self, *args, **options):
                      or upper in args[ct].free_symbols):
                 continue
 
-            if (args[ct] < lower) is True or (args[ct] > upper) is True:
+            if (args[ct] < lower) == True or (args[ct] > upper) == True:
                 return 0
 
             ct += 1
diff --git a/sympy/physics/tests/test_clebsch_gordan.py b/sympy/physics/tests/test_clebsch_gordan.py
index 2e026c4e7b23..b3e718dcd82c 100644
--- a/sympy/physics/tests/test_clebsch_gordan.py
+++ b/sympy/physics/tests/test_clebsch_gordan.py
@@ -106,7 +106,7 @@ def test_clebsch_gordan2():
 
 def test_wigner():
     def tn(a, b):
-        return abs((a - b).n(64) < S('1e-64'))
+        return (a - b).n(64) < S('1e-64')
     assert tn(wigner_9j(1, 1, 1, 1, 1, 1, 1, 1, 0, prec=64), S(1)/18)
     assert wigner_9j(3, 3, 2, 3, 3, 2, 3, 3, 2) == 3221*sqrt(
         70)/(246960*sqrt(105)) - 365/(3528*sqrt(70)*sqrt(105))
@@ -116,7 +116,7 @@ def tn(a, b):
 
 def test_gaunt():
     def tn(a, b):
-        return abs((a - b).n(64) < S('1e-64'))
+        return (a - b).n(64) < S('1e-64')
     assert gaunt(1, 0, 1, 1, 0, -1) == -1/(2*sqrt(pi))
     assert tn(gaunt(
         10, 10, 12, 9, 3, -12, prec=64), (-S(98)/62031) * sqrt(6279)/sqrt(pi))
diff --git a/sympy/polys/polyroots.py b/sympy/polys/polyroots.py
index 15e09833d001..a88671ff60a5 100644
--- a/sympy/polys/polyroots.py
+++ b/sympy/polys/polyroots.py
@@ -106,7 +106,7 @@ def roots_cubic(f):
             return [-aon3]*3
         else:
             if q.is_real:
-                if (q > 0) is True:
+                if (q > 0) == True:
                     u1 = -q**Rational(1, 3)
                 else:
                     u1 = (-q)**Rational(1, 3)
diff --git a/sympy/polys/tests/test_polytools.py b/sympy/polys/tests/test_polytools.py
index f7b49a419c3d..7b4919c22323 100644
--- a/sympy/polys/tests/test_polytools.py
+++ b/sympy/polys/tests/test_polytools.py
@@ -2675,29 +2675,29 @@ def test_nroots():
     roots = nroots(x**5 + x + 1, n=5)
     eps = Float("1e-5")
 
-    assert re(roots[0]).epsilon_eq(-0.75487, eps) is True
+    assert re(roots[0]).epsilon_eq(-0.75487, eps) is S.true
     assert im(roots[0]) == 0.0
     assert re(roots[1]) == -0.5
-    assert im(roots[1]).epsilon_eq(-0.86602, eps) is True
+    assert im(roots[1]).epsilon_eq(-0.86602, eps) is S.true
     assert re(roots[2]) == -0.5
-    assert im(roots[2]).epsilon_eq(+0.86602, eps) is True
-    assert re(roots[3]).epsilon_eq(+0.87743, eps) is True
-    assert im(roots[3]).epsilon_eq(-0.74486, eps) is True
-    assert re(roots[4]).epsilon_eq(+0.87743, eps) is True
-    assert im(roots[4]).epsilon_eq(+0.74486, eps) is True
+    assert im(roots[2]).epsilon_eq(+0.86602, eps) is S.true
+    assert re(roots[3]).epsilon_eq(+0.87743, eps) is S.true
+    assert im(roots[3]).epsilon_eq(-0.74486, eps) is S.true
+    assert re(roots[4]).epsilon_eq(+0.87743, eps) is S.true
+    assert im(roots[4]).epsilon_eq(+0.74486, eps) is S.true
 
     eps = Float("1e-6")
 
-    assert re(roots[0]).epsilon_eq(-0.75487, eps) is False
+    assert re(roots[0]).epsilon_eq(-0.75487, eps) is S.false
     assert im(roots[0]) == 0.0
     assert re(roots[1]) == -0.5
-    assert im(roots[1]).epsilon_eq(-0.86602, eps) is False
+    assert im(roots[1]).epsilon_eq(-0.86602, eps) is S.false
     assert re(roots[2]) == -0.5
-    assert im(roots[2]).epsilon_eq(+0.86602, eps) is False
-    assert re(roots[3]).epsilon_eq(+0.87743, eps) is False
-    assert im(roots[3]).epsilon_eq(-0.74486, eps) is False
-    assert re(roots[4]).epsilon_eq(+0.87743, eps) is False
-    assert im(roots[4]).epsilon_eq(+0.74486, eps) is False
+    assert im(roots[2]).epsilon_eq(+0.86602, eps) is S.false
+    assert re(roots[3]).epsilon_eq(+0.87743, eps) is S.false
+    assert im(roots[3]).epsilon_eq(-0.74486, eps) is S.false
+    assert re(roots[4]).epsilon_eq(+0.87743, eps) is S.false
+    assert im(roots[4]).epsilon_eq(+0.74486, eps) is S.false
 
     raises(DomainError, lambda: Poly(x + y, x).nroots())
     raises(MultivariatePolynomialError, lambda: Poly(x + y).nroots())
diff --git a/sympy/polys/tests/test_rootoftools.py b/sympy/polys/tests/test_rootoftools.py
index 8eb8bad340fe..b354e55bebcd 100644
--- a/sympy/polys/tests/test_rootoftools.py
+++ b/sympy/polys/tests/test_rootoftools.py
@@ -317,9 +317,9 @@ def test_RootSum_evalf():
     rs = RootSum(x**2 + 1, exp)
 
     assert rs.evalf(n=20, chop=True).epsilon_eq(
-        Float("1.0806046117362794348", 20), Float("1e-20")) is True
+        Float("1.0806046117362794348", 20), Float("1e-20")) is S.true
     assert rs.evalf(n=15, chop=True).epsilon_eq(
-        Float("1.08060461173628", 15), Float("1e-15")) is True
+        Float("1.08060461173628", 15), Float("1e-15")) is S.true
 
     rs = RootSum(x**2 + a, exp, x)
 
diff --git a/sympy/simplify/fu.py b/sympy/simplify/fu.py
index 692d2cfec360..52c2fe956175 100644
--- a/sympy/simplify/fu.py
+++ b/sympy/simplify/fu.py
@@ -496,9 +496,9 @@ def _f(rv):
         if not (rv.is_Pow and rv.base.func == f):
             return rv
 
-        if (rv.exp < 0) is True:
+        if (rv.exp < 0) == True:
             return rv
-        if (rv.exp > max) is True:
+        if (rv.exp > max) == True:
             return rv
         if rv.exp == 2:
             return h(g(rv.base.args[0])**2)
diff --git a/sympy/simplify/hyperexpand.py b/sympy/simplify/hyperexpand.py
index 4f340fe8bee6..b99006f19285 100644
--- a/sympy/simplify/hyperexpand.py
+++ b/sympy/simplify/hyperexpand.py
@@ -586,7 +586,7 @@ def _is_suitable_origin(self):
         """
         for a in self.ap:
             for b in self.bq:
-                if (a - b).is_integer and (a < b) is False:
+                if (a - b).is_integer and (a < b) == False:
                     return False
         for a in self.ap:
             if a == 0:
diff --git a/sympy/simplify/simplify.py b/sympy/simplify/simplify.py
index 5aff8ab654d9..7b2d716cbb49 100644
--- a/sympy/simplify/simplify.py
+++ b/sympy/simplify/simplify.py
@@ -2325,7 +2325,7 @@ def nc_gcd(aa, bb):
     if glogb.func is C.log or not glogb.is_Mul:
         if glogb.args[0].is_Pow or glogb.args[0].func is exp:
             glogb = _denest_pow(glogb.args[0])
-            if (abs(glogb.exp) < 1) is True:
+            if (abs(glogb.exp) < 1) == True:
                 return Pow(glogb.base, glogb.exp*e)
         return eq
 
diff --git a/sympy/solvers/inequalities.py b/sympy/solvers/inequalities.py
index baae7bca2f38..c40ff7e81ac5 100644
--- a/sympy/solvers/inequalities.py
+++ b/sympy/solvers/inequalities.py
@@ -461,11 +461,10 @@ def reduce_inequalities(inequalities, assume=True, symbols=[]):
     poly_part, abs_part, extra_assume = {}, {}, []
 
     for inequality in inequalities:
-        if isinstance(inequality, bool):
-            if inequality is False:
-                return False
-            else:
-                continue
+        if inequality == True:
+            continue
+        elif inequality == False:
+            return False
 
         if isinstance(inequality, AppliedPredicate):
             extra_assume.append(inequality)
diff --git a/sympy/solvers/recurr.py b/sympy/solvers/recurr.py
index 8e3b787a76fd..9c38fef45a42 100644
--- a/sympy/solvers/recurr.py
+++ b/sympy/solvers/recurr.py
@@ -629,7 +629,7 @@ def rsolve_hyper(coeffs, f, n, **hints):
                 # start the product with the term y(n_root + 1).
                 n0 = 0
                 for n_root in roots(ratio.as_numer_denom()[1], n).keys():
-                    if (n0 < (n_root + 1)) is True:
+                    if (n0 < (n_root + 1)) == True:
                         n0 = n_root + 1
                 K = product(ratio, (n, n0, n - 1))
                 if K.has(factorial, FallingFactorial, RisingFactorial):
diff --git a/sympy/solvers/solvers.py b/sympy/solvers/solvers.py
index 1af6bb2abfdb..008c0909d3bf 100644
--- a/sympy/solvers/solvers.py
+++ b/sympy/solvers/solvers.py
@@ -282,7 +282,7 @@ def checksol(f, symbol, sol=None, **flags):
         elif val.is_Rational:
             return val == 0
         if numerical and not val.free_symbols:
-            return abs(val.n(18).n(12, chop=True)) < 1e-9
+            return bool(abs(val.n(18).n(12, chop=True)) < 1e-9)
         was = val
 
     if flags.get('warn', False):
@@ -2714,7 +2714,7 @@ def _canonical(eq):
         # make sign canonical
         free = eq.free_symbols
         if len(free) == 1:
-            if (eq.coeff(free.pop()**degree(eq)) < 0) is True:
+            if (eq.coeff(free.pop()**degree(eq)) < 0) == True:
                 eq = -eq
         elif eq.could_extract_minus_sign():
             eq = -eq
diff --git a/sympy/solvers/tests/test_inequalities.py b/sympy/solvers/tests/test_inequalities.py
index ca60e707587d..0fd89d25fd49 100644
--- a/sympy/solvers/tests/test_inequalities.py
+++ b/sympy/solvers/tests/test_inequalities.py
@@ -70,9 +70,9 @@ def test_reduce_poly_inequalities_real_relational():
         assert reduce_rational_inequalities(
             [[Le(x**2, 0)]], x, relational=True) == Eq(x, 0)
         assert reduce_rational_inequalities(
-            [[Lt(x**2, 0)]], x, relational=True) is False
+            [[Lt(x**2, 0)]], x, relational=True) == False
         assert reduce_rational_inequalities(
-            [[Ge(x**2, 0)]], x, relational=True) is True
+            [[Ge(x**2, 0)]], x, relational=True) == True
         assert reduce_rational_inequalities(
             [[Gt(x**2, 0)]], x, relational=True) == Or(Lt(x, 0), Gt(x, 0))
         assert reduce_rational_inequalities(