Repository: sympy/sympy

Issue:
Incorrect inequality solving: False returned instead of answer
I discovered pretty weird behaviour of `solve` with rather simple inequality: it returns `False` instead of actual answer, and, which is even more strange, this behaviour fixes itself after several other inequalities are solved.

```
import sys
import sympy
from sympy import solve
from sympy.abc import x, t
f = (t / ((1 + t) ** 2))
print(f"{sys.version=}")
print(f"{sympy.__version__=}")
print(f"{solve(f.subs({t: 3 * x + 2}).diff(x) > 0, x)=}")
print(f"{solve(f.subs({t: 3 * x + 3}).diff(x) > 0, x)=}")
print(f"{solve(f.subs({t: 3 * x + 4}).diff(x) > 0, x)=}")
print(f"{solve(f.subs({t: 3 * x + 2}).diff(x) > 0, x)=}")
```

Output:

```
sys.version='3.9.0 (v3.9.0:9cf6752276, Oct  5 2020, 11:29:23) \n[Clang 6.0 (clang-600.0.57)]'
sympy.__version__='1.8.dev'
solve(f.subs({t: 3 * x + 2}).diff(x) > 0, x)=False
solve(f.subs({t: 3 * x + 3}).diff(x) > 0, x)=(-4/3 < x) & (x < -2/3)
solve(f.subs({t: 3 * x + 4}).diff(x) > 0, x)=(-5/3 < x) & (x < -1)
solve(f.subs({t: 3 * x + 2}).diff(x) > 0, x)=(-1 < x) & (x < -1/3)
```

Note that outputs of the third and last lines are different despite the fact that expressions are equal.

Expected output:

```
sys.version='3.9.0 (v3.9.0:9cf6752276, Oct  5 2020, 11:29:23) \n[Clang 6.0 (clang-600.0.57)]'
sympy.__version__='1.8.dev'
solve(f.subs({t: 3 * x + 2}).diff(x) > 0, x)=(-1 < x) & (x < -1/3)
solve(f.subs({t: 3 * x + 3}).diff(x) > 0, x)=(-4/3 < x) & (x < -2/3)
solve(f.subs({t: 3 * x + 4}).diff(x) > 0, x)=(-5/3 < x) & (x < -1)
solve(f.subs({t: 3 * x + 2}).diff(x) > 0, x)=(-1 < x) & (x < -1/3)
```

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:
* polys
  * sort real roots so that solve() returns consistent results ([#20906](https://github.com/sympy/sympy/pull/20906) by [@AaronStiff](https://github.com/AaronStiff))

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

<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 #20902

    #### Brief description of what is fixed or changed
    Amended `Poly.real_roots` so that it sorts its roots. In some cases (see linked issue), the roots were being returned in differing orders if computed at different times, leading to incorrect answers in `solve`.

    #### Other comments
    This fix uses Python's built-in `sorted` function. However, if for some reason this is not computationally efficient or internally consistent, let me know what I should use instead.

    #### 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 -->
    * polys
      * sort real roots so that solve() returns consistent results
    <!-- END RELEASE NOTES -->

</details><p>

AaronStiff: Closing due to duplicate pull request.
AaronStiff: Mistakenly closed. :)
asmeurer: I think sorted is fine. If there are computational inefficiencies those can be improved by modifying the comparison methods (e.g., perhaps a more efficient `__lt__` could be added to CRootOf).
asmeurer: Although it's also possible that the underlying algorithms should be giving roots in the correct order and there is some bug there. 
AaronStiff: @asmeurer So while @skirpichev and I both apparently decided to sort the roots independently, @skirpichev discovered a bug much deeper, in `rootisolation.is_disjoint` (https://github.com/diofant/diofant/pull/1115/commits/d8c2efe41df35d7babacbf84ec2c1cd708252775). It successfully clears up the issue, even though I have no idea why! (I'd love to know how you arrived at this @skirpichev!) If this is a preferred fix, I'll update this PR with it.
asmeurer: It would be nice to know why that is the fix, and if the same problem can't occur elsewhere. I really wish we had `hypothesis` testing in SymPy.
asmeurer: I played with hypothesis a little bit. I seems to require a lot of examples to find bugs like this (I did find https://github.com/sympy/sympy/issues/20913, however). 
oscarbenjamin: The change from @skirpichev was:
```diff
     def is_disjoint(self, other):
         """Return ``True`` if two isolation intervals are disjoint."""
-        return self.b <= other.a or other.b <= self.a
+        return self.b < other.a or other.b < self.a
```
It's easy to see why that is correct since intervals with matching start and end are not disjoint e.g.: `[1, 2], [2, 3]` have the intersection `{2}`. My guess then is that the root isolation code ends up choosing intervals that have the root at a boundary somehow. The current implementation would incorrectly report that `[1, 2]` and `[2, 3]` are disjoint which would seem to imply that they bracket distinct roots with the root in `[1, 2]` being less than the root in `[2, 3]`.


asmeurer: This is what I'm using

``py

from hypothesis import given, settings
from hypothesis.strategies import lists, integers

@given(lists(integers(), max_size=5))
@settings(max_examples=1000)
def test_real_roots_sorted(l):
    p = Poly(l, x)
    roots = p.real_roots()
    assert sorted(roots) == roots
```

But it hasn't even found this bug in master yet (you can confirm it would find it with `test_real_roots_sorted.hypothesis.inner_test([-27, -90, -108, -54, -9])`). So we might need a smarter way to generate the strategy. Or just run it for a really long time (it's also hard because it keeps running into https://github.com/sympy/sympy/issues/20913).
oscarbenjamin: That was the problem with hypothesis before. It's difficult to get it to explore examples that aren't just really slow.
asmeurer: Well https://github.com/sympy/sympy/issues/20913 is a blatent performance bug. Perhaps Poly does need to factor large integers to find the roots of larger degree polynomials, but it definitely doesn't for linear ones! You can easily limit the sizes of the integers in the `integers` strategy, though, so I can try that. 
oscarbenjamin: There is a more direct way to find the poly that causes this rather than using hypothesis:
```python
diff --git a/sympy/polys/polytools.py b/sympy/polys/polytools.py
index 39670bf874..bcd18c87f7 100644
--- a/sympy/polys/polytools.py
+++ b/sympy/polys/polytools.py
@@ -3569,6 +3569,7 @@ def real_roots(f, multiple=True, radicals=True):
         [CRootOf(x**3 + x + 1, 0)]
 
         """
+        print(f)
         reals = sympy.polys.rootoftools.CRootOf.real_roots(f, radicals=radicals)
 
         if multiple:
```
That shows that the example with the problem is:
```python
In [11]: Poly(-27*x**4 - 90*x**3 - 108*x**2 - 54*x - 9, x, domain='ZZ').real_roots()
Out[11]: [-1/3, -1, -1, -1]
```
If we print out the intervals that are tested in `is_disjoint` we get: `(-1, 0) (-1, -1)`. Here the `(-1, -1)` interval represents the root `-1` which is known exactly as a rational factor. For some reason the `-1/3` root is considered to be bracketed by `(-1, 0)`. Since the start point of one is equal to the end of the other `is_disjoint` concludes that they are disjoint and I guess it then just sorts by the start points potentially leading to non-deterministic behaviour since the start points are equal for the two intervals.

I'm not sure whether the isolation algorithms treat intervals like `(-1, 0)` as closed or open. The printed form of the interval uses round brackets suggesting an open interval. Obviously `(-1, -1)` as a bracket for `-1` would need to be a closed interval.
asmeurer: I found an example that sometimes gives the wrong order and sometimes gives the right order:

```py
>>> Poly([887, -691, 0, 0, 0], x).real_roots()
[691/887, 0, 0, 0]
```

I think this branch might fix it, but it's flaky so it's hard to tell for sure. 

I've been using this by the way, which generates a product of degree <= 5 polynomials with coefficients between -1000 and 1000:

```py
from hypothesis import given, settings
from hypothesis.strategies import lists, integers

@given(lists(lists(integers(min_value=-2**10, max_value=2**10), max_size=5)))
@settings(max_examples=100000, deadline=None)
def test_real_roots_sorted(L):
    P = [Poly(l, x) for l in L]
    p = Poly(prod(P), x)
    roots = p.real_roots()
    assert sorted(roots) == roots

test_real_roots_sorted()
```
asmeurer: > There is a more direct way to find the poly that causes this rather than using hypothesis:

The idea behind hypothesis is to see if this is a complete fix the problem, or if there are other issues that could result in roots not being sorted. Hypothesis is not a debugger. The point isn't to find out why the bug is happening, it's to find out if there are bugs in the first place. 

But since I'm having a hard time even finding the bug in master, it's not so promising that I would find other issues as well. Probably it would if ran for long enough. Simply having a test in the test suite that used hypothesis to generate a random polynomial and check if its roots were sorted would turn up any bug on CI eventually. I think it's also possible to use some tricks to share the example database between CI runs so that failures don't "disappear" in some random run.
oscarbenjamin: > The idea behind hypothesis is to see if this is a complete fix the problem, or if there are other issues that could result in roots not being sorted.

I decided hypothesis wasn't so useful when I realised that we already have 3000 open issues: we don't have a problem finding bugs!

I generally find that it's not hard to poke any part of the sympy codebase and find bugs. Fixing them is the hard part.

Here you've devised a hypothesis test that is specifically about the sorting of real roots. You have only devised that test though because someone already reported a bug resulting from this and @AaronStiff used the debugger to identify the problem and find an example polynomial.
asmeurer: > I decided hypothesis wasn't so useful when I realised that we already have 3000 open issues: we don't have a problem finding bugs!

This is the sort of thing that to me shows that we need hypothesis. Badly. We wouldn't have so many bugs if our tests used hypothesis, because hypothesis prevents you from having latent bugs. In a sense, we are using our users as our test fuzzer. They eventually find some bug, and (sometimes) report them back to us. With hypothesis, these sorts of bugs would be found before the code was ever released, ideally before it is ever merged or even committed. The alternative is the code author coming up with "good enough" test cases to cover all the corner cases, and reviewers noticing any bugs from reading through the code. I'm not saying either of those are bad or that we shouldn't do them, but they clearly have limitations. 

It's hard to have confidence that SymPy is correct when there are so many bugs found all the time. What fraction of bugs that are found are reported? And what fraction of bugs that exist are found? A single hour of hypothesis test runs can run more examples on a function than the entirety of all SymPy users ever. And unlike SymPy users, it will always recognize that there is an error (assuming the test is written well), and it will always report it.

Hypothesis is very rigorous, which can make it hard to use. It will also find bugs that are unrelated to what you are actually testing, like the performance issue I found. Ultimately this is a good thing, because all bugs are bugs. But it does make it a challenge to use it in practice.

> Here you've devised a hypothesis test that is specifically about the sorting of real roots. You have only devised that test though because someone already reported a bug resulting from this and @AaronStiff used the debugger to identify the problem and find an example polynomial.

Sure, but the next obvious question is, does this really fix the issue? Or are there still other underlying problems that make it so that roots can be unsorted. It's impossible to answer this question except by doing something like a hypothesis test. 

Anyway, writing a new test in response to a bug fix is pretty standard. That's the whole point of regression testing. If we actually used hypothesis, virtually every test would use it. Anything that asserts that some higher level fact about a function is true, like "real_roots() produces a result that is mathematically sorted" can be tested with hypothesis. That's a significant fraction of the tests in SymPy, due to its nature.
asmeurer: Hypothesis found another example in master, `(Poly(x + 1)**2*Poly(-5*x - 1)).real_roots()` (which is also flaky). I think it's fixed here too. I am now running this on this branch to see if it finds anything. I'm guessing it won't, but I'll let you know if it turns anything up. 
oscarbenjamin: > It's hard to have confidence that SymPy is correct when there are so many bugs found all the time.

Define "correct". As I said there are 3000 open issues and many of them clearly show that in some cases sympy is not correct. Inner procedures like this are usually more reliable than but there are still plenty of problems listed.

> What fraction of bugs that are found are reported?

More than we can fix or even find the time to close because they're already fixed?

> The alternative is the code author coming up with "good enough" test cases to cover all the corner cases, and reviewers noticing any bugs from reading through the code.

This is the underlying problem: much of the codebase is of low quality. You're referring to the quality of the tests but the principal problem is the quality of the code being tested. High quality code does not need to be tested so exhaustively. I think that sympy tries to use an ever-expanding regression test suite as an alternative to having well-defined, well-documented, well-maintained, high-quality code in the first place.

What we're short on when it comes to fixing issues like this is good solid contributions from careful contributors like @AaronStiff who investigate and debug the problem, identify a fix and apply careful judgement to consider how appropriate the fix is and communicate their reasoning clearly.

What we're not short on is bug reports!
AaronStiff: > Hypothesis found another example in master, `(Poly(x + 1)**2*Poly(-5*x - 1)).real_roots()` (which is also flaky). I think it's fixed here too. I am now running this on this branch to see if it finds anything. I'm guessing it won't, but I'll let you know if it turns anything up.

@asmeurer Did you want me to update this PR with the `is_disjoint` fix, or is it the `sorted` amendment that you're testing?
asmeurer: Oh, I had thought that you already pushed that fix here. The sorted change is obviously going to always pass the test. Let me try running against that commit. I would suggest pushing it here. 
AaronStiff: Oops, sorry about that. Changes have been pushed.
asmeurer: > Define "correct". As I said there are 3000 open issues and many of them clearly show that in some cases sympy is not correct. Inner procedures like this are usually more reliable than but there are still plenty of problems listed.

Correct means what you'd expect. If something is mathematically wrong, it's incorrect. There are more ambiguous instances like here where the answer isn't really wrong, but it doesn't satisfy some condition that we would want it to. But in general, if we want some property to hold for a function, it is incorrect if that property is false for some input.

> What we're short on when it comes to fixing issues like this is good solid contributions from careful contributors like @AaronStiff who investigate and debug the problem, identify a fix and apply careful judgement to consider how appropriate the fix is and communicate their reasoning clearly.

The question is, could this bug have been prevented in the first place. Someone originally wrote the code and wrote it incorrectly. If the code is right for most inputs but wrong for some corner cases, they might have been fooled by the tests that they wrote. But if they had instead written a hypothesis test, it would be much less likely. Hypothesis is really good at finding errors in corner cases.

I've written an entire library using hypothesis tests, and I can assure you that it absolutely helps with correctness. I don't know how else to convince you how useful it is, other than to suggest that you try it yourself the next time you are writing some code that is easy to check for correctness (i.e., you can easily write property based tests for it). For SymPy it's often hard because the strategies that generate the inputs can be hard to write in a good way. But there are still instances even in SymPy where it would be straightforward to add a hypothesis test, for instance, code that works purely on numerical inputs. What you'll find is that hypothesis consistently finds bugs that you would have never found yourself. It finds bugs so much that when it finally stops finding bugs and the tests pass, you gain a confidence that your code is actually right. 

Think of it another way. Imagine you could write a test that literally checked every possible input to a function. If you could do that, you would know that the function is correct. Unfortunately, except in some very restricted cases, that's impossible, because of combinatorial explosion on the input possibilities, or because the set of inputs is actually infinite. But hypothesis magically lets you do this anyway. You specify the input space, and it searches the input space in an intelligent way. Even in the cases where it isn't so intelligent, because the patterns on the inputs are nontrivial, it is at worst random testing, but with lots of nice tooling around it (it remembers failing inputs, it shrinks failing inputs to the smallest example that doesn't work). Even just random testing has a basic advantage over simple testing: every time the test is run, it becomes more robust. Normal tests don't do anything if you run them 1 time vs. 1000 times. But a random test is 1000 times better if you run it a thousand times.

In the short term, yes, adding hypothesis tests would find a lot of bugs. It would be a good idea for whoever adds them try to fix those bugs when finding them, rather than just spamming the issue tracker. But be clear, it's not *introducing* bugs. These are bugs that are already there. Bugs that users will eventually hit too. 

In the long term, the existence of a hypothesis test makes it hard for regressions to happen, because the test things in a high level way. Simple tests only test a fixed set of inputs and outputs, so any regression that isn't tested by them can slip in. Hypothesis tests all inputs, so it catches regressions even if the person who wrote the tests hadn't thought of them. The result is fewer bugs. 

I'm not claiming that it's a simple process. There are a lot of challenges involved. Hypothesis tends to find a lot of bugs, including bugs that are unrelated to what you are trying to test (like the performance issue I found here). And it's really challenging to write good strategies for SymPy expressions. But it's a worthy program in my opinion.

Anyway, I'm going to stop hijacking this issue with metadiscussions about hypothesis. We can move those to another issue if you want to discuss it more. I'm running my test above on the most recent patch that was pushed. I'll report back if anything comes up (as I said, I don't expect it to). 
oscarbenjamin: If we're going to move this discussion somewhere else then we should agree where that is. Otherwise this is the only place I can respond :)

I don't disagree that hypothesis testing *can* be very useful. The problem we have in sympy is that we already have a codebase with lots of bugs. In that situation hypothesis is less useful because we *already* know what the most important bugs are. The main thing needed is for people to take the time to debug and fix those bugs.

As you say it is possible to use something like hypothesis as part of producing the code to avoid the bugs in the first place. However if significant skill is required to produce useful hypothesis tests then it doesn't actually solve our problem: contributors with significant skill will also be better at writing ordinary tests and will produce less buggy code in the first instance.

When I tried hypothesis on sympy it quickly found some bugs that I already knew about. I asked the contributor who was suggesting to introduce hypothesis how I should move on from there because it kept showing me the same bug and the response was something like "ideally you should fix the bug" but that's not always easy in sympy. I think that hypothesis testing can be useful when you have a codebase that is close to meeting a formal spec but sympy is still some way off in most areas of interest.
asmeurer: I opened https://github.com/sympy/sympy/issues/20914, and responded to some of your points there. 

Patch:
diff --git a/doc/src/modules/evalf.rst b/doc/src/modules/evalf.rst
index 42482116f5..c1897548e3 100644
--- a/doc/src/modules/evalf.rst
+++ b/doc/src/modules/evalf.rst
@@ -392,7 +392,7 @@ Here are several more advanced examples:
         1
     ----------
     5*pi
-    ---- + 2*e
+    ---- + 2*E
      7
     >>> nsimplify(cos(atan('1/3')))
         ____
diff --git a/sympy/core/numbers.py b/sympy/core/numbers.py
index 3138c95c4b..9adb3362f5 100644
--- a/sympy/core/numbers.py
+++ b/sympy/core/numbers.py
@@ -3476,91 +3476,7 @@ def approximation_interval(self, number_cls):
 
     def _eval_power(self, expt):
         from sympy import exp
-        if global_parameters.exp_is_pow:
-            return self._eval_power_exp_is_pow(expt)
-        else:
-            return exp(expt)
-
-    def _eval_power_exp_is_pow(self, arg):
-        from ..functions.elementary.exponential import log
-        from . import Add, Mul, Pow
-        if arg.is_Number:
-            if arg is oo:
-                return oo
-            elif arg == -oo:
-                return S.Zero
-        elif isinstance(arg, log):
-            return arg.args[0]
-
-        # don't autoexpand Pow or Mul (see the issue 3351):
-        elif not arg.is_Add:
-            Ioo = I*oo
-            if arg in [Ioo, -Ioo]:
-                return nan
-
-            coeff = arg.coeff(pi*I)
-            if coeff:
-                if (2*coeff).is_integer:
-                    if coeff.is_even:
-                        return S.One
-                    elif coeff.is_odd:
-                        return S.NegativeOne
-                    elif (coeff + S.Half).is_even:
-                        return -I
-                    elif (coeff + S.Half).is_odd:
-                        return I
-                elif coeff.is_Rational:
-                    ncoeff = coeff % 2 # restrict to [0, 2pi)
-                    if ncoeff > 1: # restrict to (-pi, pi]
-                        ncoeff -= 2
-                    if ncoeff != coeff:
-                        return S.Exp1**(ncoeff*S.Pi*S.ImaginaryUnit)
-
-            # Warning: code in risch.py will be very sensitive to changes
-            # in this (see DifferentialExtension).
-
-            # look for a single log factor
-
-            coeff, terms = arg.as_coeff_Mul()
-
-            # but it can't be multiplied by oo
-            if coeff in (oo, -oo):
-                return
-
-            coeffs, log_term = [coeff], None
-            for term in Mul.make_args(terms):
-                if isinstance(term, log):
-                    if log_term is None:
-                        log_term = term.args[0]
-                    else:
-                        return
-                elif term.is_comparable:
-                    coeffs.append(term)
-                else:
-                    return
-
-            return log_term**Mul(*coeffs) if log_term else None
-        elif arg.is_Add:
-            out = []
-            add = []
-            argchanged = False
-            for a in arg.args:
-                if a is S.One:
-                    add.append(a)
-                    continue
-                newa = self**a
-                if isinstance(newa, Pow) and newa.base is self:
-                    if newa.exp != a:
-                        add.append(newa.exp)
-                        argchanged = True
-                    else:
-                        add.append(a)
-                else:
-                    out.append(newa)
-            if out or argchanged:
-                return Mul(*out)*Pow(self, Add(*add), evaluate=False)
-        elif arg.is_Matrix:
-            return arg.exp()
+        return exp(expt)
 
     def _eval_rewrite_as_sin(self, **kwargs):
         from sympy import sin
diff --git a/sympy/core/parameters.py b/sympy/core/parameters.py
index fd35cd1dc5..eff9456047 100644
--- a/sympy/core/parameters.py
+++ b/sympy/core/parameters.py
@@ -65,7 +65,7 @@ def __setattr__(self, name, value):
             clear_cache()
         return super().__setattr__(name, value)
 
-global_parameters = _global_parameters(evaluate=True, distribute=True, exp_is_pow=False)
+global_parameters = _global_parameters(evaluate=True, distribute=True)
 
 @contextmanager
 def evaluate(x):
@@ -132,30 +132,3 @@ def distribute(x):
         yield
     finally:
         global_parameters.distribute = old
-
-
-@contextmanager
-def _exp_is_pow(x):
-    """
-    Control whether `e^x` should be represented as ``exp(x)`` or a ``Pow(E, x)``.
-
-    Examples
-    ========
-
-    >>> from sympy import exp
-    >>> from sympy.abc import x
-    >>> from sympy.core.parameters import _exp_is_pow
-    >>> with _exp_is_pow(True): print(type(exp(x)))
-    <class 'sympy.core.power.Pow'>
-    >>> with _exp_is_pow(False): print(type(exp(x)))
-    exp
-    """
-    old = global_parameters.exp_is_pow
-
-    clear_cache()
-    try:
-        global_parameters.exp_is_pow = x
-        yield
-    finally:
-        clear_cache()
-        global_parameters.exp_is_pow = old
diff --git a/sympy/core/power.py b/sympy/core/power.py
index fae5590898..8d5dddf8b4 100644
--- a/sympy/core/power.py
+++ b/sympy/core/power.py
@@ -610,21 +610,11 @@ def _eval_is_integer(self):
             return False
 
     def _eval_is_extended_real(self):
-        from ..functions import arg, log, exp
-        from .mul import Mul
-
-        if self.base is S.Exp1:
-            if self.exp.is_extended_real:
-                return True
-            elif self.exp.is_imaginary:
-                return (2*S.ImaginaryUnit*self.exp/S.Pi).is_even
-
+        from sympy import arg, exp, log, Mul
         real_b = self.base.is_extended_real
         if real_b is None:
             if self.base.func == exp and self.base.args[0].is_imaginary:
                 return self.exp.is_imaginary
-            if self.base.func == Pow and self.base.base is S.Exp1 and self.base.exp.is_imaginary:
-                return self.exp.is_imaginary
             return
         real_e = self.exp.is_extended_real
         if real_e is None:
@@ -866,7 +856,7 @@ def _check(ct1, ct2, old):
                     new_l.append(Pow(self.base, expo, evaluate=False) if expo != 1 else self.base)
                     return Mul(*new_l)
 
-        if (isinstance(old, exp) or (old.is_Pow and old.base is S.Exp1)) and self.exp.is_extended_real and self.base.is_positive:
+        if isinstance(old, exp) and self.exp.is_extended_real and self.base.is_positive:
             ct1 = old.args[0].as_independent(Symbol, as_Add=False)
             ct2 = (self.exp*log(self.base)).as_independent(
                 Symbol, as_Add=False)
@@ -1257,15 +1247,7 @@ def as_real_imag(self, deep=True, **hints):
 
             rp, tp = self.func(r, self.exp), t*self.exp
 
-            return rp*cos(tp), rp*sin(tp)
-        elif self.base is S.Exp1:
-            from ..functions import exp
-            re_e, im_e = self.exp.as_real_imag()
-            if deep:
-                re_e = re_e.expand(deep, **hints)
-                im_e = im_e.expand(deep, **hints)
-            c, s = cos(im_e), sin(im_e)
-            return exp(re_e)*c, exp(re_e)*s
+            return (rp*cos(tp), rp*sin(tp))
         else:
 
             if deep:
@@ -1277,7 +1259,7 @@ def as_real_imag(self, deep=True, **hints):
                 else:
                     return (re(expanded), im(expanded))
             else:
-                return re(self), im(self)
+                return (re(self), im(self))
 
     def _eval_derivative(self, s):
         from sympy import log
@@ -1329,9 +1311,6 @@ def _eval_is_rational(self):
                     return True
             elif b.is_irrational:
                 return e.is_zero
-        if b is S.Exp1:
-            if e.is_rational and e.is_nonzero:
-                return False
 
     def _eval_is_algebraic(self):
         def _is_one(expr):
@@ -1343,18 +1322,6 @@ def _is_one(expr):
 
         if self.base.is_zero or _is_one(self.base):
             return True
-        elif self.base is S.Exp1:
-            s = self.func(*self.args)
-            if s.func == self.func:
-                if self.exp.is_nonzero:
-                    if self.exp.is_algebraic:
-                        return False
-                    elif (self.exp/S.Pi).is_rational:
-                        return False
-                    elif (self.exp/(S.ImaginaryUnit*S.Pi)).is_rational:
-                        return True
-            else:
-                return s.is_algebraic
         elif self.exp.is_rational:
             if self.base.is_algebraic is False:
                 return self.exp.is_zero
@@ -1524,29 +1491,11 @@ def _eval_nseries(self, x, n, logx, cdir=0):
         #    g has order O(x**d) where d is strictly positive.
         # 2) Then b**e = (f**e)*((1 + g)**e).
         #    (1 + g)**e is computed using binomial series.
-        from sympy import im, I, ceiling, polygamma, logcombine, EulerGamma, nan, zoo, factorial, ff, PoleError, O, powdenest, Wild
+        from sympy import im, I, ceiling, polygamma, limit, logcombine, EulerGamma, exp, nan, zoo, log, factorial, ff, PoleError, O, powdenest, Wild
         from itertools import product
-        from ..functions import exp, log
-        from ..series import Order, limit
-        from ..simplify import powsimp
-        if self.base is S.Exp1:
-            e_series = self.exp.nseries(x, n=n, logx=logx)
-            if e_series.is_Order:
-                return 1 + e_series
-            e0 = limit(e_series.removeO(), x, 0)
-            if e0 in (-S.NegativeInfinity, S.Infinity):
-                return self
-            t = e_series - e0
-            exp_series = term = exp(e0)
-            # series of exp(e0 + t) in t
-            for i in range(1, n):
-                term *= t/i
-                term = term.nseries(x, n=n, logx=logx)
-                exp_series += term
-            exp_series += Order(t**n, x)
-            return powsimp(exp_series, deep=True, combine='exp')
         self = powdenest(self, force=True).trigsimp()
         b, e = self.as_base_exp()
+
         if e.has(S.Infinity, S.NegativeInfinity, S.ComplexInfinity, S.NaN):
             raise PoleError()
 
@@ -1657,87 +1606,25 @@ def mul(d1, d2):
         return res
 
     def _eval_as_leading_term(self, x, cdir=0):
-        from ..series import Order
         from sympy import exp, I, im, log
         e = self.exp
         b = self.base
-        if self.base is S.Exp1:
-            arg = self.exp
-            if arg.is_Add:
-                return Mul(*[(S.Exp1**f).as_leading_term(x) for f in arg.args])
-            arg_1 = arg.as_leading_term(x)
-            if Order(x, x).contains(arg_1):
-                return S.One
-            if Order(1, x).contains(arg_1):
-                return S.Exp1**arg_1
-            ####################################################
-            # The correct result here should be 'None'.        #
-            # Indeed arg in not bounded as x tends to 0.       #
-            # Consequently the series expansion does not admit #
-            # the leading term.                                #
-            # For compatibility reasons, the return value here #
-            # is the original function, i.e. exp(arg),         #
-            # instead of None.                                 #
-            ####################################################
-            return S.Exp1**arg
-        elif e.has(x):
+        if e.has(x):
             return exp(e * log(b)).as_leading_term(x, cdir=cdir)
-        else:
-            f = b.as_leading_term(x, cdir=cdir)
-            if (not e.is_integer and f.is_constant() and f.is_real
-                and f.is_negative and im((b - f).dir(x, cdir)) < 0):
-                return self.func(f, e) * exp(-2 * e * S.Pi * I)
-            return self.func(f, e)
+        f = b.as_leading_term(x, cdir=cdir)
+        if (not e.is_integer and f.is_constant() and f.is_real
+            and f.is_negative and im((b - f).dir(x, cdir)) < 0):
+            return self.func(f, e)*exp(-2*e*S.Pi*I)
+        return self.func(f, e)
 
     @cacheit
     def _taylor_term(self, n, x, *previous_terms): # of (1 + x)**e
         from sympy import binomial
         return binomial(self.exp, n) * self.func(x, n)
 
-    def taylor_term(self, n, x, *previous_terms):
-        if self.base is not S.Exp1:
-            return super().taylor_term(n, x, *previous_terms)
-        from sympy import sympify, factorial
-        if n < 0:
-            return S.Zero
-        if n == 0:
-            return S.One
-        x = sympify(x)
-        if previous_terms:
-            p = previous_terms[-1]
-            if p is not None:
-                return p * x / n
-        return x**n/factorial(n)
-
     def _sage_(self):
         return self.args[0]._sage_()**self.args[1]._sage_()
 
-    def _eval_rewrite_as_sin(self, base, exp):
-        from ..functions import sin
-        if self.base is S.Exp1:
-            return sin(S.ImaginaryUnit*self.exp + S.Pi/2) - S.ImaginaryUnit*sin(S.ImaginaryUnit*self.exp)
-
-    def _eval_rewrite_as_cos(self, base, exp):
-        from ..functions import cos
-        if self.base is S.Exp1:
-            return cos(S.ImaginaryUnit*self.exp) + S.ImaginaryUnit*cos(S.ImaginaryUnit*self.exp + S.Pi/2)
-
-    def _eval_rewrite_as_tanh(self, base, exp):
-        from ..functions import tanh
-        if self.base is S.Exp1:
-            return (1 + tanh(self.exp/2))/(1 - tanh(self.exp/2))
-
-    def _eval_rewrite_as_sqrt(self, base, exp, **kwargs):
-        from sympy.functions.elementary.trigonometric import sin, cos
-        if base is not S.Exp1:
-            return None
-        if exp.is_Mul:
-            coeff = exp.coeff(S.Pi * S.ImaginaryUnit)
-            if coeff and coeff.is_number:
-                cosine, sine = cos(S.Pi*coeff), sin(S.Pi*coeff)
-                if not isinstance(cosine, cos) and not isinstance (sine, sin):
-                    return cosine + S.ImaginaryUnit*sine
-
     def as_content_primitive(self, radical=False, clear=True):
         """Return the tuple (R, self/R) where R is the positive Rational
         extracted from self.
diff --git a/sympy/functions/elementary/exponential.py b/sympy/functions/elementary/exponential.py
index 77f6336847..c4b0231395 100644
--- a/sympy/functions/elementary/exponential.py
+++ b/sympy/functions/elementary/exponential.py
@@ -6,7 +6,6 @@
 from sympy.core.logic import fuzzy_and, fuzzy_not, fuzzy_or
 from sympy.core.mul import Mul
 from sympy.core.numbers import Integer, Rational
-from sympy.core.parameters import global_parameters
 from sympy.core.power import Pow
 from sympy.core.singleton import S
 from sympy.core.symbol import Wild, Dummy
@@ -261,8 +260,6 @@ def eval(cls, arg):
         from sympy.sets.setexpr import SetExpr
         from sympy.matrices.matrices import MatrixBase
         from sympy import im, logcombine, re
-        if global_parameters.exp_is_pow:
-            return Pow(S.Exp1, arg)
         if arg.is_Number:
             if arg is S.NaN:
                 return S.NaN
@@ -672,8 +669,6 @@ def eval(cls, arg, base=None):
             elif arg.is_Rational and arg.p == 1:
                 return -cls(arg.q)
 
-        if arg.is_Pow and arg.base is S.Exp1 and arg.exp.is_extended_real:
-            return arg.exp
         I = S.ImaginaryUnit
         if isinstance(arg, exp) and arg.args[0].is_extended_real:
             return arg.args[0]
diff --git a/sympy/polys/rootisolation.py b/sympy/polys/rootisolation.py
index c2b44ef0a8..f5035dd5e2 100644
--- a/sympy/polys/rootisolation.py
+++ b/sympy/polys/rootisolation.py
@@ -1742,9 +1742,9 @@ def __repr__(self):
     def is_disjoint(self, other):
         """Return ``True`` if two isolation intervals are disjoint. """
         if isinstance(other, RealInterval):
-            return (self.b <= other.a or other.b <= self.a)
+            return (self.b < other.a or other.b < self.a)
         assert isinstance(other, ComplexInterval)
-        return (self.b <= other.ax or other.bx <= self.a
+        return (self.b < other.ax or other.bx < self.a
             or other.ay*other.by > 0)
 
     def _inner_refine(self):
@@ -2008,10 +2008,10 @@ def is_disjoint(self, other):
             return other.is_disjoint(self)
         if self.conj != other.conj:  # above and below real axis
             return True
-        re_distinct = (self.bx <= other.ax or other.bx <= self.ax)
+        re_distinct = (self.bx < other.ax or other.bx < self.ax)
         if re_distinct:
             return True
-        im_distinct = (self.by <= other.ay or other.by <= self.ay)
+        im_distinct = (self.by < other.ay or other.by < self.ay)
         return im_distinct
 
     def _inner_refine(self):
diff --git a/sympy/printing/latex.py b/sympy/printing/latex.py
index bd722de1dd..673acd9cec 100644
--- a/sympy/printing/latex.py
+++ b/sympy/printing/latex.py
@@ -1109,9 +1109,6 @@ def _print_ExpBase(self, expr, exp=None):
         tex = r"e^{%s}" % self._print(expr.args[0])
         return self._do_exponent(tex, exp)
 
-    def _print_Exp1(self, expr, exp=None):
-        return "e"
-
     def _print_elliptic_k(self, expr, exp=None):
         tex = r"\left(%s\right)" % self._print(expr.args[0])
         if exp is not None:
diff --git a/sympy/printing/pretty/pretty.py b/sympy/printing/pretty/pretty.py
index d41513e564..4f8f85215e 100644
--- a/sympy/printing/pretty/pretty.py
+++ b/sympy/printing/pretty/pretty.py
@@ -1473,9 +1473,6 @@ def _print_ExpBase(self, e):
         base = prettyForm(pretty_atom('Exp1', 'e'))
         return base ** self._print(e.args[0])
 
-    def _print_Exp1(self, e):
-        return prettyForm(pretty_atom('Exp1', 'e'))
-
     def _print_Function(self, e, sort=False, func_name=None):
         # optional argument func_name for supplying custom names
         # XXX works only for applied functions


Test patch:
diff --git a/sympy/core/tests/test_numbers.py b/sympy/core/tests/test_numbers.py
index 7e2ed5d94a..9faed1d131 100644
--- a/sympy/core/tests/test_numbers.py
+++ b/sympy/core/tests/test_numbers.py
@@ -13,7 +13,7 @@
 from sympy.polys.domains.groundtypes import PythonRational
 from sympy.utilities.decorator import conserve_mpmath_dps
 from sympy.utilities.iterables import permutations
-from sympy.testing.pytest import XFAIL, raises, _both_exp_pow
+from sympy.testing.pytest import XFAIL, raises
 
 from mpmath import mpf
 from mpmath.rational import mpq
@@ -854,7 +854,6 @@ def test_Div_By_Zero():
     assert -1/Float(0) is zoo
 
 
-@_both_exp_pow
 def test_Infinity_inequations():
     assert oo > pi
     assert not (oo < pi)
diff --git a/sympy/functions/elementary/tests/test_exponential.py b/sympy/functions/elementary/tests/test_exponential.py
index ce9aa4fc11..d2718c70fd 100644
--- a/sympy/functions/elementary/tests/test_exponential.py
+++ b/sympy/functions/elementary/tests/test_exponential.py
@@ -4,21 +4,14 @@
     adjoint, conjugate, transpose, O, refine,
     sin, cos, sinh, cosh, tanh, exp_polar, re, simplify,
     AccumBounds, MatrixSymbol, Pow, gcd, Sum, Product)
-from sympy.core.parameters import global_parameters
 from sympy.functions.elementary.exponential import match_real_imag
 from sympy.abc import x, y, z
 from sympy.core.expr import unchanged
 from sympy.core.function import ArgumentIndexError
-from sympy.testing.pytest import raises, XFAIL, _both_exp_pow
+from sympy.testing.pytest import raises, XFAIL
 
 
-@_both_exp_pow
 def test_exp_values():
-    if global_parameters.exp_is_pow:
-        assert type(exp(x)) is Pow
-    else:
-        assert type(exp(x)) is exp
-
     k = Symbol('k', integer=True)
 
     assert exp(nan) is nan
@@ -57,7 +50,6 @@ def test_exp_values():
     assert exp(oo, evaluate=False).is_finite is False
 
 
-@_both_exp_pow
 def test_exp_period():
     assert exp(I*pi*Rational(9, 4)) == exp(I*pi/4)
     assert exp(I*pi*Rational(46, 18)) == exp(I*pi*Rational(5, 9))
@@ -77,22 +69,18 @@ def test_exp_period():
     assert exp((-1 + 4*n)*I*pi/2) == -I
 
 
-@_both_exp_pow
 def test_exp_log():
     x = Symbol("x", real=True)
     assert log(exp(x)) == x
     assert exp(log(x)) == x
-
-    if not global_parameters.exp_is_pow:
-        assert log(x).inverse() == exp
-        assert exp(x).inverse() == log
+    assert log(x).inverse() == exp
+    assert exp(x).inverse() == log
 
     y = Symbol("y", polar=True)
     assert log(exp_polar(z)) == z
     assert exp(log(y)) == y
 
 
-@_both_exp_pow
 def test_exp_expand():
     e = exp(log(Rational(2))*(1 + x) - log(Rational(2))*x)
     assert e.expand() == 2
@@ -100,7 +88,6 @@ def test_exp_expand():
     assert exp(x + y).expand() == exp(x)*exp(y)
 
 
-@_both_exp_pow
 def test_exp__as_base_exp():
     assert exp(x).as_base_exp() == (E, x)
     assert exp(2*x).as_base_exp() == (E, 2*x)
@@ -157,7 +144,6 @@ def test_exp_transpose():
     assert transpose(exp(x)) == exp(transpose(x))
 
 
-@_both_exp_pow
 def test_exp_rewrite():
     from sympy.concrete.summations import Sum
     assert exp(x).rewrite(sin) == sinh(x) + cosh(x)
@@ -168,10 +154,9 @@ def test_exp_rewrite():
     assert exp(x).rewrite(tanh) == (1 + tanh(x/2))/(1 - tanh(x/2))
     assert exp(pi*I/4).rewrite(sqrt) == sqrt(2)/2 + sqrt(2)*I/2
     assert exp(pi*I/3).rewrite(sqrt) == S.Half + sqrt(3)*I/2
-    if not global_parameters.exp_is_pow:
-        assert exp(x*log(y)).rewrite(Pow) == y**x
-        assert exp(log(x)*log(y)).rewrite(Pow) in [x**log(y), y**log(x)]
-        assert exp(log(log(x))*y).rewrite(Pow) == log(x)**y
+    assert exp(x*log(y)).rewrite(Pow) == y**x
+    assert exp(log(x)*log(y)).rewrite(Pow) in [x**log(y), y**log(x)]
+    assert exp(log(log(x))*y).rewrite(Pow) == log(x)**y
 
     n = Symbol('n', integer=True)
 
@@ -181,12 +166,10 @@ def test_exp_rewrite():
             == 4/(3 - sqrt(3)*I))
 
 
-@_both_exp_pow
 def test_exp_leading_term():
     assert exp(x).as_leading_term(x) == 1
     assert exp(2 + x).as_leading_term(x) == exp(2)
     assert exp((2*x + 3) / (x+1)).as_leading_term(x) == exp(3)
-
     # The following tests are commented, since now SymPy returns the
     # original function when the leading term in the series expansion does
     # not exist.
@@ -195,7 +178,6 @@ def test_exp_leading_term():
     # raises(NotImplementedError, lambda: exp(x + 1/x).as_leading_term(x))
 
 
-@_both_exp_pow
 def test_exp_taylor_term():
     x = symbols('x')
     assert exp(x).taylor_term(1, x) == x
@@ -590,7 +572,6 @@ def test_exp_expand_NC():
     assert exp(x + y + z).expand() == exp(x)*exp(y)*exp(z)
 
 
-@_both_exp_pow
 def test_as_numer_denom():
     n = symbols('n', negative=True)
     assert exp(x).as_numer_denom() == (exp(x), 1)
diff --git a/sympy/ntheory/primetest.py b/sympy/ntheory/primetest.py
index 797fbd32f9..106bc749cb 100644
--- a/sympy/ntheory/primetest.py
+++ b/sympy/ntheory/primetest.py
@@ -653,8 +653,8 @@ def is_gaussian_prime(num):
     from sympy import sympify
     num = sympify(num)
     a, b = num.as_real_imag()
-    a = as_int(a, strict=False)
-    b = as_int(b, strict=False)
+    a = as_int(a)
+    b = as_int(b)
     if a == 0:
         b = abs(b)
         return isprime(b) and b % 4 == 3
diff --git a/sympy/polys/tests/test_rootoftools.py b/sympy/polys/tests/test_rootoftools.py
index 5ef5fd96fb..f527e8b842 100644
--- a/sympy/polys/tests/test_rootoftools.py
+++ b/sympy/polys/tests/test_rootoftools.py
@@ -284,6 +284,10 @@ def test_CRootOf_real_roots():
     assert Poly(x**5 + x + 1).real_roots(radicals=False) == [rootof(
         x**3 - x**2 + 1, 0)]
 
+    # https://github.com/sympy/sympy/issues/20902
+    p = Poly(-3*x**4 - 10*x**3 - 12*x**2 - 6*x - 1, x, domain='ZZ')
+    assert CRootOf.real_roots(p) == [S(-1), S(-1), S(-1), S(-1)/3]
+
 
 def test_CRootOf_all_roots():
     assert Poly(x**5 + x + 1).all_roots() == [
diff --git a/sympy/printing/pretty/tests/test_pretty.py b/sympy/printing/pretty/tests/test_pretty.py
index 42fae775d3..9b65539b81 100644
--- a/sympy/printing/pretty/tests/test_pretty.py
+++ b/sympy/printing/pretty/tests/test_pretty.py
@@ -42,7 +42,7 @@
 from sympy.tensor.tensor import (TensorIndexType, tensor_indices, TensorHead,
                                  TensorElement, tensor_heads)
 
-from sympy.testing.pytest import raises, _both_exp_pow
+from sympy.testing.pytest import raises
 
 from sympy.vector import CoordSys3D, Gradient, Curl, Divergence, Dot, Cross, Laplacian
 
@@ -1298,7 +1298,6 @@ def test_pretty_rational():
     assert upretty(expr) == ucode_str
 
 
-@_both_exp_pow
 def test_pretty_functions():
     """Tests for Abs, conjugate, exp, function braces, and factorial."""
     expr = (2*x + exp(x))
@@ -1321,14 +1320,9 @@ def test_pretty_functions():
 """\
  x     \n\
 ℯ + 2⋅x\
-"""
-    ucode_str_3 = \
-"""\
- x      \n\
-ℯ  + 2⋅x\
 """
     assert pretty(expr) in [ascii_str_1, ascii_str_2]
-    assert upretty(expr) in [ucode_str_1, ucode_str_2, ucode_str_3]
+    assert upretty(expr) in [ucode_str_1, ucode_str_2]
 
     expr = Abs(x)
     ascii_str = \
diff --git a/sympy/printing/tests/test_latex.py b/sympy/printing/tests/test_latex.py
index c3f183d4a8..764b18a461 100644
--- a/sympy/printing/tests/test_latex.py
+++ b/sympy/printing/tests/test_latex.py
@@ -33,7 +33,7 @@
                                 MutableSparseNDimArray,
                                 MutableDenseNDimArray,
                                 tensorproduct)
-from sympy.testing.pytest import XFAIL, raises, _both_exp_pow
+from sympy.testing.pytest import XFAIL, raises
 from sympy.functions import DiracDelta, Heaviside, KroneckerDelta, LeviCivita
 from sympy.functions.combinatorial.numbers import bernoulli, bell, lucas, \
     fibonacci, tribonacci
@@ -339,7 +339,6 @@ def test_latex_symbols_failing():
         r"{\mathrm{mass}}^{3} \cdot {\mathrm{volume}}^{3}"
 
 
-@_both_exp_pow
 def test_latex_functions():
     assert latex(exp(x)) == "e^{x}"
     assert latex(exp(1) + exp(2)) == "e + e^{2}"
diff --git a/sympy/printing/tests/test_str.py b/sympy/printing/tests/test_str.py
index 947a389b4e..77be2df38f 100644
--- a/sympy/printing/tests/test_str.py
+++ b/sympy/printing/tests/test_str.py
@@ -7,7 +7,6 @@
     AccumBounds, UnevaluatedExpr, Eq, Ne, Quaternion, Subs, MatrixSymbol, MatrixSlice,
     Q)
 from sympy.core import Expr, Mul
-from sympy.core.parameters import _exp_is_pow
 from sympy.external import import_module
 from sympy.physics.control.lti import TransferFunction, Series, Parallel, Feedback
 from sympy.physics.units import second, joule
@@ -101,8 +100,6 @@ def test_EulerGamma():
 
 def test_Exp():
     assert str(E) == "E"
-    with _exp_is_pow(True):
-        assert str(exp(x)) == "E**x"
 
 
 def test_factorial():
diff --git a/sympy/solvers/tests/test_solvers.py b/sympy/solvers/tests/test_solvers.py
index be7a567d29..7545daffcc 100644
--- a/sympy/solvers/tests/test_solvers.py
+++ b/sympy/solvers/tests/test_solvers.py
@@ -2300,3 +2300,10 @@ def test_issue_20747():
             / (1 - exp(c0/(DBH*c3 + THT*c4 + c2))))
     sol = [THT*term**(1/c1) - term**(1/c1) + 1]
     assert solve(eq, HT) == sol
+
+def test_issue_20902():
+    f = (t / ((1 + t) ** 2))
+    assert solve(f.subs({t: 3 * x + 2}).diff(x) > 0, x) == (S(-1) < x) & (x < S(-1)/3)
+    assert solve(f.subs({t: 3 * x + 3}).diff(x) > 0, x) == (S(-4)/3 < x) & (x < S(-2)/3)
+    assert solve(f.subs({t: 3 * x + 4}).diff(x) > 0, x) == (S(-5)/3 < x) & (x < S(-1))
+    assert solve(f.subs({t: 3 * x + 2}).diff(x) > 0, x) == (S(-1) < x) & (x < S(-1)/3)
diff --git a/sympy/testing/pytest.py b/sympy/testing/pytest.py
index 216d556cd1..84f0436ae0 100644
--- a/sympy/testing/pytest.py
+++ b/sympy/testing/pytest.py
@@ -214,27 +214,6 @@ def warns(warningcls, *, match=''):
             raise Failed(msg)
 
 
-def _both_exp_pow(func):
-    """
-    Decorator used to run the test twice: the first time `e^x` is represented
-    as ``Pow(E, x)``, the second time as ``exp(x)`` (exponential object is not
-    a power).
-
-    This is a temporary trick helping to manage the elimination of the class
-    ``exp`` in favor of a replacement by ``Pow(E, ...)``.
-    """
-    from sympy.core.parameters import _exp_is_pow
-
-    def func_wrap():
-        with _exp_is_pow(True):
-            func()
-        with _exp_is_pow(False):
-            func()
-
-    wrapper = functools.update_wrapper(func_wrap, func)
-    return wrapper
-
-
 @contextlib.contextmanager
 def warns_deprecated_sympy():
     '''Shorthand for ``warns(SymPyDeprecationWarning)``