Repository: sympy/sympy

Issue:
integrate 1/sqrt(x**2) is wrong
```py
>>> x=Symbol('x')
>>> f=1/sqrt(x**2)
>>> integrate(f, risch=True)
log(sqrt(x**2))
>>> integrate(f, heurisch=True)
Integral(1/sqrt(x**2), x)
>>> integrate(f, manual=True)
log(2*x + 2*sqrt(x**2))
>>> integrate(f, meijerg=True)
Integral(1/sqrt(x**2), x)
>>> F=x*log(x)/sqrt(x**2)
>>> F.diff(x)
1/sqrt(x**2)
```
Result from risch and manual are both wrong. Risch should give `Integral` after #23689.
Result from Wolfram Alpha is correct.

Discussion:
sympy-bot: :white_check_mark:

Hi, I am the [SymPy bot](https://github.com/sympy/sympy-bot) (v167). 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:
* integrals
  * Support integral steps for nested pow like (c*(a+b*x)**d)**e ([#23716](https://github.com/sympy/sympy/pull/23716) by [@eagleoflqj](https://github.com/eagleoflqj))

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

<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 #23712
    Fixes #23505
    #### Brief description of what is fixed or changed

    This is an extension to power rule.
    Nested pow like `sqrt(x**(S(5)/3))` can be treated as `x**(S(5)/6)`, but the structure should be kept in the end.

    $$
    \begin{align}
    &\int \sqrt{x^\frac{5}{3}}\mathrm{d}x\\
    =&x\sqrt{x^\frac{5}{3}} - \int x\mathrm{d}\sqrt{x^\frac{5}{3}}\\
    =&x\sqrt{x^\frac{5}{3}} - \int x\cdot\frac{1}{2}(x^\frac{5}{3})^{-\frac{1}{2}}\cdot\frac{5}{3}x^\frac{2}{3}\mathrm{d}x\\
    =&x\sqrt{x^\frac{5}{3}} - \frac{5}{6}\int \sqrt{x^\frac{5}{3}}\mathrm{d}x\\
    =&\frac{6}{11}x\sqrt{x^\frac{5}{3}}
    \end{align}
    $$

    However, if the exponent is equivalent to -1, the integrand should be multiplied by `x*log(x)` instead of `x/(exponent+1)`
    #### Other comments


    #### 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 -->
    * integrals
      * Support integral steps for nested pow like (c*(a+b*x)**d)**e
    <!-- END RELEASE NOTES -->


</details><p>

github-actions[bot]: 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)
```diff

```
Significantly changed benchmark results (master vs previous release)
```diff
       before           after         ratio
     [77f1d79c]       [b64cfcdb]
     <sympy-1.10.1^0>                 
+      96.6±0.3ms        176±0.8ms     1.83  sum.TimeSum.time_doit

```
Full benchmark results can be found as artifacts in GitHub Actions
(click on checks at the top of the PR).
eagleoflqj: @oscarbenjamin is this good to go?

Patch:
diff --git a/sympy/integrals/manualintegrate.py b/sympy/integrals/manualintegrate.py
index 06c15ca6d2..802cd0dc84 100644
--- a/sympy/integrals/manualintegrate.py
+++ b/sympy/integrals/manualintegrate.py
@@ -77,6 +77,7 @@ def __eq__(self, other):
 ConstantRule = Rule("ConstantRule", "constant")
 ConstantTimesRule = Rule("ConstantTimesRule", "constant other substep")
 PowerRule = Rule("PowerRule", "base exp")
+NestedPowRule = Rule("NestedPowRule", "base exp")
 AddRule = Rule("AddRule", "substeps")
 URule = Rule("URule", "u_var u_func constant substep")
 PartsRule = Rule("PartsRule", "u dv v_step second_step")
@@ -225,7 +226,7 @@ def test_subterm(u, u_diff):
         debug("substituted: {}, u: {}, u_var: {}".format(substituted, u, u_var))
         substituted = manual_subs(substituted, u, u_var).cancel()
 
-        if symbol in substituted.free_symbols:
+        if substituted.has_free(symbol):
             return False
         # avoid increasing the degree of a rational function
         if integrand.is_rational_function(symbol) and substituted.is_rational_function(u_var):
@@ -486,6 +487,59 @@ def _add_degenerate_step(generic_cond, generic_step, degenerate_step):
     return PiecewiseRule(subfunctions, generic_step.context, generic_step.symbol)
 
 
+def nested_pow_rule(integral: IntegralInfo):
+    # nested (c*(a+b*x)**d)**e
+    integrand, x = integral
+
+    a_ = Wild('a', exclude=[x])
+    b_ = Wild('b', exclude=[x, 0])
+    pattern = a_+b_*x
+    generic_cond = S.true
+
+    class NoMatch(Exception):
+        pass
+
+    def _get_base_exp(expr: Expr) -> tuple[Expr, Expr]:
+        if not expr.has_free(x):
+            return S.One, S.Zero
+        if expr.is_Mul:
+            _, terms = expr.as_coeff_mul()
+            if not terms:
+                return S.One, S.Zero
+            results = [_get_base_exp(term) for term in terms]
+            bases = set(b for b, _ in results)
+            bases.discard(S.One)
+            if len(bases) == 1:
+                return bases.pop(), Add(*(e for _, e in results))
+            raise NoMatch
+        if expr.is_Pow:
+            b, e = expr.base, expr.exp  # type: ignore
+            if e.has_free(x):
+                raise NoMatch
+            base_, sub_exp = _get_base_exp(b)
+            return base_, sub_exp * e
+        match = expr.match(pattern)
+        if match:
+            a, b = match[a_], match[b_]
+            base_ = x + a/b
+            nonlocal generic_cond
+            generic_cond = Ne(b, 0)
+            return base_, S.One
+        raise NoMatch
+
+    try:
+        base, exp_ = _get_base_exp(integrand)
+    except NoMatch:
+        return
+    if generic_cond is S.true:
+        degenerate_step = None
+    else:
+        # equivalent with subs(b, 0) but no need to find b
+        degenerate_step = ConstantRule(integrand.subs(x, 0), integrand, x)
+    generic_step = NestedPowRule(base, exp_, integrand, x)
+    return _add_degenerate_step(generic_cond, generic_step, degenerate_step)
+
+
 def inverse_trig_rule(integral: IntegralInfo, degenerate=True):
     """
     Set degenerate=False on recursive call where coefficient of quadratic term
@@ -531,10 +585,17 @@ def make_inverse_trig(RuleClass, a, sign_a, c, sign_c, h):
 
     if simplify(2*exp + 1) == 0:
         h, k = -b/(2*c), a - b**2/(4*c)  # rewrite base to k + c*(symbol-h)**2
-        step = general_rule = ReciprocalSqrtQuadraticRule(a, b, c, integrand, symbol)
+        non_square_cond = Ne(k, 0)
+        square_step = None
+        if non_square_cond is not S.true:
+            square_step = NestedPowRule(symbol-h, -1, 1/sqrt(c*(symbol-h)**2), symbol)
+        if non_square_cond is S.false:
+            return square_step
+        generic_step = ReciprocalSqrtQuadraticRule(a, b, c, integrand, symbol)
+        step = _add_degenerate_step(non_square_cond, generic_step, square_step)
         if k.is_real and c.is_real:
             # list of ((rule, base_exp, a, sign_a, b, sign_b), condition)
-            possibilities = []
+            rules = []
             for args, cond in (  # don't apply ArccoshRule to x**2-1
                 ((ArcsinRule, k, 1, -c, -1, h), And(k > 0, c < 0)),  # 1-x**2
                 ((ArcsinhRule, k, 1, c, 1, h), And(k > 0, c > 0)),  # 1+x**2
@@ -542,14 +603,13 @@ def make_inverse_trig(RuleClass, a, sign_a, c, sign_c, h):
                 if cond is S.true:
                     return make_inverse_trig(*args)
                 if cond is not S.false:
-                    possibilities.append((args, cond))
-            if possibilities:
-                rules = [(make_inverse_trig(*args), cond) for args, cond in possibilities]
+                    rules.append((make_inverse_trig(*args), cond))
+            if rules:
                 if not k.is_positive:  # conditions are not thorough, need fall back rule
-                    rules.append((general_rule, S.true))
+                    rules.append((generic_step, S.true))
                 step = PiecewiseRule(rules, integrand, symbol)
             else:
-                step = general_rule
+                step = generic_step
         return _add_degenerate_step(generic_cond, step, degenerate_step)
     if exp == S.Half:
         step = SqrtQuadraticRule(a, b, c, integrand, symbol)
@@ -1522,6 +1582,7 @@ def _integral_is_subclass(integral):
                 trig_powers_products_rule,
                 trig_expand_rule
             )),
+            null_safe(condition(integral_is_subclass(Mul, Pow), nested_pow_rule)),
             null_safe(trig_substitution_rule)
         ),
         fallback_rule)(integral)
@@ -1543,6 +1604,14 @@ def eval_power(base, exp, integrand, symbol):
         (log(base), True),
         )
 
+
+@evaluates(NestedPowRule)
+def eval_nested_pow(base, exp_, integrand, x):
+    m = base*integrand
+    return Piecewise((m/(exp_+1), Ne(exp_, -1)),
+                     (m*log(base), True))
+
+
 @evaluates(ExpRule)
 def eval_exp(base, exp, integrand, symbol):
     return integrand / log(base)


Test patch:
diff --git a/sympy/integrals/tests/test_integrals.py b/sympy/integrals/tests/test_integrals.py
index b8f02f0c0c..1ead89618c 100644
--- a/sympy/integrals/tests/test_integrals.py
+++ b/sympy/integrals/tests/test_integrals.py
@@ -2019,14 +2019,19 @@ def test_hyperbolic():
 
 def test_nested_pow():
     assert integrate(sqrt(x**2)) == x*sqrt(x**2)/2
+    assert integrate(sqrt(x**(S(5)/3))) == 6*x*sqrt(x**(S(5)/3))/11
+    assert integrate(1/sqrt(x**2)) == x*log(x)/sqrt(x**2)
+    assert integrate(x*sqrt(x**(-4))) == x**2*sqrt(x**-4)*log(x)
 
 
 def test_sqrt_quadratic():
     assert integrate(1/sqrt(3*x**2+4*x+5)) == sqrt(3)*asinh(3*sqrt(11)*(x + S(2)/3)/11)/3
     assert integrate(1/sqrt(-3*x**2+4*x+5)) == sqrt(3)*asin(3*sqrt(19)*(x - S(2)/3)/19)/3
     assert integrate(1/sqrt(3*x**2+4*x-5)) == sqrt(3)*log(6*x + 2*sqrt(3)*sqrt(3*x**2 + 4*x - 5) + 4)/3
+    assert integrate(1/sqrt(4*x**2-4*x+1)) == (x - S.Half)*log(x - S.Half)/(2*sqrt((x - S.Half)**2))
     assert integrate(1/sqrt(a+b*x+c*x**2), x) == \
-        Piecewise((log(b + 2*sqrt(c)*sqrt(a + b*x + c*x**2) + 2*c*x)/sqrt(c), Ne(c, 0)),
+        Piecewise((log(b + 2*sqrt(c)*sqrt(a + b*x + c*x**2) + 2*c*x)/sqrt(c), Ne(c, 0) & Ne(a - b**2/(4*c), 0)),
+                  ((b/(2*c) + x)*log(b/(2*c) + x)/sqrt(c*(b/(2*c) + x)**2), Ne(c, 0)),
                   (2*sqrt(a + b*x)/b, Ne(b, 0)), (x/sqrt(a), True))
 
     assert integrate((7*x+6)/sqrt(3*x**2+4*x+5)) == \
@@ -2036,8 +2041,10 @@ def test_sqrt_quadratic():
     assert integrate((7*x+6)/sqrt(3*x**2+4*x-5)) == \
            7*sqrt(3*x**2 + 4*x - 5)/3 + 4*sqrt(3)*log(6*x + 2*sqrt(3)*sqrt(3*x**2 + 4*x - 5) + 4)/9
     assert integrate((d+e*x)/sqrt(a+b*x+c*x**2), x) == \
-        Piecewise((e*sqrt(a + b*x + c*x**2)/c +
-                   (-b*e/(2*c) + d)*log(b + 2*sqrt(c)*sqrt(a + b*x + c*x**2) + 2*c*x)/sqrt(c), Ne(c, 0)),
+        Piecewise(((-b*e/(2*c) + d) *
+                   Piecewise((log(b + 2*sqrt(c)*sqrt(a + b*x + c*x**2) + 2*c*x)/sqrt(c), Ne(a - b**2/(4*c), 0)),
+                             ((b/(2*c) + x)*log(b/(2*c) + x)/sqrt(c*(b/(2*c) + x)**2), True)) +
+                   e*sqrt(a + b*x + c*x**2)/c, Ne(c, 0)),
                   ((2*d*sqrt(a + b*x) + 2*e*(-a*sqrt(a + b*x) + (a + b*x)**(S(3)/2)/3)/b)/b, Ne(b, 0)),
                   ((d*x + e*x**2/2)/sqrt(a), True))
 
@@ -2048,8 +2055,10 @@ def test_sqrt_quadratic():
            (x/2 - S(16683)/53225)*sqrt(53225*x**2 - 66732*x + 23013) + \
            111576969*sqrt(2129)*asinh(53225*x/10563 - S(11122)/3521)/1133160250
     assert integrate(sqrt(a+b*x+c*x**2), x) == \
-        Piecewise(((x/2 + b/(4*c))*sqrt(a + b*x + c*x**2) +
-                   (a/2 - b**2/(8*c))*log(b + 2*sqrt(c)*sqrt(a + b*x + c*x**2) + 2*c*x)/sqrt(c), Ne(c, 0)),
+        Piecewise(((a/2 - b**2/(8*c)) *
+                   Piecewise((log(b + 2*sqrt(c)*sqrt(a + b*x + c*x**2) + 2*c*x)/sqrt(c), Ne(a - b**2/(4*c), 0)),
+                             ((b/(2*c) + x)*log(b/(2*c) + x)/sqrt(c*(b/(2*c) + x)**2), True)) +
+                   (b/(4*c) + x/2)*sqrt(a + b*x + c*x**2), Ne(c, 0)),
                   (2*(a + b*x)**(S(3)/2)/(3*b), Ne(b, 0)),
                   (sqrt(a)*x, True))
 
diff --git a/sympy/integrals/tests/test_manual.py b/sympy/integrals/tests/test_manual.py
index b4109e2753..87659ffed5 100644
--- a/sympy/integrals/tests/test_manual.py
+++ b/sympy/integrals/tests/test_manual.py
@@ -37,6 +37,7 @@ def test_find_substitutions():
     assert find_substitutions((sec(x)**2 + tan(x) * sec(x)) / (sec(x) + tan(x)),
                               x, u) == [(sec(x) + tan(x), 1, 1/u)]
     assert (-x**2, Rational(-1, 2), exp(u)) in find_substitutions(x * exp(-x**2), x, u)
+    assert not find_substitutions(Derivative(f(x), x)**2, x, u)
 
 
 def test_manualintegrate_polynomials():
@@ -600,6 +601,20 @@ def test_issue_23566():
     assert str(i.n()) == '1.31695789692482'
 
 
+def test_nested_pow():
+    assert_is_integral_of(sqrt(x**2), x*sqrt(x**2)/2)
+    assert_is_integral_of(sqrt(x**(S(5)/3)), 6*x*sqrt(x**(S(5)/3))/11)
+    assert_is_integral_of(1/sqrt(x**2), x*log(x)/sqrt(x**2))
+    assert_is_integral_of(x*sqrt(x**(-4)), x**2*sqrt(x**-4)*log(x))
+    f = (c*(a+b*x)**d)**e
+    F1 = (c*(a + b*x)**d)**e*(a/b + x)/(d*e + 1)
+    F2 = (c*(a + b*x)**d)**e*(a/b + x)*log(a/b + x)
+    assert manualintegrate(f, x) == \
+        Piecewise((Piecewise((F1, Ne(d*e, -1)), (F2, True)), Ne(b, 0)), (x*(a**d*c)**e, True))
+    assert F1.diff(x).equals(f)
+    assert F2.diff(x).subs(d*e, -1).equals(f)
+
+
 def test_manualintegrate_sqrt_linear():
     assert_is_integral_of((5*x**3+4)/sqrt(2+3*x),
                           10*(3*x + 2)**(S(7)/2)/567 - 4*(3*x + 2)**(S(5)/2)/27 +
@@ -623,8 +638,10 @@ def test_manualintegrate_sqrt_quadratic():
     assert_is_integral_of(1/sqrt(3*x**2+4*x+5), sqrt(3)*asinh(3*sqrt(11)*(x + S(2)/3)/11)/3)
     assert_is_integral_of(1/sqrt(-3*x**2+4*x+5), sqrt(3)*asin(3*sqrt(19)*(x - S(2)/3)/19)/3)
     assert_is_integral_of(1/sqrt(3*x**2+4*x-5), sqrt(3)*log(6*x + 2*sqrt(3)*sqrt(3*x**2 + 4*x - 5) + 4)/3)
+    assert_is_integral_of(1/sqrt(4*x**2-4*x+1), (x - S.Half)*log(x - S.Half)/(2*sqrt((x - S.Half)**2)))
     assert manualintegrate(1/sqrt(a+b*x+c*x**2), x) == \
-        Piecewise((log(b + 2*sqrt(c)*sqrt(a + b*x + c*x**2) + 2*c*x)/sqrt(c), Ne(c, 0)),
+        Piecewise((log(b + 2*sqrt(c)*sqrt(a + b*x + c*x**2) + 2*c*x)/sqrt(c), Ne(c, 0) & Ne(a - b**2/(4*c), 0)),
+                  ((b/(2*c) + x)*log(b/(2*c) + x)/sqrt(c*(b/(2*c) + x)**2), Ne(c, 0)),
                   (2*sqrt(a + b*x)/b, Ne(b, 0)), (x/sqrt(a), True))
 
     assert_is_integral_of((7*x+6)/sqrt(3*x**2+4*x+5),
@@ -634,8 +651,10 @@ def test_manualintegrate_sqrt_quadratic():
     assert_is_integral_of((7*x+6)/sqrt(3*x**2+4*x-5),
                           7*sqrt(3*x**2 + 4*x - 5)/3 + 4*sqrt(3)*log(6*x + 2*sqrt(3)*sqrt(3*x**2 + 4*x - 5) + 4)/9)
     assert manualintegrate((d+e*x)/sqrt(a+b*x+c*x**2), x) == \
-        Piecewise((e*sqrt(a + b*x + c*x**2)/c +
-                   (-b*e/(2*c) + d)*log(b + 2*sqrt(c)*sqrt(a + b*x + c*x**2) + 2*c*x)/sqrt(c), Ne(c, 0)),
+        Piecewise(((-b*e/(2*c) + d) *
+                   Piecewise((log(b + 2*sqrt(c)*sqrt(a + b*x + c*x**2) + 2*c*x)/sqrt(c), Ne(a - b**2/(4*c), 0)),
+                             ((b/(2*c) + x)*log(b/(2*c) + x)/sqrt(c*(b/(2*c) + x)**2), True)) +
+                   e*sqrt(a + b*x + c*x**2)/c, Ne(c, 0)),
                   ((2*d*sqrt(a + b*x) + 2*e*(-a*sqrt(a + b*x) + (a + b*x)**(S(3)/2)/3)/b)/b, Ne(b, 0)),
                   ((d*x + e*x**2/2)/sqrt(a), True))
 
@@ -646,11 +665,14 @@ def test_manualintegrate_sqrt_quadratic():
                           (x/2 - S(16683)/53225)*sqrt(53225*x**2 - 66732*x + 23013) +
                           111576969*sqrt(2129)*asinh(53225*x/10563 - S(11122)/3521)/1133160250)
     assert manualintegrate(sqrt(a+c*x**2), x) == \
-        Piecewise((a*log(2*sqrt(c)*sqrt(a + c*x**2) + 2*c*x)/(2*sqrt(c)) + x*sqrt(a + c*x**2)/2, Ne(c, 0)),
+        Piecewise((a*Piecewise((log(2*sqrt(c)*sqrt(a + c*x**2) + 2*c*x)/sqrt(c), Ne(a, 0)),
+                               (x*log(x)/sqrt(c*x**2), True))/2 + x*sqrt(a + c*x**2)/2, Ne(c, 0)),
                   (sqrt(a)*x, True))
     assert manualintegrate(sqrt(a+b*x+c*x**2), x) == \
-        Piecewise(((x/2 + b/(4*c))*sqrt(a + b*x + c*x**2) +
-                   (a/2 - b**2/(8*c))*log(b + 2*sqrt(c)*sqrt(a + b*x + c*x**2) + 2*c*x)/sqrt(c), Ne(c, 0)),
+        Piecewise(((a/2 - b**2/(8*c)) *
+                   Piecewise((log(b + 2*sqrt(c)*sqrt(a + b*x + c*x**2) + 2*c*x)/sqrt(c), Ne(a - b**2/(4*c), 0)),
+                             ((b/(2*c) + x)*log(b/(2*c) + x)/sqrt(c*(b/(2*c) + x)**2), True)) +
+                   (b/(4*c) + x/2)*sqrt(a + b*x + c*x**2), Ne(c, 0)),
                   (2*(a + b*x)**(S(3)/2)/(3*b), Ne(b, 0)),
                   (sqrt(a)*x, True))