Repository: sympy/sympy

Issue:
nested Max and Min give incorrect result
While working with different types of nested Max and Min expressions, I have encountered a problem in the following expression:
```
from sympy import sympify 
sympify( "Max( Min(a, b), Min(a, b, c))" )
```

output: `Min(a, b, c)`

This simplification is not correct for all possible values of a, b, and c.
I have not found any problems with other similar expressions, and I cannot find the specific reason why this one evaluates incorrectly.



Discussion:

Patch:
diff --git a/sympy/functions/elementary/miscellaneous.py b/sympy/functions/elementary/miscellaneous.py
index 85b606e754..b458bd2af7 100644
--- a/sympy/functions/elementary/miscellaneous.py
+++ b/sympy/functions/elementary/miscellaneous.py
@@ -1,4 +1,5 @@
 from sympy.core import Function, S, sympify
+from sympy.utilities.iterables import sift
 from sympy.core.add import Add
 from sympy.core.containers import Tuple
 from sympy.core.compatibility import ordered
@@ -432,7 +433,6 @@ def _collapse_arguments(cls, args, **assumptions):
 
         >>> Min(a, Max(b, Min(c, d, Max(a, e))))
         Min(a, Max(b, Min(a, c, d)))
-
         """
         from sympy.utilities.iterables import ordered
         from sympy.simplify.simplify import walk
@@ -518,28 +518,33 @@ def do(ai, a):
         # easy case where all functions contain something in common;
         # trying to find some optimal subset of args to modify takes
         # too long
+
+        def factor_minmax(args):
+            is_other = lambda arg: isinstance(arg, other)
+            other_args, remaining_args = sift(args, is_other, binary=True)
+            if not other_args:
+                return args
+
+            # Min(Max(x, y, z), Max(x, y, u, v)) -> {x,y}, ({z}, {u,v})
+            arg_sets = [set(arg.args) for arg in other_args]
+            common = set.intersection(*arg_sets)
+            if not common:
+                return args
+
+            new_other_args = list(common)
+            arg_sets_diff = [arg_set - common for arg_set in arg_sets]
+
+            # If any set is empty after removing common then all can be
+            # discarded e.g. Min(Max(a, b, c), Max(a, b)) -> Max(a, b)
+            if all(arg_sets_diff):
+                other_args_diff = [other(*s, evaluate=False) for s in arg_sets_diff]
+                new_other_args.append(cls(*other_args_diff, evaluate=False))
+
+            other_args_factored = other(*new_other_args, evaluate=False)
+            return remaining_args + [other_args_factored]
+
         if len(args) > 1:
-            common = None
-            remove = []
-            sets = []
-            for i in range(len(args)):
-                a = args[i]
-                if not isinstance(a, other):
-                    continue
-                s = set(a.args)
-                common = s if common is None else (common & s)
-                if not common:
-                    break
-                sets.append(s)
-                remove.append(i)
-            if common:
-                sets = filter(None, [s - common for s in sets])
-                sets = [other(*s, evaluate=False) for s in sets]
-                for i in reversed(remove):
-                    args.pop(i)
-                oargs = [cls(*sets)] if sets else []
-                oargs.extend(common)
-                args.append(other(*oargs, evaluate=False))
+            args = factor_minmax(args)
 
         return args
 
diff --git a/sympy/polys/constructor.py b/sympy/polys/constructor.py
index 3ae1babce9..adcfd13dc3 100644
--- a/sympy/polys/constructor.py
+++ b/sympy/polys/constructor.py
@@ -48,7 +48,7 @@ def _construct_simple(coeffs, opt):
                         float_numbers.append(x)
                     if y.is_Float:
                         float_numbers.append(y)
-            elif is_algebraic(coeff):
+            if is_algebraic(coeff):
                 if floats:
                     # there are both algebraics and reals -> EX
                     return False
diff --git a/sympy/polys/matrices/ddm.py b/sympy/polys/matrices/ddm.py
index d1770b21b8..7b98cc16f4 100644
--- a/sympy/polys/matrices/ddm.py
+++ b/sympy/polys/matrices/ddm.py
@@ -284,9 +284,7 @@ def applyfunc(self, func, domain):
     def rref(a):
         """Reduced-row echelon form of a and list of pivots"""
         b = a.copy()
-        K = a.domain
-        partial_pivot = K.is_RealField or K.is_ComplexField
-        pivots = ddm_irref(b, _partial_pivot=partial_pivot)
+        pivots = ddm_irref(b)
         return b, pivots
 
     def nullspace(a):
diff --git a/sympy/polys/matrices/dense.py b/sympy/polys/matrices/dense.py
index dfaff42a29..bf6908805f 100644
--- a/sympy/polys/matrices/dense.py
+++ b/sympy/polys/matrices/dense.py
@@ -85,7 +85,7 @@ def ddm_imatmul(a, b, c):
             ai[j] = sum(map(mul, bi, cTj), ai[j])
 
 
-def ddm_irref(a, _partial_pivot=False):
+def ddm_irref(a):
     """a  <--  rref(a)"""
     # a is (m x n)
     m = len(a)
@@ -97,15 +97,6 @@ def ddm_irref(a, _partial_pivot=False):
     pivots = []
 
     for j in range(n):
-        # Proper pivoting should be used for all domains for performance
-        # reasons but it is only strictly needed for RR and CC (and possibly
-        # other domains like RR(x)). This path is used by DDM.rref() if the
-        # domain is RR or CC. It uses partial (row) pivoting based on the
-        # absolute value of the pivot candidates.
-        if _partial_pivot:
-            ip = max(range(i, m), key=lambda ip: abs(a[ip][j]))
-            a[i], a[ip] = a[ip], a[i]
-
         # pivot
         aij = a[i][j]
 
diff --git a/sympy/polys/matrices/linsolve.py b/sympy/polys/matrices/linsolve.py
index 337a9939ae..62e6d90aae 100644
--- a/sympy/polys/matrices/linsolve.py
+++ b/sympy/polys/matrices/linsolve.py
@@ -73,12 +73,6 @@ def _linsolve(eqs, syms):
     Aaug = sympy_dict_to_dm(eqsdict, rhs, syms)
     K = Aaug.domain
 
-    # sdm_irref has issues with float matrices. This uses the ddm_rref()
-    # function. When sdm_rref() can handle float matrices reasonably this
-    # should be removed...
-    if K.is_RealField or K.is_ComplexField:
-        Aaug = Aaug.to_ddm().rref()[0].to_sdm()
-
     # Compute reduced-row echelon form (RREF)
     Arref, pivots, nzcols = sdm_irref(Aaug)
 
diff --git a/sympy/polys/matrices/sdm.py b/sympy/polys/matrices/sdm.py
index 7c4ad43660..cfa624185a 100644
--- a/sympy/polys/matrices/sdm.py
+++ b/sympy/polys/matrices/sdm.py
@@ -904,8 +904,6 @@ def sdm_irref(A):
             Ajnz = set(Aj)
             for k in Ajnz - Ainz:
                 Ai[k] = - Aij * Aj[k]
-            Ai.pop(j)
-            Ainz.remove(j)
             for k in Ajnz & Ainz:
                 Aik = Ai[k] - Aij * Aj[k]
                 if Aik:
@@ -940,8 +938,6 @@ def sdm_irref(A):
             for l in Ainz - Aknz:
                 Ak[l] = - Akj * Ai[l]
                 nonzero_columns[l].add(k)
-            Ak.pop(j)
-            Aknz.remove(j)
             for l in Ainz & Aknz:
                 Akl = Ak[l] - Akj * Ai[l]
                 if Akl:


Test patch:
diff --git a/sympy/functions/elementary/tests/test_miscellaneous.py b/sympy/functions/elementary/tests/test_miscellaneous.py
index 0e1a5b1d08..bc6f2e2762 100644
--- a/sympy/functions/elementary/tests/test_miscellaneous.py
+++ b/sympy/functions/elementary/tests/test_miscellaneous.py
@@ -410,6 +410,10 @@ def test_issue_12638():
     assert Min(a, b, Max(a, b, c)) == Min(a, b)
     assert Min(a, b, Max(a, c)) == Min(a, b)
 
+def test_issue_21399():
+    from sympy.abc import a, b, c
+    assert Max(Min(a, b), Min(a, b, c)) == Min(a, b)
+
 
 def test_instantiation_evaluation():
     from sympy.abc import v, w, x, y, z
diff --git a/sympy/functions/elementary/tests/test_piecewise.py b/sympy/functions/elementary/tests/test_piecewise.py
index 1833939f5e..cb977261be 100644
--- a/sympy/functions/elementary/tests/test_piecewise.py
+++ b/sympy/functions/elementary/tests/test_piecewise.py
@@ -233,38 +233,40 @@ def test_piecewise_integrate1ca():
     assert g.integrate((x, 1, -2)) == g1y.subs(y, -2)
     assert g.integrate((x, 0, 1)) == gy1.subs(y, 0)
     assert g.integrate((x, 1, 0)) == g1y.subs(y, 0)
-    # XXX Make test pass without simplify
-    assert g.integrate((x, 2, 1)) == gy1.subs(y, 2).simplify()
-    assert g.integrate((x, 1, 2)) == g1y.subs(y, 2).simplify()
-
+    # since xreplace is followed by doit with deep=False
+    # to Min/Max don't get resolved so doit is needed
+    assert g.integrate((x, 2, 1)) == gy1.subs(y, 2).doit()
+    assert g.integrate((x, 1, 2)) == g1y.subs(y, 2).doit()
     assert piecewise_fold(gy1.rewrite(Piecewise)) == \
         Piecewise(
-            (1, y <= -1),
-            (-y**2/2 - y + S.Half, y <= 0),
-            (y**2/2 - y + S.Half, y < 1),
+            (-y**2/2 - y + S.Half, (y >= -1) & (y <= 0)),
+            (y**2/2 - y + S.Half, (y >= -1) & (y < 1)),
+            (1, y <= 0),
+            (y**2 + 1, y < 1),
             (0, True))
     assert piecewise_fold(g1y.rewrite(Piecewise)) == \
         Piecewise(
-            (-1, y <= -1),
-            (y**2/2 + y - S.Half, y <= 0),
-            (-y**2/2 + y - S.Half, y < 1),
+            (y**2/2 + y - S.Half, (y >= -1) & (y <= 0)),
+            (-y**2/2 + y - S.Half, (y >= -1) & (y < 1)),
+            (-1, y <= 0),
+            (-y**2 - 1, y < 1),
             (0, True))
 
-    # g1y and gy1 should simplify if the condition that y < 1
+    # g1y and gy1 simplify if the condition that y < 1
     # is applied, e.g. Min(1, Max(-1, y)) --> Max(-1, y)
-    # XXX Make test pass without simplify
-    assert gy1.simplify() == Piecewise(
+    assert gy1.doit() == Piecewise(
         (
             -Min(1, Max(-1, y))**2/2 - Min(1, Max(-1, y)) +
             Min(1, Max(0, y))**2 + S.Half, y < 1),
         (0, True)
         )
-    assert g1y.simplify() == Piecewise(
+    assert g1y.doit() == Piecewise(
         (
             Min(1, Max(-1, y))**2/2 + Min(1, Max(-1, y)) -
             Min(1, Max(0, y))**2 - S.Half, y < 1),
         (0, True))
 
+
 @slow
 def test_piecewise_integrate1cb():
     y = symbols('y', real=True)
diff --git a/sympy/polys/matrices/tests/test_linsolve.py b/sympy/polys/matrices/tests/test_linsolve.py
index a77beebdba..eda4cdbdf3 100644
--- a/sympy/polys/matrices/tests/test_linsolve.py
+++ b/sympy/polys/matrices/tests/test_linsolve.py
@@ -7,7 +7,7 @@
 from sympy.testing.pytest import raises
 
 from sympy import S, Eq, I
-from sympy.abc import x, y, z
+from sympy.abc import x, y
 
 from sympy.polys.matrices.linsolve import _linsolve
 from sympy.polys.solvers import PolyNonlinearError
@@ -23,83 +23,6 @@ def test__linsolve():
     raises(PolyNonlinearError, lambda: _linsolve([x*(1 + x)], [x]))
 
 
-def test__linsolve_float():
-
-    # This should give the exact answer:
-    eqs = [
-        y - x,
-        y - 0.0216 * x
-    ]
-    sol = {x:0.0, y:0.0}
-    assert _linsolve(eqs, (x, y)) == sol
-
-    # Other cases should be close to eps
-
-    def all_close(sol1, sol2, eps=1e-15):
-        close = lambda a, b: abs(a - b) < eps
-        assert sol1.keys() == sol2.keys()
-        return all(close(sol1[s], sol2[s]) for s in sol1)
-
-    eqs = [
-        0.8*x +         0.8*z + 0.2,
-        0.9*x + 0.7*y + 0.2*z + 0.9,
-        0.7*x + 0.2*y + 0.2*z + 0.5
-    ]
-    sol_exact = {x:-29/42, y:-11/21, z:37/84}
-    sol_linsolve = _linsolve(eqs, [x,y,z])
-    assert all_close(sol_exact, sol_linsolve)
-
-    eqs = [
-        0.9*x + 0.3*y + 0.4*z + 0.6,
-        0.6*x + 0.9*y + 0.1*z + 0.7,
-        0.4*x + 0.6*y + 0.9*z + 0.5
-    ]
-    sol_exact = {x:-88/175, y:-46/105, z:-1/25}
-    sol_linsolve = _linsolve(eqs, [x,y,z])
-    assert all_close(sol_exact, sol_linsolve)
-
-    eqs = [
-        0.4*x + 0.3*y + 0.6*z + 0.7,
-        0.4*x + 0.3*y + 0.9*z + 0.9,
-        0.7*x + 0.9*y,
-    ]
-    sol_exact = {x:-9/5, y:7/5, z:-2/3}
-    sol_linsolve = _linsolve(eqs, [x,y,z])
-    assert all_close(sol_exact, sol_linsolve)
-
-    eqs = [
-        x*(0.7 + 0.6*I) + y*(0.4 + 0.7*I) + z*(0.9 + 0.1*I) + 0.5,
-        0.2*I*x + 0.2*I*y + z*(0.9 + 0.2*I) + 0.1,
-        x*(0.9 + 0.7*I) + y*(0.9 + 0.7*I) + z*(0.9 + 0.4*I) + 0.4,
-    ]
-    sol_exact = {
-        x:-6157/7995 - 411/5330*I,
-        y:8519/15990 + 1784/7995*I,
-        z:-34/533 + 107/1599*I,
-    }
-    sol_linsolve = _linsolve(eqs, [x,y,z])
-    assert all_close(sol_exact, sol_linsolve)
-
-    # XXX: This system for x and y over RR(z) is problematic.
-    #
-    # eqs = [
-    #     x*(0.2*z + 0.9) + y*(0.5*z + 0.8) + 0.6,
-    #     0.1*x*z + y*(0.1*z + 0.6) + 0.9,
-    # ]
-    #
-    # linsolve(eqs, [x, y])
-    # The solution for x comes out as
-    #
-    #       -3.9e-5*z**2 - 3.6e-5*z - 8.67361737988404e-20
-    #  x =  ----------------------------------------------
-    #           3.0e-6*z**3 - 1.3e-5*z**2 - 5.4e-5*z
-    #
-    # The 8e-20 in the numerator should be zero which would allow z to cancel
-    # from top and bottom. It should be possible to avoid this somehow because
-    # the inverse of the matrix only has a quadratic factor (the determinant)
-    # in the denominator.
-
-
 def test__linsolve_deprecated():
     assert _linsolve([Eq(x**2, x**2+y)], [x, y]) == {x:x, y:S.Zero}
     assert _linsolve([(x+y)**2-x**2], [x]) == {x:-y/2}
diff --git a/sympy/polys/tests/test_constructor.py b/sympy/polys/tests/test_constructor.py
index 0feb681d26..0f2c747202 100644
--- a/sympy/polys/tests/test_constructor.py
+++ b/sympy/polys/tests/test_constructor.py
@@ -27,9 +27,6 @@ def test_construct_domain():
     assert isinstance(result[0], ComplexField)
     assert result[1] == [CC(3.14), CC(1.0j), CC(0.5)]
 
-    assert construct_domain([1.0+I]) == (CC, [CC(1.0, 1.0)])
-    assert construct_domain([2.0+3.0*I]) == (CC, [CC(2.0, 3.0)])
-
     assert construct_domain([1, I]) == (ZZ_I, [ZZ_I(1, 0), ZZ_I(0, 1)])
     assert construct_domain([1, I/2]) == (QQ_I, [QQ_I(1, 0), QQ_I(0, S.Half)])
 
diff --git a/sympy/polys/tests/test_polytools.py b/sympy/polys/tests/test_polytools.py
index 762e88f720..77a153c5a2 100644
--- a/sympy/polys/tests/test_polytools.py
+++ b/sympy/polys/tests/test_polytools.py
@@ -51,7 +51,6 @@
 from sympy.polys.fields import field
 from sympy.polys.domains import FF, ZZ, QQ, ZZ_I, QQ_I, RR, EX
 from sympy.polys.domains.realfield import RealField
-from sympy.polys.domains.complexfield import ComplexField
 from sympy.polys.orderings import lex, grlex, grevlex
 
 from sympy import (
@@ -390,7 +389,6 @@ def test_Poly__new__():
              modulus=65537, symmetric=False)
 
     assert isinstance(Poly(x**2 + x + 1.0).get_domain(), RealField)
-    assert isinstance(Poly(x**2 + x + I + 1.0).get_domain(), ComplexField)
 
 
 def test_Poly__args():