diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs index d97311210..0eb7f59da 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs @@ -23,6 +23,30 @@ internal static class IndefiniteIntegralSolver }); } + /// + /// A quotient of polynomials whose denominator has a rational root, split at that + /// root and integrated in two parts. + /// + /// + /// 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. + /// + 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) => diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/IntegralPatterns.cs b/Sources/AngouriMath/Functions/Continuous/Integration/IntegralPatterns.cs index 6211b028c..b08b071d8 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/IntegralPatterns.cs +++ b/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( ]); } + /// + /// ∫ (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. + /// + 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) diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs b/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs index 0a120f41d..e3cd35982 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs +++ b/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; diff --git a/Sources/AngouriMath/Functions/Simplification/PolynomialFactoring.cs b/Sources/AngouriMath/Functions/Simplification/PolynomialFactoring.cs index 583215126..159551cd7 100644 --- a/Sources/AngouriMath/Functions/Simplification/PolynomialFactoring.cs +++ b/Sources/AngouriMath/Functions/Simplification/PolynomialFactoring.cs @@ -105,21 +105,112 @@ internal static bool TryFactor(Entity expr, Variable x, [NotNullWhen(true)] out return degree > 1; } + /// + /// One step of a partial fraction decomposition, at a rational root of the + /// denominator: N/D becomes A/(x - r) + R/Q, where D = (x - r)Q. + /// + /// + /// + /// 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 + /// 1/(x^3 + 1) splits into (1/3)/(x + 1) and (2 - x)/(3(x^2 - x + 1)), + /// and the second of those is a quotient the integrator can read as it stands. + /// + /// + /// The coefficient is Heaviside's: at the root every other term of the + /// decomposition is finite, so A = N(r)/Q(r). What is left, N - A*Q, + /// then has r for a root by construction and divides exactly. + /// + /// + /// Rational roots only, and simple ones: a repeated root needs a term over + /// (x - r)^2 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. + /// + /// + 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; + } + + /// + /// 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. + /// + private static IEnumerable RootCandidates(ERational[] coefficients) + => coefficients[0].IsZero + ? new[] { ERational.Zero }.Concat(RationalRootCandidates(coefficients)) + : RationalRootCandidates(coefficients); + + /// x - r, or x where r is zero. + 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; + } + /// /// The coefficients of in , lowest /// power first, or if it is not a polynomial in /// with rational coefficients alone. /// 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(); var monomials = PolynomialSolver.GatherMonomialInformation( 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]; diff --git a/Sources/Tests/UnitTests/Calculus/PartialFractionsTest.cs b/Sources/Tests/UnitTests/Calculus/PartialFractionsTest.cs new file mode 100644 index 000000000..a98a49514 --- /dev/null +++ b/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 +{ + /// + /// 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. + /// + 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); + + /// + /// 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. + /// + [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()); + } +} diff --git a/Sources/Tests/UnitTests/Calculus/RationalIntegralsTest.cs b/Sources/Tests/UnitTests/Calculus/RationalIntegralsTest.cs new file mode 100644 index 000000000..d1baf7c6a --- /dev/null +++ b/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 +{ + /// + /// 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. + /// + 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); + } +}