Split a sum of two powers of even quadratics over the roots of its binomial

BREAKING-CHANGES.md

@@ -596,6 +596,20 @@ unevaluated, is answered by this route as well.
 | `"(a*x+2*b*n*ln(c*x^n))/(a*x^2+b*x*ln(c*x^n)^2)".ToEntity().Integrate("x")` | `integral(...)` | 37 characters |
 | `"tanh(a+2*ln(x))^2/x^3".ToEntity().Integrate("x")` | `integral(...)` | 7,796 characters |
 
+### A quotient over a sum of two powers of even quadratics is split over the roots of its binomial
+
+**Answers where there were none.** `(1 + x^2)^2/(a x^6 + b (1 + x^2)^3)` was left unevaluated. The rule that splits
+`w^m/(a + b w^n)` over the `n` roots of the binomial, for an even quadratic `w`, took a constant beside the power
+only. A second power of an even quadratic `v` is read as well. `w^m v^(n - 1 - m)/(a v^n + b w^n)` is the same
+sum in `t = w/v`, and each root's `1/(w - r v)` is again an even quadratic
+([#718](https://github.com/asc-community/AngouriMath/issues/718)).
+
+| Input | Was (2.5.0) | Now |
+|---|---|---|
+| `"(1+x^2)^2/(a*x^6+b*(1+x^2)^3)".ToEntity().Integrate("x")` | `integral(...)` | 878 characters |
+| `"x^2*(1+x^2)/(a*x^6+b*(1+x^2)^3)".ToEntity().Integrate("x")` | `integral(...)` | 896 characters |
+| `"(1+x^2)/(a*x^4+b*(1+x^2)^2)".ToEntity().Integrate("x")` | `integral(...)` | 289 characters |
+
 ### A symbolic cubic with a repeated root is written over its roots
 
 **Answers where there were none.** `1/(x^3 - 3b x + 2b^(3/2))` was left unevaluated: its coefficients are symbolic, so

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

@@ -20066,6 +20066,14 @@ Piecewise piecewise when piecewise.Cases.Any() => piecewise.Cases.First().Expres
         /// <c>P</c> named, and splits on the sign of <c>P q</c>, which means nothing for a complex
         /// <c>P</c>: its first arm, the arctangent, differentiates back for any <c>P q != 0</c>, and
         /// the answer is checked below. Rubi's 1.3.1 <c>(1 - x^2)/(a + b (1 - x^2)^4)</c> was declined.
+        /// <para>
+        /// Two powers of even quadratics, <c>w^m v^(n - 1 - m)/(a v^n + b w^n)</c>, are the same sum in
+        /// <c>t = w/v</c>: the integrand is <c>t^m/(v (a + b t^n))</c>, and <c>1/(v (t - r_j))</c> is
+        /// <c>1/(w - r_j v)</c>, again an even quadratic <c>P + Q x^2</c>, whose integral
+        /// <c>arctan(sqrt(Q) x/sqrt(P))/(sqrt(P) sqrt(Q))</c> differentiates back for any <c>P</c> and
+        /// <c>Q</c> off zero. Rubi's 1.3.1 <c>(1 + x^2)^2/(a x^6 + b (1 + x^2)^3)</c>, a sum of two cubes,
+        /// was declined.
+        /// </para>
         /// https://github.com/asc-community/AngouriMath/issues/718
         /// </remarks>
         internal static Entity? SolveAPowerOfAnEvenQuadraticOverABinomialInIt(Entity expr, Entity.Variable x, bool integrateByParts)
@@ -20081,58 +20089,98 @@ Piecewise piecewise when piecewise.Cases.Any() => piecewise.Cases.First().Expres
                 {
                     if (!factor.ContainsNode(x))
                         constant = under ? constant / factor : constant * factor;
-                    else if (!under && above is null)
-                        above = factor;
-                    else if (under && below is null)
+                    else if (!under)
+                        above = above is null ? factor : above * factor;
+                    else if (below is null)
                         below = factor;
                     else
                         return null;
                 }
             if (below is null)
                 return null;
-            // a + b w^n: one term in x, a whole power of an even quadratic, beside terms free of it.
-            Entity? inner = null;
-            Entity scale = Number.Integer.One;
+            // a + b w^n: one term in x, a whole power of an even quadratic, beside terms free of it; or
+            // a v^n + b w^n, two such powers, a power of x^2 among them written as one of x.
+            var powers = new List<(Entity Scale, Entity Base, int Degree)>();
             Entity free = Number.Integer.Zero;
-            var degree = 0;
             foreach (var term in Sumf.LinearChildren(below))
             {
                 if (!term.ContainsNode(x))
                 {
                     free = free + term;
                     continue;
                 }
-                if (inner is not null)
+                if (powers.Count == 2)
                     return null;
+                Entity termScale = Number.Integer.One;
+                Entity? termBase = null;
+                var termDegree = 0;
                 foreach (var factor in Mulf.LinearChildren(term))
                 {
                     if (!factor.ContainsNode(x))
-                        scale = scale * factor;
-                    else if (inner is null && factor is Powf(var w, Number.Integer { EInteger: var n }) && n.CompareTo(EInteger.FromInt32(2)) >= 0
-                        && n.CompareTo(EInteger.FromInt32(8)) <= 0)
+                        termScale = termScale * factor;
+                    else if (termBase is null && factor is Powf(var w, Number.Integer { EInteger: var n }) && n.CompareTo(EInteger.FromInt32(2)) >= 0
+                        && n.CompareTo(EInteger.FromInt32(16)) <= 0)
                     {
-                        inner = w;
-                        degree = n.ToInt32Unchecked();
+                        termBase = w;
+                        termDegree = n.ToInt32Unchecked();
                     }
                     else
                         return null;
                 }
+                if (termBase is null)
+                    return null;
+                powers.Add((termScale, termBase, termDegree));
+            }
+            if (powers.Count == 2)
+            {
+                if (!VanishesIdentically(free))
+                    return null;
+                // x^(2 n) beside v^n is (x^2)^n.
+                for (var k = 0; k < 2; k++)
+                    if (powers[k].Base == x && powers[k].Degree % 2 == 0 && powers[1 - k].Base != x)
+                        powers[k] = (powers[k].Scale, MathS.Pow(x, Number.Integer.Create(2)), powers[k].Degree / 2);
             }
-            if (inner is null || VanishesIdentically(free) || Functions.PartialFractions.HoldsTheImaginaryUnit(expr)
+            if (powers.Count == 0 || powers.Any(power => power.Degree > 8) || powers.Count == 2 && powers[0].Degree != powers[1].Degree)
+                return null;
+            var (scale, inner, degree) = powers[0];
+            Entity? other = powers.Count == 2 ? powers[1].Base : null;
+            if (other is not null)
+                free = powers[1].Scale;
+            Entity otherP = Number.Integer.One, otherQ = Number.Integer.Zero;
+            if (VanishesIdentically(free) || Functions.PartialFractions.HoldsTheImaginaryUnit(expr)
                 || !TreeAnalyzer.TryGetPolyQuadratic(inner, x, out var q, out var linearTerm, out var p)
                 || !VanishesIdentically(linearTerm) || VanishesIdentically(q))
                 return null;
-            // The numerator a constant multiple of w^m, m below n.
+            if (other is not null)
+            {
+                if (!TreeAnalyzer.TryGetPolyQuadratic(other, x, out var readQ, out var otherLinear, out var readP)
+                    || !VanishesIdentically(otherLinear) || VanishesIdentically(readQ))
+                    return null;
+                (otherP, otherQ) = (readP, readQ);
+            }
+            // The numerator a constant multiple of w^m v^(n - 1 - m), m below n, which is w^m alone
+            // for a binomial: then t^m/(V (a + b t^n)) for t = w/v, split over the roots in t, and
+            // 1/(V (t - r)) is 1/(w - r v).
             var power = -1;
             Entity multiple = Number.Integer.One;
-            if (above is null)
+            if (above is null && other is null)
                 power = 0;
-            else
-                for (var m = 1; m < degree && power < 0; m++)
+            else if (above is not null)
+                for (var m = other is null ? 1 : 0; m < degree && power < 0; m++)
                 {
                     var wm = MathS.Pow(inner, Number.Integer.Create(m));
-                    if (!TreeAnalyzer.TryGetPolynomial(above, x, out var aboveRead) || !TreeAnalyzer.TryGetPolynomial(wm, x, out var wmRead)
-                        || !aboveRead.TryGetValue(EInteger.FromInt32(2 * m), out var lead) || !wmRead.TryGetValue(EInteger.FromInt32(2 * m), out var wLead))
+                    if (other is not null)
+                    {
+                        // No zeroth powers: InnerSimplified writes w^0 as 1 provided w is not 0.
+                        var rest = degree - 1 - m;
+                        wm = m == 0 ? MathS.Pow(other, Number.Integer.Create(rest))
+                            : rest == 0 ? wm : wm * MathS.Pow(other, Number.Integer.Create(rest));
+                    }
+                    var top = other is null ? 2 * m : 2 * (degree - 1);
+                    // Expanded for two powers: the reader takes a sum of monomials, and v^k is not one.
+                    var (aboveAsRead, wmAsRead) = other is null ? (above, wm) : (above.Expand(), wm.Expand());
+                    if (!TreeAnalyzer.TryGetPolynomial(aboveAsRead, x, out var aboveRead) || !TreeAnalyzer.TryGetPolynomial(wmAsRead, x, out var wmRead)
+                        || !aboveRead.TryGetValue(EInteger.FromInt32(top), out var lead) || !wmRead.TryGetValue(EInteger.FromInt32(top), out var wLead))
                         continue;
                     var ratio = (lead / wLead).InnerSimplified;
                     if (!VanishesIdentically((above - ratio * wm).Expand()))
@@ -20144,16 +20192,20 @@ Piecewise piecewise when piecewise.Cases.Any() => piecewise.Cases.First().Expres
                 return null;
             var a = free.InnerSimplified;
             var b = scale.InnerSimplified;
+            Entity? each = null;
             var namedP = Variable.CreateUnique(expr, "P_b");
-            if (Integration.ComputeAsAQuestionOfItsOwn(1 / (namedP + q * MathS.Pow(x, Number.Integer.Create(2))), x, integrateByParts) is not { } each)
-                return null;
-            // See the remarks: the arctangent arm for every root.
-            each = each.Replace(node => node switch
+            if (other is null)
             {
-                Piecewise piecewise when piecewise.Cases.Any() => piecewise.Cases.First().Expression,
-                Providedf provided => provided.Expression,
-                _ => node,
-            });
+                if (Integration.ComputeAsAQuestionOfItsOwn(1 / (namedP + q * MathS.Pow(x, Number.Integer.Create(2))), x, integrateByParts) is not { } asked)
+                    return null;
+                // See the remarks: the arctangent arm for every root.
+                each = asked.Replace(node => node switch
+                {
+                    Piecewise piecewise when piecewise.Cases.Any() => piecewise.Cases.First().Expression,
+                    Providedf provided => provided.Expression,
+                    _ => node,
+                });
+            }
             var radius = MathS.Pow(Functions.PartialFractions.InLowestTermsOverTheSymbols(-a / b), Number.Rational.Create(1, degree));
             Entity sum = Number.Integer.Zero;
             for (var j = 0; j < degree; j++)
@@ -20162,7 +20214,20 @@ Piecewise piecewise when piecewise.Cases.Any() => piecewise.Cases.First().Expres
                 // the real root, and (-1)^(2/3) is 1 (#204).
                 var turn = 2 * MathS.pi * Number.Rational.Create(j, degree);
                 Entity root = j == 0 ? radius : radius * (MathS.Cos(turn) + MathS.i * MathS.Sin(turn)).InnerSimplified;
-                sum = sum + MathS.Pow(root, Number.Integer.Create(power + 1)) * each.Substitute(namedP, p - root);
+                Entity over;
+                if (each is not null)
+                    over = each.Substitute(namedP, p - root);
+                else
+                {
+                    // 1/(P + Q x^2) for w - r v: arctan(sqrt(Q) x/sqrt(P))/(sqrt(P) sqrt(Q)), whose
+                    // derivative is 1/(P + Q x^2) for any P, Q off zero, the two roots squaring back.
+                    var atP = Functions.PartialFractions.InLowestTermsOverTheSymbols(p - root * otherP);
+                    var atQ = Functions.PartialFractions.InLowestTermsOverTheSymbols(q - root * otherQ);
+                    if (VanishesIdentically(atP) || VanishesIdentically(atQ))
+                        return null;
+                    over = MathS.Arctan(MathS.Sqrt(atQ) * x / MathS.Sqrt(atP)) / (MathS.Sqrt(atP) * MathS.Sqrt(atQ));
+                }
+                sum = sum + MathS.Pow(root, Number.Integer.Create(power + 1)) * over;
             }
             var answer = (-constant * multiple / (degree * a) * sum).InnerSimplified;
             return Functions.PartialFractions.DerivativeHoldsAtSampledPoints(answer, expr, x) ? answer : null;

Sources/Tests/UnitTests/Calculus/ASumOfTwoPowersOfEvenQuadraticsIntegralTest.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 quotient over a sum of two whole powers of even quadratics, <c>a v^n + b w^n</c>, split over the
+    /// roots of <c>a + b t^n</c> in <c>t = w/v</c>: Rubi's 1.3.1 <c>(1 + x^2)^2/(a x^6 + b (1 + x^2)^3)</c>,
+    /// which was left unevaluated.
+    /// <a href="https://github.com/asc-community/AngouriMath/issues/718">#718</a>
+    /// </summary>
+    [Trait("Area", "Calculus")]
+    public sealed class ASumOfTwoPowersOfEvenQuadraticsIntegralTest
+    {
+        private static void DifferentiatesBack(string integrand, params (string, double)[] pins)
+        {
+            var integral = integrand.ToEntity().Integrate("x");
+            Assert.DoesNotContain("integral(", integral.Stringize());
+
+            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("(1+x^2)^2/(a*x^6+b*(1+x^2)^3)")]
+        [InlineData("x^2*(1+x^2)/(a*x^6+b*(1+x^2)^3)")]
+        [InlineData("x^4/(a*x^6+b*(1+x^2)^3)")]
+        [InlineData("(2+3*x^2)^2/(a*x^6+b*(2+3*x^2)^3)")]
+        [InlineData("(1+x^2)/(a*x^4+b*(1+x^2)^2)")]
+        public void IsAnsweredAndDifferentiatesBack(string integrand)
+        {
+            DifferentiatesBack(integrand, ("a", 2), ("b", 3));
+            DifferentiatesBack(integrand, ("a", -5), ("b", 0.5));
+        }
+    }
+}