Split a power of x over a sum of two cubes over the cube roots of unity

BREAKING-CHANGES.md

@@ -447,6 +447,19 @@ 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 x over a sum of two cubes is split over the cube roots of unity
+
+**Answers where there were none.** `x^k/((3 a + b x^2)^3 + 27 a^2 c x^3)` was declined: the integrator's resultant needs
+numbers for coefficients. A sextic `d6 (x^2 + s)^3 + d3 x^3` is now split over the cube roots of unity into three
+quadratics, `x^2 + w m x + s` with `m^3 = d3/d6`. For two of them `w` is complex, so the answer carries `i` the way
+Rubi's does ([#718](https://github.com/asc-community/AngouriMath/issues/718)).
+
+| Input | Was (2.5.0) | Now |
+|---|---|---|
+| `"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(...)` | 849 characters |
+| `"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 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

@@ -19588,6 +19588,95 @@ void Add(int j, int k, Entity[] poly)
             return Functions.PartialFractions.DerivativeHoldsAtSampledPoints(answer, expr, x) ? answer : null;
         }
 
+        /// <summary>
+        /// A power of <c>x</c> over a sextic that is a sum of two cubes,
+        /// <c>d6 ((x^2 + s)^3 + q x^3)</c>, split over the cube roots of unity:
+        /// <c>1/(U^3 + V^3) = sum_j w_j/(3 V^2 (U + w_j V))</c> for <c>U = x^2 + s</c>,
+        /// <c>V = m x</c>, <c>m^3 = q</c>, so that <c>x^k/(U^3 + V^3)</c> is
+        /// <c>sum_j w_j x^(k - 2)/(3 m^2 (x^2 + w_j m x + s))</c>, three quadratics.
+        /// </summary>
+        /// <remarks>
+        /// The sextic is read as <c>d6 x^6 + d4 x^4 + d3 x^3 + d2 x^2 + d0</c> with
+        /// <c>d4^2 = 3 d6 d2</c> and <c>d2^2 = 3 d4 d0</c>, which is what makes
+        /// <c>d6 x^6 + d4 x^4 + d2 x^2 + d0</c> the cube <c>d6 (x^2 + s)^3</c> with
+        /// <c>s = d4/(3 d6)</c>. The identity is the partial fractions of <c>1/(U^3 + V^3)</c> in
+        /// <c>U</c>, whose poles <c>-w_j V</c> have residues <c>1/(3 w_j^2 V^2) = w_j/(3 V^2)</c>; any
+        /// cube root <c>m</c> of <c>q</c> gives it. The three quadratics are one integral
+        /// <c>x^n/(x^2 + B x + S)</c> asked with <c>B</c> and <c>S</c> named. Rubi's 1.3.1
+        /// <c>x^k/(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)</c> was declined,
+        /// the integrator's resultant needing numbers for coefficients.
+        /// https://github.com/asc-community/AngouriMath/issues/718
+        /// </remarks>
+        internal static Entity? SolveAPowerOfXOverASumOfTwoCubes(Entity expr, Entity.Variable x, bool integrateByParts)
+        {
+            if (!Integration.AnsweringTheQuestionAsked || !expr.Vars.Any(symbol => symbol != x))
+                return null;
+            var (numerator, denominator) = Functions.SingleQuotient.Of(expr);
+            Entity? sextic = null;
+            var power = 0;
+            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;
+                    }
+                    if (factor == x || factor is Powf(var bx, Number.Integer { EInteger: var wx }) && bx == x && wx.CanFitInInt32())
+                    {
+                        var whole = factor == x ? 1 : ((Number.Integer)((Powf)factor).Exponent).EInteger.ToInt32Unchecked();
+                        power += below ? -whole : whole;
+                        continue;
+                    }
+                    if (below && sextic is null)
+                    {
+                        sextic = factor;
+                        continue;
+                    }
+                    return null;
+                }
+            if (sextic is null || Functions.PartialFractions.HoldsTheImaginaryUnit(sextic)
+                || !TreeAnalyzer.TryGetPolynomial(sextic, x, out var read) || read.Keys.Any(key => key.Sign < 0 || key.CompareTo(EInteger.FromInt32(6)) > 0))
+                return null;
+            Entity At(int degree) => read.TryGetValue(EInteger.FromInt32(degree), out var coefficient) ? coefficient.InnerSimplified : Number.Integer.Zero;
+            var (d6, d5, d4, d3, d2, d1, d0) = (At(6), At(5), At(4), At(3), At(2), At(1), At(0));
+            if (!VanishesIdentically(d5) || !VanishesIdentically(d1)
+                || VanishesIdentically(d6) || VanishesIdentically(d4) || VanishesIdentically(d3) || VanishesIdentically(d0)
+                || !VanishesIdentically((d4 * d4 - 3 * d6 * d2).Expand()) || !VanishesIdentically((d2 * d2 - 3 * d4 * d0).Expand()))
+                return null;
+            var shift = Functions.PartialFractions.InLowestTermsOverTheSymbols(d4 / (3 * d6));
+            var cube = Functions.PartialFractions.InLowestTermsOverTheSymbols(d3 / d6);
+            var root = MathS.Pow(cube, Number.Rational.Create(1, 3));
+            var namedB = Variable.CreateUnique(expr, "B_c");
+            var namedS = Variable.CreateUnique(expr + namedB, "S_c");
+            var term = MathS.Pow(x, Number.Integer.Create(power - 2)) / (MathS.Pow(x, Number.Integer.Create(2)) + namedB * x + namedS);
+            if (Integration.ComputeAsAQuestionOfItsOwn(term, x, integrateByParts) is not { } each)
+                return null;
+            // The integral over a real quadratic splits on the sign of 4 S - B^2, which means nothing
+            // for the complex B of two of the three: its first arm, the arctangent over
+            // sqrt(4 S - B^2), differentiates back to the integrand for any B with 4 S != B^2, and
+            // the answer is checked below.
+            each = each.Replace(node => node switch
+            {
+                Piecewise piecewise when piecewise.Cases.Any() => piecewise.Cases.First().Expression,
+                Providedf provided => provided.Expression,
+                _ => node,
+            });
+            var half = Number.Rational.Create(1, 2);
+            var unities = new Entity[]
+            {
+                Number.Integer.One,
+                -half + half * MathS.Sqrt(3) * MathS.i,
+                -half - half * MathS.Sqrt(3) * MathS.i,
+            };
+            Entity answer = Number.Integer.Zero;
+            foreach (var unity in unities)
+                answer = answer + unity * each.Substitute(namedB, unity * root).Substitute(namedS, shift);
+            answer = (constant / (3 * d6 * MathS.Pow(root, Number.Integer.Create(2))) * answer).InnerSimplified;
+            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

@@ -678,6 +678,9 @@ private static Entity WithItsPiecewiseSumsCombined(Entity expr)
             // 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;
+            // A power of x over a sextic that is a sum of two cubes, split over the cube roots of
+            // unity into three quadratics before the partial fractions try the symbolic sextic.
+            if ((answer = IndefiniteIntegralSolver.SolveAPowerOfXOverASumOfTwoCubes(expr, x, integrateByParts)) is { }) return answer;
             // A root or a power of a product beside a power of one of its factors, the product
             // split: `(a + b x^2)^(5/2)/sqrt(a^2 - b^2 x^4)` is `K (a + b x^2)^2/sqrt(a - b x^2)`.
             // Early, since what reads the product whole is slower and longer: Rubi's 1.2.1.2:2218,

Sources/Tests/UnitTests/Calculus/ASumOfTwoCubesIntegralTest.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 power of <c>x</c> over a sextic that is a sum of two cubes,
+    /// <c>(3 a + b x^2)^3 + 27 a^2 c x^3</c>, split over the cube roots of unity into three
+    /// quadratics. Rubi's 1.3.1 rows; the cube root of <c>27 a^2 c/b^3</c> is the principal one,
+    /// complex for a negative <c>c</c>, and the three quadratics are still the factors.
+    /// <a href="https://github.com/asc-community/AngouriMath/issues/718">#718</a>
+    /// </summary>
+    [Trait("Area", "Calculus")]
+    public sealed class ASumOfTwoCubesIntegralTest
+    {
+        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("x^4/(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)")]
+        [InlineData("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)")]
+        [InlineData("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)")]
+        [InlineData("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))")]
+        public void IsAnsweredAndDifferentiatesBack(string integrand)
+        {
+            DifferentiatesBack(integrand, ("a", 2), ("b", 3), ("c", 5));
+            DifferentiatesBack(integrand, ("a", -1), ("b", 2), ("c", -3));
+        }
+    }
+}