Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

expand(x**(a+b)) iff x != 0 or a & b both nonneg or nonpos #23962

Merged
merged 2 commits into from Sep 2, 2022

Conversation

smichr
Copy link
Member

@smichr smichr commented Aug 24, 2022

References to other Issues or PRs

Brief description of what is fixed or changed

fixes #23952 but required treating zoo/x as non-rational.

Other comments

This fixes the root issue (of the expansion to 0) but better handling needs to be implemented to get the original issue to be solved:

>>> import sympy as sp
>>> p, q = sp.symbols("p q", real=True, nonnegative=True)
>>> k1,k2 = sp.symbols("k1 k2", integer=True, nonnegative=True)
>>> n = sp.Symbol("n", integer=True, positive=True)
>>> expr = sp.Sum(sp.Abs(k1 - k2)*p**k1 *(1-q)**(n-k2), (k1,0,n), (k2,0,n))
>>> expr.subs(p,0).subs(q,1)
Sum(0**k1*0**(-k2 + n)*Abs(k1 - k2), (k1, 0, n), (k2, 0, n))
>>> _.doit()
...
sympy.polys.polyerrors.CoercionFailed: zoo is not in any domain

If the assumptions are omitted then

>>> S('''Sum(p**k1*(1 - q)**(-k2 + n)*Abs(k1 - k2), (k1, 0, n), (k2, 0, n))''')
Sum(p**k1*(1 - q)**(-k2 + n)*Abs(k1 - k2), (k1, 0, n), (k2, 0, n))
>>> ok=_;var('p q');ok.subs(p,0).subs(q,1)
(p, q)
Sum(0**k1*0**(-k2 + n)*Abs(k1 - k2), (k1, 0, n), (k2, 0, n))
>>> _.doit()
Sum(0**k1*0**(-k2 + n)*Abs(k1 - k2), (k1, 0, n), (k2, 0, n))

Release Notes

  • core
    • x**(a + b) will not expand unless a and b are both nonnegative or nonpositive or x is not zero
    • expressions containing any illegal character as defined in numbers now gives expr.is_rational_function()->False

@sympy-bot
Copy link

sympy-bot commented Aug 24, 2022

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.

Your release notes are in good order.

Here is what the release notes will look like:

  • core
    • x**(a + b) will not expand unless a and b are both nonnegative or nonpositive or x is not zero (#23962 by @smichr)

    • expressions containing any illegal character as defined in numbers now gives expr.is_rational_function()->False (#23962 by @smichr)

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

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


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

fixes #23952 but required treating zoo/x as non-rational.

#### Other comments

This fixes the root issue (of the expansion to 0) but better handling needs to be implemented to get the original issue to be solved:
```python
>>> import sympy as sp
>>> p, q = sp.symbols("p q", real=True, nonnegative=True)
>>> k1,k2 = sp.symbols("k1 k2", integer=True, nonnegative=True)
>>> n = sp.Symbol("n", integer=True, positive=True)
>>> expr = sp.Sum(sp.Abs(k1 - k2)*p**k1 *(1-q)**(n-k2), (k1,0,n), (k2,0,n))
>>> expr.subs(p,0).subs(q,1)
Sum(0**k1*0**(-k2 + n)*Abs(k1 - k2), (k1, 0, n), (k2, 0, n))
>>> _.doit()
...
sympy.polys.polyerrors.CoercionFailed: zoo is not in any domain
```
If the assumptions are omitted then
```python
>>> S('''Sum(p**k1*(1 - q)**(-k2 + n)*Abs(k1 - k2), (k1, 0, n), (k2, 0, n))''')
Sum(p**k1*(1 - q)**(-k2 + n)*Abs(k1 - k2), (k1, 0, n), (k2, 0, n))
>>> ok=_;var('p q');ok.subs(p,0).subs(q,1)
(p, q)
Sum(0**k1*0**(-k2 + n)*Abs(k1 - k2), (k1, 0, n), (k2, 0, n))
>>> _.doit()
Sum(0**k1*0**(-k2 + n)*Abs(k1 - k2), (k1, 0, n), (k2, 0, n))
```
#### 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 -->
* core
    * `x**(a + b)` will not expand unless `a` and `b` are both nonnegative or nonpositive or `x` is not zero
    * expressions containing any illegal character as defined in numbers now gives `expr.is_rational_function()->False`
<!-- END RELEASE NOTES -->

Update

The release notes on the wiki have been updated.

@sympy sympy deleted a comment from github-actions bot Aug 24, 2022
@smichr smichr marked this pull request as draft August 25, 2022 07:45
@smichr smichr marked this pull request as ready for review August 26, 2022 16:07
@smichr smichr changed the title restrict expand(0**(a+b)) expand(x**(a+b)) iff x != 0 or a + b != 0 Aug 26, 2022
@github-actions
Copy link

github-actions bot commented Aug 26, 2022

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
     [41d90958]       [ed600afb]
     <sympy-1.11.1^0>                 
-         995±2μs          628±3μs     0.63  solve.TimeSparseSystem.time_linear_eq_to_matrix(10)
-        2.86±0ms         1.17±0ms     0.41  solve.TimeSparseSystem.time_linear_eq_to_matrix(20)
-     5.71±0.02ms         1.71±0ms     0.30  solve.TimeSparseSystem.time_linear_eq_to_matrix(30)

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

@smichr
Copy link
Member Author

smichr commented Aug 29, 2022

Is there any reason we shouldn't do this?

@smichr
Copy link
Member Author

smichr commented Aug 29, 2022

I didn't think this would be so easy to fix. It is similar to not allowing x**2/x to become x unless x != 0 (something we don't prohibit). The power expanding is easier to control because it does not happen automatically. ping @asmeurer

sympy/core/expr.py Outdated Show resolved Hide resolved
@oscarbenjamin
Copy link
Contributor

The title of this PR says expand(x**(a+b)) iff x != 0 or a + b != 0 but I don't understand why it's okay to expand if a + b != 0 even though x might be zero:

In [30]: e1 = x**(y + z)

In [31]: e2 = e1.expand()

In [32]: e1
Out[32]: 
 y + z
x     

In [33]: e2
Out[33]: 
 y  z
xx 

In [34]: e1.subs(x, 0).subs({y:1,z:-2})
Out[34]: zoo

In [35]: e2.subs(x, 0).subs({y:1,z:-2})
Out[35]: 0

@smichr
Copy link
Member Author

smichr commented Aug 30, 2022

@oscarbenjamin , to get your example to work you would have to have assumptions on symbols in a sum that would allow the assumptions system to determine that the exponent wasn't zero. The example you gave relies on one number being positive and the other negative; that will not evaluate as having a non-zero assumption. I think the nearest we could come is with nonnegative + positive which must be positive. When expanded this would give 0**nonnegative*0**positive which would correctly give 0. Can you think of an example with assumptions that would behave like the expression you gave?

@smichr
Copy link
Member Author

smichr commented Aug 30, 2022

answer:

>>> one = Symbol('1',integer=1,prime=0,odd=1,positive=1)
>>> two = Symbol('2',integer=1,prime=1,even=1)
>>> eq = two-one
>>> 0**(2-1)
0
>>> 0**eq
0**(-1 + 2)
>>> expand(_)
nan
>>> (x**eq).expand().subs(x,0)
nan

>>> S.Zero**-1
zoo
>>> 0**-eq
0**(1 - 2)
>>> expand(_)
nan
>>> (x**-eq).expand().subs(x,0)
nan

@smichr smichr force-pushed the 0base branch 2 times, most recently from 2108f71 to 243b4e9 Compare August 30, 2022 05:38
@oscarbenjamin
Copy link
Contributor

to get your example to work you would have to have assumptions on symbols in a sum that would allow the assumptions system to determine that the exponent wasn't zero.

The question is under what conditions does the identity x**(a + b) == x**a * x**b hold?

The docs here claim that the identity is "always true":
https://docs.sympy.org/latest/tutorials/intro-tutorial/simplification.html#powers

I think that misses the possibility that x could be zero or infinite. If x is a nonzero finite complex number and a and b are finite complex numbers then all of x**a, x**b and x**(a + b) are well defined nonzero complex numbers and the identity holds. However if x is zero then we have
$$0^c = \begin{cases} 0 & c > 0 \\ 1 & c = 0 \\ \infty & c < 0 \\ \mathrm{undefined} & c \notin \mathbb{R} \end{cases}$$
Any situation where a + b, a and b are all nonnegative will be fine but otherwise we can have problems e.g.:

In [59]: e1 = x**(a + b)

In [60]: e2 = e1.expand()

In [61]: e1.subs(x, 0).subs({a:I, b:-I})
Out[61]: 1

In [62]: e2.subs(x, 0).subs({a:I, b:-I})
Out[62]: nan

In [67]: e1.subs(x, 0).subs({a:-1, b:2})
Out[67]: 0

In [68]: e2.subs(x, 0).subs({a:-1, b:2})
Out[68]: nan

In [69]: e1.subs(x, 0).subs({a:1, b:-2})
Out[69]: zoo

In [70]: e2.subs(x, 0).subs({a:1, b:-2})
Out[70]: 0

The last case is problematic because of the bogus evaluation of 0*x (#17224):

In [73]: e2.subs(x, 0).subs(b, -2)
Out[73]: 
 a    
0zoo

In [74]: e2.subs(x, 0).subs(b, -2).subs(a, 1)
Out[74]: nan

In [75]: e2.subs(x, 0).subs(a, 1).subs(b, -2)
Out[75]: 0

sympy/core/function.py Outdated Show resolved Hide resolved
@smichr
Copy link
Member Author

smichr commented Sep 1, 2022

I am not seeing why this is passing. I checked the Add._all_nonneg_or_nonpos and it gives None for y + 2 so why does the power expand?

    expand_power_exp(x**(y + 2))
Expected:
    x**(y + 2)
Got:
    x**2*x**y

sympy/core/power.py Outdated Show resolved Hide resolved
@smichr smichr changed the title expand(x**(a+b)) iff x != 0 or a + b != 0 expand(x**(a+b)) iff x != 0 or a & b both nonneg or nonpos Sep 1, 2022
sympy/core/power.py Outdated Show resolved Hide resolved
sympy/core/power.py Outdated Show resolved Hide resolved
sympy/core/power.py Outdated Show resolved Hide resolved
sympy/core/power.py Outdated Show resolved Hide resolved
sympy/core/power.py Outdated Show resolved Hide resolved
@smichr
Copy link
Member Author

smichr commented Sep 1, 2022

I will squash to two commits and merge if there are no other concerns.

@smichr smichr merged commit 7158ec4 into sympy:master Sep 2, 2022
@smichr smichr deleted the 0base branch September 2, 2022 02:19
Comment on lines 417 to +424
f = x**(n - 2)*cos(x)
assert fps(f, x).truncate() == \
(x**(n - 2) - x**n/2 + x**(n + 2)/24 - x**(n + 4)/720 + O(x**(n + 6), x))
(x**(n - 2) - x**n/2 + x**(n + 2)/24 + O(x**(n + 4), x))

f = x**(n - 2)*sin(x) + x**n*exp(x)
assert fps(f, x).truncate() == \
(x**(n - 1) + x**n + 5*x**(n + 1)/6 + x**(n + 2)/2 + 7*x**(n + 3)/40 +
x**(n + 4)/24 + 41*x**(n + 5)/5040 + O(x**(n + 6), x))
(x**(n - 1) + x**(n + 1) + x**(n + 2)/2 + x**n +
x**(n + 4)/24 + x**(n + 5)/60 + O(x**(n + 6), x))
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why have these changed?

@oscarbenjamin
Copy link
Contributor

I think this mostly looks good but please don't self merge PRs.

@oscarbenjamin
Copy link
Contributor

@mkoeppe after this PR the SAGE doctest is failing:

File "src/sage/calculus/calculus.py", line 1693, in sage.calculus.calculus.laplace
Failed example:
    laplace(t^n, t, s, algorithm='sympy')
Expected:
    (gamma(n + 1)/(s*s^n), 0, re(n) > -1)
Got:
    (s^(-n - 1)*gamma(n + 1), 0, re(n) > -1)

Is there a new SAGE ticket to track updating to SymPy 1.12 where a doctest like this could be updated?

@mkoeppe
Copy link
Contributor

mkoeppe commented Sep 3, 2022

I've opened https://trac.sagemath.org/ticket/34489 for this

@oscarbenjamin
Copy link
Contributor

So I guess the process is to update this line:

SAGE_TICKET: 34118

and then push a commit to the trac ticket?

@oscarbenjamin
Copy link
Contributor

I've opened a PR to update the ticket number: #24015

@mkoeppe
Copy link
Contributor

mkoeppe commented Sep 3, 2022

Yes, that's right.

Costor added a commit to Costor/sympy-contribution that referenced this pull request Mar 20, 2023
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Projects
None yet
Development

Successfully merging this pull request may close these issues.

expand(0**(1 - x)) -> 0
4 participants