Repository: sympy/sympy
Pull request: Fix wrong result in create_expand_pow_optimization(). Partially fixes #15335
<!-- 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 .-->
Fixes #15335 (partially)

#### Brief description of what is fixed or changed


#### 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 -->
* codegen
  * Fix wrong result in create_expand_pow_optimization()
<!-- END RELEASE NOTES -->


Problem statement:
create_expand_pow_optimization() is not robust
Using SymPy 1.3:

```
import sympy as sp
from sympy.codegen.rewriting import create_expand_pow_optimization
expand_opt = create_expand_pow_optimization(7)
a = sp.symbols("a")
expr = a**(5) - a**2
print("should be a**5 - a**2:",expr)
print("should be a*a*a*a*a - a*a:",expand_opt(expr))
expr = a**(-5)
print("should be a**(-5):",expr)
print("should be equal to a**(-5):", expand_opt(expr))
```

Output is:
```
should be a**5 - a**2: a**5 - a**2
should be a*a*a*a*a - a*a: -a**2 + a*a*a*a*a
should be a**(-5): a**(-5)
should be equal to a**(-5): 1
```

Notice that it failed to optimize `-a**2` and more disturbingly `a**(-5)` was evaluated to `1`.

Hints from issue discussion:
CC @bjodah 
Replacing

```
return ReplaceOptim(
        lambda e: e.is_Pow and e.base.is_symbol and e.exp.is_integer and e.exp <= limit,
        lambda p: Mul(*([p.base]*p.exp), evaluate=False)
    )
```

with
```
return ReplaceOptim(
        lambda e: e.is_Pow and e.base.is_symbol and e.exp.is_integer and e.exp.is_nonnegative and e.exp <= limit,
        lambda p: Mul(*([p.base]*p.exp), evaluate=False)
    )
```
in line 224 of sympy/codegen/rewriting.py 

fixes the most egregious problem -- i.e., `a**(-5) = 1`. I'm still confused as to why e.g., multiplying `a**2` by `-1` results in the optimization not being performed.
Looks like `replace` (via `xreplace`) doesn't respect unevaluated Muls. The issue is basically:

```py
>>> (-x).xreplace({x: Mul(a, a, evaluate=False)})
-a**2
```

because

```py
>>> Mul(-1, Mul(a, a, evaluate=False))
-a**2
```

I don't know how easy it is to fix. Unevaluated expressions are a bit fragile. It probably would be smarter to create a custom unevaluated Mul object for the code generator and replace it with that instead. 
Thank you for reporting this @zachetienne, would you mind opening a PR with your suggested fix, so that we at least aren't returning incorrect results?
Whoops, clicked the wrong button and closed the ticket. Reopening.
Here is the pull request: https://github.com/sympy/sympy/pull/15338

Commit message: Add new test of create_expand_pow_optimization() to validate against issue #15335

Patch to review:
diff --git a/sympy/codegen/rewriting.py b/sympy/codegen/rewriting.py
index c005717c56bf..f5742d2fb62d 100644
--- a/sympy/codegen/rewriting.py
+++ b/sympy/codegen/rewriting.py
@@ -221,7 +221,7 @@ def create_expand_pow_optimization(limit):
 
     """
     return ReplaceOptim(
-        lambda e: e.is_Pow and e.base.is_symbol and e.exp.is_integer and e.exp <= limit,
+        lambda e: e.is_Pow and e.base.is_symbol and e.exp.is_integer and e.exp.is_nonnegative and e.exp <= limit,
         lambda p: Mul(*([p.base]*p.exp), evaluate=False)
     )
 
diff --git a/sympy/codegen/tests/test_rewriting.py b/sympy/codegen/tests/test_rewriting.py
index 2001f95e728b..41fe7ae1daeb 100644
--- a/sympy/codegen/tests/test_rewriting.py
+++ b/sympy/codegen/tests/test_rewriting.py
@@ -169,3 +169,5 @@ def test_create_expand_pow_optimization():
 
     sin4x = sin(x)**4
     assert ccode(optimize(sin4x, [my_opt])) == 'pow(sin(x), 4)'
+
+    assert ccode(optimize((x**(-4)), [my_opt])) == 'x**(-4)'


Code review comments:
diff_hunk: @@ -169,3 +169,5 @@ def test_create_expand_pow_optimization():
 
     sin4x = sin(x)**4
     assert ccode(optimize(sin4x, [my_opt])) == 'pow(sin(x), 4)'
+
+    assert ccode(optimize((x**(-4)), [my_opt])) == 'x**(-4)'
path: sympy/codegen/tests/test_rewriting.py
text: @user1:
`x**(-4)` is not valid C code...

@author:
Whoops! Thanks for pointing that out! Here is a new patch, please take a look: https://github.com/zachetienne/sympy/commit/169d2b919ec874642dfa33172a92a3bd63fe7d57