The special functions are substitutions, and powers of one base gather across the bar

BREAKING-CHANGES.md

@@ -2278,6 +2278,23 @@ elementary integrand whose antiderivative is reached through one is answered whe
 |---|---|---|
 | `"(a+b*ln(c*x^n))/(x^2*(d+e*ln(f*x^m)))".Integrate("x")`, Rubi's 3.1.5 row 213 | `integral(...)` | an antiderivative in `Ei`, provided `f > 0` and `e^d f^e > 0` |
 
+### A special function beside its derivative is a substitution
+
+`e^(c - b^2 x^2) erf(b x)^n`, `Ei(b x) e^(b x)/x` and `Si(b x) sin(b x)/x` are each a power of a
+special function beside its derivative, and `u = erf(b x)`, `u = Ei(b x)` or `u = Si(b x)` writes
+them as a power of `u`. None of the nine special functions was a candidate for that substitution.
+Each is one now, as the logarithm is, wherever its derivative can be the differential
+([#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 |
+|---|---|---|
+| `"e^(c - b^2*x^2)*erf(b*x)".Integrate("x")` | `UnrecognizedFunctionParseException`: there is no function `erf` | `e ^ c * pi ^ (1/2) * erf(b * x) ^ 2 / 2 / (2 * b) + C` |
+| `"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` |
+
 ### 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

@@ -1610,8 +1610,7 @@ private static bool IsDifferentiatedBeforeAPolynomial(Entity factor)
                 Logf or Entity.Arcsinf or Entity.Arccosf
                     or Entity.Arctanf or Entity.Arccotanf
                     or Entity.Arcsecantf or Entity.Arccosecantf => true,
-                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 => true,
+                _ when IsASpecialFunction(factor) => true,
                 // And a whole power of one, which is the same function for this purpose:
                 // differentiating ln(x)^2 gives 2ln(x)/x, whose x cancels against the integrated
                 // polynomial exactly as ln(x)'s does, leaving x*ln(x) -- one step simpler, and
@@ -19624,6 +19623,33 @@ private static bool HoldsAnEvenRootBesides(Entity expr, Entity radicand, Entity.
                             integrandInU = SimplifiedWithoutTheImaginaryUnit(cancelled, expr);
                         }
                     }
+                    // Powers of one base on the two sides of the bar: `e^(c - b^2 x^2)` over the
+                    // `e^(-b^2 x^2)` of du/dx is `e^c`, which the one-level simplification leaves as
+                    // written, and `e^(c - b^2 x^2) erf(b x)` was refused under `u = erf(b x)` for
+                    // the x in it, where `e^(-b^2 x^2) erf(b x)` was answered.
+                    // https://github.com/asc-community/AngouriMath/issues/1501
+                    // The divisor's factors are spread first, each to the power -1: gathered as it
+                    // stands, the whole divisor `2 e^(-(b x)^2) b` is one factor below the bar, with
+                    // a base of its own.
+                    // For a special function only, which is where the constant in the exponent comes
+                    // from: done for every candidate that left an x, it gathered the exponentials of
+                    // `1/((c + d x)^3 (a + a tanh(e + f x)))` and simplified them for every one, and
+                    // a decline in a third of a second became a timeout.
+                    if (integrandInU.ContainsNode(x) && IsASpecialFunction(u))
+                    {
+                        var (above, below) = Functions.SingleQuotient.Of(Functions.SingleQuotient.Combine(integrandInU));
+                        var spread = above;
+                        foreach (var factor in Mulf.LinearChildren(below))
+                            spread *= MathS.Pow(factor, -1);
+                        // Simplified a level, not only inner-simplified, since the gathered exponent
+                        // `c - (b x)^2 + (b x)^2` is left standing by the inner simplification -- and
+                        // only where a base was gathered at all. Simplified for every candidate that
+                        // left an x, it took `x cosh(a + b x) Shi(a + b x)` from a second to a timeout.
+                        var gatheredRaw = Patterns.GatherPowersOfOneBase(spread);
+                        if (!ReferenceEquals(gatheredRaw, spread)
+                            && SimplifiedWithoutTheImaginaryUnit(gatheredRaw, expr) is var gathered && !gathered.ContainsNode(x))
+                            integrandInU = gathered;
+                    }
                     // A polynomial in x left over under a candidate that is itself a polynomial
                     // is written in the candidate where it is one in it: `(1 - x)^2` is
                     // `1 - 2x + x^2`, and that is `u` for Apostol's `(1 - 2x + x^2)^(1/5)/(1 - x)`,
@@ -20124,6 +20150,34 @@ Entity Term(Entity? sum, Entity coefficient, int power)
         /// 2. f(x^n) * x^(n-1)  ->  u = x^n
         /// 3. f(g(x)) * g'(x)  ->  u = g(x)
         /// </summary>
+        /// <summary>
+        /// The special functions of
+        /// <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>
+        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>
+        /// Whether <paramref name="expr"/> has what the derivative of the special function
+        /// <paramref name="special"/> is made of: an exponential of a quadratic in
+        /// <paramref name="x"/> for the error functions, something in <paramref name="x"/> below the
+        /// bar for the exponential, trigonometric and hyperbolic integrals, whose derivatives
+        /// divide by their argument, and a logarithm below the bar for the logarithmic integral.
+        /// </summary>
+        private static bool TheDifferentialCanBeThere(Entity special, Entity expr, Entity.Variable x)
+        {
+            var below = Functions.SingleQuotient.Of(Functions.SingleQuotient.Combine(expr)).Denominator;
+            return special switch
+            {
+                Entity.Erff or Entity.Erfcf or Entity.Erfif => expr.Nodes.Any(node => node is Powf(var @base, var exponent)
+                    && !@base.ContainsNode(x) && TreeAnalyzer.TryGetPolyQuadratic(exponent, x, out var square, out _, out _) && !TreeAnalyzer.IsZero(square)),
+                Entity.Lif => below.Nodes.Any(node => node is Logf(_, var antilogarithm) && antilogarithm.ContainsNode(x)),
+                _ => below.ContainsNode(x),
+            };
+        }
+
         private static IEnumerable<Entity> FindSubstitutionCandidates(Entity expr, Entity.Variable x)
         {
             var candidates = new List<Entity>();
@@ -20234,6 +20288,18 @@ bool APowerCanBeExact(int k)
                         candidates.Add(node); // Logarithm itself (for cases like 1/(x*ln(x)))
                         if (antilog != x && antilog.ContainsNode(x)) candidates.Add(antilog); // Also add the argument if it's not just x
                         break;
+                    // A special function itself, as the logarithm is: each has an elementary
+                    // derivative, so beside a power of it that derivative is the differential.
+                    // `e^(c - b^2 x^2) erf(b x)^n` is `sqrt(pi) e^c/(2b) u^n` under `u = erf(b x)`.
+                    // https://github.com/asc-community/AngouriMath/issues/1501
+                    // Only where that derivative can be there: an exponential of a quadratic for
+                    // the error functions, a divisor in x for the u below the bar of e^u/u, sin(u)/u
+                    // and the others, a logarithm below the bar for li. Offered beside anything,
+                    // `x cosh(a + b x) Shi(a + b x)` went from a second to ten, simplifying quotients
+                    // of exponentials that were never going to lose their x.
+                    case var special when IsASpecialFunction(special) && TheDifferentialCanBeThere(special, expr, x):
+                        candidates.Add(node);
+                        break;
                     case Sumf(var aug, var add) when !rational && !large && node.Complexity <= LargestSumOffered:
                         if (aug.ContainsNode(x) || add.ContainsNode(x)) candidates.Add(node); // Linear expressions ax + b
                         break;

Sources/Tests/UnitTests/Calculus/SpecialFunctionSubstitutionTest.cs

@@ -0,0 +1,66 @@
+//
+// 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 System.Linq;
+using AngouriMath.Extensions;
+using Xunit;
+
+namespace AngouriMath.Tests.Calculus
+{
+    /// <summary>
+    /// A special function as the substitution: beside a power of it, its derivative is the
+    /// differential, so <c>e^(c - b^2 x^2) erf(b x)^n</c> is a power of <c>u = erf(b x)</c>, and
+    /// <c>Si(b x) sin(b x)/x</c> is <c>u du</c> under <c>u = Si(b x)</c>. The rows are Rubi's, from
+    /// its files 8.1, 8.3, 8.4 and 8.5.
+    /// https://github.com/asc-community/AngouriMath/issues/1501
+    /// </summary>
+    [Trait("Area", "Calculus")]
+    public sealed class SpecialFunctionSubstitutionTest
+    {
+        /// <summary>Off 0, where the reciprocals are undefined.</summary>
+        private static readonly double[] Points = { -1.7, -0.6, 0.35, 0.9, 1.45 };
+
+        private static readonly (string, string)[] Parameters = { ("b", "13/10"), ("c", "7/10") };
+
+        /// <summary>
+        /// Integrates, pins the parameters, and compares the derivative of the answer with the
+        /// integrand at <see cref="Points"/>. The parameters are pinned after integrating, so the
+        /// rule is asked the symbolic question.
+        /// </summary>
+        private static void DifferentiatesBack(string integrand)
+        {
+            var integral = integrand.ToEntity().Integrate("x");
+            Assert.DoesNotContain("integral(", integral.Stringize());
+            Entity Pinned(Entity e) => Parameters.Aggregate(e, (current, pin) => current.Substitute(pin.Item1, pin.Item2.ToEntity()));
+            var derivative = Pinned(integral.Substitute("C", 0)).Differentiate("x");
+            var original = Pinned(integrand.ToEntity());
+            foreach (var at in Points)
+            {
+                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-9,
+                    $"d/dx of the antiderivative of {integrand} is {got} at x = {at}, where the integrand is {want}");
+            }
+        }
+
+        [Theory]
+        [InlineData("e^(c - b^2*x^2)*erf(b*x)")]
+        [InlineData("e^(c - b^2*x^2)*erf(b*x)^2")]
+        [InlineData("e^(c - b^2*x^2)/erf(b*x)")]
+        [InlineData("e^(c - b^2*x^2)/erf(b*x)^2")]
+        [InlineData("e^(c - b^2*x^2)*erfc(b*x)^3")]
+        [InlineData("e^(c + b^2*x^2)*erfi(b*x)")]
+        [InlineData("e^(b*x)*Ei(b*x)/x")]
+        [InlineData("Si(b*x)*sin(b*x)/x")]
+        [InlineData("cos(b*x)*Ci(b*x)/x")]
+        public void ThePowerBesideTheDerivativeIsASubstitution(string integrand)
+            => DifferentiatesBack(integrand);
+    }
+}