Repository: sympy/sympy

Issue:
Incorrect relational simplification with squares 
Current master:
```
from sympy import *
x = Symbol('x', real=True)
y = Symbol('y', real=True)
expr = And(Eq(x**2, 4), Ge(y, x))
expr.simplify()
```
returns
`(y >= 4) & Eq(x**2, 4)`
so the square term in the equality is not accounted for (the code seems to just handle linear relations, but doesn't detect that this is not linear). The code is here:
https://github.com/sympy/sympy/blob/ec0fe8c5f3e59840e8aa5d3d6a7c976e40f76b64/sympy/logic/boolalg.py#L760-L761

I think it is related to the recently(?) introduced exception `NonlinearError` which is thrown, but not caught here. A bit strange that it is not passed on though.

Discussion:

Patch:
diff --git a/sympy/logic/boolalg.py b/sympy/logic/boolalg.py
index fde747b538..a5ae35cb90 100644
--- a/sympy/logic/boolalg.py
+++ b/sympy/logic/boolalg.py
@@ -752,10 +752,11 @@ def _eval_simplify(self, **kwargs):
                 (i.free_symbols, i) for i in eqs]),
                 lambda x: len(x[0]))
             eqs = []
+            nonlineqs = []
             while 1 in sifted:
                 for free, e in sifted.pop(1):
                     x = free.pop()
-                    if e.lhs != x or x in e.rhs.free_symbols:
+                    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)
@@ -768,21 +769,25 @@ def _eval_simplify(self, **kwargs):
                         except ValueError:
                             pass
                     if x in reps:
-                        eqs.append(e.func(e.rhs, reps[x]))
-                    else:
+                        eqs.append(e.subs(x, reps[x]))
+                    elif e.lhs == x and x not in e.rhs.free_symbols:
                         reps[x] = e.rhs
                         eqs.append(e)
+                    else:
+                        # x is not yet identified, but may be later
+                        nonlineqs.append(e)
                 resifted = defaultdict(list)
                 for k in sifted:
                     for f, e in sifted[k]:
-                        e = e.subs(reps)
+                        e = e.xreplace(reps)
                         f = e.free_symbols
                         resifted[len(f)].append((f, e))
                 sifted = resifted
         for k in sifted:
             eqs.extend([e for f, e in sifted[k]])
+        nonlineqs = [ei.subs(reps) for ei in nonlineqs]
         other = [ei.subs(reps) for ei in other]
-        rv = rv.func(*([i.canonical for i in (eqs + other)] + nonRel))
+        rv = rv.func(*([i.canonical for i in (eqs + nonlineqs + other)] + nonRel))
         patterns = simplify_patterns_and()
         return self._apply_patternbased_simplification(rv, patterns,
                                                        measure, False)
diff --git a/sympy/utilities/iterables.py b/sympy/utilities/iterables.py
index 77ccd453bd..83fc2f48d2 100644
--- a/sympy/utilities/iterables.py
+++ b/sympy/utilities/iterables.py
@@ -1419,7 +1419,7 @@ def multiset_permutations(m, size=None, g=None):
     do = [gi for gi in g if gi[1] > 0]
     SUM = sum([gi[1] for gi in do])
     if not do or size is not None and (size > SUM or size < 1):
-        if not do and size is None or size == 0:
+        if size < 1:
             yield []
         return
     elif size == 1:


Test patch:
diff --git a/sympy/simplify/tests/test_simplify.py b/sympy/simplify/tests/test_simplify.py
index 1b3b47386b..02bea1389c 100644
--- a/sympy/simplify/tests/test_simplify.py
+++ b/sympy/simplify/tests/test_simplify.py
@@ -7,7 +7,7 @@
     Matrix, MatrixSymbol, Mul, nsimplify, oo, pi, Piecewise, Poly, posify, rad,
     Rational, S, separatevars, signsimp, simplify, sign, sin,
     sinc, sinh, solve, sqrt, Sum, Symbol, symbols, sympify, tan,
-    zoo)
+    zoo, And, Le)
 from sympy.core.mul import _keep_coeff
 from sympy.core.expr import unchanged
 from sympy.simplify.simplify import nthroot, inversecombine
@@ -864,6 +864,37 @@ def test_issue_17137():
     assert simplify(cos(x)**(2 + 3*I)) == cos(x)**(2 + 3*I)
 
 
+def test_issue_21869():
+    x = Symbol('x', real=True)
+    y = Symbol('y', real=True)
+    expr = And(Eq(x**2, 4), Le(x, y))
+    assert expr.simplify() == expr
+
+    expr = And(Eq(x**2, 4), Eq(x, 2))
+    assert expr.simplify() == Eq(x, 2)
+
+    expr = And(Eq(x**3, x**2), Eq(x, 1))
+    assert expr.simplify() == Eq(x, 1)
+
+    expr = And(Eq(sin(x), x**2), Eq(x, 0))
+    assert expr.simplify() == Eq(x, 0)
+
+    expr = And(Eq(x**3, x**2), Eq(x, 2))
+    assert expr.simplify() == S.false
+
+    expr = And(Eq(y, x**2), Eq(x, 1))
+    assert expr.simplify() == And(Eq(y,1), Eq(x, 1))
+
+    expr = And(Eq(y**2, 1), Eq(y, x**2), Eq(x, 1))
+    assert expr.simplify() == And(Eq(y,1), Eq(x, 1))
+
+    expr = And(Eq(y**2, 4), Eq(y, 2*x**2), Eq(x, 1))
+    assert expr.simplify() == And(Eq(y,2), Eq(x, 1))
+
+    expr = And(Eq(y**2, 4), Eq(y, x**2), Eq(x, 1))
+    assert expr.simplify() == S.false
+
+
 def test_issue_7971():
     z = Integral(x, (x, 1, 1))
     assert z != 0
diff --git a/sympy/utilities/tests/test_iterables.py b/sympy/utilities/tests/test_iterables.py
index 77dadbf429..221b03f618 100644
--- a/sympy/utilities/tests/test_iterables.py
+++ b/sympy/utilities/tests/test_iterables.py
@@ -423,12 +423,6 @@ def test_multiset_permutations():
         [0, 1], [0, 2], [1, 0], [1, 2], [2, 0], [2, 1]]
     assert len(list(multiset_permutations('a', 0))) == 1
     assert len(list(multiset_permutations('a', 3))) == 0
-    for nul in ([], {}, ''):
-        assert list(multiset_permutations(nul)) == [[]]
-    assert list(multiset_permutations(nul, 0)) == [[]]
-    # impossible requests give no result
-    assert list(multiset_permutations(nul, 1)) == []
-    assert list(multiset_permutations(nul, -1)) == []
 
     def test():
         for i in range(1, 7):