Repository: sympy/sympy
Pull request: NotImplementedError if nonlinsolve can't solve equations in real domain
<!-- 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://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 #15888 

#### Brief description of what is fixed or changed
If the equations are in real domain, then nonlinsolve() returns EmptySet() for most of the cases, because solveset applies `Gröbner bases` which do not work on Real domain. So instead of an EmptySet() now it will return a `NotImplementedError` with proper detail.

#### 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 -->
- solvers
    - Raise NonImplementedError when solveset fails with polynomials over inexact domain.
<!-- END RELEASE NOTES -->


Problem statement:
SolveSet: nonlinsolve return EmptySet() - Wrong Result
I was trying to solve two non-linear equations

```
eq1 = (x+y)**2 - 9
eq2 = x**2 - y**2 -3/4
nonlinsolve((eq1, eq2), (c1, c2))
```

But it returns

```
EmptySet()
```

I tried figuring it out and found that 

```
is_zero_dimensional(polys, symbols):
```

this function return False. Can someone explain me what this function does ?

Hints from issue discussion:
That function is expected to return True if the solution set is finite, as it is in this example. There must be a bug in the implementation.
You mean it just checks if the solution set is finite without actually finding the solutions.
Note that this works as expected if you use a rational value for `3/4` instead of a float.
yes you are right, I converted 3 to a sympy object using S(3) and the answer is correct.
So is that a bug or it was meant to work so.
It's still a bug since SymPy is generally expected to work with floats.
This is probably a case where floats have to be converted to rationals during computation.
using solve() gives the result even if I don't take rational(3,4) but gives answer in decimals.
Do you think same should be applied in solveset also.

Commit message: NotImplementedError if nonlinsolve can't solve

Patch to review:
diff --git a/sympy/solvers/solveset.py b/sympy/solvers/solveset.py
index 5da01c2e6261..e50fb6003421 100644
--- a/sympy/solvers/solveset.py
+++ b/sympy/solvers/solveset.py
@@ -2803,6 +2803,7 @@ def _handle_positive_dimensional(polys, symbols, denominators):
         new_system, symbols, result, [],
         denominators)
     return result
+
 # end of def _handle_positive_dimensional()
 
 
@@ -3007,6 +3008,7 @@ def nonlinsolve(system, *symbols):
 
     """
     from sympy.polys.polytools import is_zero_dimensional
+    from sympy.polys import RR
 
     if not system:
         return S.EmptySet
@@ -3051,7 +3053,11 @@ def nonlinsolve(system, *symbols):
                 return result
 
         # positive dimensional system
-        return _handle_positive_dimensional(polys, symbols, denominators)
+        res = _handle_positive_dimensional(polys, symbols, denominators)
+        if isinstance(res, EmptySet) and [True for p in polys if p.domain == RR]:
+            raise NotImplementedError("Equation not in exact domain. Try converting to rational")
+        else:
+            return res
 
     else:
         # If all the equations are not polynomial.


Code review comments:
diff_hunk: @@ -3051,7 +3053,11 @@ def nonlinsolve(system, *symbols):
                 return result
 
         # positive dimensional system
-        return _handle_positive_dimensional(polys, symbols, denominators)
+        res = _handle_positive_dimensional(polys, symbols, denominators)
+        if isinstance(res, EmptySet) and [True for p in polys if p.domain == RR]:
path: sympy/solvers/solveset.py
text: @user1:
How about `any(p.domain.is_ RR for p in polys)`?

More generally, there could be `not p.domain.is_Exact`.

@author:
Honestly I was also looking for something else instead of `True` but didn't knew about these functions. Let me quickly make another PR.