Split a rational function at a rational root of its denominator (#233)

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

@@ -23,6 +23,30 @@ internal static class IndefiniteIntegralSolver
             });
         }
 
+        /// <summary>
+        /// A quotient of polynomials whose denominator has a rational root, split at that
+        /// root and integrated in two parts.
+        /// </summary>
+        /// <remarks>
+        /// The rules for a linear or a quadratic denominator answer those in one piece, so
+        /// what is left over are the denominators of degree three and up, which nothing
+        /// read at all: 1/(x^3 + 1) had no antiderivative. Splitting at the root -1 leaves
+        /// (1/3)/(x + 1) and (2 - x)/(3(x^2 - x + 1)), and both of those are already
+        /// integrable, the second by the rule for a linear numerator over a quadratic.
+        /// Each step takes a degree off the denominator, so this ends.
+        /// </remarks>
+        internal static Entity? SolveByPartialFractions(Entity expr, Entity.Variable x, bool integrateByParts)
+        {
+            if (expr is not Entity.Divf(var numerator, var denominator)
+                || !Functions.PolynomialFactoring.TrySplitOffRationalRoot(
+                    numerator, denominator, x, out var simple, out var restNumerator, out var restDenominator))
+                return null;
+            return Integration.ComputeIndefiniteIntegral(simple, x, integrateByParts) is { } first
+                && Integration.ComputeIndefiniteIntegral(restNumerator / restDenominator, x, integrateByParts) is { } rest
+                ? first + rest
+                : null;
+        }
+
         internal static Entity? SolveAsPolynomialTerm(Entity expr, Entity.Variable x, bool integrateByParts = true) => expr switch
         {
             Entity.Mulf(var m1, var m2) =>

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/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs

@@ -58,6 +58,7 @@ internal static partial class Integration
             if ((answer = IndefiniteIntegralSolver.SolveAsPolynomialTerm(expr, x, integrateByParts)) is { }) return answer;
             if ((answer = IndefiniteIntegralSolver.SolveLogarithmic(expr, x, integrateByParts)) is { }) return answer;
             if ((answer = IndefiniteIntegralSolver.SolveBySubstitution(expr, x, integrateByParts)) is { }) return answer;
+            if ((answer = IndefiniteIntegralSolver.SolveByPartialFractions(expr, x, integrateByParts)) is { }) return answer;
             if (integrateByParts && (answer = IndefiniteIntegralSolver.SolveIntegratingByParts(expr, x)) is { }) return answer;
             // this may expand to too many terms
             if ((answer = IndefiniteIntegralSolver.SolveBySplittingSum(expr, x, integrateByParts)) is { }) return answer;

Sources/AngouriMath/Functions/Simplification/PolynomialFactoring.cs

@@ -105,21 +105,112 @@ internal static bool TryFactor(Entity expr, Variable x, [NotNullWhen(true)] out
             return degree > 1;
         }
 
+        /// <summary>
+        /// One step of a partial fraction decomposition, at a rational root of the
+        /// denominator: <c>N/D</c> becomes <c>A/(x - r) + R/Q</c>, where <c>D = (x - r)Q</c>.
+        /// </summary>
+        /// <remarks>
+        /// <para>
+        /// One step, and not the whole decomposition, because a step is all that is needed:
+        /// what is left over is a smaller problem of the same kind, and by the time its
+        /// denominator is a quadratic there is already a rule for it. So
+        /// <c>1/(x^3 + 1)</c> splits into <c>(1/3)/(x + 1)</c> and <c>(2 - x)/(3(x^2 - x + 1))</c>,
+        /// and the second of those is a quotient the integrator can read as it stands.
+        /// </para>
+        /// <para>
+        /// The coefficient is Heaviside's: at the root every other term of the
+        /// decomposition is finite, so <c>A = N(r)/Q(r)</c>. What is left, <c>N - A*Q</c>,
+        /// then has <c>r</c> for a root by construction and divides exactly.
+        /// </para>
+        /// <para>
+        /// Rational roots only, and simple ones: a repeated root needs a term over
+        /// <c>(x - r)^2</c> as well, which this does not produce, so it declines instead.
+        /// Denominators of degree below three are left alone, since the rules for a linear
+        /// or quadratic denominator already answer those and answer them in one piece.
+        /// </para>
+        /// </remarks>
+        internal static bool TrySplitOffRationalRoot(
+            Entity numerator, Entity denominator, Variable x,
+            [NotNullWhen(true)] out Entity? simplePart,
+            [NotNullWhen(true)] out Entity? restNumerator,
+            [NotNullWhen(true)] out Entity? restDenominator)
+        {
+            simplePart = restNumerator = restDenominator = null;
+            if (!TryGetRationalCoefficients(denominator, x, leastTerms: 1, leastDegree: 3, out var d)
+                || !TryGetRationalCoefficients(numerator, x, leastTerms: 1, leastDegree: 0, out var n))
+                return false;
+            // Only a proper fraction decomposes; an improper one is a polynomial plus a
+            // proper fraction and has to be divided out first, which is not done here.
+            if (n.Length >= d.Length)
+                return false;
+
+            foreach (var root in RootCandidates(d))
+            {
+                if (!Evaluate(d, root).IsZero)
+                    continue;
+                var q = DivideByLinear(d, root);
+                var atRoot = Evaluate(q, root);
+                if (atRoot.IsZero)
+                    continue;                       // repeated root, see above
+                var a = Evaluate(n, root).Divide(atRoot);
+                var left = Subtract(n, Scale(q, a));
+                if (!Evaluate(left, root).IsZero)
+                    continue;                       // cannot happen; checked rather than assumed
+                simplePart = Rational.Create(a) / LinearFactor(root, x);
+                restNumerator = BuildPolynomial(DivideByLinear(left, root), x);
+                restDenominator = BuildPolynomial(q, x);
+                return true;
+            }
+            return false;
+        }
+
+        /// <summary>
+        /// Zero first, then the candidates of the rational root theorem, which cannot offer
+        /// it: they are built from divisors of the constant term, and a polynomial with a
+        /// root at zero has none.
+        /// </summary>
+        private static IEnumerable<ERational> RootCandidates(ERational[] coefficients)
+            => coefficients[0].IsZero
+                ? new[] { ERational.Zero }.Concat(RationalRootCandidates(coefficients))
+                : RationalRootCandidates(coefficients);
+
+        /// <summary>x - r, or x where r is zero.</summary>
+        private static Entity LinearFactor(ERational root, Variable x)
+            => root.IsZero ? x
+                : root.Sign < 0 ? x + Rational.Create(root.Negate())
+                : x - Rational.Create(root);
+
+        private static ERational[] Scale(ERational[] coefficients, ERational by)
+            => coefficients.Select(c => c.Multiply(by)).ToArray();
+
+        private static ERational[] Subtract(ERational[] left, ERational[] right)
+        {
+            var result = new ERational[System.Math.Max(left.Length, right.Length)];
+            for (var i = 0; i < result.Length; i++)
+                result[i] = (i < left.Length ? left[i] : ERational.Zero)
+                    .Subtract(i < right.Length ? right[i] : ERational.Zero);
+            return result;
+        }
+
         /// <summary>
         /// The coefficients of <paramref name="expr"/> in <paramref name="x"/>, lowest
         /// power first, or <see langword="false"/> if it is not a polynomial in
         /// <paramref name="x"/> with rational coefficients alone.
         /// </summary>
         private static bool TryGetRationalCoefficients(Entity expr, Variable x, out ERational[] coefficients)
+            => TryGetRationalCoefficients(expr, x, leastTerms: 2, leastDegree: 2, out coefficients);
+
+        private static bool TryGetRationalCoefficients(
+            Entity expr, Variable x, int leastTerms, int leastDegree, out ERational[] coefficients)
         {
             coefficients = System.Array.Empty<ERational>();
             var monomials = PolynomialSolver.GatherMonomialInformation<EInteger, TreeAnalyzer.PrimitiveInteger>(
                 Sumf.LinearChildren(expr.Expand()), x);
-            if (monomials is null || monomials.Count < 2)
+            if (monomials is null || monomials.Count < leastTerms)
                 return false;
 
             var degree = monomials.Keys.Max();
-            if (degree is null || degree > MaxDegree || degree < 2)
+            if (degree is null || degree > MaxDegree || degree < leastDegree)
                 return false;
 
             var found = new ERational[degree.ToInt32Checked() + 1];

Sources/Tests/UnitTests/Calculus/PartialFractionsTest.cs

@@ -0,0 +1,86 @@
+//
+// 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 whose denominator is of degree three or more. A linear or
+    /// quadratic denominator was answered in one piece, but nothing read anything above
+    /// that, so 1/(x^3 + 1) had no antiderivative at all. Named in the issue's own list of
+    /// what is missing: https://github.com/asc-community/AngouriMath/issues/233.
+    /// 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 PartialFractionsTest
+    {
+        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, 7);
+            }
+        }
+
+        // The denominator is split at one rational root at a time, and what is left over is
+        // a smaller problem of the same kind until its denominator is a quadratic.
+        [Theory]
+        [InlineData("1 / (x ^ 3 + 1)", new[] { 0.3, 1.7, 3.2, -2.4 })]
+        [InlineData("x / (x ^ 3 + 1)", new[] { 0.3, 1.7, 3.2, -2.4 })]
+        [InlineData("1 / (x ^ 3 - x)", new[] { 0.3, 1.7, 3.2, -2.4 })]
+        [InlineData("1 / (x ^ 3 + x)", new[] { 0.3, 1.7, 3.2, -2.4 })]
+        [InlineData("1 / (x ^ 4 - 1)", new[] { 0.3, 1.7, 3.2, -2.4 })]
+        [InlineData("x / (x ^ 4 - 1)", new[] { 0.3, 1.7, 3.2 })]
+        [InlineData("(x ^ 2 + 1) / (x ^ 3 - x)", new[] { 0.3, 1.7, 3.2, -2.4 })]
+        public void ADenominatorOfDegreeThreeOrMore(string integrand, double[] points) =>
+            AssertIsAntiderivative(integrand, points);
+
+        // Whether the denominator arrives factored or multiplied out makes no difference:
+        // it is read as a polynomial either way.
+        [Theory]
+        [InlineData("1 / ((x - 1) * (x - 2) * (x - 3))", new[] { 0.3, 1.7, 3.6, 5.2 })]
+        [InlineData("1 / (x ^ 3 - 6 * x ^ 2 + 11 * x - 6)", new[] { 0.3, 1.7, 3.6 })]
+        [InlineData("(x + 1) / (x ^ 3 - x ^ 2 - 2 * x)", new[] { 0.3, 1.7, 3.2, -2.4 })]
+        public void FactoredOrMultipliedOutAlike(string integrand, double[] points) =>
+            AssertIsAntiderivative(integrand, points);
+
+        // The shapes this sits next to have to keep working, in particular the linear and
+        // quadratic denominators that are answered in one piece and must not be split.
+        [Theory]
+        [InlineData("1 / (x ^ 2 + 1)", new[] { 0.3, 1.7 })]
+        [InlineData("1 / (x ^ 2 - 1)", new[] { 0.3, 3.2 })]
+        [InlineData("1 / (x + 1)", new[] { 0.3, 1.7 })]
+        [InlineData("x / (x ^ 2 + 1)", new[] { 0.3, 1.7 })]
+        [InlineData("x / (x ^ 2 + 2 * x + 5)", 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("x * e ^ x", new[] { 0.3, 1.7 })]
+        public void NeighbouringFormsAreUnaffected(string integrand, double[] points) =>
+            AssertIsAntiderivative(integrand, points);
+
+        /// <summary>
+        /// A denominator with no rational root is not split, and neither is one whose root
+        /// repeats -- that needs a term over the square as well, which is not produced here.
+        /// Both are left unevaluated rather than answered wrongly.
+        /// </summary>
+        [Theory]
+        [InlineData("x ^ 2 / (x ^ 4 + 1)")]
+        [InlineData("1 / (x ^ 3 + x ^ 2 + x + 2)")]
+        public void WhatCannotBeSplitIsLeftAlone(string integrand) =>
+            Assert.Contains("integral(", integrand.ToEntity().Integrate("x").Stringize());
+    }
+}