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

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

The table matched k/(ax^2 + bx + c) only when k did not mention x, so anything with
a variable on top had no antiderivative at all:

  int x / (x^2 + 2x + 5)        int (x + 3) / (x^2 + 3x + 2)        int x / (x + 1)

A linear numerator is 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 that was already there, so this only needs
the rewrite. The same rewrite one degree down turns (px + q)/(bx + c) into the
constant p/b plus a remainder over the divisor.

Both are guarded on the leading coefficient being a non-zero number, since the
rewrite divides by it. A quadratic denominator with a vanishing quadratic term falls
through to the linear arm rather than being divided by zero.

Also 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 none of the four shapes was recognised.

Corpus 88/117 -> 91/117: int:table and int:rational-quadratic reach 100%, and
int:partial-fractions goes from 2 out of 5 to 3. 3901 unit tests and 127 F# tests
pass, no wrong answers in the corpus, and the property checker's 1320 checks are
clean.

x^2/(x^4 + 1) and 1/(x^3 + 1) are still out of reach: both need the denominator
factored first, which is a different piece of work.

* Split a rational function at a rational root of its denominator

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. It is on the
list of what is missing in
#233.

The decomposition is done a step at a time rather than in full, 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.
1/(x^3 + 1) splits at the root -1 into (1/3)/(x + 1) and
(2 - x)/(3(x^2 - x + 1)), and the second is a linear numerator over a quadratic,
which the previous commit reads as it stands. Each step takes a degree off the
denominator, so it ends.

The coefficient is Heaviside's -- at the root every other term of the
decomposition is finite, so A = N(r)/Q(r) -- and what is left, N - A*Q, then has
that root by construction and divides exactly. All of it on the exact ratios,
reusing the rational root search already here for factoring.

Rational roots only, and simple ones. A denominator with no rational root, such
as x^4 + 1, and one whose root repeats, which would need a term over the square
as well, are both left unevaluated rather than answered wrongly. A denominator
below degree three is left to the rules that answer it in one piece.

Ten integrals gained, each checked by differentiating back:
1/(x^3 + 1), x/(x^3 + 1), 1/(x^3 - x), 1/(x^3 + x), 1/(x^4 - 1), x/(x^4 - 1),
(x^2 + 1)/(x^3 - x), 1/((x - 1)(x - 2)(x - 3)) whether written factored or
multiplied out, and (x + 1)/(x^3 - x^2 - 2x).

21 tests, 10 of which fail without the change. 4038 in all, 127 F#, corpus
105/117 with the partial fraction category at 80% from 60% and nothing wrong,
in error or timing out. 1320 property checks.

Author: Rafael-SOWNet

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/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());
+    }
+}