Repository: sympy/sympy
Pull request: remove `Eq.rewrite(Add)`
<!-- Your title above should be a short description of what
was changed. Do not include the issue number in the title. -->

#### References to other Issues or PRs
<!-- If this pull request fixes an issue, write "Fixes #NNNN" in that exact
format, e.g. "Fixes #1234" (see
https://tinyurl.com/auto-closing for more information). Also, please
write a comment on that issue linking back to this pull request once it is
open. -->

closes #25701 

#### Brief description of what is fixed or changed

Remove `rewrite(Add)` which is misplaced in the rewrite machinery.

#### Other comments

From @oscarbenjamin in #25701
> I have said before that this rewrite rule is a mistake: rewrite is not expected to turn one kind of expression into a different kind and in this case a Boolean should not be rewritten to an Expr.
> 

This rewrite method was added as a compromise since the natural method `as_expr()` was frowned on since `Eq` aren't expressions. The problem is that it was added into machinery that wants to apply rewrite rules to an entire expression tree. And in the context of `rewrite(Add)` which only applies to `Eq` it doesn't make sense to apply this to a tree.

There are a few solutions:
1) add a guard so this can't be run on anything but `Eq`. The fact that we have to do so, however, indicates that this is a design flaw. So do one of the following:
2) Agree on a name of a method for Eq that return lhs - rhs and gives the option of doing so in an unevaluated manner. 
- `rhs0` could return Eq(Add(lhs, -rhs, evaluate=False), 0, evaluate=False) and the user could grab the lhs as in `eq.rhs0.lhs` or
- `lhs_rhs()` could return the diff in an evaluated manner or `lhs_rhs(evaluate=False)` could do so in an unevaluated manner.
3) make the user reconstruct the Add on their own: `Add(eq.lhs, -eq.rhs, evaluate=False)` 

This PR implements #2 and provides a pain-free way of updating the code-base since the name changes but the functionality is identically the same.

#### Release Notes

<!-- Write the release notes for this release below between the BEGIN and END
statements. The basic format is a bulleted list with the name of the subpackage
and the release note for this PR. For example:

* solvers
  * Added a new solver for logarithmic equations.

* functions
  * Fixed a bug with log of integers. Formerly, `log(-x)` incorrectly gave `-log(x)`.

* physics.units
  * Corrected a semantical error in the conversion between volt and statvolt which
    reported the volt as being larger than the statvolt.

or if no release note(s) should be included use:

NO ENTRY

See https://github.com/sympy/sympy/wiki/Writing-Release-Notes for more
information on how to write release notes. The bot will check your release
notes automatically to see if they are formatted correctly. -->

<!-- BEGIN RELEASE NOTES -->
* core
  * BREAKING CHANGE `eq.rewrite(Add)` should be replaced with `eq.lhs - eq.rhs`
<!-- END RELEASE NOTES -->


Problem statement:
TypeError on Eq(2*sign(x + 3)/(5*Abs(x + 3)**(3/5)), 0)
```python
from sympy import Symbol, Eq, sign, Abs, Rational

x = Symbol('x', real=True)

eq = Eq(2*sign(x + 3)/(5*Abs(x + 3)**Rational(3, 5)), 0)

eq.simplify()

```

Result:
```
Traceback (most recent call last):
  File "/home/satels/workenv/krcore/lib/python3.10/site-packages/IPython/core/interactiveshell.py", line 3508, in run_code
    exec(code_obj, self.user_global_ns, self.user_ns)
  File "<ipython-input-1-2fbb666c6791>", line 7, in <module>
    eq.simplify()
  File "/home/satels/workenv/krcore/lib/python3.10/site-packages/sympy/core/basic.py", line 1853, in simplify
    return simplify(self, **kwargs)
  File "/home/satels/workenv/krcore/lib/python3.10/site-packages/sympy/simplify/simplify.py", line 601, in simplify
    return _eval_simplify(**kwargs)
  File "/home/satels/workenv/krcore/lib/python3.10/site-packages/sympy/core/relational.py", line 692, in _eval_simplify
    e.rewrite(Add, evaluate=False), x)
  File "/home/satels/workenv/krcore/lib/python3.10/site-packages/sympy/core/basic.py", line 1981, in rewrite
    return self._rewrite(pattern, rule, method, **hints)
  File "/home/satels/workenv/krcore/lib/python3.10/site-packages/sympy/core/basic.py", line 1986, in _rewrite
    args = [a._rewrite(pattern, rule, method, **hints)
  File "/home/satels/workenv/krcore/lib/python3.10/site-packages/sympy/core/basic.py", line 1986, in <listcomp>
    args = [a._rewrite(pattern, rule, method, **hints)
  File "/home/satels/workenv/krcore/lib/python3.10/site-packages/sympy/core/basic.py", line 1986, in _rewrite
    args = [a._rewrite(pattern, rule, method, **hints)
  File "/home/satels/workenv/krcore/lib/python3.10/site-packages/sympy/core/basic.py", line 1986, in <listcomp>
    args = [a._rewrite(pattern, rule, method, **hints)
  File "/home/satels/workenv/krcore/lib/python3.10/site-packages/sympy/core/basic.py", line 2000, in _rewrite
    return self.func(*args)
  File "/home/satels/workenv/krcore/lib/python3.10/site-packages/sympy/functions/elementary/piecewise.py", line 33, in __new__
    raise TypeError(filldedent('''
TypeError: 
Second argument must be a Boolean, not `Add`
```

Hints from issue discussion:
This is because an `Eq` is rewritten as an `Add`:
```python
In [4]: Eq(x, 2).rewrite(Add)
Out[4]: x - 2
```
I have said before that this rewrite rule is a mistake: rewrite is not expected to turn one kind of expression into a different kind and in this case a Boolean should not be rewritten to an Expr.

This is a fix to the problem:
```diff
diff --git a/sympy/core/relational.py b/sympy/core/relational.py
index b2dbd1f395..a534926411 100644
--- a/sympy/core/relational.py
+++ b/sympy/core/relational.py
@@ -686,8 +686,7 @@ def _eval_simplify(self, **kwargs):
                 from .add import Add
                 from sympy.solvers.solveset import linear_coeffs
                 x = free.pop()
-                m, b = linear_coeffs(
-                    e.rewrite(Add, evaluate=False), x)
+                m, b = linear_coeffs(Add(e.lhs, -e.rhs, evaluate=False), x)
                 if m.is_zero is False:
                     enew = e.func(x, -b / m)
                 else:
```
With the diff above we have:
```python
In [4]: eq.simplify()
Out[4]: 
⎧      0         for x = -3    
⎪                              
⎪  2⋅(x + 3)                   
⎨──────────────  otherwise  = 0
⎪  │       8/5│                
⎪5⋅│(x + 3)   │                
⎩   
```
This is not completely correct but the basic bug here is fixed. The reason that this is not completely correct is that when `x = -3` the lhs should be undefined rather than `0`.

The `Eq.rewrite(Add)` should be removed entirely to fix that problem in general.

Commit message: e.rewrite(Add) -> e.lhs_rhs()

Patch to review:
diff --git a/sympy/calculus/util.py b/sympy/calculus/util.py
index e140396d8cd7..2e9bfab05108 100644
--- a/sympy/calculus/util.py
+++ b/sympy/calculus/util.py
@@ -422,7 +422,7 @@ def _check(orig_f, period):
     period = None
 
     if isinstance(f, Relational):
-        f = f.lhs - f.rhs
+        f = f.lhs_rhs()
 
     f = f.simplify()
 
diff --git a/sympy/core/expr.py b/sympy/core/expr.py
index ee4ae39abee1..457a2fba55eb 100644
--- a/sympy/core/expr.py
+++ b/sympy/core/expr.py
@@ -68,6 +68,9 @@ class Expr(Basic, EvalfMixin):
 
     is_scalar = True  # self derivative is 1
 
+    def lhs_rhs(self, evaluate=True):
+        return self  # only for Eq
+
     @property
     def _diff_wrt(self):
         """Return True if one can differentiate with respect to this
diff --git a/sympy/core/relational.py b/sympy/core/relational.py
index b2dbd1f39552..10b107765a0b 100644
--- a/sympy/core/relational.py
+++ b/sympy/core/relational.py
@@ -419,7 +419,7 @@ def _eval_simplify(self, **kwargs):
         if r.is_Relational:
             if not isinstance(r.lhs, Expr) or not isinstance(r.rhs, Expr):
                 return r
-            dif = r.lhs - r.rhs
+            dif = r.lhs_rhs()
             # replace dif with a valid Number that will
             # allow a definitive comparison with 0
             v = None
@@ -437,7 +437,7 @@ def _eval_simplify(self, **kwargs):
                 try:
                     from sympy.solvers.solveset import linear_coeffs
                     x = free.pop()
-                    dif = r.lhs - r.rhs
+                    dif = r.lhs_rhs()
                     m, b = linear_coeffs(dif, x)
                     if m.is_zero is False:
                         if m.is_negative:
@@ -467,7 +467,7 @@ def _eval_simplify(self, **kwargs):
                     from sympy.solvers.solveset import linear_coeffs
                     from sympy.polys.polytools import gcd
                     free = list(ordered(free))
-                    dif = r.lhs - r.rhs
+                    dif = r.lhs_rhs()
                     m = linear_coeffs(dif, *free)
                     constant = m[-1]
                     del m[-1]
@@ -594,7 +594,7 @@ class Equality(Relational):
 
     Since this object is already an expression, it does not respond to
     the method ``as_expr`` if one tries to create `x - y` from ``Eq(x, y)``.
-    This can be done with the ``rewrite(Add)`` method.
+    This can be done with the ``lhs_rhs()`` method.
 
     .. deprecated:: 1.5
 
@@ -626,7 +626,7 @@ def __new__(cls, lhs, rhs, **options):
     def _eval_relation(cls, lhs, rhs):
         return _sympify(lhs == rhs)
 
-    def _eval_rewrite_as_Add(self, L, R, evaluate=True, **kwargs):
+    def lhs_rhs(self, evaluate=True):
         """
         return Eq(L, R) as L - R. To control the evaluation of
         the result set pass `evaluate=True` to give L - R;
@@ -638,17 +638,18 @@ def _eval_rewrite_as_Add(self, L, R, evaluate=True, **kwargs):
         Examples
         ========
 
-        >>> from sympy import Eq, Add
+        >>> from sympy import Eq
         >>> from sympy.abc import b, x
         >>> eq = Eq(x + b, x - b)
-        >>> eq.rewrite(Add)
+        >>> eq.lhs_rhs()
         2*b
-        >>> eq.rewrite(Add, evaluate=None).args
+        >>> eq.lhs_rhs(evaluate=None).args
         (b, b, x, -x)
-        >>> eq.rewrite(Add, evaluate=False).args
+        >>> eq.lhs_rhs(evaluate=False).args
         (b, x, b, -x)
         """
         from .add import _unevaluated_Add, Add
+        L, R = self.args
         if L == 0:
             return R
         if R == 0:
@@ -683,11 +684,10 @@ def _eval_simplify(self, **kwargs):
         free = self.free_symbols
         if len(free) == 1:
             try:
-                from .add import Add
                 from sympy.solvers.solveset import linear_coeffs
                 x = free.pop()
                 m, b = linear_coeffs(
-                    e.rewrite(Add, evaluate=False), x)
+                    e.lhs_rhs(evaluate=False), x)
                 if m.is_zero is False:
                     enew = e.func(x, -b / m)
                 else:
diff --git a/sympy/core/tests/test_relational.py b/sympy/core/tests/test_relational.py
index a25650fce20f..3b1ffe9b5a39 100644
--- a/sympy/core/tests/test_relational.py
+++ b/sympy/core/tests/test_relational.py
@@ -12,6 +12,8 @@
 from sympy.core.power import Pow
 from sympy.core.singleton import S
 from sympy.core.symbol import (Symbol, symbols)
+from sympy.functions.elementary.complexes import sign, Abs
+from sympy.functions.elementary.piecewise import Piecewise
 from sympy.functions.elementary.exponential import (exp, exp_polar, log)
 from sympy.functions.elementary.integers import (ceiling, floor)
 from sympy.functions.elementary.miscellaneous import sqrt
@@ -558,14 +560,14 @@ def test_x_minus_y_not_same_as_x_lt_y():
 
     ineq = Lt(x, y, evaluate=False)
     raises(TypeError, lambda: ineq.doit())
-    assert ineq.lhs - ineq.rhs < 0
+    assert ineq.lhs_rhs() < 0
 
     t = Symbol('t', imaginary=True)
     x = 2 + t
     y = 3 + t
     ineq = Lt(x, y, evaluate=False)
     raises(TypeError, lambda: ineq.doit())
-    assert ineq.lhs - ineq.rhs < 0
+    assert ineq.lhs_rhs() < 0
 
     # this one should give error either way
     x = I + 2
@@ -1026,18 +1028,32 @@ def test_rel_args():
                 raises(TypeError, lambda: Relational(b, v, op))
 
 
-def test_Equality_rewrite_as_Add():
+def test_Equality_lhs_rhs():
+    # XXX we can be strict and make every routine that want to treat
+    # Eq like Eqn for purpose of solving explicitly test for Eq
+    #     if isinstance(eq, Eq):
+    #        eq = eq.lhs_rhs()
+    # or allow Expr to pass
+    #     eq = eq.lhs_rhs() # valid for Eq or Expr
+    assert (x + 1).lhs_rhs() == x + 1
+
     eq = Eq(x + y, y - x)
-    assert eq.rewrite(Add) == 2*x
-    assert eq.rewrite(Add, evaluate=None).args == (x, x, y, -y)
-    assert eq.rewrite(Add, evaluate=False).args == (x, y, x, -y)
+    assert eq.lhs_rhs() == 2*x
+    assert eq.lhs_rhs(evaluate=None).args == (x, x, y, -y)
+    assert eq.lhs_rhs(evaluate=False).args == (x, y, x, -y)
     for e in (True, False, None):
-        assert Eq(x, 0, evaluate=e).rewrite(Add) == x
-        assert Eq(0, x, evaluate=e).rewrite(Add) == x
+        assert Eq(x, 0, evaluate=e).lhs_rhs() == x
+        assert Eq(0, x, evaluate=e).lhs_rhs() == x
+
+    # issue 25701
+    r = symbols('r', real=True)
+    assert Eq(2*sign(r + 3)/(5*Abs(r + 3)**Rational(3, 5)), 0
+        ).simplify() == Eq(Piecewise(
+        (0, Eq(r, -3)), ((r + 3)/(5*Abs((r + 3)**Rational(8, 5)))*2, True)), 0)
 
 
 def test_issue_15847():
-    a = Ne(x*(x+y), x**2 + x*y)
+    a = Ne(x*(x + y), x**2 + x*y)
     assert simplify(a) == False
 
 
diff --git a/sympy/functions/elementary/piecewise.py b/sympy/functions/elementary/piecewise.py
index 016fe3f4e75a..67ae79e6195b 100644
--- a/sympy/functions/elementary/piecewise.py
+++ b/sympy/functions/elementary/piecewise.py
@@ -830,7 +830,7 @@ def __eval_cond(cls, cond):
             return True
         if isinstance(cond, Eq):
             try:
-                diff = cond.lhs - cond.rhs
+                diff = cond.lhs_rhs()
                 if diff.is_commutative:
                     return diff.is_zero
             except TypeError:
diff --git a/sympy/geometry/ellipse.py b/sympy/geometry/ellipse.py
index f096556f4bd7..5a76fb526370 100644
--- a/sympy/geometry/ellipse.py
+++ b/sympy/geometry/ellipse.py
@@ -1557,7 +1557,7 @@ def __new__(cls, *args, **kwargs):
             y = kwargs.get('y', 'y')
             equation = args[0].expand()
             if isinstance(equation, Eq):
-                equation = equation.lhs - equation.rhs
+                equation = equation.lhs_rhs()
             x = find(x, equation)
             y = find(y, equation)
 
diff --git a/sympy/geometry/line.py b/sympy/geometry/line.py
index e9b29b5766bd..08920623ea77 100644
--- a/sympy/geometry/line.py
+++ b/sympy/geometry/line.py
@@ -1192,7 +1192,7 @@ def __new__(cls, *args, **kwargs):
 
             equation = args[0]
             if isinstance(equation, Eq):
-                equation = equation.lhs - equation.rhs
+                equation = equation.lhs_rhs()
 
             def find_or_missing(x):
                 try:
diff --git a/sympy/logic/boolalg.py b/sympy/logic/boolalg.py
index ac190ac50883..429c0c4b676b 100644
--- a/sympy/logic/boolalg.py
+++ b/sympy/logic/boolalg.py
@@ -4,7 +4,6 @@
 
 from collections import defaultdict
 from itertools import chain, combinations, product, permutations
-from sympy.core.add import Add
 from sympy.core.basic import Basic
 from sympy.core.cache import cacheit
 from sympy.core.containers import Tuple
@@ -662,7 +661,7 @@ def _eval_simplify(self, **kwargs):
                     if (e.lhs != x or x in e.rhs.free_symbols) and x not in reps:
                         try:
                             m, b = linear_coeffs(
-                                e.rewrite(Add, evaluate=False), x)
+                                e.lhs_rhs(evaluate=False), x)
                             enew = e.func(x, -b/m)
                             if measure(enew) <= ratio*measure(e):
                                 e = enew
diff --git a/sympy/plotting/series.py b/sympy/plotting/series.py
index 731d1df89154..dd08748330e6 100644
--- a/sympy/plotting/series.py
+++ b/sympy/plotting/series.py
@@ -2485,7 +2485,7 @@ def _preprocess_meshgrid_expression(expr, adaptive):
         """
         equality = False
         if isinstance(expr, Equality):
-            expr = expr.lhs - expr.rhs
+            expr = expr.lhs_rhs()
             equality = True
         elif isinstance(expr, Relational):
             expr = expr.gts - expr.lts
diff --git a/sympy/polys/polyutils.py b/sympy/polys/polyutils.py
index 82c8c836191e..2a632f7736ff 100644
--- a/sympy/polys/polyutils.py
+++ b/sympy/polys/polyutils.py
@@ -190,7 +190,7 @@ def _parallel_dict_from_expr_if_gens(exprs, opt):
         poly = {}
 
         if expr.is_Equality:
-            expr = expr.lhs - expr.rhs
+            expr = expr.lhs_rhs()
 
         for term in Add.make_args(expr):
             coeff, monom = [], [0]*k
@@ -249,7 +249,7 @@ def _is_coeff(factor):
         terms = []
 
         if expr.is_Equality:
-            expr = expr.lhs - expr.rhs
+            expr = expr.lhs_rhs()
 
         for term in Add.make_args(expr):
             coeff, elements = [], {}
diff --git a/sympy/solvers/deutils.py b/sympy/solvers/deutils.py
index c968b65c8d51..64ca1e9fb511 100644
--- a/sympy/solvers/deutils.py
+++ b/sympy/solvers/deutils.py
@@ -173,7 +173,7 @@ def _desolve(eq, func=None, hint="default", ics=None, simplify=True, *, prep=Tru
     classify_pde(pde.py)
     """
     if isinstance(eq, Equality):
-        eq = eq.lhs - eq.rhs
+        eq = eq.lhs_rhs()
 
     # preprocess the equation and find func if not given
     if prep or func is None:
diff --git a/sympy/solvers/diophantine/diophantine.py b/sympy/solvers/diophantine/diophantine.py
index e9137170930f..9522482ee0cc 100644
--- a/sympy/solvers/diophantine/diophantine.py
+++ b/sympy/solvers/diophantine/diophantine.py
@@ -1335,7 +1335,7 @@ def diophantine(eq, param=symbols("t", integer=True), syms=None,
     eq = _sympify(eq)
 
     if isinstance(eq, Eq):
-        eq = eq.lhs - eq.rhs
+        eq = eq.lhs_rhs()
 
     try:
         var = list(eq.expand(force=True).free_symbols)
diff --git a/sympy/solvers/inequalities.py b/sympy/solvers/inequalities.py
index 246a752c3106..07f8e1297f64 100644
--- a/sympy/solvers/inequalities.py
+++ b/sympy/solvers/inequalities.py
@@ -233,7 +233,7 @@ def reduce_rational_inequalities(exprs, gen, relational=True):
                 expr, rel = expr
             else:
                 if expr.is_Relational:
-                    expr, rel = expr.lhs - expr.rhs, expr.rel_op
+                    expr, rel = expr.lhs_rhs(), expr.rel_op
                 else:
                     expr, rel = expr, '=='
 
@@ -479,7 +479,7 @@ def solve_univariate_inequality(expr, gen, relational=True, domain=S.Reals, cont
         rv = S.EmptySet
 
     else:
-        e = expr.lhs - expr.rhs
+        e = expr.lhs_rhs()
         period = periodicity(e, gen)
         if period == S.Zero:
             e = expand_mul(e)
@@ -805,7 +805,7 @@ def classify(ie, s, i):
 
     rv = None
     oo = S.Infinity
-    expr = ie.lhs - ie.rhs
+    expr = ie.lhs_rhs()
     try:
         p = Poly(expr, s)
         if p.degree() == 0:
diff --git a/sympy/solvers/ode/lie_group.py b/sympy/solvers/ode/lie_group.py
index 329b2d5d30d4..90bfeb8fd802 100644
--- a/sympy/solvers/ode/lie_group.py
+++ b/sympy/solvers/ode/lie_group.py
@@ -230,7 +230,7 @@ def infinitesimals(eq, func=None, order=None, hint='default', match=None):
     """
 
     if isinstance(eq, Equality):
-        eq = eq.lhs - eq.rhs
+        eq = eq.lhs_rhs()
     if not func:
         eq, func = _preprocess(eq)
     variables = func.args
diff --git a/sympy/solvers/ode/ode.py b/sympy/solvers/ode/ode.py
index 75bccdf8d56b..5a79876e876e 100644
--- a/sympy/solvers/ode/ode.py
+++ b/sympy/solvers/ode/ode.py
@@ -943,7 +943,7 @@ class in it.  Note that a hint may do this anyway if
         "work with functions of one variable, not %s" % func)
 
     if isinstance(eq, Equality):
-        eq = eq.lhs - eq.rhs
+        eq = eq.lhs_rhs()
 
     # Some methods want the unprocessed equation
     eq_orig = eq
@@ -1213,7 +1213,7 @@ def _sympify(eq):
     eq, funcs = (_sympify(w) for w in [eq, funcs])
     for i, fi in enumerate(eq):
         if isinstance(fi, Equality):
-            eq[i] = fi.lhs - fi.rhs
+            eq[i] = fi.lhs_rhs()
 
     t = list(list(eq[0].atoms(Derivative))[0].atoms(Symbol))[0]
     matching_hints = {"no_of_equation":i+1}
@@ -2815,7 +2815,7 @@ def checkinfsol(eq, infinitesimals, func=None, order=None):
 
     """
     if isinstance(eq, Equality):
-        eq = eq.lhs - eq.rhs
+        eq = eq.lhs_rhs()
     if not func:
         eq, func = _preprocess(eq)
     variables = func.args
diff --git a/sympy/solvers/ode/riccati.py b/sympy/solvers/ode/riccati.py
index 2ef66ed0896d..6de7fd80e32d 100644
--- a/sympy/solvers/ode/riccati.py
+++ b/sympy/solvers/ode/riccati.py
@@ -295,7 +295,7 @@ def match_riccati(eq, f, x):
     """
     # Group terms based on f(x)
     if isinstance(eq, Eq):
-        eq = eq.lhs - eq.rhs
+        eq = eq.lhs_rhs()
     eq = eq.expand().collect(f(x))
     cf = eq.coeff(f(x).diff(x))
 
diff --git a/sympy/solvers/ode/subscheck.py b/sympy/solvers/ode/subscheck.py
index 6ac7fba7d364..9cd63189bb4b 100644
--- a/sympy/solvers/ode/subscheck.py
+++ b/sympy/solvers/ode/subscheck.py
@@ -177,7 +177,7 @@ def checkodesol(ode, sol, func=None, order='auto', solve_for_func=True):
         if testnum == 0:
             # First pass, try substituting a solved solution directly into the
             # ODE. This has the highest chance of succeeding.
-            ode_diff = ode.lhs - ode.rhs
+            ode_diff = ode.lhs_rhs()
 
             if sol.lhs == func:
                 s = sub_func_doit(ode_diff, func, sol.rhs)
@@ -221,7 +221,7 @@ def checkodesol(ode, sol, func=None, order='auto', solve_for_func=True):
                 diffsols = {0: sol.lhs}
             else:
                 diffsols = {}
-            sol = sol.lhs - sol.rhs
+            sol = sol.lhs_rhs()
             for i in range(1, order + 1):
                 # Differentiation is a linear operator, so there should always
                 # be 1 solution. Nonetheless, we test just to make sure.
@@ -345,7 +345,7 @@ def _sympify(eq):
     eqs = _sympify(eqs)
     for i in range(len(eqs)):
         if isinstance(eqs[i], Equality):
-            eqs[i] = eqs[i].lhs - eqs[i].rhs
+            eqs[i] = eqs[i].lhs_rhs()
     if func is None:
         funcs = []
         for eq in eqs:
diff --git a/sympy/solvers/ode/tests/test_riccati.py b/sympy/solvers/ode/tests/test_riccati.py
index 548a1ee5b5e8..563c15939464 100644
--- a/sympy/solvers/ode/tests/test_riccati.py
+++ b/sympy/solvers/ode/tests/test_riccati.py
@@ -717,7 +717,7 @@ def check_dummy_sol(eq, solse, dummy_sym):
     contains dummy symbols.
     """
     if isinstance(eq, Eq):
-        eq = eq.lhs - eq.rhs
+        eq = eq.lhs_rhs()
     _, funcs = match_riccati(eq, f, x)
 
     sols = solve_riccati(f(x), x, *funcs)
diff --git a/sympy/solvers/pde.py b/sympy/solvers/pde.py
index f54e1e3cb049..6109b7df8713 100644
--- a/sympy/solvers/pde.py
+++ b/sympy/solvers/pde.py
@@ -468,7 +468,7 @@ def checkpdesol(pde, sol, func=None, solve_for_func=True):
 
     # try direct substitution of the solution into the PDE and simplify
     if sol.lhs == func:
-        pde = pde.lhs - pde.rhs
+        pde = pde.lhs_rhs()
         s = simplify(pde.subs(func, sol.rhs).doit())
         return s is S.Zero, s
 
diff --git a/sympy/solvers/recurr.py b/sympy/solvers/recurr.py
index ba627bbd4cb0..c7ed7c981a71 100644
--- a/sympy/solvers/recurr.py
+++ b/sympy/solvers/recurr.py
@@ -728,7 +728,7 @@ def rsolve(f, y, init=None):
 
     """
     if isinstance(f, Equality):
-        f = f.lhs - f.rhs
+        f = f.lhs_rhs()
 
     n = y.args[0]
     k = Wild('k', exclude=(n,))
diff --git a/sympy/solvers/simplex.py b/sympy/solvers/simplex.py
index ea54fefb780e..d40ff6e8a77f 100644
--- a/sympy/solvers/simplex.py
+++ b/sympy/solvers/simplex.py
@@ -734,7 +734,7 @@ def _lp_matrices(objective, constraints):
     # change Eq(x, y) to x - y <= 0 and y - x <= 0
     for i in range(len(np)):
         if isinstance(np[i], Eq):
-            np[i] = np[i].lhs - np[i].rhs <= 0
+            np[i] = np[i].lhs_rhs() <= 0
             np.append(-np[i].lhs <= 0)
 
     # convert constraints to nonpositive expressions
diff --git a/sympy/solvers/solvers.py b/sympy/solvers/solvers.py
index b79eb737804d..c04f0a99a6cd 100644
--- a/sympy/solvers/solvers.py
+++ b/sympy/solvers/solvers.py
@@ -274,7 +274,7 @@ def checksol(f, symbol, sol=None, **flags):
             if not f.is_Boolean:
                 return
         else:
-            f = f.rewrite(Add, evaluate=False, deep=False)
+            f = f.lhs_rhs(evaluate=False, deep=False)
 
     if isinstance(f, BooleanAtom):
         return bool(f)
@@ -928,7 +928,7 @@ def _sympified_list(w):
     for i, fi in enumerate(f):
         if isinstance(fi, (Eq, Ne)):
             if 'ImmutableDenseMatrix' in [type(a).__name__ for a in fi.args]:
-                fi = fi.lhs - fi.rhs
+                fi = fi.lhs_rhs()
             else:
                 L, R = fi.args
                 if isinstance(R, BooleanAtom):
@@ -948,7 +948,7 @@ def _sympified_list(w):
                             is True or False.
                         '''))
                 else:
-                    fi = fi.rewrite(Add, evaluate=False, deep=False)
+                    fi = fi.lhs_rhs(evaluate=False, deep=False)
             f[i] = fi
 
         # *** dispatch and handle as a system of relationals
@@ -2420,7 +2420,7 @@ def solve_undetermined_coeffs(equ, coeffs, *syms, **flags):
         raise ValueError('must provide symbols for coeffs')
 
     if isinstance(equ, Eq):
-        eq = equ.lhs - equ.rhs
+        eq = equ.lhs_rhs()
     else:
         eq = equ
 
@@ -3053,14 +3053,14 @@ def nsolve(*args, dict=False, **kwargs):
         f = list(f)
         for i, fi in enumerate(f):
             if isinstance(fi, Eq):
-                f[i] = fi.lhs - fi.rhs
+                f[i] = fi.lhs_rhs()
         f = Matrix(f).T
     if iterable(x0):
         x0 = list(x0)
     if not isinstance(f, Matrix):
         # assume it's a SymPy expression
         if isinstance(f, Eq):
-            f = f.lhs - f.rhs
+            f = f.lhs_rhs()
         elif f.is_Relational:
             raise TypeError('nsolve cannot accept inequalities')
         syms = f.free_symbols
@@ -3428,7 +3428,7 @@ def _take(d):
     _take = flags.setdefault('_take', _take)
 
     if isinstance(eq, Eq):
-        eq = eq.lhs - eq.rhs  # XXX legacy Eq as Eqn support
+        eq = eq.lhs_rhs()  # XXX legacy Eq as Eqn support
     elif not isinstance(eq, Expr):
         return
 
diff --git a/sympy/solvers/solveset.py b/sympy/solvers/solveset.py
index dc95ca03d6fb..7902c0e9de65 100644
--- a/sympy/solvers/solveset.py
+++ b/sympy/solvers/solveset.py
@@ -1065,7 +1065,7 @@ def _solveset(f, symbol, domain, _check=False):
             solns = solver(expr, symbol, in_set)
             result += solns
     elif isinstance(f, Eq):
-        result = solver(Add(f.lhs, - f.rhs, evaluate=False), symbol, domain)
+        result = solver(f.lhs_rhs(evaluate=False), symbol, domain)
 
     elif f.is_Relational:
         from .inequalities import solve_univariate_inequality
@@ -3311,7 +3311,7 @@ def _solve_using_known_values(result, solver):
                             # list.
                             result.remove(res)
                     continue  # skip as it's independent of desired symbols
-                depen1, depen2 = (eq2.rewrite(Add)).as_independent(*unsolved_syms)
+                depen1, depen2 = (eq2.lhs_rhs()).as_independent(*unsolved_syms)
                 if (depen1.has(Abs) or depen2.has(Abs)) and solver == solveset_complex:
                     # Absolute values cannot be inverted in the
                     # complex domain
@@ -3528,7 +3528,7 @@ def _separate_poly_nonpoly(system, symbols):
         denominators.update(_simple_dens(eq, symbols))
         # Convert equality to expression
         if isinstance(eq, Equality):
-            eq = eq.rewrite(Add)
+            eq = eq.lhs_rhs()
         # try to remove sqrt and rational power
         without_radicals = unrad(simplify(eq), *symbols)
         if without_radicals:
diff --git a/sympy/solvers/tests/test_solvers.py b/sympy/solvers/tests/test_solvers.py
index cdf3d54d9e8b..fdd8663a36ad 100644
--- a/sympy/solvers/tests/test_solvers.py
+++ b/sympy/solvers/tests/test_solvers.py
@@ -2640,7 +2640,7 @@ def test_issue_10169():
 
 def test_solve_undetermined_coeffs_issue_23927():
     A, B, r, phi = symbols('A, B, r, phi')
-    eq = Eq(A*sin(t) + B*cos(t), r*sin(t - phi)).rewrite(Add).expand(trig=True)
+    eq = Eq(A*sin(t) + B*cos(t), r*sin(t - phi)).lhs_rhs().expand(trig=True)
     soln = solve_undetermined_coeffs(eq, (r, phi), t)
     assert soln == [{
         phi: 2*atan((A - sqrt(A**2 + B**2))/B),
diff --git a/sympy/solvers/tests/test_solveset.py b/sympy/solvers/tests/test_solveset.py
index a1f756a4aa33..d00d5faa4310 100644
--- a/sympy/solvers/tests/test_solveset.py
+++ b/sympy/solvers/tests/test_solveset.py
@@ -1,6 +1,5 @@
 from math import isclose
 
-from sympy.core.add import Add
 from sympy.core.containers import Tuple
 from sympy.core.function import (Function, Lambda, nfloat, diff)
 from sympy.core.mod import Mod
@@ -3223,7 +3222,7 @@ def test_issue_23318():
         Eq(x, 0.0015 * z),
         Eq(0.0015, 7845.32 * y / z),
     ]
-    eqs_expr = [eq.rewrite(Add) for eq in eqs_eq]
+    eqs_expr = [eq.lhs_rhs() for eq in eqs_eq]
 
     sol = {(266.97755814852, 0.0340301680681629, 177985.03876568)}
 
diff --git a/sympy/stats/crv.py b/sympy/stats/crv.py
index 36e1a26e1499..bf854acab4f4 100644
--- a/sympy/stats/crv.py
+++ b/sympy/stats/crv.py
@@ -418,7 +418,7 @@ def probability(self, condition, **kwargs):
         # by computing a density handled by density computation
         except NotImplementedError:
             from sympy.stats.rv import density
-            expr = condition.lhs - condition.rhs
+            expr = condition.lhs_rhs()
             if not is_random(expr):
                 dens = self.density
                 comp = condition.rhs
diff --git a/sympy/stats/drv.py b/sympy/stats/drv.py
index 13517e0f6dd3..4e804afa83cb 100644
--- a/sympy/stats/drv.py
+++ b/sympy/stats/drv.py
@@ -238,7 +238,7 @@ def probability(self, condition):
             prob = self.eval_prob(_domain)
         except NotImplementedError:
             from sympy.stats.rv import density
-            expr = condition.lhs - condition.rhs
+            expr = condition.lhs_rhs()
             dens = density(expr)
             if not isinstance(dens, DiscreteDistribution):
                 from sympy.stats.drv_types import DiscreteDistributionHandmade
diff --git a/sympy/stats/rv.py b/sympy/stats/rv.py
index 24d0aa8e90c6..963ec2c818c2 100644
--- a/sympy/stats/rv.py
+++ b/sympy/stats/rv.py
@@ -467,7 +467,7 @@ def probability(self, condition, **kwargs):
             return Mul(*[self.probability(arg) for arg in condition.args])
         elif isinstance(condition, Or): # they are independent
             return Add(*[self.probability(arg) for arg in condition.args])
-        expr = condition.lhs - condition.rhs
+        expr = condition.lhs_rhs()
         rvs = random_symbols(expr)
         dens = self.compute_density(expr)
         if any(pspace(rv).is_Continuous for rv in rvs):
diff --git a/sympy/vector/implicitregion.py b/sympy/vector/implicitregion.py
index fd608a07f044..5653c88feebd 100644
--- a/sympy/vector/implicitregion.py
+++ b/sympy/vector/implicitregion.py
@@ -59,7 +59,7 @@ def __new__(cls, variables, equation):
             variables = Tuple(*variables)
 
         if isinstance(equation, Eq):
-            equation = equation.lhs - equation.rhs
+            equation = equation.lhs_rhs()
 
         return super().__new__(cls, variables, equation)
 


Code review comments:
diff_hunk: @@ -422,7 +422,7 @@ def _check(orig_f, period):
     period = None
 
     if isinstance(f, Relational):
-        f = f.lhs - f.rhs
+        f = f.lhs_rhs()
path: sympy/calculus/util.py
text: @author:
```suggestion
        f = f.lhs - f.rhs
```

diff_hunk: @@ -437,7 +437,7 @@ def _eval_simplify(self, **kwargs):
                 try:
                     from sympy.solvers.solveset import linear_coeffs
                     x = free.pop()
-                    dif = r.lhs - r.rhs
+                    dif = r.lhs_rhs()
path: sympy/core/relational.py
text: @author:
```suggestion
                    dif = r.lhs - r.rhs
```

diff_hunk: @@ -467,7 +467,7 @@ def _eval_simplify(self, **kwargs):
                     from sympy.solvers.solveset import linear_coeffs
                     from sympy.polys.polytools import gcd
                     free = list(ordered(free))
-                    dif = r.lhs - r.rhs
+                    dif = r.lhs_rhs()
path: sympy/core/relational.py
text: @author:
```suggestion
                    dif = r.lhs - r.rhs
```

diff_hunk: @@ -594,7 +594,7 @@ class Equality(Relational):
 
     Since this object is already an expression, it does not respond to
     the method ``as_expr`` if one tries to create `x - y` from ``Eq(x, y)``.
-    This can be done with the ``rewrite(Add)`` method.
+    This can be done with the ``lhs_rhs()`` method.
path: sympy/core/relational.py
text: @author:
```suggestion
    This can be done for ``eq = Eq(x, y)`` with `eq.lhs - eq.rhs`.
```

diff_hunk: @@ -558,14 +560,14 @@ def test_x_minus_y_not_same_as_x_lt_y():
 
     ineq = Lt(x, y, evaluate=False)
     raises(TypeError, lambda: ineq.doit())
-    assert ineq.lhs - ineq.rhs < 0
+    assert ineq.lhs_rhs() < 0
path: sympy/core/tests/test_relational.py
text: @author:
```suggestion
    assert ineq.lhs - ineq.rhs < 0
```

diff_hunk: @@ -1026,18 +1028,32 @@ def test_rel_args():
                 raises(TypeError, lambda: Relational(b, v, op))
 
 
-def test_Equality_rewrite_as_Add():
+def test_Equality_lhs_rhs():
+    # XXX we can be strict and make every routine that want to treat
+    # Eq like Eqn for purpose of solving explicitly test for Eq
+    #     if isinstance(eq, Eq):
+    #        eq = eq.lhs_rhs()
+    # or allow Expr to pass
+    #     eq = eq.lhs_rhs() # valid for Eq or Expr
+    assert (x + 1).lhs_rhs() == x + 1
path: sympy/core/tests/test_relational.py
text: @author:
```suggestion
```

diff_hunk: @@ -1026,18 +1028,32 @@ def test_rel_args():
                 raises(TypeError, lambda: Relational(b, v, op))
 
 
-def test_Equality_rewrite_as_Add():
+def test_Equality_lhs_rhs():
+    # XXX we can be strict and make every routine that want to treat
+    # Eq like Eqn for purpose of solving explicitly test for Eq
+    #     if isinstance(eq, Eq):
+    #        eq = eq.lhs_rhs()
+    # or allow Expr to pass
+    #     eq = eq.lhs_rhs() # valid for Eq or Expr
+    assert (x + 1).lhs_rhs() == x + 1
+
     eq = Eq(x + y, y - x)
path: sympy/core/tests/test_relational.py
text: @author:
```suggestion
```

diff_hunk: @@ -1026,18 +1028,32 @@ def test_rel_args():
                 raises(TypeError, lambda: Relational(b, v, op))
 
 
-def test_Equality_rewrite_as_Add():
+def test_Equality_lhs_rhs():
+    # XXX we can be strict and make every routine that want to treat
+    # Eq like Eqn for purpose of solving explicitly test for Eq
+    #     if isinstance(eq, Eq):
+    #        eq = eq.lhs_rhs()
+    # or allow Expr to pass
+    #     eq = eq.lhs_rhs() # valid for Eq or Expr
+    assert (x + 1).lhs_rhs() == x + 1
+
     eq = Eq(x + y, y - x)
-    assert eq.rewrite(Add) == 2*x
-    assert eq.rewrite(Add, evaluate=None).args == (x, x, y, -y)
-    assert eq.rewrite(Add, evaluate=False).args == (x, y, x, -y)
+    assert eq.lhs_rhs() == 2*x
+    assert eq.lhs_rhs(evaluate=None).args == (x, x, y, -y)
+    assert eq.lhs_rhs(evaluate=False).args == (x, y, x, -y)
path: sympy/core/tests/test_relational.py
text: @author:
```suggestion
```

diff_hunk: @@ -1026,18 +1028,32 @@ def test_rel_args():
                 raises(TypeError, lambda: Relational(b, v, op))
 
 
-def test_Equality_rewrite_as_Add():
+def test_Equality_lhs_rhs():
+    # XXX we can be strict and make every routine that want to treat
+    # Eq like Eqn for purpose of solving explicitly test for Eq
+    #     if isinstance(eq, Eq):
+    #        eq = eq.lhs_rhs()
+    # or allow Expr to pass
+    #     eq = eq.lhs_rhs() # valid for Eq or Expr
+    assert (x + 1).lhs_rhs() == x + 1
+
     eq = Eq(x + y, y - x)
-    assert eq.rewrite(Add) == 2*x
-    assert eq.rewrite(Add, evaluate=None).args == (x, x, y, -y)
-    assert eq.rewrite(Add, evaluate=False).args == (x, y, x, -y)
+    assert eq.lhs_rhs() == 2*x
+    assert eq.lhs_rhs(evaluate=None).args == (x, x, y, -y)
+    assert eq.lhs_rhs(evaluate=False).args == (x, y, x, -y)
     for e in (True, False, None):
-        assert Eq(x, 0, evaluate=e).rewrite(Add) == x
-        assert Eq(0, x, evaluate=e).rewrite(Add) == x
+        assert Eq(x, 0, evaluate=e).lhs_rhs() == x
+        assert Eq(0, x, evaluate=e).lhs_rhs() == x
+
path: sympy/core/tests/test_relational.py
text: @author:
```suggestion
```

diff_hunk: @@ -830,7 +830,7 @@ def __eval_cond(cls, cond):
             return True
         if isinstance(cond, Eq):
             try:
-                diff = cond.lhs - cond.rhs
+                diff = cond.lhs_rhs()
path: sympy/functions/elementary/piecewise.py
text: @author:
```suggestion
                diff = cond.lhs - cond.rhs
```

diff_hunk: @@ -1557,7 +1557,7 @@ def __new__(cls, *args, **kwargs):
             y = kwargs.get('y', 'y')
             equation = args[0].expand()
             if isinstance(equation, Eq):
-                equation = equation.lhs - equation.rhs
+                equation = equation.lhs_rhs()
path: sympy/geometry/ellipse.py
text: @author:
```suggestion
               equation = equation.lhs - equation.rhs
```

@user1:
I don't see how
```
equation = (_:= equation).lhs - _.rhs
```
is better than the original:
```
equation = equation.lhs - equation.rhs
```

@author:
I was interrupted mid-edit

diff_hunk: @@ -1192,7 +1192,7 @@ def __new__(cls, *args, **kwargs):
 
             equation = args[0]
             if isinstance(equation, Eq):
-                equation = equation.lhs - equation.rhs
+                equation = equation.lhs_rhs()
path: sympy/geometry/line.py
text: @author:
```suggestion
                equation = equation.lhs - equation.rhs
```

diff_hunk: @@ -662,7 +661,7 @@ def _eval_simplify(self, **kwargs):
                     if (e.lhs != x or x in e.rhs.free_symbols) and x not in reps:
                         try:
                             m, b = linear_coeffs(
-                                e.rewrite(Add, evaluate=False), x)
+                                e.lhs_rhs(evaluate=False), x)
path: sympy/logic/boolalg.py
text: @author:
```suggestion
                                Add(e.lhs, -e.rhs, evaluate=False), x)
```

diff_hunk: @@ -2485,7 +2485,7 @@ def _preprocess_meshgrid_expression(expr, adaptive):
         """
         equality = False
         if isinstance(expr, Equality):
-            expr = expr.lhs - expr.rhs
+            expr = expr.lhs_rhs()
path: sympy/plotting/series.py
text: @author:
```suggestion
            expr = expr.lhs - expr.rhs
```

diff_hunk: @@ -190,7 +190,7 @@ def _parallel_dict_from_expr_if_gens(exprs, opt):
         poly = {}
 
         if expr.is_Equality:
-            expr = expr.lhs - expr.rhs
+            expr = expr.lhs_rhs()
path: sympy/polys/polyutils.py
text: @author:
```suggestion
            expr = expr.lhs - expr.rhs
```

diff_hunk: @@ -249,7 +249,7 @@ def _is_coeff(factor):
         terms = []
 
         if expr.is_Equality:
-            expr = expr.lhs - expr.rhs
+            expr = expr.lhs_rhs()
path: sympy/polys/polyutils.py
text: @author:
```suggestion
            expr = expr.lhs - expr.rhs
```

diff_hunk: @@ -173,7 +173,7 @@ def _desolve(eq, func=None, hint="default", ics=None, simplify=True, *, prep=Tru
     classify_pde(pde.py)
     """
     if isinstance(eq, Equality):
-        eq = eq.lhs - eq.rhs
+        eq = eq.lhs_rhs()
path: sympy/solvers/deutils.py
text: @author:
```suggestion
        eq = eq.lhs - eq.rhs
```

diff_hunk: @@ -1335,7 +1335,7 @@ def diophantine(eq, param=symbols("t", integer=True), syms=None,
     eq = _sympify(eq)
 
     if isinstance(eq, Eq):
-        eq = eq.lhs - eq.rhs
+        eq = eq.lhs_rhs()
path: sympy/solvers/diophantine/diophantine.py
text: @author:
```suggestion
        eq = eq.lhs - eq.rhs
```

diff_hunk: @@ -233,7 +233,7 @@ def reduce_rational_inequalities(exprs, gen, relational=True):
                 expr, rel = expr
             else:
                 if expr.is_Relational:
-                    expr, rel = expr.lhs - expr.rhs, expr.rel_op
+                    expr, rel = expr.lhs_rhs(), expr.rel_op
path: sympy/solvers/inequalities.py
text: @author:
```suggestion
                    expr, rel = expr.lhs - expr.rhs, expr.rel_op
```

diff_hunk: @@ -479,7 +479,7 @@ def solve_univariate_inequality(expr, gen, relational=True, domain=S.Reals, cont
         rv = S.EmptySet
 
     else:
-        e = expr.lhs - expr.rhs
+        e = expr.lhs_rhs()
path: sympy/solvers/inequalities.py
text: @author:
```suggestion
        e = expr.lhs - expr.rhs
```

diff_hunk: @@ -805,7 +805,7 @@ def classify(ie, s, i):
 
     rv = None
     oo = S.Infinity
-    expr = ie.lhs - ie.rhs
+    expr = ie.lhs_rhs()
path: sympy/solvers/inequalities.py
text: @author:
```suggestion
    expr = ie.lhs - ie.rhs
```

diff_hunk: @@ -230,7 +230,7 @@ def infinitesimals(eq, func=None, order=None, hint='default', match=None):
     """
 
     if isinstance(eq, Equality):
-        eq = eq.lhs - eq.rhs
+        eq = eq.lhs_rhs()
path: sympy/solvers/ode/lie_group.py
text: @author:
```suggestion
        eq = eq.lhs - eq.rhs
```

diff_hunk: @@ -943,7 +943,7 @@ class in it.  Note that a hint may do this anyway if
         "work with functions of one variable, not %s" % func)
 
     if isinstance(eq, Equality):
-        eq = eq.lhs - eq.rhs
+        eq = eq.lhs_rhs()
path: sympy/solvers/ode/ode.py
text: @author:
```suggestion
        eq = eq.lhs - eq.rhs
```

diff_hunk: @@ -1213,7 +1213,7 @@ def _sympify(eq):
     eq, funcs = (_sympify(w) for w in [eq, funcs])
     for i, fi in enumerate(eq):
         if isinstance(fi, Equality):
-            eq[i] = fi.lhs - fi.rhs
+            eq[i] = fi.lhs_rhs()
path: sympy/solvers/ode/ode.py
text: @author:
```suggestion
            eq[i] = fi.lhs - fi.rhs
```

diff_hunk: @@ -2815,7 +2815,7 @@ def checkinfsol(eq, infinitesimals, func=None, order=None):
 
     """
     if isinstance(eq, Equality):
-        eq = eq.lhs - eq.rhs
+        eq = eq.lhs_rhs()
path: sympy/solvers/ode/ode.py
text: @author:
```suggestion
        eq = eq.lhs - eq.rhs
```

diff_hunk: @@ -295,7 +295,7 @@ def match_riccati(eq, f, x):
     """
     # Group terms based on f(x)
     if isinstance(eq, Eq):
-        eq = eq.lhs - eq.rhs
+        eq = eq.lhs_rhs()
path: sympy/solvers/ode/riccati.py
text: @author:
```suggestion
        eq = eq.lhs - eq.rhs
```

diff_hunk: @@ -177,7 +177,7 @@ def checkodesol(ode, sol, func=None, order='auto', solve_for_func=True):
         if testnum == 0:
             # First pass, try substituting a solved solution directly into the
             # ODE. This has the highest chance of succeeding.
-            ode_diff = ode.lhs - ode.rhs
+            ode_diff = ode.lhs_rhs()
path: sympy/solvers/ode/subscheck.py
text: @author:
```suggestion
            ode_diff = ode.lhs - ode.rhs
```

diff_hunk: @@ -221,7 +221,7 @@ def checkodesol(ode, sol, func=None, order='auto', solve_for_func=True):
                 diffsols = {0: sol.lhs}
             else:
                 diffsols = {}
-            sol = sol.lhs - sol.rhs
+            sol = sol.lhs_rhs()
path: sympy/solvers/ode/subscheck.py
text: @author:
```suggestion
            sol = sol.lhs - sol.rhs
```

diff_hunk: @@ -345,7 +345,7 @@ def _sympify(eq):
     eqs = _sympify(eqs)
     for i in range(len(eqs)):
         if isinstance(eqs[i], Equality):
-            eqs[i] = eqs[i].lhs - eqs[i].rhs
+            eqs[i] = eqs[i].lhs_rhs()
path: sympy/solvers/ode/subscheck.py
text: @author:
```suggestion
            eqs[i] = eqs[i].lhs - eqs[i].rhs
```

diff_hunk: @@ -717,7 +717,7 @@ def check_dummy_sol(eq, solse, dummy_sym):
     contains dummy symbols.
     """
     if isinstance(eq, Eq):
-        eq = eq.lhs - eq.rhs
+        eq = eq.lhs_rhs()
path: sympy/solvers/ode/tests/test_riccati.py
text: @author:
```suggestion
        eq = eq.lhs - eq.rhs
```

diff_hunk: @@ -468,7 +468,7 @@ def checkpdesol(pde, sol, func=None, solve_for_func=True):
 
     # try direct substitution of the solution into the PDE and simplify
     if sol.lhs == func:
-        pde = pde.lhs - pde.rhs
+        pde = pde.lhs_rhs()
path: sympy/solvers/pde.py
text: @author:
```suggestion
        pde = pde.lhs - pde.rhs
```

diff_hunk: @@ -728,7 +728,7 @@ def rsolve(f, y, init=None):
 
     """
     if isinstance(f, Equality):
-        f = f.lhs - f.rhs
+        f = f.lhs_rhs()
path: sympy/solvers/recurr.py
text: @author:
```suggestion
        f = f.lhs - f.rhs
```

diff_hunk: @@ -734,7 +734,7 @@ def _lp_matrices(objective, constraints):
     # change Eq(x, y) to x - y <= 0 and y - x <= 0
     for i in range(len(np)):
         if isinstance(np[i], Eq):
-            np[i] = np[i].lhs - np[i].rhs <= 0
+            np[i] = np[i].lhs_rhs() <= 0
path: sympy/solvers/simplex.py
text: @author:
```suggestion
            np[i] = np[i].lhs - np[i].rhs <= 0
```

diff_hunk: @@ -274,7 +274,7 @@ def checksol(f, symbol, sol=None, **flags):
             if not f.is_Boolean:
                 return
         else:
-            f = f.rewrite(Add, evaluate=False, deep=False)
+            f = f.lhs_rhs(evaluate=False, deep=False)
path: sympy/solvers/solvers.py
text: @author:
```suggestion
            f = Add(f.lhs, -f.rhs, evaluate=False)
```

diff_hunk: @@ -928,7 +928,7 @@ def _sympified_list(w):
     for i, fi in enumerate(f):
         if isinstance(fi, (Eq, Ne)):
             if 'ImmutableDenseMatrix' in [type(a).__name__ for a in fi.args]:
-                fi = fi.lhs - fi.rhs
+                fi = fi.lhs_rhs()
path: sympy/solvers/solvers.py
text: @author:
```suggestion
                fi = fi.lhs - fi.rhs
```

diff_hunk: @@ -948,7 +948,7 @@ def _sympified_list(w):
                             is True or False.
                         '''))
                 else:
-                    fi = fi.rewrite(Add, evaluate=False, deep=False)
+                    fi = fi.lhs_rhs(evaluate=False, deep=False)
path: sympy/solvers/solvers.py
text: @author:
```suggestion
                    fi = Add(fi.lhs, -fi.rhs, evaluate=False)
```

diff_hunk: @@ -2420,7 +2420,7 @@ def solve_undetermined_coeffs(equ, coeffs, *syms, **flags):
         raise ValueError('must provide symbols for coeffs')
 
     if isinstance(equ, Eq):
-        eq = equ.lhs - equ.rhs
+        eq = equ.lhs_rhs()
path: sympy/solvers/solvers.py
text: @author:
```suggestion
        eq = equ.lhs - equ.rhs
```

diff_hunk: @@ -3053,14 +3053,14 @@ def nsolve(*args, dict=False, **kwargs):
         f = list(f)
         for i, fi in enumerate(f):
             if isinstance(fi, Eq):
-                f[i] = fi.lhs - fi.rhs
+                f[i] = fi.lhs_rhs()
path: sympy/solvers/solvers.py
text: @author:
```suggestion
                f[i] = fi.lhs - fi.rhs
```

diff_hunk: @@ -1335,7 +1335,7 @@ def diophantine(eq, param=symbols("t", integer=True), syms=None,
     eq = _sympify(eq)
 
     if isinstance(eq, Eq):
-        eq = eq.lhs - eq.rhs
+        eq = eq.lhs_rhs()
path: sympy/solvers/diophantine/diophantine.py
text: @user1:
I don't think that this is an improvement.

diff_hunk: @@ -3428,7 +3428,7 @@ def _take(d):
     _take = flags.setdefault('_take', _take)
 
     if isinstance(eq, Eq):
-        eq = eq.lhs - eq.rhs  # XXX legacy Eq as Eqn support
+        eq = eq.lhs_rhs()  # XXX legacy Eq as Eqn support
path: sympy/solvers/solvers.py
text: @author:
```suggestion
        eq = eq.lhs - eq.rhs  # XXX legacy Eq as Eqn support
```

diff_hunk: @@ -3053,14 +3053,14 @@ def nsolve(*args, dict=False, **kwargs):
         f = list(f)
         for i, fi in enumerate(f):
             if isinstance(fi, Eq):
-                f[i] = fi.lhs - fi.rhs
+                f[i] = fi.lhs_rhs()
         f = Matrix(f).T
     if iterable(x0):
         x0 = list(x0)
     if not isinstance(f, Matrix):
         # assume it's a SymPy expression
         if isinstance(f, Eq):
-            f = f.lhs - f.rhs
+            f = f.lhs_rhs()
path: sympy/solvers/solvers.py
text: @author:
```suggestion
            f = f.lhs - f.rhs
```

diff_hunk: @@ -3311,7 +3311,7 @@ def _solve_using_known_values(result, solver):
                             # list.
                             result.remove(res)
                     continue  # skip as it's independent of desired symbols
-                depen1, depen2 = (eq2.rewrite(Add)).as_independent(*unsolved_syms)
+                depen1, depen2 = (eq2.lhs_rhs()).as_independent(*unsolved_syms)
path: sympy/solvers/solveset.py
text: @author:
```suggestion
                depen1, depen2 = Add(eq2.lhs - eq2.rhs, evaluate=False).as_independent(*unsolved_syms)
```

diff_hunk: @@ -3528,7 +3528,7 @@ def _separate_poly_nonpoly(system, symbols):
         denominators.update(_simple_dens(eq, symbols))
         # Convert equality to expression
         if isinstance(eq, Equality):
-            eq = eq.rewrite(Add)
+            eq = eq.lhs_rhs()
path: sympy/solvers/solveset.py
text: @author:
```suggestion
            eq = eq.lhs - eq.rhs
```

diff_hunk: @@ -2640,7 +2640,7 @@ def test_issue_10169():
 
 def test_solve_undetermined_coeffs_issue_23927():
     A, B, r, phi = symbols('A, B, r, phi')
-    eq = Eq(A*sin(t) + B*cos(t), r*sin(t - phi)).rewrite(Add).expand(trig=True)
+    eq = Eq(A*sin(t) + B*cos(t), r*sin(t - phi)).lhs_rhs().expand(trig=True)
path: sympy/solvers/tests/test_solvers.py
text: @author:
```suggestion
    _ = Eq(A*sin(t) + B*cos(t), r*sin(t - phi))
    eq = (_.lhs - _.rhs).expand(trig=True)
```

diff_hunk: @@ -3223,7 +3222,7 @@ def test_issue_23318():
         Eq(x, 0.0015 * z),
         Eq(0.0015, 7845.32 * y / z),
     ]
-    eqs_expr = [eq.rewrite(Add) for eq in eqs_eq]
+    eqs_expr = [eq.lhs_rhs() for eq in eqs_eq]
path: sympy/solvers/tests/test_solveset.py
text: @author:
```suggestion
    eqs_expr = [eq.lhs - eq.rhs for eq in eqs_eq]
```

diff_hunk: @@ -418,7 +418,7 @@ def probability(self, condition, **kwargs):
         # by computing a density handled by density computation
         except NotImplementedError:
             from sympy.stats.rv import density
-            expr = condition.lhs - condition.rhs
+            expr = condition.lhs_rhs()
path: sympy/stats/crv.py
text: @author:
```suggestion
            expr = condition.lhs - condition.rhs
```

diff_hunk: @@ -1065,7 +1065,7 @@ def _solveset(f, symbol, domain, _check=False):
             solns = solver(expr, symbol, in_set)
             result += solns
     elif isinstance(f, Eq):
-        result = solver(Add(f.lhs, - f.rhs, evaluate=False), symbol, domain)
+        result = solver(f.lhs_rhs(evaluate=False), symbol, domain)
path: sympy/solvers/solveset.py
text: @author:
```suggestion
        result = solver(Add(f.lhs, -f.rhs, evaluate=False), symbol, domain)
```

diff_hunk: @@ -467,7 +467,7 @@ def probability(self, condition, **kwargs):
             return Mul(*[self.probability(arg) for arg in condition.args])
         elif isinstance(condition, Or): # they are independent
             return Add(*[self.probability(arg) for arg in condition.args])
-        expr = condition.lhs - condition.rhs
+        expr = condition.lhs_rhs()
path: sympy/stats/rv.py
text: @author:
```suggestion
        expr = condition.lhs - condition.rhs
```

diff_hunk: @@ -59,7 +59,7 @@ def __new__(cls, variables, equation):
             variables = Tuple(*variables)
 
         if isinstance(equation, Eq):
-            equation = equation.lhs - equation.rhs
+            equation = equation.lhs_rhs()
path: sympy/vector/implicitregion.py
text: @author:
```suggestion
            equation = equation.lhs - equation.rhs
```

diff_hunk: @@ -238,7 +238,7 @@ def probability(self, condition):
             prob = self.eval_prob(_domain)
         except NotImplementedError:
             from sympy.stats.rv import density
-            expr = condition.lhs - condition.rhs
+            expr = condition.lhs_rhs()
path: sympy/stats/drv.py
text: @author:
```suggestion
            expr = condition.lhs - condition.rhs
```

diff_hunk: @@ -1026,18 +1028,18 @@ def test_rel_args():
                 raises(TypeError, lambda: Relational(b, v, op))
 
 
-def test_Equality_rewrite_as_Add():
-    eq = Eq(x + y, y - x)
-    assert eq.rewrite(Add) == 2*x
-    assert eq.rewrite(Add, evaluate=None).args == (x, x, y, -y)
-    assert eq.rewrite(Add, evaluate=False).args == (x, y, x, -y)
+def test_Equality_lhs_rhs():
+
     for e in (True, False, None):
path: sympy/core/tests/test_relational.py
text: @author:
```suggestion
```

diff_hunk: @@ -565,7 +567,7 @@ def test_x_minus_y_not_same_as_x_lt_y():
     y = 3 + t
     ineq = Lt(x, y, evaluate=False)
     raises(TypeError, lambda: ineq.doit())
-    assert ineq.lhs - ineq.rhs < 0
+    assert ineq.lhs_rhs() < 0
path: sympy/core/tests/test_relational.py
text: @author:
```suggestion
    assert ineq.lhs - ineq.rhs < 0
```

diff_hunk: @@ -683,11 +684,10 @@ def _eval_simplify(self, **kwargs):
         free = self.free_symbols
         if len(free) == 1:
             try:
-                from .add import Add
                 from sympy.solvers.solveset import linear_coeffs
                 x = free.pop()
                 m, b = linear_coeffs(
-                    e.rewrite(Add, evaluate=False), x)
+                    e.lhs_rhs(evaluate=False), x)
path: sympy/core/relational.py
text: @author:
```suggestion
                    Add(e.lhs, -e.rhs, evaluate=False), x)
```

diff_hunk: @@ -3311,7 +3311,7 @@ def _solve_using_known_values(result, solver):
                             # list.
                             result.remove(res)
                     continue  # skip as it's independent of desired symbols
-                depen1, depen2 = (eq2.rewrite(Add)).as_independent(*unsolved_syms)
+                depen1, depen2 = Add(eq2.lhs - eq2.rhs, evaluate=False).as_independent(*unsolved_syms)
path: sympy/solvers/solveset.py
text: @author:
```suggestion
                depen1, depen2 = (eq2.lhs - eq2.rhs).as_independent(*unsolved_syms)
```

diff_hunk: @@ -594,7 +594,7 @@ class Equality(Relational):
 
     Since this object is already an expression, it does not respond to
     the method ``as_expr`` if one tries to create `x - y` from ``Eq(x, y)``.
-    This can be done with the ``rewrite(Add)`` method.
+    This can be done for ``eq = Eq(x, y)`` with `eq.lhs - eq.rhs`.
path: sympy/core/relational.py
text: @author:
```suggestion
    If ``eq = Eq(x, y)`` then write `eq.lhs - eq.rhs` to get ``x - y``.
```

diff_hunk: @@ -1026,18 +1028,17 @@ def test_rel_args():
                 raises(TypeError, lambda: Relational(b, v, op))
 
 
-def test_Equality_rewrite_as_Add():
-    eq = Eq(x + y, y - x)
-    assert eq.rewrite(Add) == 2*x
-    assert eq.rewrite(Add, evaluate=None).args == (x, x, y, -y)
-    assert eq.rewrite(Add, evaluate=False).args == (x, y, x, -y)
-    for e in (True, False, None):
-        assert Eq(x, 0, evaluate=e).rewrite(Add) == x
-        assert Eq(0, x, evaluate=e).rewrite(Add) == x
+def test_nothing_happens_to_Eq_condition_during_simplify():
+
path: sympy/core/tests/test_relational.py
text: @author:
```suggestion
```

diff_hunk: @@ -3311,7 +3311,7 @@ def _solve_using_known_values(result, solver):
                             # list.
                             result.remove(res)
                     continue  # skip as it's independent of desired symbols
-                depen1, depen2 = (eq2.rewrite(Add)).as_independent(*unsolved_syms)
+                depen1, depen2 = (eq2.lhs - eq2.rhs).as_independent(*unsolved_syms)
path: sympy/solvers/solveset.py
text: @author:
```suggestion
                # XXX this should have happened in the initially called routine
                # not in the private method
                _ = eq2.lhs - eq2.rhs if isinstance(eq2, Eq) else eq2
                depen1, depen2 = _.as_independent(*unsolved_syms)
```

diff_hunk: @@ -274,7 +274,7 @@ def checksol(f, symbol, sol=None, **flags):
             if not f.is_Boolean:
                 return
         else:
path: sympy/solvers/solvers.py
text: @author:
```suggestion
        elif isinstance(f, Eq):
```

diff_hunk: @@ -948,7 +948,7 @@ def _sympified_list(w):
                             is True or False.
                         '''))
                 else:
path: sympy/solvers/solvers.py
text: @author:
```suggestion
                elif isinstance(fi, Eq):
```

diff_hunk: @@ -419,7 +419,7 @@ def _eval_simplify(self, **kwargs):
         if r.is_Relational:
             if not isinstance(r.lhs, Expr) or not isinstance(r.rhs, Expr):
                 return r
-            dif = r.lhs - r.rhs
+            dif = r.lhs_rhs()
path: sympy/core/relational.py
text: @author:
```suggestion
            dif = r.lhs - r.rhs
```

diff_hunk: @@ -1026,18 +1028,32 @@ def test_rel_args():
                 raises(TypeError, lambda: Relational(b, v, op))
 
 
-def test_Equality_rewrite_as_Add():
+def test_Equality_lhs_rhs():
+    # XXX we can be strict and make every routine that want to treat
+    # Eq like Eqn for purpose of solving explicitly test for Eq
+    #     if isinstance(eq, Eq):
+    #        eq = eq.lhs_rhs()
+    # or allow Expr to pass
+    #     eq = eq.lhs_rhs() # valid for Eq or Expr
path: sympy/core/tests/test_relational.py
text: @user1:
I think it is better to be strict. Code that does not know whether `eq` is an `Eq` or an `Expr` is always fragile. In any situation where either `Eq` or `Expr` is allowed then a conversion should be done at the earliest opportunity either like:
```python
if isinstance(eq, Eq):
     eq = eq.lhs_rhs()
elif isinstance(eq, Expr):
     eq = eq
else:
    raise Error
# Now we know we have Expr
```
Or it should be like:
```python
if isinstance(eq, Eq):
    lhs, rhs = eq.lhs, eq.rhs
elif isinstance(eq, Expr):
    lhs, rhs = eq, S.Zero
else:
    raise Error
# Now we know we have a pair of Expr
```
Or
```python
if isinstance(eq, Eq):
    eq = eq
elif isinstance(eq, Expr):
    eq = Eq(eq, 0)
else:
    raise Error
# Now we know we have Eq
```
I don't think that burying this logic into a convenience method helps for clarity of code.

In general we should minimise having situations where either `Eq` or `Expr` can be used or at least try to localise very carefully where the conversions are done if needed.

diff_hunk: @@ -68,6 +68,9 @@ class Expr(Basic, EvalfMixin):
 
     is_scalar = True  # self derivative is 1
 
+    def lhs_rhs(self, evaluate=True):
+        return self  # only for Eq
+
path: sympy/core/expr.py
text: @author:
```suggestion
```