Integrate a whole power of the square of a polynomial over its root

BREAKING-CHANGES.md

@@ -488,6 +488,21 @@ Rubi's does ([#718](https://github.com/asc-community/AngouriMath/issues/718)).
 | `"1/(27*a^3+27*a^2*b*x^2+27*a^2*c*x^3+9*a*b^2*x^4+b^3*x^6)".ToEntity().Integrate("x")` | `integral(...)` | 2,091 characters |
 | `"1/(x^2*(27*a^3+27*a^2*b*x^2+27*a^2*c*x^3+9*a*b^2*x^4+b^3*x^6))".ToEntity().Integrate("x")` | `integral(...)` | 3,820 characters |
 
+### A whole power of the square of a polynomial is integrated over its root
+
+**Answers where there were none.** `1/((d x)^(3/2) (a^2 + 2 a b x^2 + b^2 x^4)^2)` ran past any budget. A whole power
+`P^n` of a polynomial with symbols in it that is the square of a polynomial `R` is now written `R^(2 n)`, which the
+rules for a power of a binomial answer ([#718](https://github.com/asc-community/AngouriMath/issues/718)).
+
+| Input | Was (2.5.0) | Now |
+|---|---|---|
+| `"1/((d*x)^(3/2)*(a^2+2*a*b*x^2+b^2*x^4)^2)".ToEntity().Integrate("x")` | `integral(...)` | 965 characters |
+| `"1/((d*x)^(5/2)*(a^2+2*a*b*x^2+b^2*x^4)^2)".ToEntity().Integrate("x")` | `integral(...)` | 892 characters |
+| `"1/((d*x)^(7/2)*(a^2+2*a*b*x^2+b^2*x^4)^2)".ToEntity().Integrate("x")` | `integral(...)` | 1,055 characters |
+| `"1/((d*x)^(3/2)*(a^2+2*a*b*x^2+b^2*x^4)^3)".ToEntity().Integrate("x")` | `integral(...)` | 1,122 characters |
+| `"1/((d*x)^(5/2)*(a^2+2*a*b*x^2+b^2*x^4)^3)".ToEntity().Integrate("x")` | `integral(...)` | 1,044 characters |
+| `"1/((d*x)^(7/2)*(a^2+2*a*b*x^2+b^2*x^4)^3)".ToEntity().Integrate("x")` | `integral(...)` | 1,197 characters |
+
 ### A polynomial of degree six with a large leading coefficient is factored
 
 **Answers where there were none.** The factoring clears a leading coefficient by rescaling the variable, and declined

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

@@ -19862,6 +19862,63 @@ Piecewise piecewise when piecewise.Cases.Any() => piecewise.Cases.First().Expres
             return null;
         }
 
+        /// <summary>
+        /// A whole power of a polynomial that is the square of a polynomial in <c>x</c>, written as
+        /// twice that power of the root: <c>(a^2 + 2 a b x^2 + b^2 x^4)^n</c> is
+        /// <c>(a + b x^2)^(2 n)</c>, whatever the signs, since the power is whole.
+        /// </summary>
+        /// <remarks>
+        /// The rules for a power of a binomial read the root; the trinomial is the partial fractions'
+        /// whole. Rubi's 1.2.2.2 <c>1/((d x)^(3/2) (a^2 + 2 a b x^2 + b^2 x^4)^2)</c> ran past the
+        /// budget, and over <c>(a + b x^2)^4</c> it takes a twentieth of a second.
+        /// https://github.com/asc-community/AngouriMath/issues/718
+        /// </remarks>
+        internal static Entity? SolveAWholePowerOfASquareOverItsRoot(Entity expr, Entity.Variable x, bool integrateByParts)
+        {
+            if (!Integration.AnsweringTheQuestionAsked)
+                return null;
+            var (numerator, denominator) = Functions.SingleQuotient.Of(expr);
+            var changed = false;
+            Entity? Rewritten(Entity side)
+            {
+                Entity? product = null;
+                foreach (var factor in Mulf.LinearChildren(side))
+                {
+                    var written = factor;
+                    var (@base, whole) = factor is Powf(var b0, Number.Integer { EInteger: var w }) && w.Sign > 0 && w.CanFitInInt32()
+                        ? (b0, w.ToInt32Unchecked()) : (factor, 1);
+                    // With numbers for coefficients the partial fractions answer the trinomial.
+                    if (@base is Sumf && @base.ContainsNode(x) && @base.Vars.Any(symbol => symbol != x)
+                        && Sumf.LinearChildren(@base).Count() >= 3
+                        && TreeAnalyzer.TryGetPolynomial(@base, x, out var read) && read.Keys.All(key => key.Sign >= 0)
+                        && read.Keys.Any(key => key.CompareTo(EInteger.FromInt32(2)) >= 0))
+                    {
+                        var symbols = @base.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(@base, indices) is { } parsed && parsed.TrySquareRoot() is { } root
+                            && root.ToEntity(symbols) is var rootEntity && rootEntity.ContainsNode(x))
+                        {
+                            written = MathS.Pow(rootEntity, Number.Integer.Create(2 * whole));
+                            changed = true;
+                        }
+                    }
+                    product = product is null ? written : product * written;
+                }
+                return product;
+            }
+            var above = Rewritten(numerator);
+            var below = Rewritten(denominator);
+            if (!changed || above is null)
+                return null;
+            var asked = below is null || below == Number.Integer.One ? above : above / below;
+            if (Integration.ComputeAsAQuestionOfItsOwn(asked, x, integrateByParts) is not { } answer)
+                return null;
+            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:

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

@@ -673,6 +673,10 @@ private static Entity WithItsPiecewiseSumsCombined(Entity expr)
             // A rational function with numbers for coefficients over a denominator past the twelfth
             // degree, by Hermite's reduction before anything searches: nothing below answered one.
             if ((answer = IndefiniteIntegralSolver.SolveAHighDegreeRationalFunctionByHermite(expr, x)) is { }) return answer;
+            // A whole power of the square of a polynomial written out, written over its root: the rules
+            // for a power of a binomial read `(a + b x^2)^4` where `(a^2 + 2 a b x^2 + b^2 x^4)^2` ran
+            // past the budget.
+            if ((answer = IndefiniteIntegralSolver.SolveAWholePowerOfASquareOverItsRoot(expr, x, integrateByParts)) is { }) return answer;
             // 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;

Sources/Tests/UnitTests/Calculus/AWholePowerOfASquareIntegralTest.cs

@@ -0,0 +1,58 @@
+//
+// 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 whole power of a trinomial that is the square of a binomial in <c>x^2</c>,
+    /// <c>(a^2 + 2 a b x^2 + b^2 x^4)^n</c>, written as <c>(a + b x^2)^(2 n)</c>. Rubi's 1.2.2.2
+    /// rows beside <c>(d x)^(k/2)</c>, sampled where <c>d x</c> is positive and on both sides of
+    /// the root of <c>a + b x^2</c> when there is one.
+    /// <a href="https://github.com/asc-community/AngouriMath/issues/718">#718</a>
+    /// </summary>
+    [Trait("Area", "Calculus")]
+    public sealed class AWholePowerOfASquareIntegralTest
+    {
+        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[] { 0.3, 0.7, 1.1, 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/((d*x)^(3/2)*(a^2+2*a*b*x^2+b^2*x^4)^2)")]
+        [InlineData("1/((d*x)^(7/2)*(a^2+2*a*b*x^2+b^2*x^4)^2)")]
+        [InlineData("1/((d*x)^(5/2)*(a^2+2*a*b*x^2+b^2*x^4)^3)")]
+        public void IsAnsweredAndDifferentiatesBack(string integrand)
+        {
+            DifferentiatesBack(integrand, ("a", 2), ("b", 3), ("d", 0.7));
+            DifferentiatesBack(integrand, ("a", -2), ("b", 3), ("d", 1.5));
+        }
+    }
+}