fix(core): don't simplify sqrt(x**2) if x is not positive

[oscarbenjamin]

References to other Issues or PRs

This is an old bug dating back to #7638.

This was discovered during #22993 (#22993 (comment)).

Brief description of what is fixed or changed

This PR makes Pow._eval_power more conservative about simplifying (z**a)**b to z**(a*b) in the case of square rooting a square.

On master we have:

In [2]: p = (-1 + I)**2

In [3]: p
Out[3]: 
        2
(-1 + ⅈ) 

In [4]: sqrt(p)
Out[4]: -1 + ⅈ

In [5]: sqrt(expand(p))
Out[5]: 1 - ⅈ

It is In [4] that is incorrect as it does not give the principal root:

In [6]: arg(sqrt(p))
Out[6]: 
3⋅π
───
 4 

In [7]: arg(sqrt(expand(p)))
Out[7]: 
-π 
───
 4 

The change here just removes a special case from Pow._eval_power:

sympy/sympy/core/power.py

Lines 473 to 474 in e614ccb

	elif fuzzy_not(im(b).is_zero) and abs(e) == 2:
	s = 1 # floor = 0

The logic can the fall back to the more general code that can correctly handle this case:

sympy/sympy/core/power.py

Lines 475 to 481 in e614ccb

	elif _half(other):
	s = exp(2*S.Pi*S.ImaginaryUnit*other*floor(
	S.Half - e*arg(b)/(2*S.Pi)))
	if s.is_extended_real and _n2(sign(s) - s) == 0:
	s = sign(s)
	else:
	s = None

See the discussion in #7638 and certik/theoretical-physics#23 to understand the formula that is being used for simplification in this case.

Other comments

I expect this will break a few tests so I'm just putting it up first to see how it fares.

Release Notes

NO ENTRY

[sympy-bot]

✅

Hi, I am the SymPy bot (v167). I'm here to help you write a release notes entry. Please read the guide on how to write release notes.

• No release notes entry will be added for this pull request.

Click here to see the pull request description that was parsed.

<!-- 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. -->

This is an old bug dating back to #7638.

This was discovered during #22993 (https://github.com/sympy/sympy/pull/22993#issuecomment-1159711004).

#### Brief description of what is fixed or changed

This PR makes `Pow._eval_power` more conservative about simplifying `(z**a)**b` to `z**(a*b)` in the case of square rooting a square.

On master we have:
```python
In [2]: p = (-1 + I)**2

In [3]: p
Out[3]: 
        2
(-1 + ⅈ) 

In [4]: sqrt(p)
Out[4]: -1 + ⅈ

In [5]: sqrt(expand(p))
Out[5]: 1 - ⅈ
```
It is `In [4]` that is incorrect as it does not give the principal root:
```python
In [6]: arg(sqrt(p))
Out[6]: 
3⋅π
───
 4 

In [7]: arg(sqrt(expand(p)))
Out[7]: 
-π 
───
 4 
```
The change here just removes a special case from `Pow._eval_power`:
https://github.com/sympy/sympy/blob/e614ccb66ad8d75e817de17f419f4b8e4469fd1f/sympy/core/power.py#L473-L474
The logic can the fall back to the more general code that can correctly handle this case:
https://github.com/sympy/sympy/blob/e614ccb66ad8d75e817de17f419f4b8e4469fd1f/sympy/core/power.py#L475-L481

See the discussion in #7638 and https://github.com/certik/theoretical-physics/issues/23 to understand the formula that is being used for simplification in this case.

#### Other comments

I expect this will break a few tests so I'm just putting it up first to see how it fares.

#### 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 -->
NO ENTRY
<!-- END RELEASE NOTES -->

[github-actions]

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)

Significantly changed benchmark results (master vs previous release)

       before           after         ratio
     [77f1d79c]       [3e51f99f]
     <sympy-1.10.1^0>                 
+         129±2ms          233±1ms     1.80  sum.TimeSum.time_doit

Full benchmark results can be found as artifacts in GitHub Actions
(click on checks at the top of the PR).

[oscarbenjamin]

Okay, tests have passed. I expected this change to be more problematic...

I need to take a little time to work through the logic of Pow._eval_power more carefully to see if there are other problems and at least add some comments to it that explain what is going on.

The two test lines I've commented out also need more attention. I think they might be mathematically correct as tests but they test an evaluation that does not necessarily need to happen. I want to make sure I understand them but I don't currently see a need to try to fix them unless it can be done in a simple (and correct!) way.

[oscarbenjamin]

CC @smichr

[oscarbenjamin]

These tests are clearly wrong. The left hand sides are obviously equal but the right hand sides are not.

[eagleoflqj]

The first case seems correct. Is there a way to keep that and correct the 2nd case?

[oscarbenjamin]

I think it can be done like this:

diff --git a/sympy/core/power.py b/sympy/core/power.py
index 018ae5d..e46715c 100644
--- a/sympy/core/power.py
+++ b/sympy/core/power.py
@@ -470,6 +470,11 @@ def _n2(e):
                     s = 1  # floor = 0
                 elif re(b).is_extended_nonnegative and (abs(e) < 2) == True:
                     s = 1  # floor = 0
+                elif fuzzy_not(im(b).is_zero) and abs(e) == 2:
+                    if re(b).is_extended_nonnegative:
+                        s = 1
+                    elif re(b).is_negative:
+                        s = -1
                 elif _half(other):
                     s = exp(2*S.Pi*S.ImaginaryUnit*other*floor(
                         S.Half - e*arg(b)/(2*S.Pi)))

With that we have:

In [2]: p = symbols('p', positive=True)

In [3]: sqrt((p**2*I - p)**2)
Out[3]: 
     2    
- ⅈ⋅p  + p

In [4]: sqrt((p - p**2*I)**2)
Out[4]: 
     2    
- ⅈ⋅p  + p

I'm not sure that it's really a good idea though. This routine already has too many cases and is very complicated for ordinary evaluation. Really it would be better to handle this sort of thing in something like powsimp.

[eagleoflqj]

I don't have strong opinion. We can apply this when others complain.

[oscarbenjamin]

I think that the other cases in Pow._eval_power are correct. The fix here is just to remove the erroneous case. I've deleted two tests but they were testing incorrect behaviour.

I started looking at improving Pow._eval_power to add comments etc but then I realised I was going to end up completely rewriting it so I decided to do that in a separate PR and just keep this one as a straight-forward bugfix.

This is ready for review.

[oscargus]

Looks good to me!

[oscarbenjamin]

Thanks both