Integrate a sum of powers of an inverse or error function by parts down the powers

BREAKING-CHANGES.md

@@ -6394,6 +6394,29 @@ no integrand without one changes: Rubi's independent suites and a sample of its
 | `"e^(b*x)*Ei(b*x)/x".Integrate("x")` | `Ei * e ^ (b * x) + C`, with `Ei` a variable | `Ei(b * x) ^ 2 / 2 + C` |
 | `"Si(b*x)*sin(b*x)/x".Integrate("x")` | `Si * -cos(b * x) + C`, with `Si` a variable | `Si(b * x) ^ 2 / 2 + C` |
 
+### A sum of powers of an inverse trigonometric or error function is integrated by parts down the powers
+
+**Answers where there were none.** `x^3/((1 + a^2 x^2) atan(a x)^3) - 3 x^2/(2 a atan(a x)^2)`,
+`x asin(x)^(3/2)/(1 - x^2)^2 - 3 x/(8 (1 - x^2) sqrt(asin(x)))` and `erf(b x) e^(-b^2 x^2) (1/x^3 + b^2/x)` with its
+`erfc` and `erfi` twins, Rubi's 5.1.4, 5.3.4 and 8.1, were declined: no term is elementary alone, and the split and
+the parts each read one term at a time. A sum of terms each a power of one function `K`, whose derivative is free of
+it, times a cofactor free of it is now integrated by parts from the top power,
+`integral A K^p = F K^p - integral p F K^(p - 1) K'` with `F` the antiderivative of `A`, and the last integral joins
+the terms a power below. Where every power cancels so, the
+parts are the answer, `(-x^3/(2 a atan(a x)^2))'` for the first; where what is left is free of `K` it is asked. The
+powers a whole number apart and at most five apart, a cancelling group decided at sampled points, and the answer
+differentiated back. Three sums of the same kind the unreleased master answered come back shorter
+([#233](https://github.com/asc-community/AngouriMath/issues/233)).
+
+| Input | Was (2.5.0) | Now |
+|---|---|---|
+| `"x^3/((1+a^2*x^2)*arctan(a*x)^3)-3/2*x^2/(a*arctan(a*x)^2)".ToEntity().Integrate("x")` | `integral(...)` | `-x^3/(2 a arctan(a x)^2)`, 51 characters |
+| `"x*arcsin(x)^(3/2)/(1-x^2)^2-3/8*x/((1-x^2)*sqrt(arcsin(x)))".ToEntity().Integrate("x")` | `integral(...)` | 99 characters |
+| `"erf(b*x)/(e^(b^2*x^2)*x^3)+b^2*erf(b*x)/(e^(b^2*x^2)*x)".ToEntity().Integrate("x")` | `UnrecognizedFunctionParseException`: there is no function `erf` | 219 characters |
+| `"erfi(b*x)/(e^(b^2*x^2)*x^3)+b^2*erfi(b*x)/(e^(b^2*x^2)*x)".ToEntity().Integrate("x")` | `integral(...)`, reading `erfi` as a symbol | 131 characters |
+| `"erf(x)^2+4*x*erf(x)*e^(-x^2)/sqrt(pi)".ToEntity().Integrate("x")` | `UnrecognizedFunctionParseException`; 352 characters on the unreleased master | `x erf(x)^2`, 18 characters |
+| `"2*x/arccos(x)+x^2/(sqrt(1-x^2)*arccos(x)^2)".ToEntity().Integrate("x")` | `integral(...)`; 121 characters on the unreleased master | `x^2/arccos(x)`, 33 characters |
+
 ### A square of a special function, and one beside a power of `x` and its elementary derivative, are integrated by parts
 
 `x Si(b x) sin(b x)` was left unevaluated where `Si(b x) sin(b x)` was not: by parts it is

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

@@ -15991,6 +15991,137 @@ static bool More(Entity difference)
             return changed ? rebuilt / bottom : null;
         }
 
+        /// <summary>
+        /// A sum of terms each a power of one function <c>K</c> whose derivative is free of it -- an
+        /// inverse trigonometric function or an error function -- times a cofactor free of it,
+        /// integrated by parts down the powers: <c>∫ A K^p = F K^p - ∫ p F K^(p - 1) K'</c> for
+        /// <c>F = ∫ A</c>, the cofactor of the top power, and the last integral joins the terms a
+        /// power below. Where every power cancels so, the parts are the antiderivative:
+        /// <c>x^3/((1 + a^2 x^2) atan(a x)^3) - 3 x^2/(2 a atan(a x)^2)</c> is
+        /// <c>(-x^3/(2 a atan(a x)^2))'</c>, in one step, and
+        /// <c>x asin(x)^(3/2)/(1 - x^2)^2 - 3 x/(8 (1 - x^2) sqrt(asin(x)))</c> in two. Where what is
+        /// left is free of <c>K</c>, it is asked: <c>erf(b x) e^(-b^2 x^2) (1/x^3 + b^2/x)</c> leaves
+        /// <c>b e^(-2 b^2 x^2)/(sqrt(pi) x^2)</c>.
+        /// </summary>
+        /// <remarks>
+        /// <para>
+        /// The terms are not elementary apart -- <c>∫ x/((1 - x^2) sqrt(asin(x)))</c> is
+        /// <c>∫ tan(t)/sqrt(t) dt</c> -- so the split declined each, and by parts reads one term
+        /// at a time. Rubi's 5.1.4, 5.3.4 and 8.1 write such sums to test that the parts meet.
+        /// </para>
+        /// <para>
+        /// One such function in the sum, every term a rational power of it times a cofactor free
+        /// of it, the powers a whole number apart, at most five apart, and each step's
+        /// <c>F</c> free of <c>K</c>. A group that cancels is decided at sampled points, and the
+        /// answer is checked by differentiating it back. A remainder free of <c>K</c> is asked
+        /// only for a top power of two terms or more with a cofactor among them that is not
+        /// algebraic, which is where the cofactors integrate together and not apart: a sum of
+        /// algebraic cofactors, <c>a + b asin(c x)</c>, is what the split answers term by term,
+        /// and this rule would only respell it. At the question asked only, and not inside its
+        /// own questions, which are free of <c>K</c> anyway.
+        /// </para>
+        /// https://github.com/asc-community/AngouriMath/issues/233
+        /// </remarks>
+        internal static Entity? SolveASumOfPowersOfAnInverseOrErrorFunctionByPartsDownThePowers(Entity expr, Entity.Variable x, bool integrateByParts)
+        {
+            if (!Integration.AnsweringTheQuestionAsked || expr is not (Sumf or Minusf) || descendingThePowers)
+                return null;
+            Entity? kernel = null;
+            foreach (var node in expr.Nodes)
+                if (node is Arcsinf or Arccosf or Arctanf or Arccotanf or Erff or Erfcf or Erfif && node.ContainsNode(x))
+                {
+                    if (kernel is null)
+                        kernel = node;
+                    else if (node != kernel)
+                        return null;
+                }
+            if (kernel is null)
+                return null;
+            var slope = kernel.Differentiate(x).InnerSimplified;
+            if (slope.ContainsNode(kernel))
+                return null;
+            Number.Rational? PowerOfTheKernel(Entity factor)
+                => factor == kernel ? Number.Integer.One
+                   : factor is Powf(var @base, Number.Rational power) && PowerOfTheKernel(@base) is { } inner ? inner * power
+                   : null;
+            var terms = new List<(Number.Rational Power, Entity Cofactor)>();
+            foreach (var term in Sumf.LinearChildren(expr))
+            {
+                Number.Rational power = Number.Integer.Zero;
+                Entity cofactor = Number.Integer.One;
+                foreach (var factor in Mulf.LinearChildren(term))
+                    if (!factor.ContainsNode(kernel))
+                        cofactor *= factor;
+                    else if (PowerOfTheKernel(factor) is { } read)
+                        power += read;
+                    else
+                        return null;
+                terms.Add((power, cofactor));
+            }
+            var top = terms.Select(term => term.Power).Aggregate((l, r) => Above(r, l) ? r : l);
+            var lowest = terms.Select(term => term.Power).Aggregate((l, r) => Above(l, r) ? r : l);
+            if (terms.Any(term => !(top - term.Power).ERational.Denominator.Equals(EInteger.One))
+                || Above(top - lowest, Number.Integer.Create(5)))
+                return null;
+            var groups = new List<(Number.Rational Power, Entity Cofactor)>();
+            foreach (var (power, cofactor) in terms)
+                Join(groups, power, cofactor);
+            var atTheTop = terms.Where(term => Same(term.Power, top)).ToList();
+            var mayAskWhatIsLeft = atTheTop.Count >= 2 && atTheTop.Any(term => !IsAlgebraicIn(term.Cofactor, x));
+            if (groups.Count < 2 && !mayAskWhatIsLeft)
+                return null;
+            descendingThePowers = true;
+            try
+            {
+                Entity parts = Number.Integer.Zero;
+                for (var step = 0; step < 8; step++)
+                {
+                    groups.RemoveAll(group => Functions.PartialFractions.HoldsAtSampledPoints(group.Cofactor, Number.Integer.Zero, x));
+                    if (groups.Count == 0)
+                        return Functions.PartialFractions.DerivativeHoldsAtSampledPoints(parts, expr, x) ? parts : null;
+                    var (power, cofactor) = groups.Aggregate((l, r) => Above(r.Power, l.Power) ? r : l);
+                    if (power.ERational.IsZero)
+                    {
+                        // Free of K: asked where the top power's cofactors integrate together only.
+                        if (groups.Count > 1 || !mayAskWhatIsLeft
+                            || Integration.ComputeAsAQuestionOfItsOwn(cofactor, x, integrateByParts) is not { } rest)
+                            return null;
+                        var answer = parts + rest;
+                        return Functions.PartialFractions.DerivativeHoldsAtSampledPoints(answer, expr, x) ? answer : null;
+                    }
+                    // Below zero the powers never reach a remainder free of K, so below the lowest
+                    // written one nothing is left to cancel.
+                    if (power.ERational.Sign < 0 && Above(lowest, power)
+                        || Integration.ComputeAsAQuestionOfItsOwn(cofactor, x, integrateByParts) is not { } integrated
+                        || integrated.ContainsNode(kernel))
+                        return null;
+                    parts += integrated * (Same(power, Number.Integer.One) ? kernel : MathS.Pow(kernel, power));
+                    groups.RemoveAll(group => Same(group.Power, power));
+                    Join(groups, power - Number.Integer.One, -power * integrated * slope);
+                    for (var i = 0; i < groups.Count; i++)
+                        groups[i] = (groups[i].Power, groups[i].Cofactor.InnerSimplified);
+                }
+                return null;
+            }
+            finally
+            {
+                descendingThePowers = false;
+            }
+            static bool Same(Number.Rational a, Number.Rational b) => a.ERational.CompareTo(b.ERational) == 0;
+            static bool Above(Number.Rational a, Number.Rational b) => a.ERational.CompareTo(b.ERational) > 0;
+            static void Join(List<(Number.Rational Power, Entity Cofactor)> groups, Number.Rational power, Entity cofactor)
+            {
+                var at = groups.FindIndex(group => Same(group.Power, power));
+                if (at < 0)
+                    groups.Add((power, cofactor));
+                else
+                    groups[at] = (power, groups[at].Cofactor + cofactor);
+            }
+        }
+
+        /// <summary>Set while <see cref="SolveASumOfPowersOfAnInverseOrErrorFunctionByPartsDownThePowers"/> asks its questions.</summary>
+        [System.ThreadStatic] private static bool descendingThePowers;
+
         /// <summary>
         /// A product of powers times a sum that is the derivative of the product with some of
         /// the powers raised by one: <c>e^x x^2 ln(x)^2 (3 + (3 + x) ln(x))</c> is

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

@@ -932,6 +932,11 @@ private static Entity WithItsPiecewiseSumsCombined(Entity expr)
             // recurrence that raises the whole power to -1: closed, where the substitution for the
             // root makes a rational function over a power of a symbolic quadratic.
             if ((answer = IndefiniteIntegralSolver.SolveANegativePowerOfALinearBesideAPowerOfAnotherByTheRecurrence(expr, x, integrateByParts)) is { }) return answer;
+            // A sum of powers of one inverse trigonometric or error function beside cofactors free
+            // of it, by parts down the powers from the top: `x^3/((1 + a^2 x^2) atan(a x)^3) -
+            // 3 x^2/(2 a atan(a x)^2)` is the derivative of one product, and neither term is
+            // elementary on its own. Before the split, which loses it.
+            if ((answer = IndefiniteIntegralSolver.SolveASumOfPowersOfAnInverseOrErrorFunctionByPartsDownThePowers(expr, x, integrateByParts)) is { }) return answer;
             // A product of powers times a sum that is the derivative of the product with some
             // 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.

Sources/Tests/UnitTests/Calculus/ASumOfPowersOfAnInverseOrErrorFunctionByPartsDownThePowersIntegralTest.cs

@@ -0,0 +1,60 @@
+//
+// 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 sum of powers of one inverse trigonometric or error function beside cofactors free of
+    /// it, by parts down the powers from the top. Its terms are not elementary apart, so the
+    /// split declined each, and the sum was declined or ran out of time. Rubi's 5.1.4, 5.3.4
+    /// and 8.1.
+    /// <a href="https://github.com/asc-community/AngouriMath/issues/233">#233</a>
+    /// </summary>
+    [Trait("Area", "Calculus")]
+    public sealed class ASumOfPowersOfAnInverseOrErrorFunctionByPartsDownThePowersIntegralTest
+    {
+        [Theory]
+        // The parts cancel every power below the top: one step, and two.
+        [InlineData("x^3/((1+a^2*x^2)*arctan(a*x)^3)-3/2*x^2/(a*arctan(a*x)^2)")]
+        [InlineData("x*arcsin(x)^(3/2)/(1-x^2)^2-3/8*x/((1-x^2)*sqrt(arcsin(x)))")]
+        [InlineData("1/((1-x^2)*arcsin(x)^2)-x/((1-x^2)^(3/2)*arcsin(x))")]
+        [InlineData("2*x/arccos(x)+x^2/(sqrt(1-x^2)*arccos(x)^2)")]
+        [InlineData("erf(x)^2+4*x*erf(x)*e^(-x^2)/sqrt(pi)")]
+        // What is left is free of the function, and the cofactors integrate together only.
+        [InlineData("erf(b*x)/(e^(b^2*x^2)*x^3)+b^2*erf(b*x)/(e^(b^2*x^2)*x)")]
+        [InlineData("erfc(b*x)/(e^(b^2*x^2)*x^3)+b^2*erfc(b*x)/(e^(b^2*x^2)*x)")]
+        [InlineData("erfi(b*x)/(e^(b^2*x^2)*x^3)+b^2*erfi(b*x)/(e^(b^2*x^2)*x)")]
+        public void ByPartsDownThePowers(string integrand)
+        {
+            var integral = integrand.ToEntity().Integrate("x");
+            var text = integral.Stringize();
+            Assert.DoesNotContain("integral(", text);
+            Assert.True(text.Length < 2000, $"{text.Length} characters of answer for {integrand}");
+            Entity Pinned(Entity e) => e.Substitute("a", 1.3).Substitute("b", 0.7);
+            var derivative = Pinned(integral.Substitute("C", 0)).Differentiate("x");
+            var original = Pinned(integrand.ToEntity());
+            var compared = 0;
+            // Inside (0, 1), where the arcsine and the arccosine are real and positive.
+            foreach (var at in new[] { 0.15, 0.3, 0.45, 0.6, 0.75, 0.9 })
+            {
+                var want = original.Substitute("x", at).EvalNumerical();
+                var got = derivative.Substitute("x", at).EvalNumerical();
+                if (want.IsNaN)
+                    continue;
+                compared++;
+                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}");
+            }
+            Assert.True(compared >= 5, $"only {compared} points could be compared for {integrand}");
+        }
+    }
+}