Reduce a power of a linear over a power of a quadratic a power at a time

BREAKING-CHANGES.md

@@ -447,6 +447,23 @@ case is most of its length ([#718](https://github.com/asc-community/AngouriMath/
 | `"x^4/((a + b*x^2 + c*x^4)*sqrt(d + k*x^2))".ToEntity().Integrate("x")` | `integral(...)` | 28,551 characters |
 | `"1/(x^2*(a + b*x^2 + c*x^4)*sqrt(d + k*x^2))".ToEntity().Integrate("x")` | `integral(...)` | 10,194 characters |
 
+### A power of a linear over a power of a quadratic is reduced a power at a time
+
+**Answers where there were none.** `sqrt(d + k x)/(a + b x + c x^2)^3` took thirty seconds for 26,437 characters on the
+unreleased master, and `(d + k x)^(7/2)/(a + b x + c x^2)^2` ran past a minute. A power of a linear that is not whole,
+times a linear, over a whole power of a quadratic is now reduced one power of the quadratic at a time: what is left each
+time is the same shape with the quadratic's power one lower, and over the first power the integrator answers it
+([#718](https://github.com/asc-community/AngouriMath/issues/718)).
+
+| Input | Was (2.5.0) | Now |
+|---|---|---|
+| `"(d + k*x)^(7/2)/(a + b*x + c*x^2)^2".ToEntity().Integrate("x")` | `integral(...)` | 4,691 characters |
+| `"1/((d + k*x)^(3/2)*(a + b*x + c*x^2)^2)".ToEntity().Integrate("x")` | `integral(...)` | 2,789 characters |
+| `"(d + k*x)^(7/2)/(a + b*x + c*x^2)^3".ToEntity().Integrate("x")` | `integral(...)` | 7,149 characters |
+| `"sqrt(d + k*x)/(a + b*x + c*x^2)^3".ToEntity().Integrate("x")` | `integral(...)` | 4,819 characters |
+| `"(A + B*x)/(x^(3/2)*(a + b*x + c*x^2)^2)".ToEntity().Integrate("x")` | `integral(...)` | 1,169 characters |
+| `"(b + 2*c*x)*(d + k*x)^(3/2)/(a + b*x + c*x^2)^3".ToEntity().Integrate("x")` | `integral(...)` | 2,505 characters |
+
 ### A power of the cosine beside a secant binomial is integrated in the cosine
 
 **Answers where there were none.** `cos(c + d x)^2/(a + b sec(c + d x))` was declined, and so were the cosine's powers

Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs

@@ -19544,6 +19544,137 @@ void Add(int j, int k, Entity[] poly)
             return Functions.PartialFractions.DerivativeHoldsAtSampledPoints(answer, expr, x) ? answer : null;
         }
 
+        /// <summary>
+        /// A power of a linear that is not whole, times a linear, over a whole power of a quadratic,
+        /// with the power of the quadratic lowered one at a time:
+        /// <c>integral L^m N Q^p = (alpha + beta x) L^(m + 1) Q^(p + 1) + integral L^m W Q^(p + 1)</c>
+        /// for <c>p &lt;= -2</c>, where <c>alpha</c> and <c>beta</c> make
+        /// <c>N - (alpha + beta x) (p + 1) L Q'</c> a multiple of <c>Q</c> and
+        /// <c>W = (N - (alpha + beta x) (p + 1) L Q')/Q - beta L - (alpha + beta x) (m + 1) e</c> is
+        /// linear again; at <c>p = -1</c> the integrator answers what is left.
+        /// </summary>
+        /// <remarks>
+        /// The derivative of the first term is <c>L^m Q^p</c> times
+        /// <c>beta L Q + (alpha + beta x) (m + 1) e Q + (alpha + beta x) L (p + 1) Q'</c>, which with
+        /// <c>L^m Q^(p + 1) W</c> added is <c>L^m Q^p N</c>. <c>(alpha + beta x) L Q'</c> modulo
+        /// <c>Q</c> is <c>(alpha + beta x) (r1 x + r0)</c> with <c>r1 = 2 c d - b e</c> and
+        /// <c>r0 = b d - 2 a e</c>, for <c>L = d + e x</c> and <c>Q = a + b x + c x^2</c>, and
+        /// matching it to <c>N</c> is two equations in two unknowns. Rubi's 1.2.1.2
+        /// <c>(d + e x)^(7/2)/(a + b x + c x^2)^2</c> ran past a minute and
+        /// <c>sqrt(d + e x)/(a + b x + c x^2)^3</c> took thirty seconds for 26,437 characters.
+        /// https://github.com/asc-community/AngouriMath/issues/718
+        /// </remarks>
+        internal static Entity? SolveAPowerOfALinearOverAPowerOfAQuadratic(Entity expr, Entity.Variable x, bool integrateByParts)
+        {
+            if (!Integration.AnsweringTheQuestionAsked)
+                return null;
+            var (numerator, denominator) = Functions.SingleQuotient.Of(expr);
+            Entity? quadratic = null;
+            var power = 0;
+            Entity? linear = null;
+            Entity exponent = Number.Integer.Zero;
+            Entity? other = null;
+            Entity constant = Number.Integer.One;
+            foreach (var (side, below) in new[] { (numerator, false), (denominator, true) })
+                foreach (var factor in Mulf.LinearChildren(side))
+                {
+                    if (!factor.ContainsNode(x))
+                    {
+                        constant = below ? constant / factor : constant * factor;
+                        continue;
+                    }
+                    var (@base, raised) = factor is Powf(var b0, var raisedTo) && !raisedTo.ContainsNode(x) ? (b0, raisedTo) : (factor, (Entity)Number.Integer.One);
+                    if (below && quadratic is null && raised is Number.Integer { EInteger: var whole } && whole.CompareTo(EInteger.FromInt32(2)) >= 0
+                        && whole.CanFitInInt32() && TreeAnalyzer.TryGetPolyQuadratic(@base, x, out var leading, out _, out _) && !VanishesIdentically(leading))
+                    {
+                        quadratic = @base;
+                        power = whole.ToInt32Unchecked();
+                        continue;
+                    }
+                    if (!TreeAnalyzer.TryGetPolyLinear(@base, x, out var slope, out _) || VanishesIdentically(slope))
+                        return null;
+                    // The linear to a power that is not whole, either side; one more linear above, once.
+                    if (raised is not Number.Integer && (linear is null || linear == @base))
+                    {
+                        linear = @base;
+                        exponent = below ? exponent - raised : exponent + raised;
+                        continue;
+                    }
+                    if (!below && other is null && raised == Number.Integer.One)
+                    {
+                        other = @base;
+                        continue;
+                    }
+                    return null;
+                }
+            exponent = exponent.InnerSimplified;
+            if (quadratic is null || linear is null || exponent is Number.Integer
+                || !TreeAnalyzer.TryGetPolyQuadratic(quadratic, x, out var c, out var b, out var a)
+                || !TreeAnalyzer.TryGetPolyLinear(linear, x, out var e, out var d)
+                || VanishesIdentically(b * b - 4 * a * c))
+                return null;
+            // Not a positive power of x itself, which u = sqrt(x) answers sooner -- the
+            // reduction took `x^(9/2)/(a + b x + c x^2)^3` from 2.6 seconds to 8.3; not a quadratic
+            // over the imaginary unit; and not one that factors over linears, b^2 - 4 a c the
+            // square of a polynomial, which the partial fractions take apart cheaply.
+            if (linear == x && exponent is Number.Rational { IsPositive: true } || Functions.PartialFractions.HoldsTheImaginaryUnit(quadratic))
+                return null;
+            var opposite = Functions.PartialFractions.Bare((b * b - 4 * a * c).Expand().InnerSimplified);
+            var symbols = opposite.Vars.Distinct().OrderBy(v => v.Name, System.StringComparer.Ordinal).ToArray();
+            var indices = new Dictionary<Variable, int>();
+            for (var i = 0; i < symbols.Length; i++)
+                indices[symbols[i]] = i;
+            if (symbols.Length <= MultivariatePolynomial.MaxVariables
+                && MultivariatePolynomial.TryParse(opposite, indices) is { } parsed && parsed.TrySquareRoot() is not null)
+                return null;
+            Entity nb = Number.Integer.Zero, na = constant;
+            if (other is not null)
+            {
+                if (!TreeAnalyzer.TryGetPolyLinear(other, x, out var ob, out var oa))
+                    return null;
+                nb = (constant * ob).InnerSimplified;
+                na = (constant * oa).InnerSimplified;
+            }
+            var r1 = (2 * c * d - b * e).InnerSimplified;
+            var r0 = (b * d - 2 * a * e).InnerSimplified;
+            Entity answer = Number.Integer.Zero;
+            for (var p = -power; p <= -2; p++)
+            {
+                // (p + 1) [r1 alpha + (r0 - r1 b/c) beta] = B, (p + 1) [r0 alpha - (r1 a/c) beta] = A.
+                var m11 = r1;
+                var m12 = r0 - r1 * b / c;
+                var m21 = r0;
+                var m22 = -r1 * a / c;
+                var determinant = Functions.PartialFractions.InLowestTermsOverTheSymbols((p + 1) * (m11 * m22 - m12 * m21));
+                if (VanishesIdentically(determinant))
+                    return null;
+                var alpha = Functions.PartialFractions.InLowestTermsOverTheSymbols((nb * m22 - m12 * na) / determinant);
+                var beta = Functions.PartialFractions.InLowestTermsOverTheSymbols((m11 * na - m21 * nb) / determinant);
+                // N - (alpha + beta x) (p + 1) L Q', a cubic, over Q: its quotient is linear, read off
+                // the two highest coefficients, the remainder being zero by the choice of alpha, beta.
+                var c3 = -(p + 1) * beta * e * 2 * c;
+                var c2 = -(p + 1) * (beta * (d * 2 * c + e * b) + alpha * e * 2 * c);
+                var q1 = (c3 / c).InnerSimplified;
+                var q0 = ((c2 - q1 * b) / c).InnerSimplified;
+                // W = quotient - beta L - (alpha + beta x)(m + 1) e.
+                var wb = Functions.PartialFractions.InLowestTermsOverTheSymbols(q1 - beta * e - beta * (exponent + 1) * e);
+                var wa = Functions.PartialFractions.InLowestTermsOverTheSymbols(q0 - beta * d - alpha * (exponent + 1) * e);
+                answer = answer + (alpha + beta * x) * MathS.Pow(linear, (exponent + 1).InnerSimplified) * MathS.Pow(quadratic, Number.Integer.Create(p + 1));
+                nb = wb;
+                na = wa;
+            }
+            // What is left asked with its two coefficients named: written out they are long, and the
+            // integrator took past a minute over them where with two symbols it takes a third of a
+            // second.
+            var namedA = Variable.CreateUnique(expr, "A_l");
+            var namedB = Variable.CreateUnique(expr + namedA, "B_l");
+            var left = (namedA + namedB * x) * MathS.Pow(linear, exponent) / quadratic;
+            if (Integration.ComputeAsAQuestionOfItsOwn(left, x, integrateByParts) is not { } rest)
+                return null;
+            answer = (answer + rest.Substitute(namedA, na).Substitute(namedB, nb)).InnerSimplified;
+            return Functions.PartialFractions.DerivativeHoldsAtSampledPoints(answer, expr, x) ? answer : null;
+        }
+
         /// <summary>
         /// A power of a linear beside the derivative of a quadratic, over a whole power of that
         /// quadratic, by parts once: <c>Q' L^m/Q^n</c> with <c>u = L^m</c> and

Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs

@@ -669,6 +669,9 @@ private static Entity WithItsPiecewiseSumsCombined(Entity expr)
             // A power of a linear beside the derivative of a quadratic over a whole power of it, by
             // parts once, before anything takes the symbolic quadratic apart.
             if ((answer = IndefiniteIntegralSolver.SolveAPowerOfALinearBesideTheDerivativeOfAQuadraticOverItsPower(expr, x, integrateByParts)) is { }) return answer;
+            // A power of a linear that is not whole, times a linear, over a whole power of a quadratic,
+            // the power of the quadratic lowered one at a time.
+            if ((answer = IndefiniteIntegralSolver.SolveAPowerOfALinearOverAPowerOfAQuadratic(expr, x, integrateByParts)) is { }) return answer;
             // A rational function over a whole power of a symbolic quadratic, the power lowered one
             // at a time, before the partial fractions take the quadratic apart over its roots.
             if ((answer = IndefiniteIntegralSolver.SolveARationalFunctionOverAPowerOfASymbolicQuadratic(expr, x, integrateByParts)) is { }) return answer;

Sources/Tests/UnitTests/Calculus/APowerOfALinearOverAPowerOfAQuadraticIntegralTest.cs

@@ -0,0 +1,59 @@
+//
+// Copyright (c) 2019-2026 Angouri.
+// AngouriMath is licensed under MIT.
+// Details: https://github.com/asc-community/AngouriMath/blob/master/LICENSE.md.
+// Website: https://am.angouri.org.
+//
+
+using System;
+using AngouriMath.Extensions;
+using Xunit;
+
+namespace AngouriMath.Tests.Calculus
+{
+    /// <summary>
+    /// A half-odd power of a linear over a whole power of a quadratic, the power of the quadratic
+    /// lowered one at a time. Rubi's 1.2.1.2 <c>sqrt(d + e x)/(a + b x + c x^2)^3</c>, written here
+    /// with <c>k</c> for <c>e</c>, which is Euler's number.
+    /// <a href="https://github.com/asc-community/AngouriMath/issues/718">#718</a>
+    /// </summary>
+    [Trait("Area", "Calculus")]
+    public sealed class APowerOfALinearOverAPowerOfAQuadraticIntegralTest
+    {
+        private static void DifferentiatesBack(string integrand, params (string, double)[] pins)
+        {
+            var integral = integrand.ToEntity().Integrate("x");
+            Assert.DoesNotContain("integral(", integral.Stringize());
+
+            // Pinned before it is differentiated: what the integrator answers over the first power
+            // is piecewise on signs read where they are numbers.
+            var pinned = integral.Substitute("C", 0);
+            Entity original = integrand.ToEntity();
+            foreach (var (name, value) in pins)
+            {
+                pinned = pinned.Substitute(name, value);
+                original = original.Substitute(name, value);
+            }
+            var derivative = pinned.Differentiate("x");
+            foreach (var at in new[] { -2.3, -0.7, 0.4, 1.9 })
+            {
+                var got = derivative.Substitute("x", at).EvalNumerical();
+                var want = original.Substitute("x", at).EvalNumerical();
+                var difference = Math.Abs((double)(got - want).RealPart) + Math.Abs((double)(got - want).ImaginaryPart);
+                var scale = Math.Max(1.0, Math.Abs((double)want.RealPart) + Math.Abs((double)want.ImaginaryPart));
+                Assert.True(difference / scale < 1e-8,
+                    $"d/dx of the antiderivative of {integrand} is {got} at x = {at}, where the integrand is {want}");
+            }
+        }
+
+        [Theory]
+        [InlineData("sqrt(d+k*x)/(a+b*x+c*x^2)^3")]
+        [InlineData("(d+k*x)^(7/2)/(a+b*x+c*x^2)^2")]
+        [InlineData("1/((d+k*x)^(3/2)*(a+b*x+c*x^2)^2)")]
+        public void IsAnsweredAndDifferentiatesBack(string integrand)
+        {
+            DifferentiatesBack(integrand, ("a", 2), ("b", 1.5), ("c", 1.25), ("d", 1.25), ("k", 2));
+            DifferentiatesBack(integrand, ("a", -1), ("b", 2), ("c", 2), ("d", -0.5), ("k", 1.5));
+        }
+    }
+}