Repository: sympy/sympy

Issue:
Complexes and ComplexRegion have incorrect .contains() logic
Below `x` is just a general symbol:
```
In [1]: i1 = Interval(1, 2)

In [2]: i2 = Interval(3, 5)

In [3]: cr = ComplexRegion(i1*i2)

In [4]: cr.contains(x)
Out[4]: True

In [5]: cr
Out[5]: {x + y⋅ⅈ | x, y ∊ [1, 2] × [3, 5]}
```
This is because a symbolic `Boolean` expression is incorrectly used in a plain `if`:
https://github.com/sympy/sympy/blob/957f468b80325e33a40c37bbb50a6ca09e9784f8/sympy/sets/fancysets.py#L1264-L1271
Since the Sympy logic expression is neither a (Python) `False` nor `None`, it is interpreted as `True`.

The same bug is a few lines further down in the polar region handling.

If I'm reading the code correctly, there's more to fix: after all product sets are checked, `_contains` returns `False`, but the actual result could just as well be unknown, i.e. `None`.
(IIRC it is `Set.contains()` that then transforms a `None` value returned by a `_contains` method into a `Contains` object.)

Discussion:
sympy-bot: :white_check_mark:

Hi, I am the [SymPy bot](https://github.com/sympy/sympy-bot) (v149). 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:
* sets
  * Fixed bugs in `ComplexRegion.contains()` that could lead to incorrect results when membership is in fact undecidable or at least unknown. ([#18308](https://github.com/sympy/sympy/pull/18308) by [@gschintgen](https://github.com/gschintgen))

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

Note: This comment will be updated with the latest check if you edit the pull request. You need to reload the page to see it. <details><summary>Click here to see the pull request description that was parsed.</summary>

    #### 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://github.com/blog/1506-closing-issues-via-pull-requests . Please also
    write a comment on that issue linking back to this pull request once it is
    open. -->
    Fixes #18280

    #### Brief description of what is fixed or changed
    The following are fixed in _contains():

    - #18280: The method had no provisions for undecidable membership;
      instead of `None` it would return `False`.
    - #18280: The method checked `arg(other)` for membership in the
      respective theta intervals but those are normalized to [0, 2*pi)
      and not (-pi, pi].
    - For PolarComplexRegions the method incorrectly set theta = 0 if
      `other` was 0. So it returned a somewhat arbitrary truth value in
      this case. Now it will exclusively check complex modulus if other==0.

    Tests are added.


    #### Other comments


    #### Release Notes

    <!-- Write the release notes for this release below. 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 -->
    * sets
      * Fixed bugs in `ComplexRegion.contains()` that could lead to incorrect results when membership is in fact undecidable or at least unknown.
    <!-- END RELEASE NOTES -->


</details><p>

oscarbenjamin: I see you have a test failure for this
```
File "/home/travis/virtualenv/python3.8.0/lib/python3.8/site-packages/sympy-1.6.dev0-py3.8.egg/sympy/solvers/solveset.py", line 3286, in sympy.solvers.solveset.nonlinsolve
F
Failed example:
    pprint(nonlinsolve(system, [a, b, c, d]), use_unicode=False)
Expected:
      -1       1               1      -1
    {(---, -d, -, {d} \ {0}), (-, -d, ---, {d} \ {0})}
       d       d               d       d
Got:
       -1                       1                         1                      -
    {({---} \ S.Complexes, -d, {-} \ S.Complexes, {d}), ({-} \ S.Complexes, -d, {-
        d                       d                         d                       
    <BLANKLINE>
    1                       
    --} \ S.Complexes, {d})}
    d  
```
I've seen this precise example come up in many PRs including my own. It seems to be very flakey.
gschintgen: > 
> 
> I see you have a test failure for this
> 
> ```
> File "/home/travis/virtualenv/python3.8.0/lib/python3.8/site-packages/sympy-1.6.dev0-py3.8.egg/sympy/solvers/solveset.py", line 3286, in sympy.solvers.solveset.nonlinsolve
> F
> Failed example:
>     pprint(nonlinsolve(system, [a, b, c, d]), use_unicode=False)
> Expected:
>       -1       1               1      -1
>     {(---, -d, -, {d} \ {0}), (-, -d, ---, {d} \ {0})}
>        d       d               d       d
> Got:
>        -1                       1                         1                      -
>     {({---} \ S.Complexes, -d, {-} \ S.Complexes, {d}), ({-} \ S.Complexes, -d, {-
>         d                       d                         d                       
>     <BLANKLINE>
>     1                       
>     --} \ S.Complexes, {d})}
>     d  
> ```
> 
> I've seen this precise example come up in many PRs including my own. It seems to be very flakey.

It's flaky indeed. It's the one I meant [here](https://github.com/sympy/sympy/issues/18280#issuecomment-573101883) and that got me to investigate.

I noticed that `nonlinsolve` has some broken logic: https://github.com/sympy/sympy/blob/8a7dc8ef786cccbd7dc0fe0b0e6beec4dd1f7b49/sympy/solvers/solveset.py#L2759-L2771
The above code block is called with both `intersection_true` and `complement_true` `True`, so that the second block overrides the first one...

The commit I pushed just now handles almost all the other test failures. But there's another remaining test failure that pointed me to a very unfortunate fact: Many tests in `test_solveset.py` assume `x` and `y` to have been created beforehand, which is true indeed, but they're actually instantiated with `real=True` assumption! I don't believe that's intentional and it could falsify a *lot* of tests...
https://github.com/sympy/sympy/blob/8a7dc8ef786cccbd7dc0fe0b0e6beec4dd1f7b49/sympy/solvers/tests/test_solveset.py#L47-L55
oscarbenjamin: > I don't believe that's intentional and it could falsify a _lot_ of tests...

Blame shows that it dates back to the creation of the solveset module.
gschintgen: Strange: With the latest commit the number of Travis test failures in `test_solveset` went up instead of down. The failures on Travis all mention `_timeout` so I thought that maybe there was a Travis hiccup, so I restarted this specific job, but the results are the same. I'm also now noticing that in general tests are started via a _timeout function so that's normal.

What's bothering me is that when testing locally I only have a single failure (`test_solve_polynomial_symbolic_param`) with `test("solveset")`. I also ran the two first failing Travis tests manually on my local setup and everything was fine.
gschintgen: >  I also ran the two first failing Travis tests manually on my local setup and everything was fine.

Never mind. I think I figured it out. I had a closer look at the Travis logs and it seems what is actually tested is not the PR branch but a "fictitious" merged branch. That's sensible I guess.
oscarbenjamin: Yes, Travis tests the PR after merging it into master (that's why the tests don't run when there's a merge conflict). That should mean that conflicting PRs can't lead to test failures on master but there is a loop-hole: if the tests have passed then other PRs can be merged to master and the tests won't re-run but the PR remains mergeable.

Patch:
diff --git a/sympy/polys/polytools.py b/sympy/polys/polytools.py
index 1f15d8d98c..e148ad1a1a 100644
--- a/sympy/polys/polytools.py
+++ b/sympy/polys/polytools.py
@@ -6625,7 +6625,7 @@ def cancel(f, *gens, **args):
         (F, G), opt = parallel_poly_from_expr((p, q), *gens, **args)
     except PolificationFailed:
         if not isinstance(f, (tuple, Tuple)):
-            return f.expand()
+            return f
         else:
             return S.One, p, q
     except PolynomialError as msg:
diff --git a/sympy/sets/fancysets.py b/sympy/sets/fancysets.py
index 8421aa0742..34a35faed2 100644
--- a/sympy/sets/fancysets.py
+++ b/sympy/sets/fancysets.py
@@ -6,7 +6,7 @@
 from sympy.core.containers import Tuple
 from sympy.core.expr import Expr
 from sympy.core.function import Lambda
-from sympy.core.logic import fuzzy_not, fuzzy_or
+from sympy.core.logic import fuzzy_not, fuzzy_or, fuzzy_and
 from sympy.core.numbers import oo, Integer
 from sympy.core.relational import Eq
 from sympy.core.singleton import Singleton, S
@@ -1257,25 +1257,28 @@ def _contains(self, other):
         # self in rectangular form
         if not self.polar:
             re, im = other if isTuple else other.as_real_imag()
-            for element in self.psets:
-                if And(element.args[0]._contains(re),
-                        element.args[1]._contains(im)):
-                    return True
-            return False
+            return fuzzy_or(fuzzy_and([
+                pset.args[0]._contains(re),
+                pset.args[1]._contains(im)])
+                for pset in self.psets)
 
         # self in polar form
         elif self.polar:
+            if other.is_zero:
+                # ignore undefined complex argument
+                return fuzzy_or(pset.args[0]._contains(S.Zero)
+                    for pset in self.psets)
             if isTuple:
                 r, theta = other
-            elif other.is_zero:
-                r, theta = S.Zero, S.Zero
             else:
                 r, theta = Abs(other), arg(other)
-            for element in self.psets:
-                if And(element.args[0]._contains(r),
-                        element.args[1]._contains(theta)):
-                    return True
-            return False
+            if theta.is_real and theta.is_number:
+                # angles in psets are normalized to [0, 2pi)
+                theta %= 2*S.Pi
+                return fuzzy_or(fuzzy_and([
+                    pset.args[0]._contains(r),
+                    pset.args[1]._contains(theta)])
+                    for pset in self.psets)
 
 
 class CartesianComplexRegion(ComplexRegion):
diff --git a/sympy/solvers/solveset.py b/sympy/solvers/solveset.py
index 3481c40b21..e7e6466890 100644
--- a/sympy/solvers/solveset.py
+++ b/sympy/solvers/solveset.py
@@ -167,7 +167,12 @@ def _invert(f_x, y, x, domain=S.Complexes):
     if not isinstance(s, FiniteSet) or x1 != x:
         return x1, s
 
-    return x1, s.intersection(domain)
+    # Avoid adding gratuitous intersections with S.Complexes. Actual
+    # conditions should be handled by the respective inverters.
+    if domain is S.Complexes:
+        return x1, s
+    else:
+        return x1, s.intersection(domain)
 
 
 invert_complex = _invert
@@ -702,7 +707,9 @@ def _solve_as_poly(f, symbol, domain=S.Complexes):
                     for s in result]):
                 s = Dummy('s')
                 result = imageset(Lambda(s, expand_complex(s)), result)
-        if isinstance(result, FiniteSet):
+        if isinstance(result, FiniteSet) and domain != S.Complexes:
+            # Avoid adding gratuitous intersections with S.Complexes. Actual
+            # conditions should be handled elsewhere.
             result = result.intersection(domain)
         return result
     else:
@@ -2708,9 +2715,7 @@ def substitution(system, symbols, result=[{}], known_symbols=[],
                'Not type %s: %s')
         raise TypeError(filldedent(msg % (type(symbols), symbols)))
 
-    sym = getattr(symbols[0], 'is_Symbol', False)
-
-    if not sym:
+    if not getattr(symbols[0], 'is_Symbol', False):
         msg = ('Iterable of symbols must be given as '
                'second argument, not type %s: %s')
         raise ValueError(filldedent(msg % (type(symbols[0]), symbols[0])))
@@ -2746,36 +2751,39 @@ def _unsolved_syms(eq, sort=False):
     eqs_in_better_order = list(
         ordered(system, lambda _: len(_unsolved_syms(_))))
 
-    def add_intersection_complement(result, sym_set, **flags):
-        # If solveset have returned some intersection/complement
-        # for any symbol. It will be added in final solution.
+    def add_intersection_complement(result, intersection_dict, complement_dict):
+        # If solveset has returned some intersection/complement
+        # for any symbol, it will be added in the final solution.
         final_result = []
         for res in result:
             res_copy = res
             for key_res, value_res in res.items():
-                # Intersection/complement is in Interval or Set.
-                intersection_true = flags.get('Intersection', True)
-                complements_true = flags.get('Complement', True)
-                for key_sym, value_sym in sym_set.items():
+                intersect_set, complement_set = None, None
+                for key_sym, value_sym in intersection_dict.items():
                     if key_sym == key_res:
-                        if intersection_true:
-                            # testcase is not added for this line(intersection)
-                            new_value = \
-                                Intersection(FiniteSet(value_res), value_sym)
-                            if new_value is not S.EmptySet:
-                                res_copy[key_res] = new_value
-                        if complements_true:
-                            new_value = \
-                                Complement(FiniteSet(value_res), value_sym)
-                            if new_value is not S.EmptySet:
-                                res_copy[key_res] = new_value
+                        intersect_set = value_sym
+                for key_sym, value_sym in complement_dict.items():
+                    if key_sym == key_res:
+                        complement_set = value_sym
+                if intersect_set or complement_set:
+                    new_value = FiniteSet(value_res)
+                    if intersect_set and intersect_set != S.Complexes:
+                        new_value = Intersection(new_value, intersect_set)
+                    if complement_set:
+                        new_value = Complement(new_value, complement_set)
+                    if new_value is S.EmptySet:
+                        res_copy = {}
+                    elif new_value.is_FiniteSet and len(new_value) == 1:
+                        res_copy[key_res] = set(new_value).pop()
+                    else:
+                        res_copy[key_res] = new_value
             final_result.append(res_copy)
         return final_result
     # end of def add_intersection_complement()
 
-    def _extract_main_soln(sol, soln_imageset):
-            """separate the Complements, Intersections, ImageSet lambda expr
-            and it's base_set.
+    def _extract_main_soln(sym, sol, soln_imageset):
+            """Separate the Complements, Intersections, ImageSet lambda expr
+            and its base_set.
             """
             # if there is union, then need to check
             # Complement, Intersection, Imageset.
@@ -2788,11 +2796,10 @@ def _extract_main_soln(sol, soln_imageset):
                 # using `add_intersection_complement` method
             if isinstance(sol, Intersection):
                 # Interval/Set will be at 0th index always
-                if sol.args[0] != Interval(-oo, oo):
-                    # sometimes solveset returns soln
-                    # with intersection `S.Reals`, to confirm that
-                    # soln is in `domain=S.Reals` or not. We don't consider
-                    # that intersection.
+                if sol.args[0] not in (S.Reals, S.Complexes):
+                    # Sometimes solveset returns soln with intersection
+                    # S.Reals or S.Complexes. We don't consider that
+                    # intersection.
                     intersections[sym] = sol.args[0]
                 sol = sol.args[1]
             # after intersection and complement Imageset should
@@ -3027,13 +3034,13 @@ def _solve_using_known_values(result, solver):
 
                     if soln is not S.EmptySet:
                         soln, soln_imageset = _extract_main_soln(
-                            soln, soln_imageset)
+                            sym, soln, soln_imageset)
 
                     for sol in soln:
                         # sol is not a `Union` since we checked it
                         # before this loop
                         sol, soln_imageset = _extract_main_soln(
-                            sol, soln_imageset)
+                            sym, sol, soln_imageset)
                         sol = set(sol).pop()
                         free = sol.free_symbols
                         if got_symbol and any([
@@ -3117,17 +3124,9 @@ def _solve_using_known_values(result, solver):
         # eg : [{x: -1, y : 1}, {x : -y , y: y}] then
         # return [{x : -y, y : y}]
         result_all_variables = result_infinite
-    if intersections and complements:
-        # no testcase is added for this block
-        result_all_variables = add_intersection_complement(
-            result_all_variables, intersections,
-            Intersection=True, Complement=True)
-    elif intersections:
-        result_all_variables = add_intersection_complement(
-            result_all_variables, intersections, Intersection=True)
-    elif complements:
+    if intersections or complements:
         result_all_variables = add_intersection_complement(
-            result_all_variables, complements, Complement=True)
+            result_all_variables, intersections, complements)
 
     # convert to ordered tuple
     result = S.EmptySet


Test patch:
diff --git a/sympy/polys/tests/test_polytools.py b/sympy/polys/tests/test_polytools.py
index 236766afc7..89360b21ee 100644
--- a/sympy/polys/tests/test_polytools.py
+++ b/sympy/polys/tests/test_polytools.py
@@ -3315,7 +3315,3 @@ def test_issue_17988():
     p = poly(x - 1)
     M = Matrix([[poly(x + 1), poly(x + 1)]])
     assert p * M == M * p == Matrix([[poly(x**2 - 1), poly(x**2 - 1)]])
-
-def test_issue_18205():
-    assert cancel((2 + I)*(3 - I)) == 7 + I
-    assert cancel((2 + I)*(2 - I)) == 5
diff --git a/sympy/sets/tests/test_fancysets.py b/sympy/sets/tests/test_fancysets.py
index 59b9ffae9b..4d45df5175 100644
--- a/sympy/sets/tests/test_fancysets.py
+++ b/sympy/sets/tests/test_fancysets.py
@@ -691,6 +691,7 @@ def test_ImageSet_contains():
 
 
 def test_ComplexRegion_contains():
+    r = Symbol('r', real=True)
     # contains in ComplexRegion
     a = Interval(2, 3)
     b = Interval(4, 6)
@@ -703,6 +704,10 @@ def test_ComplexRegion_contains():
     assert 8 + 2.5*I in c2
     assert 2.5 + 6.1*I not in c1
     assert 4.5 + 3.2*I not in c1
+    assert c1.contains(x) == Contains(x, c1, evaluate=False)
+    assert c1.contains(r) == False
+    assert c2.contains(x) == Contains(x, c2, evaluate=False)
+    assert c2.contains(r) == False
 
     r1 = Interval(0, 1)
     theta1 = Interval(0, 2*S.Pi)
@@ -716,6 +721,25 @@ def test_ComplexRegion_contains():
     assert 0 in c3
     assert 1 + I not in c3
     assert 1 - I not in c3
+    assert c3.contains(x) == Contains(x, c3, evaluate=False)
+    assert c3.contains(r + 2*I) == Contains(
+        r + 2*I, c3, evaluate=False)  # is in fact False
+    assert c3.contains(1/(1 + r**2)) == Contains(
+        1/(1 + r**2), c3, evaluate=False)  # is in fact True
+
+    r2 = Interval(0, 3)
+    theta2 = Interval(pi, 2*pi, left_open=True)
+    c4 = ComplexRegion(r2*theta2, polar=True)
+    assert c4.contains(0) == True
+    assert c4.contains(2 + I) == False
+    assert c4.contains(-2 + I) == False
+    assert c4.contains(-2 - I) == True
+    assert c4.contains(2 - I) == True
+    assert c4.contains(-2) == False
+    assert c4.contains(2) == True
+    assert c4.contains(x) == Contains(x, c4, evaluate=False)
+    assert c4.contains(3/(1 + r**2)) == Contains(
+        3/(1 + r**2), c4, evaluate=False)  # is in fact True
 
     raises(ValueError, lambda: ComplexRegion(r1*theta1, polar=2))
 
diff --git a/sympy/solvers/tests/test_solveset.py b/sympy/solvers/tests/test_solveset.py
index c479d8eacd..0f92865426 100644
--- a/sympy/solvers/tests/test_solveset.py
+++ b/sympy/solvers/tests/test_solveset.py
@@ -1,5 +1,5 @@
 from sympy.core.containers import Tuple
-from sympy.core.function import (Function, Lambda, nfloat)
+from sympy.core.function import (Function, Lambda, nfloat, diff)
 from sympy.core.mod import Mod
 from sympy.core.numbers import (E, I, Rational, oo, pi)
 from sympy.core.relational import (Eq, Gt,
@@ -179,23 +179,39 @@ def test_issue_11536():
 
 
 def test_issue_17479():
-    import sympy as sb
-    from sympy.solvers.solveset import nonlinsolve
-    x, y, z = sb.symbols("x, y, z")
+    x, y, z = symbols("x, y, z")
     f = (x**2 + y**2)**2 + (x**2 + z**2)**2 - 2*(2*x**2 + y**2 + z**2)
-    fx = sb.diff(f, x)
-    fy = sb.diff(f, y)
-    fz = sb.diff(f, z)
+    fx = diff(f, x)
+    fy = diff(f, y)
+    fz = diff(f, z)
     sol = nonlinsolve([fx, fy, fz], [x, y, z])
-    # FIXME: This previously gave 18 solutions and now gives 20 due to fixes
-    # in the handling of intersection of FiniteSets or possibly a small change
-    # to ImageSet._contains. However Using expand I can turn this into 16
-    # solutions either way:
-    #
-    #    >>> len(FiniteSet(*(Tuple(*(expand(w) for w in s)) for s in sol)))
-    #    16
-    #
-    assert len(sol) == 20
+    assert len(sol) >= 4 and len(sol) <= 20
+    # nonlinsolve has been giving a varying number of solutions
+    # (originally 18, then 20, now 19) due to various internal changes.
+    # Unfortunately not all the solutions are actually valid and some are
+    # redundant. Since the original issue was that an exception was raised,
+    # this first test only checks that nonlinsolve returns a "plausible"
+    # solution set. The next test checks the result for correctness.
+
+
+@XFAIL
+def test_issue_18449():
+    x, y, z = symbols("x, y, z")
+    f = (x**2 + y**2)**2 + (x**2 + z**2)**2 - 2*(2*x**2 + y**2 + z**2)
+    fx = diff(f, x)
+    fy = diff(f, y)
+    fz = diff(f, z)
+    sol = nonlinsolve([fx, fy, fz], [x, y, z])
+    for (xs, ys, zs) in sol:
+        d = {x: xs, y: ys, z: zs}
+        assert tuple(_.subs(d).simplify() for _ in (fx, fy, fz)) == (0, 0, 0)
+    # After simplification and removal of duplicate elements, there should
+    # only be 4 parametric solutions left:
+    # simplifiedsolutions = FiniteSet((sqrt(1 - z**2), z, z),
+    #                                 (-sqrt(1 - z**2), z, z),
+    #                                 (sqrt(1 - z**2), -z, z),
+    #                                 (-sqrt(1 - z**2), -z, z))
+    # TODO: Is the above solution set definitely complete?
 
 
 def test_is_function_class_equation():