A special function beside an elementary factor is integrated by parts, and so is its square

BREAKING-CHANGES.md

@@ -2295,6 +2295,26 @@ 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 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
+`Si(b x)` against `x sin(b x)`, whose integral needs parts of its own, and the integral of that
+factor was taken with parts off. It is integrated by the polynomial's parts now, beside a special
+function of a linear argument, and what is left, `sin(b x)/x` times sines and cosines, by parts in
+turn. A square of one of `b x` is two rounds of parts, and what the first round leaves comes back
+as one product over a sum -- `(b x Ei(b x) - e^(b x)) e^(b x)/(b x)` for `Ei(b x)^2` -- that no rule
+reads whole: its terms are asked one at a time now
+([#1501](https://github.com/asc-community/AngouriMath/issues/1501)). An integrand holding one of
+these functions had no reading in 2.5.0, which the entries for the functions themselves record, and
+no integrand without one changes: Rubi's independent suites and a sample of its families 1 to 7,
+2726 problems, are answered alone as they were.
+
+| Input | Was (2.5.0) | Now |
+|---|---|---|
+| `"x*Si(b*x)*sin(b*x)".Integrate("x")` | a polynomial in a variable `Si` | an antiderivative in `Si(b x)`, `Ci(2 b x)` and `ln(x)` |
+| `"Ei(b*x)^2".Integrate("x")` | `Ei * (b * x) ^ 3 / 3 / b + C`, with `Ei` a variable | an antiderivative in `Ei(b x)` and `Ei(2 b x)` |
+| `"x*erf(b*x)^2".Integrate("x")` | `UnrecognizedFunctionParseException`: there is no function `erf` | an antiderivative in `erf(b x)` |
+
 ### An inverse trigonometric function below the bar is integrated to the sine and cosine integrals
 
 `1/arcsin(x)` was left unintegrated. Under the substitution that undoes the inverse function,

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

@@ -1722,6 +1722,19 @@ static Entity WithTheConstantMatchedTo(Entity antiderivative, Entity derivativeO
                 // a coefficient's sign whose conditions are decided, `0 = 0`, and a remainder
                 // holding that piecewise is read by nothing.
                 var integralOfU = Integration.ComputeIndefiniteIntegral(u, x, false)?.InnerSimplified;
+                // Beside a special function of a linear argument, a polynomial times something
+                // else is integrated by the polynomial's parts, which end in as many steps as its
+                // degree: `x sin(b x)` beside `Si(b x)`, which nothing answers with parts off.
+                // https://github.com/asc-community/AngouriMath/issues/1501
+                var besideASpecialFunction = IsASpecialFunctionOfALinearOrAPowerOfOne(v, x);
+                var byThePolynomialsParts = false;
+                if (integralOfU is null && besideASpecialFunction
+                    && ThePolynomialFactor(u, x) is var (polynomialBeside, restBeside)
+                    && polynomialBeside is not null && restBeside is not null)
+                {
+                    integralOfU = IntegrateByPartsPolynomial(polynomialBeside, restBeside, x)?.InnerSimplified;
+                    byThePolynomialsParts = integralOfU is not null;
+                }
                 if (integralOfU is null) return null;
 
                 // Differentiate v
@@ -1780,7 +1793,12 @@ static Entity WithTheConstantMatchedTo(Entity antiderivative, Entity derivativeO
                 // https://github.com/asc-community/AngouriMath/issues/1501
                 var partsOnTheRemainder = remaining.Nodes.Count() < wholeSize
                     || (remainingPower >= 1 && remainingPower < wholePower)
-                    || IsASpecialFunctionOfALinear(v, x) && MathS.TryPolynomial(u, x, out _);
+                    || IsASpecialFunctionOfALinear(v, x) && MathS.TryPolynomial(u, x, out _)
+                    // And where what was integrated beside it was a polynomial times an
+                    // elementary function, whose integral the polynomial's parts wrote: the
+                    // remainder is then the special function's derivative times polynomials and
+                    // elementary functions, and parts on it descend on their degrees.
+                    || byThePolynomialsParts;
                 // Spelled as the product its own next step reads, where there is one. `Simplify`
                 // writes `2 arctan(x)/(1 + x^2) * x^2/2` as `arctan(x) x^2/(x^2 + 1)`, a
                 // quotient, and this rule runs on a product. `x*arctan(x)^2` is answered by
@@ -1804,6 +1822,46 @@ static Entity WithTheConstantMatchedTo(Entity antiderivative, Entity derivativeO
                 // search above it carry on, thirty seconds where the gate leaves it at one.
                 var remainingIntegral = (Integration.AnsweringTheQuestionAsked ? SolveByEulerSubstitution(remaining, x) : null)
                     ?? Integration.ComputeIndefiniteIntegral(remaining, x, partsOnTheRemainder);
+                // After a special function, or a power of one, the remainder is its derivative
+                // times what was integrated, written as one product over a sum: for `Ei(b x)^2`,
+                // `(b x Ei(b x) - e^(b x)) e^(b x)/(b x)`, which no rule reads, where its two
+                // terms are `e^(b x) Ei(b x)` and `e^(2 b x)/(b x)`, each answered when asked.
+                // So its terms are asked, each as a question of its own. From the top only,
+                // which is where that asking reaches, and once: a term asked so is at the top
+                // itself, and its own remainder asked the same way would lift the search again,
+                // each level with the depth reset. And of a multiple of x only: of `a + b x`, the
+                // remainder is over the linear, and its terms are the exponentials a hyperbolic
+                // function is written as, each a harder question than the whole -- the decline of
+                // `x Shi(a + b x)^2` went from a second to past twenty asking them.
+                // https://github.com/asc-community/AngouriMath/issues/1501
+                if (remainingIntegral is null && besideASpecialFunction && OfAMultipleOfTheVariable(v, x)
+                    && Integration.AnsweringTheQuestionAsked && !askingTheTermsOfARemainder)
+                {
+                    var terms = remaining is Divf(var remainingAbove, var remainingBelow)
+                        ? DistributedOverTheSum(remainingAbove).Select(term => term / remainingBelow).ToList()
+                        : DistributedOverTheSum(remaining);
+                    askingTheTermsOfARemainder = true;
+                    try
+                    {
+                        foreach (var term in terms)
+                        {
+                            if ((Integration.ComputeIndefiniteIntegral(term, x, integrateByParts: true)
+                                    ?? Integration.ComputeAsAQuestionOfItsOwn(term, x, integrateByParts: true)) is not { } termIntegral)
+                            {
+                                remainingIntegral = null;
+                                break;
+                            }
+                            remainingIntegral = remainingIntegral is null ? termIntegral : remainingIntegral + termIntegral;
+                        }
+                    }
+                    finally
+                    {
+                        askingTheTermsOfARemainder = false;
+                    }
+                }
+                else if (remainingIntegral is null && besideASpecialFunction && OfAMultipleOfTheVariable(v, x)
+                    && Integration.AnsweringTheQuestionAsked)
+                    Integration.DeclinedForItsScope();
                 if (remainingIntegral is null) return null;
 
                 return v * integralOfU - remainingIntegral;
@@ -20155,10 +20213,49 @@ Entity Term(Entity? sum, Entity coefficient, int power)
         /// <a href="https://github.com/asc-community/AngouriMath/issues/1501">#1501</a>: the error
         /// functions and the exponential, logarithmic, sine, cosine and hyperbolic integrals.
         /// </summary>
+        /// <summary>
+        /// Set while the terms of a remainder after a special function are asked as questions of
+        /// their own, so that the asking is not nested; see <c>TryIntegrateByPartsOnce</c>.
+        /// </summary>
+        [System.ThreadStatic] private static bool askingTheTermsOfARemainder;
+
         private static bool IsASpecialFunction(Entity node)
             => node is Entity.Erff or Entity.Erfcf or Entity.Erfif or Entity.Eif or Entity.Lif
                 or Entity.Sif or Entity.Cif or Entity.Shif or Entity.Chif;
 
+        /// <summary>
+        /// A special function of a linear argument, as <see cref="IsASpecialFunctionOfALinear"/>
+        /// reads one, or a positive whole power of one: <c>Si(b x)^2</c>.
+        /// </summary>
+        private static bool IsASpecialFunctionOfALinearOrAPowerOfOne(Entity factor, Variable x)
+            => IsASpecialFunctionOfALinear(factor, x)
+                || factor is Powf(var @base, Number.Integer power) && power.EInteger.Sign > 0 && IsASpecialFunctionOfALinear(@base, x);
+
+        /// <summary>
+        /// Whether the special function in <paramref name="factor"/>, or in the base of its power,
+        /// is of a multiple of <paramref name="x"/>, <c>b x</c>, with no offset.
+        /// </summary>
+        private static bool OfAMultipleOfTheVariable(Entity factor, Variable x)
+            => (factor is Powf(var @base, _) ? @base : factor).DirectChildren.FirstOrDefault() is { } argument
+                && TreeAnalyzer.TryGetPolyLinear(argument, x, out _, out var offset) && TreeAnalyzer.IsZero(offset);
+
+        /// <summary>
+        /// The factors of <paramref name="expr"/> that are polynomials in <paramref name="x"/> of
+        /// positive degree, against the rest: <c>x sin(b x)</c> is <c>x</c> and <c>sin(b x)</c>.
+        /// Both halves <see langword="null"/> where either would be empty.
+        /// </summary>
+        private static (Entity? Polynomial, Entity? Others) ThePolynomialFactor(Entity expr, Variable x)
+        {
+            Entity? polynomial = null;
+            Entity? rest = null;
+            foreach (var factor in Mulf.LinearChildren(expr))
+                if (factor.ContainsNode(x) && MathS.TryPolynomial(factor, x, out _))
+                    polynomial = polynomial is null ? factor : polynomial * factor;
+                else
+                    rest = rest is null ? factor : rest * factor;
+            return polynomial is null || rest is null || !rest.ContainsNode(x) ? (null, null) : (polynomial, rest);
+        }
+
         /// <summary>
         /// Whether <paramref name="expr"/> has what the derivative of the special function
         /// <paramref name="special"/> is made of: an exponential of a quadratic in

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

@@ -412,6 +412,13 @@ private static Entity Normalized(Entity expr, Entity.Variable x) =>
         /// </summary>
         [System.ThreadStatic] private static bool descentTruncated;
 
+        /// <summary>
+        /// Marks what is being worked out as declined for the scope it was asked in rather than
+        /// for its mathematics, as running out of descent is, so that the decline is not
+        /// remembered: the same question asked where that scope does not bind may be answered.
+        /// </summary>
+        internal static void DeclinedForItsScope() => descentTruncated = true;
+
         /// <summary>
         /// The integrals this thread is part-way through, so that asking for one again while it is
         /// still being worked out is recognised as a cycle rather than followed round again.

Sources/Tests/UnitTests/Calculus/SpecialFunctionsByPartsTest.cs

@@ -95,5 +95,43 @@ public void TheLogarithmicIntegralIsByParts(string integrand)
         [InlineData("x^3*Chi(b*x)")]
         public void APowerTimesASpecialFunctionIsByParts(string integrand)
             => DifferentiatesBack(integrand, Parameters);
+
+        /// <summary>
+        /// Beside a power of <c>x</c> times the elementary function its derivative is made of, the
+        /// special function is still the factor differentiated: <c>x sin(b x)</c> is integrated by
+        /// its polynomial's parts, and what is left, <c>sin(b x)/x</c> times that, is products of
+        /// sines and cosines over powers of <c>x</c>.
+        /// </summary>
+        [Theory]
+        [InlineData("x*Si(b*x)*sin(b*x)")]
+        [InlineData("x^3*Si(b*x)*sin(b*x)")]
+        [InlineData("x^2*Ci(b*x)*cos(b*x)")]
+        [InlineData("x*Ci(b*x)*sin(b*x)")]
+        [InlineData("x*Si(a + b*x)*sin(a + b*x)")]
+        [InlineData("x*Si(c + d*x)*sin(a + b*x)")]
+        public void BesideAPowerAndTheElementaryFactorOfItsDerivative(string integrand)
+            => DifferentiatesBack(integrand, Parameters);
+
+        /// <summary>
+        /// A square of a special function of <c>b x</c> is two rounds of parts: the first against
+        /// the power of <c>x</c> leaves the special function once, beside its derivative, and that
+        /// is the case above or a substitution. The remainder of a round is asked term by term,
+        /// since it comes back as one product over a sum -- <c>(b x Ei(b x) - e^(b x)) e^(b x)/(b x)</c>
+        /// for <c>Ei(b x)^2</c> -- that no rule reads whole.
+        /// </summary>
+        [Theory]
+        [InlineData("Ei(b*x)^2")]
+        [InlineData("x*Ei(b*x)^2")]
+        [InlineData("x^2*Ei(b*x)^2")]
+        [InlineData("x*Si(b*x)^2")]
+        [InlineData("Ci(b*x)^2")]
+        [InlineData("x*Ci(b*x)^2")]
+        [InlineData("x^2*erf(b*x)^2")]
+        [InlineData("x*erfc(b*x)^2")]
+        [InlineData("erfi(b*x)^2/x^3")]
+        [InlineData("x*Shi(b*x)^2")]
+        [InlineData("Chi(b*x)^2")]
+        public void ASquareIsTwoRoundsOfParts(string integrand)
+            => DifferentiatesBack(integrand, Parameters);
     }
 }