Add five families of standard integrals (#233)

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

@@ -5,6 +5,8 @@
 // Website: https://am.angouri.org.
 //
 
+using PeterO.Numbers;
+
 namespace AngouriMath.Functions.Algebra
 {
     internal static class IntegralPatterns
@@ -19,17 +21,12 @@ internal static class IntegralPatterns
                 TreeAnalyzer.TryGetPolyLinear(arg, x, out var a, out _) =>
                     MathS.Sin(arg) / a,
 
-            // By power reduction: sin(u)^2 = (1 - cos(2u)) / 2, so the integral is
-            // x/2 - sin(2u)/(4a). Without this, integrating sin(x)^2 fell through to
-            // integration by parts and cycled there.
-            Entity.Powf(Entity.Sinf(var arg), Entity.Number.Integer(2)) when
-                TreeAnalyzer.TryGetPolyLinear(arg, x, out var a, out _) =>
-                    x / 2 - MathS.Sin(2 * arg) / (4 * a),
-
-            // cos(u)^2 = (1 + cos(2u)) / 2
-            Entity.Powf(Entity.Cosf(var arg), Entity.Number.Integer(2)) when
-                TreeAnalyzer.TryGetPolyLinear(arg, x, out var a, out _) =>
-                    x / 2 + MathS.Sin(2 * arg) / (4 * a),
+            // sin(u)^n * cos(u)^m for whole n and m, which covers sin(x)^2 and cos(x)^2
+            // as much as sin(x)^3 or sin(x)^2 * cos(x)^2.
+            _ when TryReadSineCosinePowers(expr, x, out var trigArg, out var sinePower, out var cosinePower)
+                   && sinePower + cosinePower >= 2
+                   && TreeAnalyzer.TryGetPolyLinear(trigArg, x, out var trigRate, out _) =>
+                       IntegrateSineCosinePowers(trigArg, sinePower, cosinePower, trigRate, x),
 
             Entity.Secantf(var arg) when
                 TreeAnalyzer.TryGetPolyLinear(arg, x, out var a, out _) =>
@@ -73,9 +70,247 @@ internal static class IntegralPatterns
                 && TreeAnalyzer.TryGetPolyQuadratic(denominator, x, out var a, out var b, out var c) // ∫ k/(ax^2 + bx + c) dx
                     => IntegrateRationalQuadratic(numerator, a, b, c, x),
 
+            // The inverse trigonometric functions, each of which is integration by parts
+            // against 1 -- a shape the by-parts solver does not look for, since there is no
+            // product to split.
+            Entity.Arcsinf(var arg) when
+                TreeAnalyzer.TryGetPolyLinear(arg, x, out var a, out _) =>
+                    (arg * MathS.Arcsin(arg) + MathS.Sqrt(1 - arg * arg)) / a,
+
+            Entity.Arccosf(var arg) when
+                TreeAnalyzer.TryGetPolyLinear(arg, x, out var a, out _) =>
+                    (arg * MathS.Arccos(arg) - MathS.Sqrt(1 - arg * arg)) / a,
+
+            Entity.Arctanf(var arg) when
+                TreeAnalyzer.TryGetPolyLinear(arg, x, out var a, out _) =>
+                    (arg * MathS.Arctan(arg) - MathS.Ln(MathS.Abs(1 + arg * arg)) / 2) / a,
+
+            Entity.Arccotanf(var arg) when
+                TreeAnalyzer.TryGetPolyLinear(arg, x, out var a, out _) =>
+                    (arg * MathS.Arccotan(arg) + MathS.Ln(MathS.Abs(1 + arg * arg)) / 2) / a,
+
+            // ∫ B^(px + q) * sin(mx + n) dx and its cosine twin. Integrating by parts
+            // twice returns the integral it started from, so the usual machinery cycles
+            // rather than terminating; solving that equation for the integral once gives
+            // the closed form below, which is what goes in the table.
+            Entity.Mulf(var exponential, Entity.Sinf(var wave)) when
+                IsExponentialRate(exponential, x, out var rate)
+                && TreeAnalyzer.TryGetPolyLinear(wave, x, out var frequency, out _) =>
+                    exponential * (rate * MathS.Sin(wave) - frequency * MathS.Cos(wave))
+                        / (rate * rate + frequency * frequency),
+
+            Entity.Mulf(Entity.Sinf(var wave), var exponential) when
+                IsExponentialRate(exponential, x, out var rate)
+                && TreeAnalyzer.TryGetPolyLinear(wave, x, out var frequency, out _) =>
+                    exponential * (rate * MathS.Sin(wave) - frequency * MathS.Cos(wave))
+                        / (rate * rate + frequency * frequency),
+
+            Entity.Mulf(var exponential, Entity.Cosf(var wave)) when
+                IsExponentialRate(exponential, x, out var rate)
+                && TreeAnalyzer.TryGetPolyLinear(wave, x, out var frequency, out _) =>
+                    exponential * (rate * MathS.Cos(wave) + frequency * MathS.Sin(wave))
+                        / (rate * rate + frequency * frequency),
+
+            Entity.Mulf(Entity.Cosf(var wave), var exponential) when
+                IsExponentialRate(exponential, x, out var rate)
+                && TreeAnalyzer.TryGetPolyLinear(wave, x, out var frequency, out _) =>
+                    exponential * (rate * MathS.Cos(wave) + frequency * MathS.Sin(wave))
+                        / (rate * rate + frequency * frequency),
+
+            // ∫ sqrt(ax^2 + bx + c) dx, which is one integration by parts away from the
+            // reciprocal form below and is written in terms of it.
+            Entity.Powf(var radicand, Entity.Number.Rational(Entity.Number.Integer(1), Entity.Number.Integer(2))) when
+                TreeAnalyzer.TryGetPolyQuadratic(radicand, x, out var qa, out var qb, out var qc)
+                && qa.Evaled is Entity.Number.Complex { IsZero: false }
+                    => IntegrateRootOfQuadratic(qa, qb, qc, radicand, x),
+
+            // ∫ k / sqrt(ax^2 + bx + c) dx -- the arcsine and logarithm forms. Without
+            // these, 1/sqrt(1 - x^2) had no antiderivative at all.
+            Entity.Divf(var numerator,
+                        Entity.Powf(var radicand, Entity.Number.Rational(Entity.Number.Integer(1), Entity.Number.Integer(2)))) when
+                !numerator.ContainsNode(x)
+                && TreeAnalyzer.TryGetPolyQuadratic(radicand, x, out var ra, out var rb, out var rc)
+                    => IntegrateOverRootOfQuadratic(numerator, ra, rb, rc, radicand, x),
+
             _ => null
         };
 
+
+        /// <summary>
+        /// Reads an expression as <c>sin(arg)^n * cos(arg)^m</c>, where either power may be
+        /// zero, or reports that it is not of that shape. Both factors have to be functions
+        /// of the same argument.
+        /// </summary>
+        private static bool TryReadSineCosinePowers(
+            Entity expr, Entity.Variable x, out Entity arg, out int sinePower, out int cosinePower)
+        {
+            arg = 0;
+            sinePower = cosinePower = 0;
+            switch (expr)
+            {
+                case Entity.Sinf(var a):
+                    arg = a; sinePower = 1; return a.ContainsNode(x);
+                case Entity.Cosf(var a):
+                    arg = a; cosinePower = 1; return a.ContainsNode(x);
+                case Entity.Powf(Entity.Sinf(var a), Entity.Number.Integer power)
+                    when power.EInteger.Sign > 0 && power.EInteger.CanFitInInt32():
+                    arg = a; sinePower = power.EInteger.ToInt32Checked(); return a.ContainsNode(x);
+                case Entity.Powf(Entity.Cosf(var a), Entity.Number.Integer power)
+                    when power.EInteger.Sign > 0 && power.EInteger.CanFitInInt32():
+                    arg = a; cosinePower = power.EInteger.ToInt32Checked(); return a.ContainsNode(x);
+                case Entity.Mulf(var left, var right)
+                    when TryReadSineCosinePowers(left, x, out var leftArg, out var leftSine, out var leftCosine)
+                         && TryReadSineCosinePowers(right, x, out var rightArg, out var rightSine, out var rightCosine)
+                         && leftArg == rightArg:
+                    arg = leftArg;
+                    sinePower = leftSine + rightSine;
+                    cosinePower = leftCosine + rightCosine;
+                    return true;
+                default:
+                    return false;
+            }
+        }
+
+        /// <summary>
+        /// The antiderivative of <c>sin(u)^n * cos(u)^m</c> where <c>u</c> is linear in x
+        /// with slope <paramref name="rate"/>.
+        /// </summary>
+        /// <remarks>
+        /// With an odd power there is a substitution: peel one factor off to be the
+        /// differential and write what is left in the other function, which turns the
+        /// integral into a polynomial. With both powers even there is no such factor to
+        /// peel, so the halved-angle identities go in instead and the result is integrated
+        /// again -- the total power halves each time, so this ends.
+        /// </remarks>
+        private static Entity IntegrateSineCosinePowers(
+            Entity arg, int sinePower, int cosinePower, Entity rate, Entity.Variable x)
+        {
+            if (sinePower == 0 && cosinePower == 0)
+                return x;
+
+            // cos^(2k+1) * sin^n: let s = sin(u), ds = cos(u) du
+            //   = (1/rate) * sum_j C(k, j) (-1)^j s^(n + 2j + 1) / (n + 2j + 1)
+            //
+            // Tried before the sine substitution, so that where both powers are odd the
+            // answer comes back in sin rather than cos: sin(x)cos(x) integrates to
+            // sin(x)^2/2, which is the form everyone writes, rather than the equally
+            // correct -cos(x)^2/2 that differs from it by a constant.
+            if (cosinePower % 2 == 1)
+            {
+                var k = (cosinePower - 1) / 2;
+                Entity sum = 0;
+                for (var j = 0; j <= k; j++)
+                {
+                    var power = sinePower + 2 * j + 1;
+                    sum += Binomial(k, j) * (j % 2 == 0 ? 1 : -1) * MathS.Sin(arg).Pow(power) / power;
+                }
+                return sum / rate;
+            }
+
+            // sin^(2k+1) * cos^m: let c = cos(u), dc = -sin(u) du
+            if (sinePower % 2 == 1)
+            {
+                var k = (sinePower - 1) / 2;
+                Entity sum = 0;
+                for (var j = 0; j <= k; j++)
+                {
+                    var power = cosinePower + 2 * j + 1;
+                    sum += Binomial(k, j) * (j % 2 == 0 ? 1 : -1) * MathS.Cos(arg).Pow(power) / power;
+                }
+                return -sum / rate;
+            }
+
+            // Both even: sin^2 = (1 - cos(2u))/2 and cos^2 = (1 + cos(2u))/2, expanded into
+            // powers of cos(2u), each of which is integrated by the same rules.
+            var p = sinePower / 2;
+            var q = cosinePower / 2;
+            var doubled = 2 * arg;
+            Entity result = 0;
+            // (1 - t)^p (1 + t)^q with t = cos(2u), over 2^(p+q)
+            for (var i = 0; i <= p; i++)
+                for (var j = 0; j <= q; j++)
+                {
+                    var coefficient = Binomial(p, i) * Binomial(q, j) * (i % 2 == 0 ? 1 : -1);
+                    var term = IntegrateSineCosinePowers(doubled, 0, i + j, 2 * rate, x);
+                    result += coefficient * term;
+                }
+            return result / Entity.Number.Integer.Create(EInteger.One.ShiftLeft(p + q));
+        }
+
+        private static Entity.Number.Integer Binomial(int n, int k)
+        {
+            var result = EInteger.One;
+            for (var i = 0; i < k; i++)
+                result = result * EInteger.FromInt32(n - i) / EInteger.FromInt32(i + 1);
+            return Entity.Number.Integer.Create(result);
+        }
+
+        /// <summary>
+        /// Whether <paramref name="expr"/> is an exponential in <paramref name="x"/>, and
+        /// at what rate: <c>B^(px + q)</c> grows as <c>e^(rate * x)</c> with
+        /// <c>rate = p * ln(B)</c>. The constant factor <c>B^q</c> needs no separating out,
+        /// because the antiderivative is written in terms of the original expression.
+        /// </summary>
+        private static bool IsExponentialRate(Entity expr, Entity.Variable x, out Entity rate)
+        {
+            rate = 0;
+            if (expr is not Entity.Powf(var @base, var exponent)
+                || @base.ContainsNode(x)
+                || !TreeAnalyzer.TryGetPolyLinear(exponent, x, out var perX, out _))
+                return false;
+            rate = perX * MathS.Ln(@base);
+            return true;
+        }
+
+        /// <summary>
+        /// The antiderivative of <c>sqrt(a x^2 + b x + c)</c>:
+        /// <c>(2ax + b) sqrt(Q) / (4a) + ((4ac - b^2) / (8a))</c> times the integral of
+        /// <c>1/sqrt(Q)</c> -- integration by parts once, leaving the reciprocal form that
+        /// <see cref="IntegrateOverRootOfQuadratic"/> already knows.
+        /// </summary>
+        /// <remarks>
+        /// Only where the leading coefficient is a number other than zero. With a = 0 this
+        /// is the square root of something linear, which the ordinary power rule already
+        /// integrates, and dividing by a would not be allowed anyway.
+        /// </remarks>
+        private static Entity IntegrateRootOfQuadratic(
+            Entity a, Entity b, Entity c, Entity radicand, Entity.Variable x)
+            => (2 * a * x + b) * MathS.Sqrt(radicand) / (4 * a)
+               + (4 * a * c - b * b) / (8 * a) * IntegrateOverRootOfQuadratic(1, a, b, c, radicand, x);
+
+        /// <summary>
+        /// The antiderivative of <c>k / sqrt(a x^2 + b x + c)</c>, which takes one of two
+        /// forms depending on the sign of the leading coefficient:
+        /// <list type="bullet">
+        /// <item>a &lt; 0, an arc of a circle: <c>-k/sqrt(-a) * arcsin((2ax + b) / sqrt(b^2 - 4ac))</c></item>
+        /// <item>a &gt; 0, a hyperbolic arc: <c>k/sqrt(a) * ln|2ax + b + 2 sqrt(a) sqrt(a x^2 + b x + c)|</c></item>
+        /// </list>
+        /// Returned as a piecewise on that sign, the way the rational quadratic below is,
+        /// since which one applies is not known until a and the coefficients are.
+        /// </summary>
+        private static Entity IntegrateOverRootOfQuadratic(
+            Entity numerator, Entity a, Entity b, Entity c, Entity radicand, Entity.Variable x)
+        {
+            var twoAxPlusB = 2 * a * x + b;
+
+            // a < 0: the radicand is a downward parabola, positive between its roots
+            var arcsinCase =
+                -numerator * MathS.Arcsin(twoAxPlusB / MathS.Sqrt(b * b - 4 * a * c)) / MathS.Sqrt(-a);
+
+            // a > 0
+            var logarithmCase =
+                numerator * MathS.Ln(MathS.Abs(twoAxPlusB + 2 * MathS.Sqrt(a) * MathS.Sqrt(radicand))) / MathS.Sqrt(a);
+
+            // a = 0: sqrt(bx + c), which integrates as an ordinary power
+            var linearCase = 2 * numerator * MathS.Sqrt(b * x + c) / b;
+
+            return MathS.Piecewise([
+                new Entity.Providedf(linearCase, a.EqualTo(0)),
+                new Entity.Providedf(arcsinCase, a < 0),
+                new Entity.Providedf(logarithmCase, a > 0)
+            ]);
+        }
+
         private static Entity IntegrateRationalQuadratic(Entity numerator, Entity a, Entity b, Entity c, Entity.Variable x)
         {
             // The formula depends on whether it's linear (a = 0) or quadratic (a ≠ 0)

Sources/Tests/UnitTests/Calculus/IntegrationTest.cs

@@ -268,20 +268,33 @@ public void TestLnAbsSquared()
             }
         }
 
-        [Theory(Skip = "TODO: integration by parts multiple times")]
-        [InlineData("ln(abs(x)) ^ 3", "C + x * (ln(abs(x)) ^ 3 - ln(abs(x)) ^ 2 - ln(abs(x)) ^ 2 - ln(abs(x)) ^ 2) + 6 * (x * (ln(abs(x)) - 1) + -x)")] // Triple integration by parts
-        [InlineData("e^x * sin(x)", "-1/2 * cos(x) * e ^ x + 1/2 * sin(x) * e ^ x + C")] // Classic integration by parts
-        [InlineData("e^x * cos(x)", "1/2 * cos(x) * e ^ x + 1/2 * sin(x) * e ^ x + C")] // Classic integration by parts
-        [InlineData("arctan(x)", "x * arctan(x) - 1/2 * ln(abs(x ^ 2 + 1)) + C")] // Integration by parts with 1 * arctan(x)
-        [InlineData("arcsin(x)", "x * arcsin(x) + sqrt(1 - x ^ 2) + C")] // Integration by parts with 1 * arcsin(x)
-        [InlineData("arccos(x)", "x * arccos(x) - sqrt(1 - x ^ 2) + C")] // Integration by parts with 1 * arccos(x)
+        // These five no longer need integration by parts at all. e^x*sin(x) and its
+        // cosine twin cycle under by parts, so they are solved as the closed form that
+        // cycle resolves to; the inverse trigonometric ones are by parts against 1, which
+        // has no product for the by-parts solver to split, so they are table entries.
+        [Theory]
+        [InlineData("e^x * sin(x)", "-1/2 * cos(x) * e ^ x + 1/2 * sin(x) * e ^ x + C")]
+        [InlineData("e^x * cos(x)", "1/2 * cos(x) * e ^ x + 1/2 * sin(x) * e ^ x + C")]
+        [InlineData("arctan(x)", "x * arctan(x) - 1/2 * ln(abs(x ^ 2 + 1)) + C")]
+        [InlineData("arcsin(x)", "x * arcsin(x) + sqrt(1 - x ^ 2) + C")]
+        [InlineData("arccos(x)", "x * arccos(x) - sqrt(1 - x ^ 2) + C")]
         public void TestIntegrationByPartsNonPolynomial(string initial, string expected)
         {
             var result = initial.Integrate("x").InnerSimplified;
             var expectedResult = expected.ToEntity().InnerSimplified;
             Assert.Equal(MathS.Boolean.True, result.EqualTo(expectedResult).Simplify());
         }
 
+        // Still open: this one wants by parts applied three times over.
+        [Theory(Skip = "TODO: integration by parts multiple times")]
+        [InlineData("ln(abs(x)) ^ 3", "C + x * (ln(abs(x)) ^ 3 - ln(abs(x)) ^ 2 - ln(abs(x)) ^ 2 - ln(abs(x)) ^ 2) + 6 * (x * (ln(abs(x)) - 1) + -x)")]
+        public void TestTripleIntegrationByParts(string initial, string expected)
+        {
+            var result = initial.Integrate("x").InnerSimplified;
+            var expectedResult = expected.ToEntity().InnerSimplified;
+            Assert.Equal(MathS.Boolean.True, result.EqualTo(expectedResult).Simplify());
+        }
+
         [Theory(Skip = "TODO: integration by parts multiple times")]
         [InlineData("sin(ln(abs(x)))", "x / 2 * (sin(ln(abs(x))) - cos(ln(abs(x)))) + C")] // Integration by parts twice
         [InlineData("cos(ln(abs(x)))", "x / 2 * (sin(ln(abs(x))) + cos(ln(abs(x)))) + C")] // Integration by parts twice