Repository: sympy/sympy
Pull request: Fix irrational exponent case in solveset._invert_real
Fixes #21236

Original code would raise ValueError on 'not base.is_positive'
assuming that the only situation for this case is when the
exponent is irrational. Issue 21236 showed that this case is
reached even when exponent is a Float. Also, the 'not base.is_positive'
is invalid for three-valued logic. solveset removes assumptions before
calling invert causing 'not base.is_positive' to be True regardless.

- [x] change odd/even tests for numer and denom to be accurate for rational -- they are currently written for Rational cases. The `rational` case tests should be something like 
```
numer, denom = expo.as_numer_denom()
denom.is_odd and numer.is_even # for numer % 2 == 0 case
denom.is_even and numer.is_odd # for denom % 2 == 0 case
else:  # unknown
```
<!-- BEGIN RELEASE NOTES -->
* solvers
  * Fixed bug in solveset 
<!-- END RELEASE NOTES -->


Problem statement:
computing probability involving a weibull distribution raises ValueError
OS: Linux Ubuntu 20.04
Python: 3.9.3
SymPy: 1.7.1

```python
>>> import sympy
>>> sympy.__version__
'1.7.1'
>>> from sympy.stats import Weibull, P
>>> X = Weibull("x", 8000, 1.5)
>>> P(X <= 4000)
...
ValueError: x**w where w is irrational is not defined for negative x
```
The problem appears to be with the shape factor being a floating point value causing SymPy to interpret the value 1.5 as an irrational number. Forcing the shape factor to be Rational gives the expected result.

```python
>>> from sympy import Rational, N
>>> X = Weibull("x", 8000, Rational(3, 2))
>>> P(X <= 4000)
1 - exp(-sqrt(2)/4)
>>> N(_)
0.297811498673440
```
Working directly with the cdf also gives the expected result.

```python
>>> from sympy.stats import cdf
>>> X = Weibull("x", 8000, 1.5)
>>> cdf(X)(4000)
0.297811498673440
```

This problem does not appear in SymPy version 1.5.1.

```python
>>> import sympy
>>> sympy.__version__
'1.5.1'
>>> from sympy.stats import Weibull, P
>>> X = Weibull("x", 8000, 1.5)
>>> P(X <= 4000)
0.297811498673440
```

Hints from issue discussion:
I'd like to work on this issue.
The particular issue described here could be fixed by this diff:
```diff
diff --git a/sympy/solvers/solveset.py b/sympy/solvers/solveset.py
index 17da5ae025..786041be5e 100644
--- a/sympy/solvers/solveset.py
+++ b/sympy/solvers/solveset.py
@@ -258,7 +258,7 @@ def _invert_real(f, g_ys, symbol):
                 else:
                     return _invert_real(base, res, symbol)
             else:
-                if not base.is_positive:
+                if base.is_positive is False:
                     raise ValueError("x**w where w is irrational is not "
                                      "defined for negative x")
                 return _invert_real(base, res, symbol)
```
The problem is that in this case `base` is a `Symbol` with `real=True` and it is not known to be positive:
```python
In [1]: base = Symbol('base', real=True)

In [2]: print(base.is_positive)
None
```
The `is_positive` property is part of the core assumptions and uses three-valued logic: True, False or None with None meaning unknown. It is not usually valid to use the plain Python `not` operator with a three-valued bool because `not None` gives `True`. In this case that means that an error is raised if the `base` is not known to be positive which is probably not correct.

A simpler way to see the error is:
```python
In [13]: solveset(x**1.5, x, Reals)
---------------------------------------------------------------------------
ValueError: x**w where w is irrational is not defined for negative x
```
@oscarbenjamin , as usual things are more complicated. `solveset` strips assumptions so that `base.is_positive` is always `None` no matter what the end user provides to `solveset`. While the change you propose works for the original problem, it breaks the existing test suite and does not allow for the case where `base` is explicitly non-positive.
So we're trying to rearrange
```
f(x)**w  in  B
```
 into
```
f(x) in {real_root(b, w) where b in B}
```
 where `x`, `f(x)`, `w`, are real and `B` is a subset of the Reals and `w` is irrational.

If `f(x)` is positive (for every real value of x) then that's simple. If `f(x)`.

Now consider if `f(x)` is negative for some x. We can write `f(x)` as `|f(x)| * exp(I*pi)` and hence
```
f(x)**w = |f(x)|**w * exp(I*pi)**w
```
Here `exp(I*pi)**w` can never be real if `w` is irrational because it has argument `pi*w + 2*n*pi`. That means that there can be no solutions where `f(x) < 0`.

So maybe the solution is something like
```
f(x) in {real_root(b, w), b in B} & Interval(0, oo)
```
In other words rather than raising an error we can adjoin this condition on `f(x)` to the rhs set `res`.
I think this does it:
```diff
diff --git a/sympy/solvers/solveset.py b/sympy/solvers/solveset.py
index aecdb2c2cb..6c654cf708 100644
--- a/sympy/solvers/solveset.py
+++ b/sympy/solvers/solveset.py
@@ -257,10 +257,9 @@ def _invert_real(f, g_ys, symbol):
 
                 else:
                     return _invert_real(base, res, symbol)
-            else:
-                if not base.is_positive:
-                    raise ValueError("x**w where w is irrational is not "
-                                     "defined for negative x")
+            elif expo.is_irrational or expo.is_Float:
+                g_ys_pos = g_ys & Interval(0, oo)
+                res = imageset(Lambda(n, real_root(n, expo)), g_ys_pos)
                 return _invert_real(base, res, symbol)
 
         if not base_has_sym:
```
That gives e.g.:
```python
In [7]: solveset(x**pi-E, x, Reals)
Out[7]: 
⎧ 1⎫
⎪ ─⎪
⎨ π⎬
⎪ℯ ⎪
⎩  ⎭
```

@oscarbenjamin, well not quite. What about the case when `expo=1.0`? With your proposed solution:
```python
>>> solveset(x**1.0 + 1, x, Reals)
EmptySet
>>> solveset(x**1 + 1, x, Reals)
FiniteSet(-1)
```
@oscarbenjamin, moving the `Float` case up appears to resolve the previous issue:

```diff
diff --git a/a/sympy/solvers/solveset.py b/b/sympy/solvers/solveset.py
index 4eb82a9..40683a7 100644
--- a/a/sympy/solvers/solveset.py
+++ b/b/sympy/solvers/solveset.py
@@ -240,7 +240,7 @@ def _invert_real(f, g_ys, symbol):

         if not expo_has_sym:
             res = imageset(Lambda(n, real_root(n, expo)), g_ys)
-            if expo.is_rational:
+            if expo.is_rational or expo.is_Float:
                 numer, denom = expo.as_numer_denom()
                 if denom % 2 == 0:
                     base_positive = solveset(base >= 0, symbol, S.Reals)
@@ -257,10 +257,9 @@ def _invert_real(f, g_ys, symbol):

                 else:
                     return _invert_real(base, res, symbol)
-            else:
-                if not base.is_positive:
-                    raise ValueError("x**w where w is irrational is not "
-                                     "defined for negative x")
+            elif expo.is_irrational:
+                g_ys_pos = g_ys & Interval(0, oo)
+                res = imageset(Lambda(n, real_root(n, expo)), g_ys_pos)
                 return _invert_real(base, res, symbol)

         if not base_has_sym:
```

```python
>>> solveset(x**1.0 + 1, x, Reals)
FiniteSet(-1.0)
```
 @oscarbenjamin I have made your most recent suggested changes on my build and

solveset(x ** 1.0 + 1, x, Reals) -> FiniteSet(-1.0)
solveset(x ** pi-E, x, Reals) -> FiniteSet(exp(1/pi))
solveset(x ** 1.5, x, Reals) -> FiniteSet(0)

all of these tests appears to return the correct answer.

Is there anything I can do to make this change happen? This is my first issue and I'd like to make a contribution somehow. I could potentially add these tests and a few more into the regression test suite.
You can send a pull request with these changes. It would be good to have some tests for cases with either a non-trivial base or a non-trivial rhs e.g. if the rhs is always/sometimes negative or something like that. I expect that there are already tests for this somewhere so these can be added to those test functions along with a comment pointing to this issue e.g.:
```
   # https://github.com/sympy/sympy/issues/21236
```
Alright, was testing out the fix with some less trivial functions.
@l-johnston , I've found a test that your last fix does not pass:
```
solveset(sin(x) ** 1.0, x, Reals) -> ImageSet(Lambda(_n, 2*_n*pi), Integers)
```
If we change 1.0 to an int we get:
```
solveset(sin(x) ** 1, x, Reals) -> Union(ImageSet(Lambda(_n, 2*_n*pi + pi), Integers), ImageSet(Lambda(_n, 2*_n*pi), Integers))
```
Which is the correct answer

Not sure what is happening here though.

Commit message: Refactor and add test

Patch to review:
diff --git a/sympy/solvers/solveset.py b/sympy/solvers/solveset.py
index 17da5ae0255d..bbfbb47a5494 100644
--- a/sympy/solvers/solveset.py
+++ b/sympy/solvers/solveset.py
@@ -239,29 +239,27 @@ def _invert_real(f, g_ys, symbol):
         expo_has_sym = expo.has(symbol)
 
         if not expo_has_sym:
-            res = imageset(Lambda(n, real_root(n, expo)), g_ys)
+            root = Lambda(n, real_root(n, expo))
             if expo.is_rational:
-                numer, denom = expo.as_numer_denom()
-                if denom % 2 == 0:
+                n, d = expo.as_numer_denom()
+                if d % 2 == 0 and n % 2 == 1 and d.is_zero is False:
+                    g_ys_pos = g_ys & Interval(0, oo)
+                    res = imageset(root, g_ys_pos)
                     base_positive = solveset(base >= 0, symbol, S.Reals)
-                    res = imageset(Lambda(n, real_root(n, expo)
-                        ), g_ys.intersect(
-                        Interval.Ropen(S.Zero, S.Infinity)))
                     _inv, _set = _invert_real(base, res, symbol)
                     return (_inv, _set.intersect(base_positive))
-
-                elif numer % 2 == 0:
-                        n = Dummy('n')
-                        neg_res = imageset(Lambda(n, -n), res)
-                        return _invert_real(base, res + neg_res, symbol)
-
-                else:
-                    return _invert_real(base, res, symbol)
-            else:
-                if not base.is_positive:
-                    raise ValueError("x**w where w is irrational is not "
-                                     "defined for negative x")
+                if n % 2 == 0 and d % 2 == 1:
+                    res = imageset(root, g_ys)
+                    n = Dummy("n")
+                    neg_res = imageset(Lambda(n, -n), res)
+                    return _invert_real(base, res + neg_res, symbol)
+                return _invert_real(base, imageset(root, g_ys), symbol)
+            if expo.is_irrational:
+                g_ys_pos = g_ys & Interval(0, oo)
+                res = imageset(root, g_ys_pos)
                 return _invert_real(base, res, symbol)
+            # indeterminate exponent, e.g. Float
+            return (f, g_ys)
 
         if not base_has_sym:
             rhs = g_ys.args[0]
diff --git a/sympy/solvers/tests/test_solveset.py b/sympy/solvers/tests/test_solveset.py
index b30d370c383c..97befcded6d7 100644
--- a/sympy/solvers/tests/test_solveset.py
+++ b/sympy/solvers/tests/test_solveset.py
@@ -92,7 +92,13 @@ def ireal(x, s=S.Reals):
     assert invert_real(x**S.Half, y, x) == (x, FiniteSet(y**2))
 
     raises(ValueError, lambda: invert_real(x, x, x))
-    raises(ValueError, lambda: invert_real(x**pi, y, x))
+
+    # issue 21236
+    assert invert_real(x**pi, y, x) == (x, FiniteSet(y**(1/pi)))
+    assert invert_real(x**pi, -E, x) == (x, EmptySet())
+    assert invert_real(x**Rational(3/2), 1000, x) == (x, FiniteSet(100))
+    assert invert_real(x**1.0, 1, x) == (x**1.0, FiniteSet(1))
+
     raises(ValueError, lambda: invert_real(S.One, y, x))
 
     assert invert_real(x**31 + x, y, x) == (x**31 + x, FiniteSet(y))
@@ -2811,3 +2817,7 @@ def test_substitution_with_infeasible_solution():
         (0, Complement(FiniteSet(p01), FiniteSet(0)), 0, p11, 0, 0, 0, 0, 0, 0, 0, -l2*p11/p01, -l3*p11/p01, l2, l3),
     )
     assert sol != nonlinsolve(system, solvefor)
+
+def test_issue_21236():
+    x, y, z = symbols("x y z")
+    assert solveset(x**y - z, x, S.Reals) == ConditionSet(x, Eq(x**y - z, 0), S.Reals)


Code review comments:
diff_hunk: @@ -2811,3 +2817,7 @@ def test_substitution_with_infeasible_solution():
         (0, Complement(FiniteSet(p01), FiniteSet(0)), 0, p11, 0, 0, 0, 0, 0, 0, 0, -l2*p11/p01, -l3*p11/p01, l2, l3),
     )
     assert sol != nonlinsolve(system, solvefor)
+
+def test_issue_21236():
+    x, y, z = symbols("x y z")
path: sympy/solvers/tests/test_solveset.py
text: @user1:
```suggestion
    x, z = symbols("x z")
    y = symbols('y', rational=True)
```

diff_hunk: @@ -2811,3 +2817,7 @@ def test_substitution_with_infeasible_solution():
         (0, Complement(FiniteSet(p01), FiniteSet(0)), 0, p11, 0, 0, 0, 0, 0, 0, 0, -l2*p11/p01, -l3*p11/p01, l2, l3),
     )
     assert sol != nonlinsolve(system, solvefor)
+
+def test_issue_21236():
+    x, y, z = symbols("x y z")
+    assert solveset(x**y - z, x, S.Reals) == ConditionSet(x, Eq(x**y - z, 0), S.Reals)
path: sympy/solvers/tests/test_solveset.py
text: @user1:
```suggestion
    assert solveset(x**y - z, x, S.Reals) == ConditionSet(x, Eq(x**y - z, 0), S.Reals)
    e1, e2 = symbols('e1 e2', even=True)
    y = e1/e2  # don't know if num or den will be odd and the other even
    assert solveset(x**y - z, x, S.Reals) == ConditionSet(x, Eq(x**y - z, 0), S.Reals)
```