Substitution for powers should be like polynomial division #22249

[robinarmingaud]

References to other Issues or PRs
Fixes #22249
Brief description of what is fixed or changed
_eval_subs in power.py modified : simple division replaced by euclidean division in order to substitute even if pow is not integer
test_issue_22249() Added
Other comments
Release Notes

• power
  • Changed how substitution is performed with polynomes

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

❌ There was an issue with the release notes. Please do not close this pull request; instead edit the description after reading the guide on how to write release notes.

• power is not a valid release notes header. Release notes headers should be SymPy submodule names, like

  * core
    * made Add faster
  
  * printing
    * improve LaTeX printing of fractions
  

  or other.

  If you have added a new submodule, please add it to the list of valid release notes headers at https://github.com/sympy/sympy-bot/blob/master/sympy_bot/submodules.txt.

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

References to other Issues or PRs
Fixes #22249 
Brief description of what is fixed or changed
_eval_subs in power.py modified : simple division replaced by euclidean division in order to substitute even if pow is not integer
test_issue_22249() Added
Other comments
Release Notes
<!-- BEGIN RELEASE NOTES -->
- power
   - Changed how substitution is performed with polynomes
<!-- END RELEASE NOTES -->

[asmeurer]

Does this now apply to symbolic powers too?

[faze-geek]

Ping @asmeurer @oscarbenjamin for suggestions/help.
I was working on this PR and built upon the previous code discussed above .My changes resolve all issues faced inlimits.pyand subs.py (2 obvious tests which will be resolved by this PR also changed ) .But probably since ode module is dependent on subs() there are some failed tests in single.py , systems.py as faced in the earlier checks too .I cannot figure out the root cause of those fails as of now .I'll attach the changes below .I would request some guidance .

             coeff2, terms2 = ct2
             if terms1 == terms2:
                 if old.is_commutative:
-                    # Allow fractional powers for commutative objects
-                    pow = coeff1/coeff2
-                    try:
-                        as_int(pow, strict=False)
-                        combines = True
-                    except ValueError:
-                        b, e = old.as_base_exp()
-                        # These conditions ensure that (b**e)**f == b**(e*f) for any f
-                        combines = b.is_positive and e.is_real or b.is_nonnegative and e.is_nonnegative
-
-                    return combines, pow, None
+                    lst = any(sym.is_nonnegative for sym in list(old.free_symbols.union(new.free_symbols)))
+                    if type(coeff1) is type(coeff2) and (coeff1 > coeff2) and isinstance(coeff1, Integer) and not (lst):
+                         pow = coeff1//coeff2
+                         remainder_pow = coeff1%coeff2
+                         return True, pow, remainder_pow
+                    else:
+                        # Allow fractional powers for commutative objects
+                        pow = coeff1/coeff2
+                        try:
+                            as_int(pow, strict=False)
+                            combines = True
+                        except ValueError:
+                            b, e = old.as_base_exp()
+                            # These conditions ensure that (b**e)**f == b**(e*f) for any f
+                            combines = b.is_positive and e.is_real or b.is_nonnegative and e.is_nonnegative
+                        return combines, pow, None
                 else:
                     # With noncommutative symbols, substitute only integer powers
                     if not isinstance(terms1, tuple):
diff --git a/sympy/core/tests/test_subs.py b/sympy/core/tests/test_subs.py
index 9c4a5e90c8..a0abf20ebc 100644
--- a/sympy/core/tests/test_subs.py
+++ b/sympy/core/tests/test_subs.py
@@ -829,11 +829,11 @@ def test_issue_13333():
 def test_issue_15234():
     x, y = symbols('x y', real=True)
     p = 6*x**5 + x**4 - 4*x**3 + 4*x**2 - 2*x + 3
-    p_subbed = 6*x**5 - 4*x**3 - 2*x + y**4 + 4*y**2 + 3
+    p_subbed = 6*x*y**4 - 4*x*y**2 - 2*x + y**4 + 4*y**2 + 3
     assert p.subs([(x**i, y**i) for i in [2, 4]]) == p_subbed
     x, y = symbols('x y', complex=True)
     p = 6*x**5 + x**4 - 4*x**3 + 4*x**2 - 2*x + 3
-    p_subbed = 6*x**5 - 4*x**3 - 2*x + y**4 + 4*y**2 + 3
+    p_subbed = 6*x*y**4 - 4*x*y**2 - 2*x + y**4 + 4*y**2 + 3
     assert p.subs([(x**i, y**i) for i in [2, 4]]) == p_subbed


@@ -842,7 +842,7 @@ def test_issue_6976():
     assert (sqrt(x)**3 + sqrt(x) + x + x**2).subs(sqrt(x), y) == \
         y**4 + y**3 + y**2 + y
     assert (x**4 + x**3 + x**2 + x + sqrt(x)).subs(x**2, y) == \
-        sqrt(x) + x**3 + x + y**2 + y
+        sqrt(x) + y*x + x + y**2 + y
     assert x.subs(x**3, y) == x
     assert x.subs(x**Rational(1, 3), y) == y**3

I would also report the issues faced below

Traceback (most recent call last):
  File "c:\users\kunni\sympy\sympy\sympy\solvers\ode\tests\test_single.py", line 612, in test_nth_order_linear_euler_eq_nonhomogeneous_undetermined_coefficients
    _ode_solver_test(_get_examples_ode_sol_euler_undetermined_coeff)
  File "c:\users\kunni\sympy\sympy\sympy\solvers\ode\tests\test_single.py", line 185, in _ode_solver_test
    result = _test_particular_example(example['hint'], example, solver_flag=True)
  File "c:\users\kunni\sympy\sympy\sympy\solvers\ode\tests\test_single.py", line 271, in _test_particular_example
    raise AssertionError(message)
AssertionError: Different solution found from dsolve for example euler_undet_02.

The ODE is:
Eq(x**2*Derivative(f(x), (x, 2)) - 2*x*Derivative(f(x), x) + 2*f(x), x**3)

The expected solution was
[Eq(f(x), x*(C1 + C2*x + x**2/2))]

What dsolve returned is:
Eq(f(x), x*(C1 + C2*x + E*x)/2)

Traceback (most recent call last):
  File "c:\users\kunni\sympy\sympy\sympy\solvers\ode\tests\test_systems.py", line 535, in test_sysode_linear_neq_order1_type1
    assert dsolve(eqs2) == sol2
AssertionError

Traceback (most recent call last):
  File "c:\users\kunni\sympy\sympy\sympy\solvers\ode\tests\test_systems.py", line 1639, in test_higher_order_to_first_order
    assert dsolve(eqs15) == sol15
AssertionError

Thank You.

[smichr]

I wonder if some of the ode related issues will stem from the assumption that this partial extraction is not going to work. So the code would have to be inspected and corrected where this assumption is made by checking the substitution result to verify that the target power was removed completely, e.g. ok = x not in eq.subs(x**2, y)

[smichr]

I would use something like this...

                         combines = b.is_positive and e.is_real or b.is_nonnegative and e.is_nonnegative
-
+                        if not combines:
+                            # check for partial extraction as in (x**3).subs(x**2, y) -> x*y
+                            if coeff1.is_Integer and coeff2.is_Integer and coeff1 >= coeff2:
+                                return True, *divmod(coeff1, coeff2)
                     return combines, pow, None