Integrate a quotient whose numerator is not constant (#233)

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

@@ -66,7 +66,7 @@ _ when TryReadSineCosinePowers(expr, x, out var trigArg, out var sinePower, out
                     (arg / a) * (MathS.Ln(MathS.Abs(arg)) - 1),
 
             Entity.Divf(var numerator, var denominator) when
-                !numerator.ContainsNode(x) 
+                !numerator.ContainsNode(x)
                 && 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),
 
@@ -132,6 +132,42 @@ _ when TryReadSineCosinePowers(expr, x, out var trigArg, out var sinePower, out
                 && TreeAnalyzer.TryGetPolyQuadratic(radicand, x, out var ra, out var rb, out var rc)
                     => IntegrateOverRootOfQuadratic(numerator, ra, rb, rc, radicand, x),
 
+            // ∫ (px + q)/(ax^2 + bx + c) dx. Only the constant numerator was covered, so
+            // x/(x^2 + 2x + 5) and (x + 3)/(x^2 + 3x + 2) had no antiderivative at all.
+            Entity.Divf(var numerator, var denominator) when
+                numerator.ContainsNode(x)
+                && TreeAnalyzer.TryGetPolyLinear(numerator, x, out var p, out var q)
+                && TreeAnalyzer.TryGetPolyQuadratic(denominator, x, out var a, out var b, out var c)
+                && a.Evaled is Entity.Number.Complex { IsZero: false }
+                    => IntegrateLinearOverQuadratic(p, q, a, b, c, denominator, x),
+
+            // ∫ (px + q)/(bx + c) dx, the same rewrite one degree down: the quotient is the
+            // constant p/b plus a remainder over the divisor. x/(x + 1) had no antiderivative.
+            Entity.Divf(var numerator, var denominator) when
+                numerator.ContainsNode(x)
+                && TreeAnalyzer.TryGetPolyLinear(numerator, x, out var p, out var q)
+                && TreeAnalyzer.TryGetPolyLinear(denominator, x, out var b, out var c)
+                && b.Evaled is Entity.Number.Complex { IsZero: false }
+                    => p * x / b + (q - p * c / b) * MathS.Ln(MathS.Abs(denominator)) / b,
+
+            // 1/cos(u)^2 and 1/sin(u)^2, which are written that way at least as often as
+            // sec(u)^2 and csc(u)^2 and were not recognised in that form.
+            Entity.Divf(var numerator, Entity.Powf(Entity.Cosf(var arg), Entity.Number.Integer(2))) when
+                !numerator.ContainsNode(x) && TreeAnalyzer.TryGetPolyLinear(arg, x, out var a, out _) =>
+                    numerator * MathS.Tan(arg) / a,
+
+            Entity.Divf(var numerator, Entity.Powf(Entity.Sinf(var arg), Entity.Number.Integer(2))) when
+                !numerator.ContainsNode(x) && TreeAnalyzer.TryGetPolyLinear(arg, x, out var a, out _) =>
+                    -numerator * MathS.Cotan(arg) / a,
+
+            Entity.Powf(Entity.Secantf(var arg), Entity.Number.Integer(2)) when
+                TreeAnalyzer.TryGetPolyLinear(arg, x, out var a, out _) =>
+                    MathS.Tan(arg) / a,
+
+            Entity.Powf(Entity.Cosecantf(var arg), Entity.Number.Integer(2)) when
+                TreeAnalyzer.TryGetPolyLinear(arg, x, out var a, out _) =>
+                    -MathS.Cotan(arg) / a,
+
             _ => null
         };
 
@@ -311,6 +347,17 @@ private static Entity IntegrateOverRootOfQuadratic(
             ]);
         }
 
+        /// <summary>
+        /// ∫ (px + q)/(ax^2 + bx + c) dx, by writing the numerator as a multiple of the
+        /// denominator's derivative plus a constant:
+        /// px + q = (p/2a)(2ax + b) + (q - pb/2a). The first part integrates to a logarithm
+        /// and the second is the constant-numerator case above.
+        /// </summary>
+        private static Entity IntegrateLinearOverQuadratic(
+            Entity p, Entity q, Entity a, Entity b, Entity c, Entity denominator, Entity.Variable x)
+            => p / (2 * a) * MathS.Ln(MathS.Abs(denominator))
+               + IntegrateRationalQuadratic(q - p * b / (2 * a), a, b, c, x);
+
         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/RationalIntegralsTest.cs

@@ -0,0 +1,81 @@
+//
+// Copyright (c) 2019-2022 Angouri.
+// AngouriMath is licensed under MIT.
+// Details: https://github.com/asc-community/AngouriMath/blob/master/LICENSE.md.
+// Website: https://am.angouri.org.
+//
+
+using AngouriMath;
+using AngouriMath.Extensions;
+using Xunit;
+
+namespace AngouriMath.Tests.Calculus
+{
+    /// <summary>
+    /// Quotients of polynomials of low degree. Only a constant numerator was recognised, so
+    /// x/(x^2 + 2x + 5) and x/(x + 1) had no antiderivative at all. Each answer is checked by
+    /// differentiating it back and comparing at points, since what matters is that it is an
+    /// antiderivative and not what form it is written in.
+    /// </summary>
+    public sealed class RationalIntegralsTest
+    {
+        private static void AssertIsAntiderivative(string integrand, params double[] points)
+        {
+            var f = integrand.ToEntity();
+            var antiderivative = f.Integrate("x");
+            Assert.DoesNotContain("integral(", antiderivative.Stringize());
+            var derivative = antiderivative.Substitute("C", 0).Differentiate("x");
+            foreach (var point in points)
+            {
+                var expected = f.Substitute("x", point).EvalNumerical().RealPart.EDecimal.ToDouble();
+                var actual = derivative.Substitute("x", point).EvalNumerical().RealPart.EDecimal.ToDouble();
+                Assert.Equal(expected, actual, 8);
+            }
+        }
+
+        // (px + q) / (ax^2 + bx + c). The numerator is written as a multiple of the
+        // denominator's derivative plus a constant, which splits it into a logarithm and the
+        // constant-numerator case that was already there.
+        [Theory]
+        [InlineData("x / (x ^ 2 + 2 * x + 5)", new[] { 0.3, 1.7, -0.6 })]
+        [InlineData("(x + 3) / (x ^ 2 + 3 * x + 2)", new[] { 0.3, 1.7, 4.2 })]
+        [InlineData("(2 * x + 1) / (x ^ 2 + x + 1)", new[] { 0.3, 1.7, -0.6 })]
+        [InlineData("x / (x ^ 2 + 1)", new[] { 0.3, 1.7, -0.6 })]
+        [InlineData("(3 * x - 2) / (2 * x ^ 2 + 5)", new[] { 0.3, 1.7, -0.6 })]
+        public void ALinearNumeratorOverAQuadratic(string integrand, double[] points) =>
+            AssertIsAntiderivative(integrand, points);
+
+        // (px + q) / (bx + c), the same rewrite one degree down.
+        [Theory]
+        [InlineData("x / (x + 1)", new[] { 0.3, 1.7 })]
+        [InlineData("(2 * x + 1) / (x - 3)", new[] { 0.3, 1.7 })]
+        [InlineData("(3 * x) / (2 * x + 5)", new[] { 0.3, 1.7 })]
+        public void ALinearNumeratorOverALinearDenominator(string integrand, double[] points) =>
+            AssertIsAntiderivative(integrand, points);
+
+        // 1/cos(u)^2 and 1/sin(u)^2 are written that way at least as often as sec(u)^2 and
+        // csc(u)^2, and neither of the four shapes was recognised.
+        [Theory]
+        [InlineData("1 / cos(x) ^ 2", new[] { 0.3, 1.1, -0.6 })]
+        [InlineData("1 / sin(x) ^ 2", new[] { 0.3, 1.1, 2.2 })]
+        [InlineData("2 / cos(3 * x) ^ 2", new[] { 0.3, 0.9 })]
+        [InlineData("sec(x) ^ 2", new[] { 0.3, 1.1, -0.6 })]
+        [InlineData("cosec(x) ^ 2", new[] { 0.3, 1.1, 2.2 })]
+        public void ReciprocalSquaresOfTheWaves(string integrand, double[] points) =>
+            AssertIsAntiderivative(integrand, points);
+
+        // The shapes these sit next to in the table have to keep working. A constant over a
+        // quadratic in particular is matched by the earlier arm and must stay there.
+        [Theory]
+        [InlineData("1 / (x ^ 2 + 1)", new[] { 0.3, 1.7, -0.6 })]
+        [InlineData("1 / (x ^ 2 + 2 * x + 5)", new[] { 0.3, 1.7, -0.6 })]
+        [InlineData("1 / (2 * x + 3)", new[] { 0.3, 1.7 })]
+        [InlineData("1 / x", new[] { 0.3, 1.7 })]
+        [InlineData("x ^ 2", new[] { 0.3, 1.7 })]
+        [InlineData("sin(x)", new[] { 0.3, 1.7 })]
+        [InlineData("tan(x)", new[] { 0.3, 1.1 })]
+        [InlineData("x * e ^ x", new[] { 0.3, 1.7 })]
+        public void NeighbouringFormsAreUnaffected(string integrand, double[] points) =>
+            AssertIsAntiderivative(integrand, points);
+    }
+}