From 49a77c01cfd45142bfa0d5ad084d9f145912b0d1 Mon Sep 17 00:00:00 2001 From: Rafael Vuijk Date: Thu, 3 Sep 2026 20:25:31 +0000 Subject: [PATCH] A repeated quadratic denominator is integrated `1/(x^2 + 1)^2` had no antiderivative, while `1/(x^2 - 1)^2` and `1/(x^2 + 2x + 1)^2` both did: a denominator with real roots comes apart into linear factors and never reached the gap. What was missing was the irreducible case. Writing Q for the quadratic, u for its derivative 2ax + b and D for 4ac - b^2, the identity u^2 = 4aQ - D turns the derivative of u/Q^(m-1) into an equation for the integral of 1/Q^m: J_m = (u / Q^(m-1) + 2a(2m - 3) J_(m-1)) / ((m - 1) D) unrolled down to J_1, which the table already answered. It is an identity in a, b and c rather than a fact about signs, so one line serves both signs of the discriminant; D = 0 is the sole exclusion, and the division by D is why. The two shapes the reduction cannot speak for get their own arms -- a zero leading coefficient leaves a linear power, a zero discriminant leaves a power of u -- which is what lets a symbolic quadratic be answered as a piecewise instead of carrying a division by zero. A numerator of higher degree is divided by the quadratic first: N = qQ + r takes one power off the denominator and two degrees off the numerator, and ends. An improper fraction is declined rather than divided out here, since the recursion bottoms out at a single power with a numerator of degree two or more -- exactly the case SolveByPartialFractions already opens with a long division. The remainder is computed as N - qQ rather than read from the division, which reports it already divided by the divisor: x^2 over x^2 + 2 answers (1, -2/(x^2 + 2)), and taking that as a polynomial gives a numerator that is not one, so every polynomial-numerator case declined. A numerator already along the denominator's derivative leaves no second integral and the term is dropped rather than multiplied by zero. That is not tidiness: 0 * (2x/(x^2 + 1)) is `0 provided not 1 + x^2 = 0` and not 0, so multiplying would attach a condition to x/(x^2 + 1)^2, which had none and is owed none. A denominator with real roots now takes this route rather than partial fractions. The value is unchanged and the form is shorter -- one logarithm where there were two -- which is why the metric that ranks candidates selects it. Measured, both arms on one machine, the same way. Where the integrand is rational this is 26 to 165 times faster, because TryStandardIntegrals runs before every search and a shape answered outright never enters the candidate exploration: 1/(x^2 + 1)^2 goes from 2879 ms unevaluated to 111 ms answered, and x^2/(x^2 + 2)^2 from 92,990 ms unevaluated to 565 ms answered. Where the integrand is a transcendental function over a repeated quadratic, none of which has an elementary antiderivative, the same search takes 3.4 to 10.8 times longer to conclude nothing -- e^x/(x^2 + 1)^2 from 7944 ms to 85,683 ms. That is growth in integration by parts, which those integrands reach and this rule never matches; it is filed separately rather than bounded here, since by-parts is machinery every integral goes through. Full suite 9515 passed, 0 failed. 26 new tests, each answer differentiated back and compared with the integrand. Part of #180. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_012sonx8iAspMiwRwokT1Ura --- BREAKING-CHANGES.md | 55 ++++++ Sources/.editorconfig | 3 + .../Integration/IntegralPatterns.cs | 177 ++++++++++++++++++ .../Calculus/RepeatedQuadraticIntegralTest.cs | 155 +++++++++++++++ 4 files changed, 390 insertions(+) create mode 100644 Sources/Tests/UnitTests/Calculus/RepeatedQuadraticIntegralTest.cs diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index 415a76396..e00978a73 100644 --- a/BREAKING-CHANGES.md +++ b/BREAKING-CHANGES.md @@ -111,6 +111,61 @@ read first. | | `"sum(k, k, 1, n)".ToEntity().Simplify()`, and every summation whose body is a polynomial in the index | `sum(k, k, 1, n)` — carried | `piecewise((n + n ^ 2) / 2 provided n >= 0, 0)` | | | `"sum(k, k, 1, 100000)".ToEntity().Simplify()`, and every concrete range past a hundred terms | `sum(k, k, 1, 100000)` — carried | `5000050000` | | | `"product(k, k, 1, n)".ToEntity().Simplify()`, and every product whose body is a monomial in the index | `product(k, k, 1, n)` — carried | `piecewise(n! provided n >= 1, 1)` | +| | `"1/(x^2 + 1)^2".Integrate("x")`, and every proper rational function over a repeated irreducible quadratic | `integral(1 / (1 + x ^ 2) ^ 2, x)` — left unevaluated, after 2.9 s | `arctan(x) / 2 + C + 1/2 * x / (x ^ 2 + 1)`, in 0.11 s | +| | `"x^2/(x^2 + 2)^2".Integrate("x")`, and every polynomial numerator over one | `integral(x ^ 2 / (x ^ 2 + 2) ^ 2, x)` — left unevaluated, after 93 s | the antiderivative, in 0.57 s | +| | `"1/(x^2 - 1)^2".Integrate("x")`, and every repeated quadratic whose roots are real | `C - ln(x + -1) / 4 + ln(1 + 2 * x + x ^ 2) / 8 + -1/2 * x / (x ^ 2 + -1)` | `C - ln(1 + (-2) / (x + 1)) / 4 + -1/2 * x / (x ^ 2 + -1)` — the same value, one logarithm rather than two | + +### A repeated quadratic denominator is integrated + +`1/(x^2 + 1)^2` had no antiderivative, while `1/(x^2 - 1)^2` and `1/(x^2 + 2x + 1)^2` both did — a +denominator with real roots comes apart into linear factors and never reached the gap. What was +missing was the irreducible case +([#180](https://github.com/asc-community/AngouriMath/issues/180)). + +``` +"1/(x^2 + 1)^2".Integrate("x") was integral(1 / (1 + x ^ 2) ^ 2, x) + is arctan(x) / 2 + C + 1/2 * x / (x ^ 2 + 1) +"x^2/(x^2 + 2)^2".Integrate("x") was integral(x ^ 2 / (x ^ 2 + 2) ^ 2, x) + is (sqrt(2) * arctan(sqrt(2) * x / 2) * 2 + (-4) * x / (x ^ 2 + 2)) / 8 + + C provided not x ^ 2 + 2 = 0 +``` + +The condition on that second one is the library's own and it is true: the antiderivative has no +value where the denominator vanishes. It is attached to some of these answers and not others, and I +did not pin down what decides which — recorded as observed rather than explained. + +Writing `Q` for the quadratic, `u` for its derivative `2ax + b` and `D` for `4ac - b^2`, the +identity `u^2 = 4aQ - D` turns the derivative of `u/Q^(m-1)` into an equation for the integral of +`1/Q^m`, which is unrolled down to the single power the table already answered. It is an identity in +`a`, `b` and `c` rather than a fact about signs, so one line serves both signs of the discriminant; +`D = 0` is the sole exclusion, and the division by `D` is why. A numerator of higher degree is +divided by `Q` first. + +**A denominator with real roots now takes this route rather than partial fractions.** The value is +unchanged and the form is shorter — one logarithm where there were two — which is why the +node-count metric that ranks candidates selects it. No test pinned the old spelling. + +**The cost, measured rather than estimated.** Where the integrand is rational this is between 26 and +165 times faster, because `TryStandardIntegrals` runs before every search and a shape answered +outright never enters the candidate exploration. Where the integrand is a transcendental function +over a repeated quadratic — `sin(x)/(x^2 + 1)^2`, `e^x/(x^2 + 1)^2`, `ln(x)/(x^2 + 1)^2`, none of +which has an elementary antiderivative — the same search now takes **3.4 to 10.8 times longer to +conclude nothing**, because the sub-integral it used to fail on immediately now succeeds and +integration by parts carries on from a larger expression. Both arms measured on one machine, master +and branch, the same way: + +| | master | now | +|---|---|---| +| `1/(x^2 + 1)^2` | 2879 ms, unevaluated | **111 ms, answered** | +| `x^2/(x^2 + 2)^2` | 92 990 ms, unevaluated | **565 ms, answered** | +| `1/(x^4 + 1)^2` | 3811 ms | 3870 ms — the base is not quadratic, and nothing changes | +| `ln(x)/(x^2 + 1)^2` | 15 034 ms, unevaluated | 51 072 ms, unevaluated | +| `sin(x)/(x^2 + 1)^2` | 26 492 ms, unevaluated | 114 851 ms, unevaluated | +| `e^x/(x^2 + 1)^2` | 7944 ms, unevaluated | 85 683 ms, unevaluated | + +That last group is a growth in integration by parts rather than in this rule, which those integrands +never match; it is filed separately rather than fixed here, since bounding by-parts is a change to +machinery every integral goes through. ### A nested radical comes apart diff --git a/Sources/.editorconfig b/Sources/.editorconfig index febbf95b7..35cd5c5f6 100644 --- a/Sources/.editorconfig +++ b/Sources/.editorconfig @@ -189,3 +189,6 @@ file_header_template=\nCopyright (c) 2019-2026 Angouri.\nAngouriMath is licensed [Tests/UnitTests/PatternsTest/DenestRadicalTest.cs] file_header_template=\nCopyright (c) 2019-2026 Angouri.\nAngouriMath is licensed under MIT.\nDetails: https://github.com/asc-community/AngouriMath/blob/master/LICENSE.md.\nWebsite: https://am.angouri.org.\n + +[Tests/UnitTests/Calculus/RepeatedQuadraticIntegralTest.cs] +file_header_template=\nCopyright (c) 2019-2026 Angouri.\nAngouriMath is licensed under MIT.\nDetails: https://github.com/asc-community/AngouriMath/blob/master/LICENSE.md.\nWebsite: https://am.angouri.org.\n diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/IntegralPatterns.cs b/Sources/AngouriMath/Functions/Continuous/Integration/IntegralPatterns.cs index a5fd1dff5..a3ad55db2 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/IntegralPatterns.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/IntegralPatterns.cs @@ -223,6 +223,18 @@ _ when TryReadOverSquareTimesRoot(expr, x, out var overFactor, out var overRadic && a.Evaled is Entity.Number.Complex { IsZero: false } => IntegrateLinearOverQuadratic(p, q, a, b, c, denominator, x), + // ∫ N(x)/(ax^2 + bx + c)^n dx for a whole n of two or more and a proper fraction. + // After the single-power arms above, which already read a denominator like + // (x + 1)^2 -- expanded, it is a quadratic, so the perfect-square arm there answers + // it and nothing about that changes. + Entity.Divf(var numerator, Entity.Powf(var repeated, Entity.Number.Integer repetitions)) when + repetitions.EInteger.CanFitInInt32() + && repetitions.EInteger.ToInt32Checked() is var n and >= 2 + && TreeAnalyzer.TryGetPolyQuadratic(repeated, x, out var qa, out var qb, out var qc) + && (!numerator.ContainsNode(x) || qa.Evaled is Entity.Number.Complex { IsZero: false }) + && IntegrateRationalOverPowerOfQuadratic(numerator, qa, qb, qc, repeated, n, x) is { } repeatedAnswer + => repeatedAnswer, + // ∫ (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 @@ -529,5 +541,170 @@ private static Entity IntegrateRationalQuadratic(Entity numerator, Entity a, Ent new Entity.Providedf(lnCase, discriminant < 0) ]).InnerSimplified; } + + /// + /// ∫ k / (a x^2 + b x + c)^n dx for a whole n of two or more, by the reduction that + /// takes one power off the denominator at a time. + /// + /// + /// + /// Write Q for the quadratic, u for its derivative 2ax + b and + /// D for 4ac - b^2. Then u^2 is 4aQ - D, so differentiating + /// u / Q^(m-1) gives 2a(3 - 2m)/Q^(m-1) + (m-1)D/Q^m, and reading that as + /// an equation for the integral of 1/Q^m leaves + /// + /// J_m = (u / Q^(m-1) + 2a(2m - 3) J_(m-1)) / ((m - 1) D) + /// + /// which is unrolled down to J_1, the case the table above already answers. It is + /// an identity in a, b and c rather than a fact about signs, so one line serves a + /// positive and a negative discriminant alike; D = 0 is the only exclusion, and + /// the division by D is why. + /// + /// + /// The two shapes the reduction cannot speak for get their own arms, and each is + /// elementary: with a = 0 the denominator is a linear power, and with + /// D = 0 the quadratic is u^2/(4a) so the whole integrand is a power of + /// u. Both are already answered for a numeric quadratic — by the power rule and + /// by partial fractions — and the arms are here for a symbolic one, where the + /// discriminant's sign is not known and dropping an arm is not available. + /// + /// + /// What was missing was the irreducible case. 1/(x^2 - 1)^2 had an antiderivative + /// and 1/(x^2 + 2x + 1)^2 had one, because a denominator with real roots comes + /// apart into linear factors and never reaches here; 1/(x^2 + 1)^2 had none. + /// https://github.com/asc-community/AngouriMath/issues/180 + /// + /// + private static Entity IntegrateOverPowerOfQuadratic( + Entity numerator, Entity a, Entity b, Entity c, int power, Entity.Variable x) + { + var quadratic = a * x * x + b * x + c; + var derivative = 2 * a * x + b; + var discriminant = 4 * a * c - b * b; + + // a = 0: k/(bx + c)^n, an ordinary power -- and the constant k/c^n when b is zero + // too, where writing the power would divide by it. Same reason as the rule above. + var linearCase = TreeAnalyzer.IsZero(b) + ? numerator * x / MathS.Pow(c, power) + : numerator * MathS.Pow(b * x + c, 1 - power) / (b * (1 - power)); + + // D = 0: the quadratic is u^2/(4a), so the integrand is k(4a)^n / u^(2n). + var perfectSquareCase = + numerator * MathS.Pow(4 * a, power) * MathS.Pow(derivative, 1 - 2 * power) + / (2 * a * (1 - 2 * power)); + + var sqrtDiscriminant = MathS.Sqrt(discriminant); + var sqrtNegDiscriminant = MathS.Sqrt(-discriminant); + + return MathS.Piecewise([ + new Entity.Providedf(linearCase, a.EqualTo(0)), + new Entity.Providedf(perfectSquareCase, discriminant.EqualTo(0)), + new Entity.Providedf( + Reduce(2 * MathS.Arctan(derivative / sqrtDiscriminant) / sqrtDiscriminant), + discriminant > 0), + new Entity.Providedf( + Reduce(AntiderivativeLog((derivative - sqrtNegDiscriminant) / (derivative + sqrtNegDiscriminant)) + / sqrtNegDiscriminant), + discriminant < 0) + ]).InnerSimplified; + + Entity Reduce(Entity firstPower) + { + var accumulated = firstPower; + for (var m = 2; m <= power; m++) + accumulated = + (derivative / MathS.Pow(quadratic, m - 1) + 2 * a * (2 * m - 3) * accumulated) + / ((m - 1) * discriminant); + return numerator * accumulated; + } + } + + /// + /// ∫ (px + q)/(a x^2 + b x + c)^n dx, by the same rewrite of the numerator the single + /// power uses: px + q = (p/2a)(2ax + b) + (q - pb/2a). The first part is the + /// denominator's own derivative over a power of it, which integrates as a power; the + /// second is the constant-numerator case above. + /// + /// + /// A numerator that is already a multiple of that derivative leaves nothing of the + /// second part, and it is dropped rather than multiplied by zero. Multiplying is not + /// harmless: the reduction is a quotient by the quadratic, so 0 * (2x/(x^2 + 1)) + /// is 0 provided not 1 + x^2 = 0 and not 0 — correctly, since the factor + /// has no value there. Carrying that condition would attach it to the answer for + /// x/(x^2 + 1)^2, which had none before and is owed none. + /// + private static Entity IntegrateLinearOverPowerOfQuadratic( + Entity p, Entity q, Entity a, Entity b, Entity c, Entity quadratic, int power, Entity.Variable x) + { + var alongTheDerivative = p / (2 * a) * MathS.Pow(quadratic, 1 - power) / (1 - power); + var constantPart = q - p * b / (2 * a); + return TreeAnalyzer.IsZero(constantPart) + ? alongTheDerivative + : alongTheDerivative + IntegrateOverPowerOfQuadratic(constantPart, a, b, c, power, x); + } + + /// + /// ∫ N(x)/(a x^2 + b x + c)^n dx for a numerator of any degree, or + /// where the fraction is improper and this is not the rule for it. + /// + /// + /// + /// Dividing the numerator by the quadratic writes N = qQ + r with r linear + /// at worst, so N/Q^n is q/Q^(n-1) + r/Q^n — one power off the denominator + /// and two degrees off the numerator at each step, which ends. The remainder is answered + /// by the two helpers above and the quotient goes round again. + /// + /// + /// Improper fractions are declined rather than divided out here. The recursion + /// bottoms out at a single power, and a numerator still of degree two or more there is + /// exactly the improper case, which + /// already opens with + /// . Declining hands it there instead of + /// answering it twice, in two shapes, from two places. + /// + /// + /// This is not only reach. TryStandardIntegrals runs before every search, so a + /// shape answered here never reaches the candidate exploration that re-integrates + /// rewritten forms — which is where this family's cost has always been. Both + /// x^2/(x^2 + 2)^2 and (x^2 + 1)/(x^2 + 2)^2 were tens of seconds spent to + /// return the integral unevaluated. + /// + /// + private static Entity? IntegrateRationalOverPowerOfQuadratic( + Entity numerator, Entity a, Entity b, Entity c, Entity quadratic, int power, Entity.Variable x) + { + if (!numerator.ContainsNode(x)) + return power == 1 + ? IntegrateRationalQuadratic(numerator, a, b, c, x) + : IntegrateOverPowerOfQuadratic(numerator, a, b, c, power, x); + + if (TreeAnalyzer.TryGetPolyLinear(numerator, x, out var p, out var q)) + return power == 1 + ? IntegrateLinearOverQuadratic(p, q, a, b, c, quadratic, x) + : IntegrateLinearOverPowerOfQuadratic(p, q, a, b, c, quadratic, power, x); + + if (power < 2) + return null; + + var division = TreeAnalyzer.PolynomialLongDivision(numerator, quadratic); + if (division is null) + return null; + + // The quotient is what is wanted; the remainder that comes back with it is already + // divided by the divisor -- x^2 over x^2 + 2 answers (1, -2/(x^2 + 2)) -- so taking + // it as a polynomial and dividing it again is a numerator that is not one, and the + // next division declines. Subtracting gives the polynomial remainder whatever shape + // the division chose to report. + // Reduced before either is used: the division reports its quotient built up term by + // term, so x^2 over x^2 + 2 comes back as `0 + 1` rather than `1`, and the arms below + // ask whether a numerator is constant. + var quotient = division.Value.Divided.InnerSimplified; + var remainder = (numerator - quotient * quadratic).InnerSimplified; + + return IntegrateRationalOverPowerOfQuadratic(remainder, a, b, c, quadratic, power, x) is { } head + && IntegrateRationalOverPowerOfQuadratic(quotient, a, b, c, quadratic, power - 1, x) is { } rest + ? head + rest + : null; + } } } diff --git a/Sources/Tests/UnitTests/Calculus/RepeatedQuadraticIntegralTest.cs b/Sources/Tests/UnitTests/Calculus/RepeatedQuadraticIntegralTest.cs new file mode 100644 index 000000000..31c443c45 --- /dev/null +++ b/Sources/Tests/UnitTests/Calculus/RepeatedQuadraticIntegralTest.cs @@ -0,0 +1,155 @@ +// +// Copyright (c) 2019-2026 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 +{ + /// + /// A quadratic denominator raised to a power. 1/(x^2 + 1)^2 had no antiderivative, + /// while 1/(x^2 - 1)^2 and 1/(x^2 + 2x + 1)^2 both did — a denominator with + /// real roots comes apart into linear factors and never reached the gap. + /// https://github.com/asc-community/AngouriMath/issues/180 + /// + [Trait("Area", "Calculus")] + public sealed class RepeatedQuadraticIntegralTest + { + /// + /// Differentiates the answer back and compares it with the integrand at points. The + /// guard on integral( is the whole point of the helper: differentiating an + /// unevaluated integral(f, x) gives back f, so without it every case + /// the rule declines would pass. + /// + private static Entity 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); + } + return antiderivative; + } + + // k/(ax^2 + bx + c)^n with no real roots -- the shape that had no antiderivative at all. + [Theory] + [InlineData("1 / (x ^ 2 + 1) ^ 2", new[] { 0.37, 1.4, -0.83 })] + [InlineData("1 / (x ^ 2 + 4) ^ 2", new[] { 0.37, 1.4, -0.83 })] + [InlineData("1 / (3 * x ^ 2 + 5) ^ 2", new[] { 0.37, 1.4, -0.83 })] + [InlineData("1 / (2 * x ^ 2 + 3 * x + 7) ^ 2", new[] { 0.37, 1.4, -0.83 })] + [InlineData("5 / (x ^ 2 + 2 * x + 5) ^ 2", new[] { 0.37, 1.4, -0.83 })] + public void AConstantOverARepeatedIrreducibleQuadratic(string integrand, double[] points) => + AssertIsAntiderivative(integrand, points); + + // The reduction takes one power off at a time, so a third and a fourth power are the + // same line applied again rather than new cases. + [Theory] + [InlineData("1 / (x ^ 2 + 1) ^ 3", new[] { 0.37, 1.4, -0.83 })] + [InlineData("1 / (x ^ 2 + 1) ^ 4", new[] { 0.37, 1.4, -0.83 })] + [InlineData("1 / (x ^ 2 + 2) ^ 5", new[] { 0.37, 1.4, -0.83 })] + public void AHigherPowerIsTheSameReductionAppliedAgain(string integrand, double[] points) => + AssertIsAntiderivative(integrand, points); + + // (px + q)/(ax^2 + bx + c)^n, by writing the numerator as a multiple of the + // denominator's derivative plus a constant. + [Theory] + [InlineData("(x + 3) / (x ^ 2 + 2 * x + 5) ^ 2", new[] { 0.37, 1.4, -0.83 })] + [InlineData("x / (x ^ 2 + 1) ^ 2", new[] { 0.37, 1.4, -0.83 })] + [InlineData("x / (x ^ 2 + 1) ^ 3", new[] { 0.37, 1.4, -0.83 })] + [InlineData("(2 * x + 1) / (x ^ 2 + x + 1) ^ 2", new[] { 0.37, 1.4, -0.83 })] + public void ALinearNumeratorOverARepeatedQuadratic(string integrand, double[] points) => + AssertIsAntiderivative(integrand, points); + + // A numerator of higher degree is divided by the quadratic: N = qQ + r takes one power + // off the denominator and two degrees off the numerator, and ends. + [Theory] + [InlineData("x ^ 2 / (x ^ 2 + 2) ^ 2", new[] { 0.37, 1.4, -0.83 })] + [InlineData("(x ^ 2 + 1) / (x ^ 2 + 2) ^ 2", new[] { 0.37, 1.4, -0.83 })] + [InlineData("x ^ 3 / (x ^ 2 + 1) ^ 2", new[] { 0.37, 1.4, -0.83 })] + [InlineData("(x ^ 3 + x + 1) / (x ^ 2 + 1) ^ 3", new[] { 0.37, 1.4, -0.83 })] + [InlineData("(x ^ 2 + x) / (x ^ 2 + 3) ^ 2", new[] { 0.37, 1.4, -0.83 })] + public void APolynomialNumeratorIsDividedByTheQuadratic(string integrand, double[] points) => + AssertIsAntiderivative(integrand, points); + + /// + /// The one case checked exactly rather than at points: the derivative of the answer is + /// the integrand, with a residual of zero, which is the identity itself rather than a + /// consequence of it holding at the places sampled. + /// + /// + /// The zero comes back conditional, and that is right rather than a shortfall. + /// Both the integrand and its antiderivative are undefined where 1 + x^2 = 0 — + /// off the real line, at +-i — and "these two agree everywhere" would be a + /// claim about points where neither has a value. Asserted on the node rather than on + /// the printed form, so that a change of spelling does not read as a change of answer. + /// + [Fact] + public void TheAnswerDifferentiatesBackToAConditionalZero() + { + var antiderivative = "1 / (x ^ 2 + 1) ^ 2".ToEntity().Integrate("x").Substitute("C", 0); + var residual = (antiderivative.Differentiate("x") - "1 / (x ^ 2 + 1) ^ 2".ToEntity()).Simplify(); + + var zero = residual is Entity.Providedf(var expression, _) ? expression : residual; + Assert.Equal(Entity.Number.Integer.Create(0), zero); + } + + /// + /// A numerator that is already a multiple of the denominator's derivative leaves no + /// second integral, and the answer carries no condition. Multiplying the reduction by + /// zero instead of dropping it would attach one, since 0 * (2x/(x^2 + 1)) is + /// 0 provided not 1 + x^2 = 0 and not 0. + /// + [Fact] + public void ANumeratorAlongTheDerivativeIsOwedNoCondition() + { + var answer = "x / (x ^ 2 + 1) ^ 2".ToEntity().Integrate("x").Stringize(); + Assert.DoesNotContain("provided", answer); + } + + // A symbolic quadratic cannot have its discriminant's sign decided, so every arm + // survives and the answer is the piecewise. Without the two arms the reduction cannot + // speak for -- a zero leading coefficient and a zero discriminant -- those inputs would + // carry a division by zero into the answer. + [Theory] + [InlineData("1 / (a * x ^ 2 + c) ^ 2")] + [InlineData("1 / (a * x ^ 2 + b * x + c) ^ 2")] + public void ASymbolicQuadraticIsAnsweredAsAPiecewise(string integrand) + { + var answer = integrand.ToEntity().Integrate("x").Stringize(); + Assert.DoesNotContain("integral(", answer); + Assert.Contains("piecewise", answer); + Assert.Contains("arctan", answer); + } + + // Already answered before this rule, by the power rule and by partial fractions, and + // still answered: a denominator with real roots comes apart into linear factors, and a + // perfect square is a linear factor to twice the power. + [Theory] + [InlineData("1 / (x ^ 2 - 1) ^ 2", new[] { 1.4, 2.6, -3.1 })] + [InlineData("1 / (x + 1) ^ 2", new[] { 0.37, 1.4, 2.6 })] + [InlineData("1 / (x ^ 2 + 2 * x + 1) ^ 2", new[] { 0.37, 1.4, 2.6 })] + [InlineData("1 / (x ^ 2 - 2) ^ 2", new[] { 0.37, 2.6, -3.1 })] + public void ADenominatorThatFactorsIsStillAnswered(string integrand, double[] points) => + AssertIsAntiderivative(integrand, points); + + /// + /// An improper fraction is declined by this rule rather than divided out in it — the + /// recursion bottoms out at a single power with a numerator still of degree two, which + /// is the case SolveByPartialFractions already opens with a long division. It is + /// answered; the point is that it is answered there. + /// + [Fact] + public void AnImproperFractionIsAnsweredByTheDivisionThatAlreadyExisted() => + AssertIsAntiderivative("x ^ 4 / (x ^ 2 + 1) ^ 2", 0.37, 1.4, -0.83); + } +}