Powers of two linears whose exponents sum to a whole number below -2 are integrated

BREAKING-CHANGES.md

@@ -1168,6 +1168,24 @@ in both quadratics, the cube's coefficients grow past what the reduction takes,
 | `"(7 + 13*x)/((5 + x + 2*x^2)^3*sqrt(2 + x + 3*x^2))".ToEntity().Integrate("x")` | `integral(...)` | an antiderivative |
 | `"sqrt(a + a*sec(x))/(c + d*sec(x))^2".ToEntity().Integrate("x")` | `integral(...)` | an antiderivative in `tan(x/2)` |
 
+### Powers of two linears whose exponents sum to a whole number below `-2` are integrated
+
+**Answers where there were none.** `(a + b x)^m (c + d x)^(-3 - m)`, `x^(n - 4)/(a + b x)^n` and
+the rest of Rubi's `(a + b x)^m (c + d x)^n` problems whose exponents sum to a whole number `-k`
+below `-2`, beside a polynomial of degree at most `k - 2`, were declined; the sum `-2` alone was
+answered. Under `t = (a + b x)/(c + d x)` each is a power of `t` beside a polynomial in `t`, so the
+antiderivative is the two powers as written times a polynomial in `x`, exact wherever the
+integrand is defined and taking each `m + i` it divides by to be nonzero, as the integrator does
+everywhere. Numeric exponents that sum so go the same way, so an answer master already gave for
+one, `(a + b x)^(5/2)/(c + d x)^(11/2)`, is written so now, the same function
+([#718](https://github.com/asc-community/AngouriMath/issues/718)).
+
+| Input | Was (2.5.0) | Now |
+|---|---|---|
+| `"(a + b*x)^m*(c + d*x)^(-3 - m)".ToEntity().Integrate("x")` | `integral(...)` | the two powers times `-(d (a + b x)^2 (c + d x)/(m + 2) - b (a + b x) (c + d x)^2/(m + 1))/(a d - b c)^2` |
+| `"x^(n - 4)/(a + b*x)^n".ToEntity().Integrate("x")` | `integral(...)` | the two powers times a cubic in `x` and `a + b x`, over `a^3` |
+| `"(a + b*x)^(5/2)/(c + d*x)^(11/2)".ToEntity().Integrate("x")` | `integral(...)` | the two powers times `-(2 d (a + b x)^2 (c + d x)/9 - 2 b (a + b x) (c + d x)^2/7)/(a d - b c)^2` |
+
 ### A polynomial beside the square roots of two linears with symbols in them is integrated
 
 `P(x) sqrt(c + d x) sqrt(e + f x)`, a polynomial beside the roots of two linears with symbols in

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

@@ -11936,6 +11936,129 @@ static bool IsRationalInBoth(Entity e, Entity.Variable x, Entity.Variable L)
             return null;
         }
 
+        /// <summary>
+        /// Powers of two linears whose exponents sum to a whole number <c>-k</c>, <c>k &gt;= 2</c>,
+        /// beside a polynomial of degree at most <c>k - 2</c>: <c>(a + b x)^m (c + d x)^(-3 - m)</c>,
+        /// <c>x^(n - 4)/(a + b x)^n</c>, <c>(e + f x)^2 (a + b x)^m (c + d x)^(-4 - m)</c>. Under
+        /// <c>t = L1/L2</c> each is <c>t^A</c> beside a polynomial in <c>t</c>, so the answer is
+        /// the two powers as written times a polynomial in <c>x</c>.
+        /// </summary>
+        /// <remarks>
+        /// <para>
+        /// With <c>L1 = a1 + b1 x</c>, <c>L2 = a2 + b2 x</c> and <c>D = a1 b2 - a2 b1</c>, the
+        /// substitution makes <c>x = (a1 - a2 t)/(b2 t - b1)</c>, <c>L2 = D/(b2 t - b1)</c> and
+        /// <c>dx = -D dt/(b2 t - b1)^2</c>, so <c>L1^A L2^B P(x) dx</c> is, up to a factor constant
+        /// on each interval, <c>-D^(1 - k) t^A S(t) dt</c> for the polynomial
+        /// <c>S(t) = sum_j p_j (a1 - a2 t)^j (b2 t - b1)^(k - 2 - j)</c>. Its antiderivative is
+        /// <c>sum_i s_i t^(A + i + 1)/(A + i + 1)</c>, and since <c>t^(A + i + 1)</c> is
+        /// <c>t^A t^(i + 1)</c> for a whole <c>i + 1</c>, the constant factor and <c>t^A</c> come
+        /// back as the two powers as written: the answer is
+        /// <c>L1^A L2^B sum_i -D^(1 - k) s_i L1^(i + 1) L2^(k - 1 - i)/(A + i + 1)</c>, exact
+        /// wherever the integrand is defined. The generic case, as the integrator answers
+        /// everywhere: <c>A + i + 1</c> is taken to be nonzero.
+        /// </para>
+        /// <para>
+        /// Only the sum <c>-2</c> was answered, as the derivative of the product of the two powers
+        /// raised by one; Rubi's <c>(a + b x)^m (c + d x)^n</c> files hold the rest, and under
+        /// <c>u = sin(y)</c> so do the powers of <c>1 ± sin(y)</c> beside one another. A closed
+        /// rule, so at any depth. https://github.com/asc-community/AngouriMath/issues/718
+        /// </para>
+        /// </remarks>
+        internal static Entity? SolveTwoLinearPowersWhoseExponentsSumToAWholeNumber(Entity expr, Entity.Variable x)
+        {
+            Entity constant = Number.Integer.One;
+            var powers = new List<(Entity Base, Entity Exponent)>();
+            var polynomial = new List<(Entity Factor, bool Underneath)>();
+            foreach (var (factor, underneath) in FactorsOfTheIntegrand(expr))
+            {
+                if (!factor.ContainsNode(x))
+                {
+                    constant = underneath ? constant / factor : constant * factor;
+                    continue;
+                }
+                var (@base, exponent) = factor is Powf(var powerBase, var power) ? (powerBase, power) : (factor, (Entity)Number.Integer.One);
+                if (exponent.ContainsNode(x))
+                    return null;
+                if (underneath)
+                    exponent = (-exponent).InnerSimplified;
+                var index = powers.FindIndex(pair => pair.Base == @base);
+                if (index >= 0)
+                    powers[index] = (@base, (powers[index].Exponent + exponent).InnerSimplified);
+                else if (exponent is Number.Integer)
+                    polynomial.Add((factor, underneath));
+                else
+                    powers.Add((@base, exponent));
+            }
+            if (powers.Count != 2)
+                return null;
+            // A whole power of a base that turned out to be one of the two goes into its exponent.
+            for (var i = polynomial.Count - 1; i >= 0; i--)
+            {
+                var (factor, underneath) = polynomial[i];
+                var (@base, exponent) = factor is Powf(var powerBase, Number.Integer power) ? (powerBase, (Entity)power) : (factor, (Entity)Number.Integer.One);
+                var index = powers.FindIndex(pair => pair.Base == @base);
+                if (index < 0)
+                {
+                    if (underneath)
+                        return null;
+                    continue;
+                }
+                powers[index] = (@base, (powers[index].Exponent + (underneath ? -exponent : exponent)).InnerSimplified);
+                polynomial.RemoveAt(i);
+            }
+            var ((first, a), (second, b)) = (powers[0], powers[1]);
+            if (a is Number.Integer || b is Number.Integer)
+                return null;
+            if ((a + b).Simplify() is not Number.Integer { EInteger: var sum } || sum.CompareTo(EInteger.FromInt32(-2)) > 0
+                || sum.CompareTo(EInteger.FromInt32(-MaximumLinearPowerSum)) < 0)
+                return null;
+            var k = -sum.ToInt32Checked();
+            if (!TreeAnalyzer.TryGetPolyLinear(first, x, out var b1, out var a1) || !TreeAnalyzer.TryGetPolyLinear(second, x, out var b2, out var a2))
+                return null;
+            var d = (a1 * b2 - a2 * b1).InnerSimplified;
+            if (d.Evaled is Number.Complex { IsZero: true })
+                return null;
+            Entity product = Number.Integer.One;
+            foreach (var (factor, _) in polynomial)
+                product *= factor;
+            if (!TreeAnalyzer.TryGetPolynomial(product, x, out var coefficients)
+                || coefficients.Keys.Any(power => power.Sign < 0 || power.CompareTo(EInteger.FromInt32(k - 2)) > 0))
+                return null;
+            var t = Variable.CreateUnique(expr, "t_quotient");
+            // Whole powers written out, so that what is read as a polynomial in t is one.
+            static Entity Raised(Entity @base, int power) => power switch
+            {
+                0 => Number.Integer.One,
+                1 => @base,
+                _ => MathS.Pow(@base, power),
+            };
+            Entity inT = Number.Integer.Zero;
+            foreach (var pair in coefficients)
+            {
+                var j = pair.Key.ToInt32Checked();
+                inT += pair.Value * Raised(a1 - a2 * t, j) * Raised(b2 * t - b1, k - 2 - j);
+            }
+            if (!TreeAnalyzer.TryGetPolynomial(inT.InnerSimplified, t, out var ofT))
+                return null;
+            Entity polynomialInX = Number.Integer.Zero;
+            foreach (var pair in ofT)
+            {
+                var (i, s) = (pair.Key.ToInt32Checked(), pair.Value);
+                if (i < 0 || i > k - 2)
+                    return null;
+                polynomialInX += s / (a + (i + 1)) * Raised(first, i + 1) * Raised(second, k - 1 - i);
+            }
+            var answer = (-constant * MathS.Pow(d, 1 - k)).InnerSimplified * MathS.Pow(first, a) * MathS.Pow(second, b) * polynomialInX.InnerSimplified;
+            return answer.Nodes.Any(node => node == MathS.NaN) ? null : answer;
+        }
+
+        /// <summary>
+        /// The most negative sum of the two exponents
+        /// <see cref="SolveTwoLinearPowersWhoseExponentsSumToAWholeNumber"/> takes: the polynomial
+        /// it writes has <c>k - 1</c> terms, each a product of powers summing to <c>k</c>.
+        /// </summary>
+        private const int MaximumLinearPowerSum = 12;
+
         /// <summary>
         /// The constant <c>k</c> with <paramref name="left"/> equal to <c>k</c> times
         /// <paramref name="right"/>, read off one monomial with both written as polynomials

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

@@ -762,6 +762,9 @@ private static Entity Normalized(Entity expr, Entity.Variable x) =>
             // of the powers raised by one -- `e^x x^2 ln(x)^2 (3 + (3 + x) ln(x))` -- read off
             // the sum, for symbolic exponents too. Before the split, which loses it.
             if ((answer = IndefiniteIntegralSolver.SolveAsTheDerivativeOfAProductOfPowers(expr, x)) is { }) return answer;
+            // Powers of two linears whose exponents sum to a whole number below -2, beside a polynomial of
+            // low degree, by t = L1/L2: a power of t beside a polynomial. The rule above reads the sum -2.
+            if ((answer = IndefiniteIntegralSolver.SolveTwoLinearPowersWhoseExponentsSumToAWholeNumber(expr, x)) is { }) return answer;
             // A rational function of x and of a tower of exponentials and logarithms over it,
             // by the Risch-Norman ansatz: a rational function of the same plus logarithms of
             // the denominator's factors, the coefficients unknown and solved for. Before the

Sources/Tests/UnitTests/Calculus/TwoLinearPowersIntegralTest.cs

@@ -0,0 +1,50 @@
+//
+// 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>
+    /// Powers of two linears whose exponents sum to a whole number below <c>-2</c>, beside a
+    /// polynomial of degree at most two less than its negative: under <c>t = L1/L2</c> a power of
+    /// <c>t</c> beside a polynomial, so the antiderivative is the two powers times a polynomial.
+    /// Rubi's <c>(a + b x)^m (c + d x)^n</c> files.
+    /// <a href="https://github.com/asc-community/AngouriMath/issues/718">#718</a>
+    /// </summary>
+    [Trait("Area", "Calculus")]
+    public sealed class TwoLinearPowersIntegralTest
+    {
+        [Theory]
+        [InlineData("(a + b*x)^m*(c + d*x)^(-3 - m)")]
+        [InlineData("x^(n - 4)/(a + b*x)^n")]
+        [InlineData("(a + b*x)^m*(c + d*x)^(-5 - m)*(p + q*x)^3")]
+        [InlineData("(a + b*x)*(c + d*x)^(n - 4)/(p + q*x)^n")]
+        [InlineData("(a + b*x)^m*(c + d*x)^(-4 - m)*(p + q*x)*(g + h*x)")]
+        [InlineData("(a + b*x)^(5/2)/(c + d*x)^(11/2)")]
+        public void UnderTheQuotientOfTheLinears(string integrand)
+        {
+            var integral = integrand.ToEntity().Integrate("x");
+            Assert.DoesNotContain("integral(", integral.Stringize());
+            Entity Pinned(Entity e) => e.Substitute("a", 1.3).Substitute("b", 0.9).Substitute("c", 0.7)
+                .Substitute("d", 1.1).Substitute("g", 0.8).Substitute("h", 0.6).Substitute("m", 0.37)
+                .Substitute("n", 1.21).Substitute("p", 0.4).Substitute("q", 1.3);
+            var derivative = Pinned(integral.Substitute("C", 0)).Differentiate("x");
+            var original = Pinned(integrand.ToEntity());
+            foreach (var at in new[] { -2.4, -1.5, 0.3, 0.8, 1.5, 2.4 })
+            {
+                var want = original.Substitute("x", at).EvalNumerical();
+                var got = derivative.Substitute("x", at).EvalNumerical();
+                Assert.True(Math.Abs((double)(got - want).RealPart) + Math.Abs((double)(got - want).ImaginaryPart)
+                            < 1e-9 * Math.Max(1, Math.Abs((double)want.RealPart) + Math.Abs((double)want.ImaginaryPart)),
+                    $"d/dx of the antiderivative of {integrand} is {got} at x = {at}, where the integrand is {want}");
+            }
+        }
+    }
+}