diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md
index 5e39bfc78..90c8c587b 100644
--- a/BREAKING-CHANGES.md
+++ b/BREAKING-CHANGES.md
@@ -2335,6 +2335,24 @@ beside a power of `x`, the answer is written with the same constant: `x Ei(a + b
 | `"x*Shi(a+b*x)^2".Integrate("x")` | `a ^ 2 * Shi * x ^ 2 / 2 + a * b * 2 * Shi * x ^ 3 / 3 + b ^ 2 * Shi * x ^ 4 / 4 + C`, with `Shi` a variable | an antiderivative in `Shi(a + b x)` and `Ei(2 (a + b x))`, `Ei(-2 (a + b x))` |
 | `"x*Ei(a+b*x)^2".Integrate("x")` | `a ^ 2 * Ei * x ^ 2 / 2 + a * b * 2 * Ei * x ^ 3 / 3 + b ^ 2 * Ei * x ^ 4 / 4 + C`, with `Ei` a variable | an antiderivative in `Ei(a + b x)` and `Ei(2 (a + b x))` |
 
+### A square of a special function of a shifted argument is integrated by parts
+
+The remainder of a special function's square was asked term by term only for an argument `b x`.
+For `a + b x`, the remainder divided by the linear until the constant of integration was matched to
+it; with that matched, the terms are the case without the offset, and the asking is offered for any
+linear argument ([#1501](https://github.com/asc-community/AngouriMath/issues/1501)). The exception
+is an error function beside anything but itself, whose derivative is a Gaussian of the shifted
+argument, so a polynomial beside it stays in every term: `(c + d x) erf(a + b x)^2` is still
+declined. 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 are answered alone as they were.
+
+| Input | Was (2.5.0) | Now |
+|---|---|---|
+| `"Shi(a+b*x)^2".Integrate("x")` | `Shi * (a + b * x) ^ 3 / 3 / b + C`, with `Shi` a variable | an antiderivative in `Shi(a + b x)`, `Ei(2 (a + b x))` and `Ei(-2 (a + b x))` |
+| `"x^2*Ei(a+b*x)^2".Integrate("x")` | `a ^ 2 * Ei * x ^ 3 / 3 + a * b * 2 * Ei * x ^ 4 / 4 + b ^ 2 * Ei * x ^ 5 / 5 + C`, with `Ei` a variable | an antiderivative in `Ei(a + b x)` and `Ei(2 (a + b x))` |
+| `"erf(a+b*x)^2".Integrate("x")` | `UnrecognizedFunctionParseException`: there is no function `erf` | an antiderivative in `erf` |
+
 ### An exponential or a hyperbolic function over several linears is split into partial fractions over them
 
 `e^x/(x (x + 1))` was left unevaluated, where `e^x/x` and `e^x/(x + 1)` were each answered with
diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs
index 8694cc5d6..c4354b9a5 100644
--- a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs
+++ b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs
@@ -1842,12 +1842,10 @@ static Entity WithTheConstantMatchedTo(Entity antiderivative, Entity derivativeO
                 // 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.
+                // each level with the depth reset. And not for every argument: see
+                // TheRemainderIsAskedTermByTerm.
                 // https://github.com/asc-community/AngouriMath/issues/1501
-                if (remainingIntegral is null && besideASpecialFunction && OfAMultipleOfTheVariable(v, x)
+                if (remainingIntegral is null && besideASpecialFunction && TheRemainderIsAskedTermByTerm(v, u, x)
                     && Integration.AnsweringTheQuestionAsked && !askingTheTermsOfARemainder)
                 {
                     var terms = remaining is Divf(var remainingAbove, var remainingBelow)
@@ -1872,7 +1870,7 @@ static Entity WithTheConstantMatchedTo(Entity antiderivative, Entity derivativeO
                         askingTheTermsOfARemainder = false;
                     }
                 }
-                else if (remainingIntegral is null && besideASpecialFunction && OfAMultipleOfTheVariable(v, x)
+                else if (remainingIntegral is null && besideASpecialFunction && TheRemainderIsAskedTermByTerm(v, u, x)
                     && Integration.AnsweringTheQuestionAsked)
                     Integration.DeclinedForItsScope();
                 if (remainingIntegral is null) return null;
@@ -20345,12 +20343,29 @@ private static bool IsASpecialFunctionOfALinearOrAPowerOfOne(Entity factor, Vari
                 || 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.
+        /// Whether the remainder after differentiating <paramref name="v"/>, a special function of
+        /// a linear argument or a power of one, against <paramref name="u"/> is asked term by term
+        /// when nothing answers it whole. Of a multiple of <paramref name="x"/>, always. Of
+        /// <c>a + b x</c>, where the function's derivative divides by its argument -- <c>Ei</c>,
+        /// <c>Si</c>, <c>Ci</c>, <c>Shi</c>, <c>Chi</c> -- since the constant of the factor
+        /// integrated beside it is then matched to the linear and the terms are the case without
+        /// the offset; and where the function is integrated against itself, the square alone.
+        /// Not an error function of <c>a + b x</c> beside anything else: its derivative is a
+        /// Gaussian of the shifted argument rather than a quotient by it, a polynomial beside it
+        /// stays in every term, and asking them took the decline of <c>(c + d x) erf(a + b x)^2</c>
+        /// from six seconds to twenty-four, answering nothing.
+        /// https://github.com/asc-community/AngouriMath/issues/1501
         /// </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);
+        private static bool TheRemainderIsAskedTermByTerm(Entity v, Entity u, Variable x)
+        {
+            var special = v is Powf(var @base, _) ? @base : v;
+            if (special.DirectChildren.FirstOrDefault() is not { } argument
+                || !TreeAnalyzer.TryGetPolyLinear(argument, x, out _, out var offset))
+                return false;
+            return TreeAnalyzer.IsZero(offset)
+                || special is not (Entity.Erff or Entity.Erfcf or Entity.Erfif)
+                || u == special;
+        }
 
         /// <summary>
         /// The factors of <paramref name="expr"/> that are polynomials in <paramref name="x"/> of
diff --git a/Sources/Tests/UnitTests/Calculus/SpecialFunctionsByPartsTest.cs b/Sources/Tests/UnitTests/Calculus/SpecialFunctionsByPartsTest.cs
index 6df826679..c83e8806d 100644
--- a/Sources/Tests/UnitTests/Calculus/SpecialFunctionsByPartsTest.cs
+++ b/Sources/Tests/UnitTests/Calculus/SpecialFunctionsByPartsTest.cs
@@ -142,6 +142,10 @@ public void ASquareIsTwoRoundsOfParts(string integrand)
         [Theory]
         [InlineData("x*Shi(a + b*x)^2")]
         [InlineData("x*Ei(a + b*x)^2")]
+        [InlineData("Shi(a + b*x)^2")]
+        [InlineData("x*Si(a + b*x)^2")]
+        [InlineData("x^2*Ei(a + b*x)^2")]
+        [InlineData("erf(a + b*x)^2")]
         public void ASquareOfAShiftedArgumentIsTwoRoundsOfParts(string integrand)
             => DifferentiatesBack(integrand, Parameters);
     }