Repository: sympy/sympy

Issue:
Simplification of Relational and BooleanFunction
Consider the following:

``` python
In [6]: x = Symbol('x', real=True)

In [7]: simplify(Or(x<1, x>= 1))
Out[7]: Or(x < 1, x >= 1)
```

I'd expect the result to be `True` or somehow else to designate "all `x`".

Even if I force the symbol to be a positive integer, this does not simplify:

``` python
In [3]: r = Symbol('r', integer=True, positive=True)

In [4]: Or(r<1, r>= 1)
Out[4]: Or(r < 1, r >= 1)

In [5]: simplify(Or(r<1, r>= 1))
Out[5]: Or(r < 1, r >= 1)
```

which is a bit odd. All of this causes over-complicated outputs sometimes. For example, trying to sum the PMF of a negative binomial distribution with `p=1/2`, I get

``` python
In [3]: k = Symbol('k', integer=True, nonnegative=True)

In [4]: r = Symbol('r', integer=True, positive=True)

In [5]: Sum(binomial(k+r-1,k)*(S.Half)**r*S.Half**k, (k, 0, oo)).doit()
Out[5]: Piecewise((1, Or(r < 1, r >= 1)), (Sum(2**(-k)*2**(-r)*binomial(k + r - 1, k), (k, 0, oo)), True))
```

simplify And(Le|Eq(x, 3), Ge(x, 3)) to Eq(x, 3)
```
x = Symbol("x")
```
With
```
expr = And(Ge(x, 3), Le(x, 3))
```
Or:
```
expr = And(Eq(x, 3), Le(x, 3))
```
I would expect that `expr.simplify()` yields `Eq(x, 3)` in both cases.
Maybe this is related to #18642 or #18389 but I'm not sure.


Discussion:
sympy-bot: :white_check_mark:

Hi, I am the [SymPy bot](https://github.com/sympy/sympy-bot) (v161). I'm here to help you write a release notes entry. Please read the [guide on how to write release notes](https://github.com/sympy/sympy/wiki/Writing-Release-Notes).



Your release notes are in good order.

Here is what the release notes will look like:
* logic
  * Pattern based logic simplification is actually used when doing `simplify` with relations. ([#21854](https://github.com/sympy/sympy/pull/21854) by [@oscargus](https://github.com/oscargus))

  * Improved logic simplification, dealing better with some symmetric cases. ([#21854](https://github.com/sympy/sympy/pull/21854) by [@oscargus](https://github.com/oscargus))

This will be added to https://github.com/sympy/sympy/wiki/Release-Notes-for-1.9.

<details><summary>Click here to see the pull request description that was parsed.</summary>

    <!-- 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. -->
    Fixes #18645
    Fixes #19114
    Fixes #19882
    Fixes #20870
    Added (XFAIL) tests for #18389 #18642

    #### Brief description of what is fixed or changed
    Fixed logic simplification so that the pattern based simplification is used even if there are no `Eq`-constrains involved (removed early exit)

    Improved set cover algorithm which only used to work in case one term fully covered another term. Now, a greedy heuristic is used in that case.

    #### Other comments


    #### 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.

    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 -->
    * logic
      * Pattern based logic simplification is actually used when doing `simplify` with relations.
      * Improved logic simplification, dealing better with some symmetric cases.
    <!-- END RELEASE NOTES -->


</details><p>

oscargus: @oscarbenjamin There is a test that fails using optional dependencies (currently XFAILed). I have no idea why and which dependency might affect it. Any ideas? Anything being used to simplify dnf-forms now?
oscargus: Ahh, OK! I guess I need to sort at some stages then to get consistent results.
oscargus: The results should now be consistent. Also improved the run time a bit for `simplify_logic` by not testing for cnf or dnf forms unless required.
oscargus: For now, I just check that the resulting expression is equal and minimal. Not sure how to deal with it to get results that are consistent. 
oscarbenjamin: > Is there are way to access the test server to see what is happening? If I recall correctly simply printing intermediate results doesn't work.

It's not hard to reproduce the setup outside of CI:
https://github.com/sympy/sympy/blob/12827227e1d71aef6967724df3f8dd019fc18a4e/.github/workflows/runtests.yml#L80-L99
oscargus: You need to try the penultimate commit for that... But you are probably correct in that it cannot produce consistent results (I think because all boolean functions use a set for the args, so you cannot even guarantee that the order you type it in is retained.)

github-actions[bot]: Benchmark results from GitHub Actions

Lower numbers are good, higher numbers are bad. A ratio less than 1
means a speed up and greater than 1 means a slowdown. Green lines
beginning with `+` are slowdowns (the PR is slower then master or
master is slower than the previous release). Red lines beginning
with `-` are speedups.

Significantly changed benchmark results (PR vs master)
```diff

```
Significantly changed benchmark results (master vs previous release)
```diff
       before           after         ratio
     [ed9a550f]       [12827227]
     <sympy-1.8^0>                 
-      1.12±0.01s          157±9ms     0.14  dsolve.TimeDsolve01.time_dsolve
+         120±2ms          184±7ms     1.54  logic.LogicSuite.time_dpll2
-        82.1±7μs         32.2±1μs     0.39  matrices.TimeDiagonalEigenvals.time_eigenvals
-        94.8±5μs         60.9±4μs     0.64  matrices.TimeMatrixGetItem.time_ImmutableDenseMatrix_getitem
-         102±8μs         64.4±3μs     0.63  matrices.TimeMatrixGetItem.time_MutableSparseMatrix_getitem
-        90.1±4μs         56.6±2μs     0.63  solve.TimeMatrixArithmetic.time_dense_add(10, 0)
+        97.2±6μs          157±4μs     1.62  solve.TimeMatrixArithmetic.time_dense_add(10, 5)
+      12.2±0.4μs       24.8±0.2μs     2.03  solve.TimeMatrixArithmetic.time_dense_add(3, 0)
+      16.3±0.4μs         37.4±1μs     2.29  solve.TimeMatrixArithmetic.time_dense_add(3, 5)
+        21.1±1μs         44.2±2μs     2.10  solve.TimeMatrixArithmetic.time_dense_add(4, 5)
+        41.3±2μs         69.1±2μs     1.67  solve.TimeMatrixArithmetic.time_dense_add(6, 5)
-     1.33±0.06ms         275±10μs     0.21  solve.TimeMatrixArithmetic.time_dense_multiply(10, 0)
-        57.1±2μs         27.7±2μs     0.48  solve.TimeMatrixArithmetic.time_dense_multiply(3, 0)
+        88.4±2μs          151±7μs     1.70  solve.TimeMatrixArithmetic.time_dense_multiply(3, 5)
-         115±6μs         34.2±2μs     0.30  solve.TimeMatrixArithmetic.time_dense_multiply(4, 0)
-        352±20μs       75.7±0.8μs     0.22  solve.TimeMatrixArithmetic.time_dense_multiply(6, 0)
-        96.8±4μs       45.0±0.9μs     0.47  solve.TimeMatrixOperations.time_det(3, 0)
-        176±10μs          113±3μs     0.64  solve.TimeMatrixOperations.time_det(3, 2)
-         185±9μs         98.0±3μs     0.53  solve.TimeMatrixOperations.time_det(3, 5)
-         112±3μs         46.7±2μs     0.42  solve.TimeMatrixOperations.time_det_bareiss(3, 0)
-         178±7μs          114±4μs     0.64  solve.TimeMatrixOperations.time_det_bareiss(3, 2)
-         184±6μs          108±9μs     0.58  solve.TimeMatrixOperations.time_det_bareiss(3, 5)
-         107±3μs         48.6±2μs     0.45  solve.TimeMatrixOperations.time_det_berkowitz(3, 0)
-         168±3μs          106±6μs     0.63  solve.TimeMatrixOperations.time_det_berkowitz(3, 5)
+        814±70μs      1.33±0.06ms     1.64  solve.TimeMatrixOperations.time_det_berkowitz(4, 0)
+     1.04±0.07ms      2.07±0.05ms     2.00  solve.TimeMatrixOperations.time_det_berkowitz(4, 2)
+     1.09±0.07ms      2.16±0.08ms     1.98  solve.TimeMatrixOperations.time_det_berkowitz(4, 5)
+        313±20μs         574±20μs     1.83  solve.TimeMatrixOperations.time_rank(3, 0)
-      5.44±0.3ms       2.94±0.1ms     0.54  solve.TimeRationalSystem.time_linsolve(10)
-     1.18±0.08ms          707±9μs     0.60  solve.TimeRationalSystem.time_linsolve(5)
+     3.19±0.06ms       5.36±0.1ms     1.68  solve.TimeSparseSystem.time_linear_eq_to_matrix(20)
+      6.37±0.5ms       10.5±0.6ms     1.64  solve.TimeSparseSystem.time_linear_eq_to_matrix(30)

```
Full benchmark results can be found as artifacts in GitHub Actions
(click on checks at the top of the PR).
oscargus: Is there anything here to be sorted out? I am quite convinced that the inconsistency is not caused by the current fixes.
smichr: > I still see nondeterministic results:

With this diff I get `B & C` always appearing first:

```diff
diff --git a/sympy/logic/boolalg.py b/sympy/logic/boolalg.py
index 5021354dd5..d59141528d 100644
--- a/sympy/logic/boolalg.py
+++ b/sympy/logic/boolalg.py
@@ -2865,7 +2865,7 @@ def simplify_logic(expr, form=None, deep=True, force=False):
     if not force and len(variables) > 8:
         return expr
     # group into constants and variable values
-    c, v = sift(variables, lambda x: x in (True, False), binary=True)
+    c, v = sift(ordered(variables), lambda x: x in (True, False), binary=True)
     variables = c + v
     truthtable = []
     # standardize constants to be 1 or 0 in keeping with truthtable
```
oscargus: Thanks @smichr ! I've added that one now.

Interesting that this has never manifested itself earlier. I didn't have to change any of the other tests (in `test_boolalg.py` and `test_simplify.py` at least, let's see how the rest goes...).
oscarbenjamin: > With this diff I get `B & C` always appearing first:

I also get deterministic results now with the PR.
oscargus: Nice! And the tests pass! Good to go?
oscarbenjamin: It looks good to me. My only question would be about performance since I know that some of these routines can be very slow in large cases.

Do we have any benchmarks for this sort of stuff?

Could you maybe suggest an example for how to time this in a potentially slow case?
smichr: > My only question would be about performance

Does this concern block this merge? With the number of issues being closed, it seems like this would be good to merge. And if it can be slow then we would still want the issues to be fixed, but in a better way.
oscarbenjamin: > Does this concern block this merge? With the number of issues being closed, it seems like this would be good to merge. And if it can be slow then we would still want the issues to be fixed, but in a better way.

I think this is a good fix but I'd like to at least consider the performance implications. If there are slowdowns it doesn't mean that it should not be merged but at least we can consider whether there is a better place to implement this.

I ask about this because the tests are all for relatively simple examples but I can see several quadruply nested loops which is the kind of thing that can be slow so it would be good to consider under what conditions those loops could have large numbers of iterations e.g.:
```python
    while anythingchanged:
        ...
        for n, rowi in enumerate(oldndterms):
                ...
                for n2, row2i in enumerate(oldndterms):
                        ...
                        if all(a >= b for (a, b) in zip(filteredrows[n2], row)):
```
That's a loop in a loop in a loop in a loop. Maybe all loops only ever have a couple of iterations but I don't know. Maybe this is the best algorithm even if it's `O(n^4)` since many logic algorithms are `O(2^n)` or worse.

Testing with master it's not hard to see that some of these things are slow:
```python
In [10]: xs = symbols('x:10')

In [11]: fact = And(*(x1 < x2 for x1 , x2 in zip(xs[:-1], xs[1:])))

In [12]: %time simplify(fact) == fact
CPU times: user 8.91 s, sys: 49.8 ms, total: 8.96 s
Wall time: 9.01 s
Out[12]: True
```
That's 9 seconds to simplify an `And` with 10 terms and the simplified result as still the same as the input.

Increasing the size of the input leads to a nonlinear increase in time:
```python
In [18]: xs = symbols('x:30')

In [19]: fact = And(*(x1 < x2 for x1 , x2 in zip(xs[:-1], xs[1:])))

In [20]: %time ok = simplify(fact)
CPU times: user 1min 10s, sys: 870 ms, total: 1min 10s
Wall time: 1min 15s
```
(For some reason with 30 terms the result does not come back identical...)

Profiling suggests that the code here isn't involved:
```python
In [21]: xs = symbols('x:10')

In [22]: fact = And(*(x1 < x2 for x1 , x2 in zip(xs[:-1], xs[1:])))

In [24]: %prun -s cumulative simplify(fact)

         11091774 function calls (10736229 primitive calls) in 14.515 seconds

   Ordered by: cumulative time

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
      2/1    0.000    0.000   14.515   14.515 {built-in method builtins.exec}
        1    0.000    0.000   14.515   14.515 <string>:1(<module>)
     91/1    0.004    0.000   14.515   14.515 simplify.py:411(simplify)
        1    0.000    0.000   14.515   14.515 boolalg.py:730(_eval_simplify)
        1    0.024    0.024   11.914   11.914 boolalg.py:543(_apply_patternbased_simplification)
     3050    0.063    0.000   10.945    0.004 operations.py:481(__new__)
     3050    0.095    0.000    9.841    0.003 boolalg.py:692(_new_args_filter)
    12319    0.197    0.000    7.970    0.001 relational.py:202(canonical)
24796/24688    0.119    0.000    5.309    0.000 expr.py:2452(could_extract_minus_sign)
74558/49708    0.863    0.000    4.741    0.000 expr.py:2182(extract_multiplicatively)
235829/233547    0.562    0.000    4.262    0.000 cache.py:69(wrapper)
172874/62092    0.947    0.000    4.256    0.000 compatibility.py:501(ordered)
185857/160614    0.930    0.000    3.371    0.000 expr.py:143(__eq__)
        1    0.000    0.000    2.504    2.504 boolalg.py:444(_eval_simplify)
       18    0.002    0.000    2.481    0.138 relational.py:303(_eval_simplify)
    54/18    0.002    0.000    2.421    0.135 expr.py:710(equals)
51223/51200    0.088    0.000    2.340    0.000 decorators.py:254(_func)
51189/51166    0.078    0.000    2.215    0.000 decorators.py:129(binary_op_wrapper)
       18    0.001    0.000    2.049    0.114 expr.py:538(is_constant)
50142/50139    0.149    0.000    2.005    0.000 expr.py:259(__truediv__)
    36966    0.548    0.000    1.947    0.000 compatibility.py:315(default_sort_key)
       18    0.001    0.000    1.895    0.105 expr.py:457(_random)
       18    0.000    0.000    1.879    0.104 expr.py:498(<listcomp>)
       36    0.001    0.000    1.879    0.052 randtest.py:12(random_complex_number)
   144/72    0.004    0.000    1.878    0.026 simplify.py:1439(nsimplify)
       72    0.003    0.000    1.864    0.026 simplify.py:1578(_real_to_rational)
       36    0.001    0.000    1.722    0.048 simplify.py:1533(nsimplify_real)
    59/36    0.015    0.000    1.689    0.047 identification.py:533(identify)
   257942    0.499    0.000    1.664    0.000 symbol.py:306(_hashable_content)
      898    1.242    0.001    1.620    0.002 identification.py:17(pslq)
    45/27    0.000    0.000    1.325    0.049 basic.py:1673(simplify)
        1    0.000    0.000    1.325    1.325 boolalg.py:445(<listcomp>)
374293/373977    0.566    0.000    1.201    0.000 sympify.py:92(sympify)
        1    0.000    0.000    1.168    1.168 boolalg.py:2796(simplify_logic)
```
I don't know if that is the right sort of example to test this code but it does show the impact that some of these algorithms have as soon as you scale up beyond tiny inputs. The input I used here is not what I would call "large":
```python
In [25]: fact
Out[25]: x₀ < x₁ ∧ x₁ < x₂ ∧ x₂ < x₃ ∧ x₃ < x₄ ∧ x₄ < x₅ ∧ x₅ < x₆ ∧ x₆ < x₇ ∧ x₇ < x₈ ∧ x₈ < x₉
```
@oscargus can you suggest what would be a good way to time this code in a case where it succeeds in simplifying and also a case where it does not?
oscargus: Well, it is a set cover problem that is solved, one of Karp's original 21 NP-complete problems, https://en.wikipedia.org/wiki/Set_cover_problem. So even now it is a heuristic, but slightly better than the previous (and I have updated it at least twice earlier, first to get better results/smaller expressions, second to get better performance). The new code is not really affecting the performance as it is "only" n^3 as far as I can tell. Another issue is that n is here (up to) 2^terms, hence, the `force` keyword argument which, unless set to `True`, disables this part if there are more than eight variables.

It would indeed be interesting to get some benchmark numbers in though. I think I've looked into that earlier, but didn't really come up with anything (and at that time, the current benchmark action didn't exist).

In your example, you have nine variables (each inequality becomes a variable in this respect), so not really sure if this part of the code is called. On the other hand, there is only one expression for the boolean simplification, so it wouldn't really become that time consuming. Basically, for this part of the code, your problem is just simplifying A & B & C & D ..., so that is very easy (already simplified as there is a single AND-expression, cannot really be any smaller).
oscargus: I think that there are probably slightly more efficient ways to obtain a similar heuristic as the current one. Not sure if it is really the purpose of this PR though. No worse asymptotic performance and better results in some cases (not worse compares to the previous).
oscarbenjamin: Okay that makes sense. Then this looks good to merge.

I do think that the performance of logic simplification needs to be improved in general though.
oscargus: Sounds good!

I'm rewriting it a bit now, but not sure if there is any major changes in performance. Some things will be better, some may be worse, so probably depends on case if it improves. A bit better structure though compared to the current one...
oscargus: The thing is that for small problems the current algorithm is enough and I think it may give close to optimal results. For larger problems, it is not. One could possibly use things like https://en.wikipedia.org/wiki/Espresso_heuristic_logic_minimizer but then we are talking external dependencies etc. Fundamentally hard problem when the problem size grows.
oscargus: Nah, I didn't manage to get the improved version to work and will need to start all over at some later stage. So if this is merged, I can rebase other related PRs.
oscargus: Gah! Found the problem, so will push a new version in a few minutes...
oscarbenjamin: > Fundamentally hard problem when the problem size grows.

Yes, but then it shouldn't be attempted even by `simplify` without some additional `force` flag or something. The flipside of any algorithm like this in a system like SymPy is that if the algorithm is being invoked indirectly then users don't actually control when that happens. Here this is called by `simplify` so in many cases the cost of the algorithm will be paid for no benefit e.g. the example above where it takes 10 seconds for simplify to achieve nothing.

We have to remember that performance is not orthogonal to capability: often in SymPy something is so slow that effectively the result can not be computed (e.g. "simplify hangs for this expression"). When that happens the impact can be much wider than the type of problem that the slow algorithm was designed for.
oscargus: I just bench marked a bit for the problem in #14700 which has eight variables.

Of the 2.5 seconds it takes, 5 ms is in the currently discussed code (n^4). In your code example above, simplify_logic takes 1.168 seconds out of the 14. I would expect that the time for the n^4-code is about the same there. Most of the time in simplify_logic is spent doing `xreplace` for small number of variables.

Yes, hanging is one thing, getting "why isn't this Boolean expression fully simplified?" is another. Yes, the time may increase further for larger problems, but currently it is repeated calls to small pieces of code that takes quite a bit of time. For example, in my test, 47193 calls to `sympify`, requiring a total of 127 ms local time, for simplifying an expression with plain symbols and operators.

One should also distinguish between `simplify` and `simplify_logic`. Two quite different functions, where the second is only a small part of the first.



Patch:
diff --git a/sympy/logic/boolalg.py b/sympy/logic/boolalg.py
index a5ae35cb90..220e7289b4 100644
--- a/sympy/logic/boolalg.py
+++ b/sympy/logic/boolalg.py
@@ -573,6 +573,8 @@ def _apply_patternbased_simplification(self, rv, patterns, measure,
 
         """
         from sympy.core.relational import Relational, _canonical
+        from sympy.functions.elementary.miscellaneous import Min, Max
+
         if replacementvalue is None and dominatingvalue is not None:
             replacementvalue = dominatingvalue
         # Use replacement patterns for Relationals
@@ -627,9 +629,9 @@ def _apply_patternbased_simplification(self, rv, patterns, measure,
                                 # will be replacementvalue
                                 return replacementvalue
                             # add replacement
-                            if not isinstance(np, ITE):
-                                # We only want to use ITE replacements if
-                                # they simplify to a relational
+                            if not isinstance(np, ITE) and not np.has(Min, Max):
+                                # We only want to use ITE and Min/Max
+                                # replacements if they simplify away
                                 costsaving = measure(oldexpr) - measure(np)
                                 if costsaving > 0:
                                     results.append((costsaving, (i, j, np)))
@@ -740,13 +742,12 @@ def _eval_simplify(self, **kwargs):
         if not Rel:
             return rv
         eqs, other = sift(Rel, lambda i: isinstance(i, Equality), binary=True)
-        if not eqs:
-            return rv
 
-        measure, ratio = kwargs['measure'], kwargs['ratio']
-        reps = {}
-        sifted = {}
+        measure = kwargs['measure']
         if eqs:
+            ratio = kwargs['ratio']
+            reps = {}
+            sifted = {}
             # group by length of free symbols
             sifted = sift(ordered([
                 (i.free_symbols, i) for i in eqs]),
@@ -783,11 +784,11 @@ def _eval_simplify(self, **kwargs):
                         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 + nonlineqs + other)] + nonRel))
+            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 + nonlineqs + other)] + nonRel))
         patterns = simplify_patterns_and()
         return self._apply_patternbased_simplification(rv, patterns,
                                                        measure, False)
@@ -2154,7 +2155,7 @@ def truth_table(expr, variables, input=True):
     if not isinstance(expr, BooleanFunction) and not is_literal(expr):
         return
 
-    table = product([0, 1], repeat=len(variables))
+    table = product((0, 1), repeat=len(variables))
     for term in table:
         term = list(term)
         value = expr.xreplace(dict(zip(variables, term)))
@@ -2251,11 +2252,9 @@ def _simplified_pairs(terms):
 
     # Count number of ones as _check_pair can only potentially match if there
     # is at most a difference of a single one
-    def ones_count(term):
-        return sum([1 for t in term if t == 1])
     termdict = defaultdict(list)
     for n, term in enumerate(terms):
-        ones = ones_count(term)
+        ones = sum([1 for t in term if t == 1])
         termdict[ones].append(n)
 
     variables = len(terms[0])
@@ -2279,22 +2278,10 @@ def ones_count(term):
         simplified_terms = _simplified_pairs(simplified_terms)
 
     # Add remaining, non-simplified, terms
-    simplified_terms.extend(
-        [terms[i] for i in [_ for _ in todo if _ is not None]])
+    simplified_terms.extend([terms[i] for i in todo if i is not None])
     return simplified_terms
 
 
-def _compare_term(minterm, term):
-    """
-    Return True if a binary term is satisfied by the given term. Used
-    for recognizing prime implicants.
-    """
-    for m, t in zip(minterm, term):
-        if t != 3 and m != t:
-            return False
-    return True
-
-
 def _rem_redundancy(l1, terms):
     """
     After the truth table has been sufficiently simplified, use the prime
@@ -2307,18 +2294,18 @@ def _rem_redundancy(l1, terms):
 
     nterms = len(terms)
     nl1 = len(l1)
+
     # Create dominating matrix
     dommatrix = [[0]*nl1 for n in range(nterms)]
+    colcount = [0]*nl1
+    rowcount = [0]*nterms
     for primei, prime in enumerate(l1):
         for termi, term in enumerate(terms):
-            if _compare_term(term, prime):
+            # Check prime implicant covering term
+            if all(t == 3 or t == mt for t, mt in zip(prime, term)):
                 dommatrix[termi][primei] = 1
-
-    # Non-dominated prime implicants, dominated to be removed
-    ndprimeimplicants = set(range(nl1))
-
-    # Non-dominated terms, dominated to be removed
-    ndterms = set(range(nterms))
+                colcount[primei] += 1
+                rowcount[termi] += 1
 
     # Keep track if anything changed
     anythingchanged = True
@@ -2326,41 +2313,75 @@ def _rem_redundancy(l1, terms):
     while anythingchanged:
         anythingchanged = False
 
-        # Make copy for iteration
-        oldndterms = ndterms.copy()
-        # Filter matrix to only get non-dominated items
-        filteredrows = [[dommatrix[rowi][i] for i in list(ndprimeimplicants)]
-                        for rowi in oldndterms]
-        for n, rowi in enumerate(oldndterms):
+        for rowi in range(nterms):
             # Still non-dominated?
-            if rowi in ndterms:
-                row = filteredrows[n]
-                for n2, row2i in enumerate(oldndterms):
+            if rowcount[rowi]:
+                row = dommatrix[rowi]
+                for row2i in range(nterms):
                     # Still non-dominated?
-                    if n != n2 and row2i in ndterms:
-                        if all(a >= b for (a, b) in zip(filteredrows[n2], row)):
+                    if rowi != row2i and rowcount[rowi] and (rowcount[rowi] <= rowcount[row2i]):
+                        row2 = dommatrix[row2i]
+                        if all(row2[n] >= row[n] for n in range(nl1)):
                             # row2 dominating row, remove row2
-                            ndterms.remove(row2i)
+                            rowcount[row2i] = 0
                             anythingchanged = True
+                            for primei, prime in enumerate(row2):
+                                if prime:
+                                    # Make corresponding entry 0
+                                    dommatrix[row2i][primei] = 0
+                                    colcount[primei] -= 1
+
+        colcache = dict()
 
-        # Make copy for iteration
-        oldndprimeimplicants = ndprimeimplicants.copy()
-        # Filter matrix to only get non-dominated items
-        filteredcols = [[dommatrix[i][coli] for i in list(ndterms)]
-                        for coli in oldndprimeimplicants]
-        for n, coli in enumerate(oldndprimeimplicants):
+        for coli in range(nl1):
             # Still non-dominated?
-            if coli in ndprimeimplicants:
-                col = filteredcols[n]
-                for n2, col2i in enumerate(oldndprimeimplicants):
+            if colcount[coli]:
+                if coli in colcache:
+                    col = colcache[coli]
+                else:
+                    col = [dommatrix[i][coli] for i in range(nterms)]
+                    colcache[coli] = col
+                for col2i in range(nl1):
                     # Still non-dominated?
-                    if coli != col2i and col2i in ndprimeimplicants:
-                        if all(a >= b for (a, b) in zip(col, filteredcols[n2])):
+                    if coli != col2i and colcount[col2i] and (colcount[coli] >= colcount[col2i]):
+                        if col2i in colcache:
+                            col2 = colcache[col2i]
+                        else:
+                            col2 = [dommatrix[i][col2i] for i in range(nterms)]
+                            colcache[col2i] = col2
+                        if all(col[n] >= col2[n] for n in range(nterms)):
                             # col dominating col2, remove col2
-                            ndprimeimplicants.remove(col2i)
+                            colcount[col2i] = 0
                             anythingchanged = True
-
-    return [l1[i] for i in ndprimeimplicants]
+                            for termi, term in enumerate(col2):
+                                if term and dommatrix[termi][col2i]:
+                                    # Make corresponding entry 0
+                                    dommatrix[termi][col2i] = 0
+                                    rowcount[termi] -= 1
+
+        if not anythingchanged:
+            # Heuristically select the prime implicant covering most terms
+            maxterms = 0
+            bestcolidx = -1
+            for coli in range(nl1):
+                s = colcount[coli]
+                if s > maxterms:
+                    bestcolidx = coli
+                    maxterms = s
+
+            # In case we found a prime implicant covering at least two terms
+            if bestcolidx != -1 and maxterms > 1:
+                for primei, prime in enumerate(l1):
+                    if primei != bestcolidx:
+                        for termi, term in enumerate(colcache[bestcolidx]):
+                            if term and dommatrix[termi][primei]:
+                                # Make corresponding entry 0
+                                dommatrix[termi][primei] = 0
+                                anythingchanged = True
+                                rowcount[termi] -= 1
+                                colcount[primei] -= 1
+
+    return [l1[i] for i in range(nl1) if colcount[i]]
 
 
 def _input_to_binlist(inputlist, variables):
@@ -2373,7 +2394,7 @@ def _input_to_binlist(inputlist, variables):
             nonspecvars = list(variables)
             for key in val.keys():
                 nonspecvars.remove(key)
-            for t in product([0, 1], repeat=len(nonspecvars)):
+            for t in product((0, 1), repeat=len(nonspecvars)):
                 d = dict(zip(nonspecvars, t))
                 d.update(val)
                 binlist.append([d[v] for v in variables])
@@ -2443,10 +2464,12 @@ def SOPform(variables, minterms, dontcares=None):
     .. [1] https://en.wikipedia.org/wiki/Quine-McCluskey_algorithm
 
     """
-    variables = [sympify(v) for v in variables]
     if minterms == []:
         return false
 
+    variables = tuple(map(sympify, variables))
+
+
     minterms = _input_to_binlist(minterms, variables)
     dontcares = _input_to_binlist((dontcares or []), variables)
     for d in dontcares:
@@ -2513,10 +2536,10 @@ def POSform(variables, minterms, dontcares=None):
     .. [1] https://en.wikipedia.org/wiki/Quine-McCluskey_algorithm
 
     """
-    variables = [sympify(v) for v in variables]
     if minterms == []:
         return false
 
+    variables = tuple(map(sympify, variables))
     minterms = _input_to_binlist(minterms, variables)
     dontcares = _input_to_binlist((dontcares or []), variables)
     for d in dontcares:
@@ -2524,7 +2547,7 @@ def POSform(variables, minterms, dontcares=None):
             raise ValueError('%s in minterms is also in dontcares' % d)
 
     maxterms = []
-    for t in product([0, 1], repeat=len(variables)):
+    for t in product((0, 1), repeat=len(variables)):
         t = list(t)
         if (t not in minterms) and (t not in dontcares):
             maxterms.append(t)
@@ -2583,12 +2606,12 @@ def ANFform(variables, truthvalues):
         raise ValueError("The number of truth values must be equal to 2^%d, "
                          "got %d" % (n_vars, n_values))
 
-    variables = [sympify(v) for v in variables]
+    variables = tuple(map(sympify, variables))
 
     coeffs = anf_coeffs(truthvalues)
     terms = []
 
-    for i, t in enumerate(product([0, 1], repeat=n_vars)):
+    for i, t in enumerate(product((0, 1), repeat=n_vars)):
         if coeffs[i] == 1:
             terms.append(t)
 
@@ -2682,7 +2705,7 @@ def bool_minterm(k, variables):
     """
     if isinstance(k, int):
         k = integer_to_term(k, len(variables))
-    variables = list(map(sympify, variables))
+    variables = tuple(map(sympify, variables))
     return _convert_to_varsSOP(k, variables)
 
 
@@ -2718,7 +2741,7 @@ def bool_maxterm(k, variables):
     """
     if isinstance(k, int):
         k = integer_to_term(k, len(variables))
-    variables = list(map(sympify, variables))
+    variables = tuple(map(sympify, variables))
     return _convert_to_varsPOS(k, variables)
 
 
@@ -2759,7 +2782,7 @@ def bool_monomial(k, variables):
     """
     if isinstance(k, int):
         k = integer_to_term(k, len(variables))
-    variables = list(map(sympify, variables))
+    variables = tuple(map(sympify, variables))
     return _convert_to_varsANF(k, variables)
 
 
@@ -2772,7 +2795,7 @@ def _find_predicates(expr):
     """
     if not isinstance(expr, BooleanFunction):
         return {expr}
-    return set().union(*(_find_predicates(i) for i in expr.args))
+    return set().union(*(map(_find_predicates, expr.args)))
 
 
 def simplify_logic(expr, form=None, deep=True, force=False):
@@ -2822,20 +2845,22 @@ def simplify_logic(expr, form=None, deep=True, force=False):
     if form not in (None, 'cnf', 'dnf'):
         raise ValueError("form can be cnf or dnf only")
     expr = sympify(expr)
-    # check for quick exit: right form and all args are
+    # check for quick exit if form is given: right form and all args are
     # literal and do not involve Not
-    isc = is_cnf(expr)
-    isd = is_dnf(expr)
-    form_ok = (
-        isc and form == 'cnf' or
-        isd and form == 'dnf')
-    if form_ok and all(is_literal(a)
-            for a in expr.args):
-        return expr
+    if form:
+        form_ok = False
+        if form == 'cnf':
+            form_ok = is_cnf(expr)
+        elif form == 'dnf':
+            form_ok = is_dnf(expr)
+
+        if form_ok and all(is_literal(a)
+                for a in expr.args):
+            return expr
     if deep:
         variables = _find_predicates(expr)
         from sympy.simplify.simplify import simplify
-        s = [simplify(v) for v in variables]
+        s = tuple(map(simplify, variables))
         expr = expr.xreplace(dict(zip(variables, s)))
     if not isinstance(expr, BooleanFunction):
         return expr
@@ -2845,12 +2870,12 @@ def simplify_logic(expr, form=None, deep=True, force=False):
     if not force and len(variables) > 8:
         return expr
     # group into constants and variable values
-    c, v = sift(variables, lambda x: x in (True, False), binary=True)
+    c, v = sift(ordered(variables), lambda x: x in (True, False), binary=True)
     variables = c + v
     truthtable = []
     # standardize constants to be 1 or 0 in keeping with truthtable
     c = [1 if i == True else 0 for i in c]
-    for t in product([0, 1], repeat=len(v)):
+    for t in product((0, 1), repeat=len(v)):
         if expr.xreplace(dict(zip(v, t))) == True:
             truthtable.append(c + list(t))
     big = len(truthtable) >= (2 ** (len(variables) - 1))


Test patch:
diff --git a/sympy/logic/tests/test_boolalg.py b/sympy/logic/tests/test_boolalg.py
index 9013af9e5b..a5d4598e45 100644
--- a/sympy/logic/tests/test_boolalg.py
+++ b/sympy/logic/tests/test_boolalg.py
@@ -5,7 +5,6 @@
 from sympy.core.singleton import S
 from sympy.core.symbol import (Dummy, symbols)
 from sympy.functions import Piecewise
-from sympy.functions.elementary.miscellaneous import Max, Min
 from sympy.functions.elementary.trigonometric import sin
 from sympy.sets.sets import (EmptySet, Interval, Union)
 from sympy.simplify.simplify import simplify
@@ -1080,12 +1079,6 @@ def test_relational_simplification():
     assert And(x >= y, x < y).simplify() == S.false
     assert Or(x >= y, Eq(y, x)).simplify() == (x >= y)
     assert And(x >= y, Eq(y, x)).simplify() == Eq(x, y)
-    assert Or(Eq(x, y), x >= y, w < y, z < y).simplify() == \
-        Or(x >= y, y > Min(w, z))
-    assert And(Eq(x, y), x >= y, w < y, y >= z, z < y).simplify() == \
-        And(Eq(x, y), y > Max(w, z))
-    assert Or(Eq(x, y), x >= 1, 2 < y, y >= 5, z < y).simplify() == \
-        (Eq(x, y) | (x >= 1) | (y > Min(2, z)))
     assert And(Eq(x, y), x >= 1, 2 < y, y >= 5, z < y).simplify() == \
         (Eq(x, y) & (x >= 1) & (y >= 5) & (y > z))
     assert (Eq(x, y) & Eq(d, e) & (x >= y) & (d >= e)).simplify() == \
@@ -1207,6 +1200,19 @@ def test_check_pair():
     assert _check_pair([0, 1, 0], [1, 1, 1]) == -1
 
 
+def test_issue_19114():
+    expr = (B & C) | (A & ~C) | (~A & ~B)
+    result = to_dnf(expr, simplify=True)
+    assert result == expr
+
+
+def test_issue_20870():
+    result = SOPform([a, b, c, d], [1, 2, 3, 4, 5, 6, 8, 9, 11, 12, 14, 15])
+    expected = ((d & ~b) | (a & b & c) | (a & ~c & ~d) |
+                (b & ~a & ~c) | (c & ~a & ~d))
+    assert result == expected
+
+
 def test_convert_to_varsSOP():
     assert _convert_to_varsSOP([0, 1, 0], [x, y, z]) ==  And(Not(x), y, Not(z))
     assert _convert_to_varsSOP([3, 1, 0], [x, y, z]) ==  And(y, Not(z))
diff --git a/sympy/simplify/tests/test_simplify.py b/sympy/simplify/tests/test_simplify.py
index 02bea1389c..f787592c3c 100644
--- a/sympy/simplify/tests/test_simplify.py
+++ b/sympy/simplify/tests/test_simplify.py
@@ -7,13 +7,13 @@
     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, And, Le)
+    zoo, And, Gt, Ge, Le, Or)
 from sympy.core.mul import _keep_coeff
 from sympy.core.expr import unchanged
 from sympy.simplify.simplify import nthroot, inversecombine
 from sympy.testing.pytest import XFAIL, slow, _both_exp_pow
 
-from sympy.abc import x, y, z, t, a, b, c, d, e, f, g, h, i
+from sympy.abc import x, y, z, t, a, b, c, d, e, f, g, h, i, n
 
 
 def test_issue_7263():
@@ -947,6 +947,43 @@ def test_issue_17292():
     # this is bigger than the issue: check that deep processing works
     assert simplify(5*abs((x**2 - 1)/(x - 1))) == 5*Abs(x + 1)
 
+
+def test_issue_19822():
+    expr = And(Gt(n-2, 1), Gt(n, 1))
+    assert simplify(expr) == Gt(n, 3)
+
+
+def test_issue_18645():
+    expr = And(Ge(x, 3), Le(x, 3))
+    assert simplify(expr) == Eq(x, 3)
+    expr = And(Eq(x, 3), Le(x, 3))
+    assert simplify(expr) == Eq(x, 3)
+
+
+@XFAIL
+def test_issue_18642():
+    i = Symbol("i", integer=True)
+    n = Symbol("n", integer=True)
+    expr = And(Eq(i, 2 * n), Le(i, 2*n -1))
+    assert simplify(expr) == S.false
+
+
+@XFAIL
+def test_issue_18389():
+    n = Symbol("n", integer=True)
+    expr = Eq(n, 0) | (n >= 1)
+    assert simplify(expr) == Ge(n, 0)
+
+
+def test_issue_8373():
+    x = Symbol('x', real=True)
+    assert simplify(Or(x < 1, x >= 1)) == S.true
+
+
+def test_issue_7950():
+    expr = And(Eq(x, 1), Eq(x, 2))
+    assert simplify(expr) == S.false
+
 def test_issue_19484():
     assert simplify(sign(x) * Abs(x)) == x