diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index 09c9c3b8a..cdc3c4a0a 100644 --- a/BREAKING-CHANGES.md +++ b/BREAKING-CHANGES.md @@ -232,6 +232,20 @@ Rubi's `x^m (a + b x^n)^p` and `(a + b x^n)^p (c + d x^n)^q` | `"1/(2+3/x^2)^3".Integrate("x")`, `1/(a + b/x^3)` | left unevaluated | an antiderivative | | `"(a+b*x^n)*(c+d*x^n)^3".Integrate("x")` | left unevaluated | written out, eight powers of `x` integrated: `a c^3 x + ... + b d^3 x^(4 n + 1)/(4 n + 1)` | +### A polynomial with symbols in it that is a binomial or an even quartic in a shifted variable is integrated in it + +**Answers where there were none.** `1/(c^2 x^3 + 3 b c x^2 + 3 b^2 x + 3 a b)` was left unevaluated, +although its denominator is `((c x + b)^3 + 3 a b c - b^3)/c`, which the rule for a binomial reads at +once: nothing factors a polynomial with a symbol among its coefficients. A polynomial below the bar +of the third to the sixth degree, with a symbol in it, that written in `y = x + s` with +`s = a_(n-1)/(n a_n)` is a binomial, or a quartic even in `y`, is integrated in `y`. Rubi's 1.3.1 +([#718](https://github.com/asc-community/AngouriMath/issues/718)). + +| Input | Was (2.5.0) | Now | +|---|---|---| +| `"1/(3*a*b + 3*b^2*x + 3*b*c*x^2 + c^2*x^3)".ToEntity().Integrate("x")`, and its square | `integral(...)` | two logarithms and an arctangent in `x + b/c` | +| `"x/(a + 8*x - 8*x^2 + 4*x^3 - x^4)".ToEntity().Integrate("x")`, and `1` over it | `integral(...)` | the antiderivative in `x - 1` | + ### The hyperbolic functions have antiderivatives, and so does anything rational in `e^(k x)` An integrand rational in `e^(k x)` becomes a rational function of one variable under diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs index 01e0779d5..60538086d 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs @@ -745,6 +745,17 @@ is var (multiple, leftover) ?? Integration.ComputeIndefiniteIntegral(numerator / respelled, x, integrateByParts)) is { } overPolynomials) return overPolynomials; + // A polynomial below the bar that is a binomial, or a quartic even in its variable, once + // written in `y = x + s` with `s = a_(n-1)/(n a_n)`: `c^2 x^3 + 3 b c x^2 + 3 b^2 x + 3 a b` + // is `((c x + b)^3 + 3 a b c - b^3)/c`, and `a + 8 x - 8 x^2 + 4 x^3 - x^4` is + // `a + 3 - 2 y^2 - y^4` in `y = x - 1`. Nothing above factors a polynomial with a symbol + // in it, and the rules for a binomial and for an even quartic read it in y at once. + // Once: in y the term that s takes away is gone. + if (InTheVariableThatDepressesIt(numerator, denominator, x) is var (inY, y, shift) + && (SolveByPartialFractions(inY, y, integrateByParts) + ?? Integration.ComputeIndefiniteIntegral(inY, y, integrateByParts)) is { } inTheShiftedVariable) + return inTheShiftedVariable.Substitute(y, (x + shift).InnerSimplified); + // A denominator with a written repeated factor takes the Hermite reduction first: // the rational part of the answer in one linear solve, and what is left is a proper // fraction over a squarefree denominator for the splits below. `(1 + x^2)/(x (1 + x^3)^2)` @@ -850,6 +861,79 @@ is var (multiple, leftover) return null; } + /// + /// over written in y = x + s, + /// where the denominator is a constant times a power of one polynomial in x of the third to + /// the sixth degree, written as a sum of its monomials, that in y is a binomial or, of the + /// fourth degree, even; s = a_(n-1)/(n a_n) is what takes its term of degree + /// n - 1 away. Null otherwise, and where that term is not there to take away. + /// + private static (Entity InY, Entity.Variable Y, Entity Shift)? InTheVariableThatDepressesIt(Entity numerator, Entity denominator, Entity.Variable x) + { + Entity? polynomial = null; + var power = 0; + Entity constant = Number.Integer.One; + foreach (var factor in Mulf.LinearChildren(denominator)) + { + if (!factor.ContainsNode(x)) + { + constant *= factor; + continue; + } + if (polynomial is not null) + return null; + (polynomial, power) = factor is Powf(var raised, Number.Integer { EInteger.Sign: > 0 } exponent) && exponent.EInteger.CanFitInInt32() + ? (raised, exponent.EInteger.ToInt32Unchecked()) : (factor, 1); + } + // With a symbol among its coefficients, since a polynomial over the rationals is split + // over its factors above; and written out, since one written in a linear, + // `a + (b + c x)^3`, is read as it stands. + if (polynomial is null || !polynomial.Vars.Any(symbol => symbol != x) + || polynomial.Nodes.Any(node => node is Powf(Sumf or Minusf, _) && node.ContainsNode(x)) + || !TreeAnalyzer.TryGetPolynomial(polynomial, x, out var terms) + || terms.Any(term => term.Key.Sign < 0 || !term.Key.CanFitInInt32() || term.Value.ContainsNode(x))) + return null; + if (numerator.ContainsNode(x) + && (!TreeAnalyzer.TryGetPolynomial(numerator, x, out var above) || above.Any(term => term.Key.Sign < 0 || term.Value.ContainsNode(x)))) + return null; + var degree = terms.Keys.Max()!.ToInt32Unchecked(); + if (degree < 3 || degree > 6) + return null; + Entity Coefficient(int k) => terms.TryGetValue(EInteger.FromInt32(k), out var coefficient) ? coefficient : Number.Integer.Zero; + if (VanishesIdentically(Coefficient(degree)) || VanishesIdentically(Coefficient(degree - 1))) + return null; + var shift = Functions.PartialFractions.InLowestTermsOverTheSymbols(Coefficient(degree - 1) / (degree * Coefficient(degree))); + // P(y - s) has at y^i the sum over m >= i of a_m C(m, i) (-s)^(m - i). + var inY = new Entity[degree + 1]; + for (var i = 0; i <= degree; i++) + { + Entity sum = Number.Integer.Zero; + var choose = EInteger.One; + for (var m = i; m <= degree; m++) + { + if (m > i) + choose = choose.Multiply(EInteger.FromInt32(m)).Divide(EInteger.FromInt32(m - i)); + sum += Coefficient(m) * Number.Integer.Create(choose) * (m == i ? Number.Integer.One : MathS.Pow(-shift, m - i)); + } + inY[i] = Functions.PartialFractions.InLowestTermsOverTheSymbols(sum); + } + var vanishing = inY.Select(VanishesIdentically).ToArray(); + var aBinomial = Enumerable.Range(1, degree - 1).All(i => vanishing[i]); + var anEvenQuartic = degree == 4 && vanishing[1] && vanishing[3]; + if (!aBinomial && !anEvenQuartic) + return null; + var y = Variable.CreateUnique(numerator + denominator, "y_shift"); + Entity depressed = Number.Integer.Zero; + for (var i = degree; i >= 0; i--) + if (!vanishing[i]) + { + var monomial = i == 0 ? inY[i] : inY[i] * (i == 1 ? y : MathS.Pow(y, i)); + depressed = depressed == Number.Integer.Zero ? monomial : depressed + monomial; + } + var below = power == 1 ? depressed : MathS.Pow(depressed, power); + return (numerator.Substitute(x, y - shift) / (constant == Number.Integer.One ? below : constant * below), y, shift); + } + /// /// A polynomial over a power of a linear, P(x)/(a + b x)^k with k >= 2, /// under t = a + b x: P((t - a)/b) t^(-k) / b is a sum of powers of diff --git a/Sources/Tests/UnitTests/Calculus/ShiftedPolynomialIntegralTest.cs b/Sources/Tests/UnitTests/Calculus/ShiftedPolynomialIntegralTest.cs new file mode 100644 index 000000000..7c98be466 --- /dev/null +++ b/Sources/Tests/UnitTests/Calculus/ShiftedPolynomialIntegralTest.cs @@ -0,0 +1,61 @@ +// +// 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 System; +using AngouriMath.Extensions; +using Xunit; + +namespace AngouriMath.Tests.Calculus +{ + /// + /// A polynomial below the bar with a symbol among its coefficients that is a binomial, or a + /// quartic even in its variable, once written in y = x + s with + /// s = a_(n-1)/(n a_n): Rubi's 1.3.1. Nothing factors a polynomial with a symbol in it, + /// and in y the rules for a binomial and for an even quartic read it at once. + /// #718 + /// + /// + /// Checked by differentiating back with the symbols pinned after the integral is taken. + /// + [Trait("Area", "Calculus")] + public sealed class ShiftedPolynomialIntegralTest + { + private static readonly double[] Points = { 0.3, 0.6, 0.9 }; + + private static void DifferentiatesBack(string integrand) + { + var integral = integrand.ToEntity().Integrate("x"); + Assert.DoesNotContain("integral(", integral.Stringize()); + Entity Pin(Entity e) => e.Substitute("a", 0.7).Substitute("b", 1.3).Substitute("c", 2.1); + var derivative = Pin(integral.Substitute("C", 0).Differentiate("x")); + var original = Pin(integrand.ToEntity()); + foreach (var point in Points) + { + var expected = original.Substitute("x", point).EvalNumerical().RealPart.EDecimal.ToDouble(); + var actual = derivative.Substitute("x", point).EvalNumerical().RealPart.EDecimal.ToDouble(); + Assert.True(Math.Abs(expected - actual) < 1e-8 * Math.Max(1, Math.Abs(expected)), + $"d/dx of the antiderivative of {integrand} is {actual} at x = {point}, where the integrand is {expected}"); + } + } + + /// + /// c^2 x^3 + 3 b c x^2 + 3 b^2 x + 3 a b is ((c x + b)^3 + 3 a b c - b^3)/c. + /// + [Theory] + [InlineData("1/(3*a*b + 3*b^2*x + 3*b*c*x^2 + c^2*x^3)")] + [InlineData("1/(3*a*b + 3*b^2*x + 3*b*c*x^2 + c^2*x^3)^2")] + public void ACubicThatIsABinomialInALinear(string integrand) => DifferentiatesBack(integrand); + + /// + /// a + 8 x - 8 x^2 + 4 x^3 - x^4 is a + 3 - 2 y^2 - y^4 in y = x - 1. + /// + [Theory] + [InlineData("x/(a + 8*x - 8*x^2 + 4*x^3 - x^4)")] + [InlineData("1/(a + 8*x - 8*x^2 + 4*x^3 - x^4)")] + public void AQuarticThatIsEvenInALinear(string integrand) => DifferentiatesBack(integrand); + } +}