diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index e6e6df82f..8052a83ad 100644 --- a/BREAKING-CHANGES.md +++ b/BREAKING-CHANGES.md @@ -1168,6 +1168,24 @@ in both quadratics, the cube's coefficients grow past what the reduction takes, | `"(7 + 13*x)/((5 + x + 2*x^2)^3*sqrt(2 + x + 3*x^2))".ToEntity().Integrate("x")` | `integral(...)` | an antiderivative | | `"sqrt(a + a*sec(x))/(c + d*sec(x))^2".ToEntity().Integrate("x")` | `integral(...)` | an antiderivative in `tan(x/2)` | +### Powers of two linears whose exponents sum to a whole number below `-2` are integrated + +**Answers where there were none.** `(a + b x)^m (c + d x)^(-3 - m)`, `x^(n - 4)/(a + b x)^n` and +the rest of Rubi's `(a + b x)^m (c + d x)^n` problems whose exponents sum to a whole number `-k` +below `-2`, beside a polynomial of degree at most `k - 2`, were declined; the sum `-2` alone was +answered. Under `t = (a + b x)/(c + d x)` each is a power of `t` beside a polynomial in `t`, so the +antiderivative is the two powers as written times a polynomial in `x`, exact wherever the +integrand is defined and taking each `m + i` it divides by to be nonzero, as the integrator does +everywhere. Numeric exponents that sum so go the same way, so an answer master already gave for +one, `(a + b x)^(5/2)/(c + d x)^(11/2)`, is written so now, the same function +([#718](https://github.com/asc-community/AngouriMath/issues/718)). + +| Input | Was (2.5.0) | Now | +|---|---|---| +| `"(a + b*x)^m*(c + d*x)^(-3 - m)".ToEntity().Integrate("x")` | `integral(...)` | the two powers times `-(d (a + b x)^2 (c + d x)/(m + 2) - b (a + b x) (c + d x)^2/(m + 1))/(a d - b c)^2` | +| `"x^(n - 4)/(a + b*x)^n".ToEntity().Integrate("x")` | `integral(...)` | the two powers times a cubic in `x` and `a + b x`, over `a^3` | +| `"(a + b*x)^(5/2)/(c + d*x)^(11/2)".ToEntity().Integrate("x")` | `integral(...)` | the two powers times `-(2 d (a + b x)^2 (c + d x)/9 - 2 b (a + b x) (c + d x)^2/7)/(a d - b c)^2` | + ### A polynomial beside the square roots of two linears with symbols in them is integrated `P(x) sqrt(c + d x) sqrt(e + f x)`, a polynomial beside the roots of two linears with symbols in diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs index c367f55ed..0ef9bca59 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs @@ -11936,6 +11936,129 @@ static bool IsRationalInBoth(Entity e, Entity.Variable x, Entity.Variable L) return null; } + /// + /// Powers of two linears whose exponents sum to a whole number -k, k >= 2, + /// beside a polynomial of degree at most k - 2: (a + b x)^m (c + d x)^(-3 - m), + /// x^(n - 4)/(a + b x)^n, (e + f x)^2 (a + b x)^m (c + d x)^(-4 - m). Under + /// t = L1/L2 each is t^A beside a polynomial in t, so the answer is + /// the two powers as written times a polynomial in x. + /// + /// + /// + /// With L1 = a1 + b1 x, L2 = a2 + b2 x and D = a1 b2 - a2 b1, the + /// substitution makes x = (a1 - a2 t)/(b2 t - b1), L2 = D/(b2 t - b1) and + /// dx = -D dt/(b2 t - b1)^2, so L1^A L2^B P(x) dx is, up to a factor constant + /// on each interval, -D^(1 - k) t^A S(t) dt for the polynomial + /// S(t) = sum_j p_j (a1 - a2 t)^j (b2 t - b1)^(k - 2 - j). Its antiderivative is + /// sum_i s_i t^(A + i + 1)/(A + i + 1), and since t^(A + i + 1) is + /// t^A t^(i + 1) for a whole i + 1, the constant factor and t^A come + /// back as the two powers as written: the answer is + /// L1^A L2^B sum_i -D^(1 - k) s_i L1^(i + 1) L2^(k - 1 - i)/(A + i + 1), exact + /// wherever the integrand is defined. The generic case, as the integrator answers + /// everywhere: A + i + 1 is taken to be nonzero. + /// + /// + /// Only the sum -2 was answered, as the derivative of the product of the two powers + /// raised by one; Rubi's (a + b x)^m (c + d x)^n files hold the rest, and under + /// u = sin(y) so do the powers of 1 ± sin(y) beside one another. A closed + /// rule, so at any depth. https://github.com/asc-community/AngouriMath/issues/718 + /// + /// + internal static Entity? SolveTwoLinearPowersWhoseExponentsSumToAWholeNumber(Entity expr, Entity.Variable x) + { + Entity constant = Number.Integer.One; + var powers = new List<(Entity Base, Entity Exponent)>(); + var polynomial = new List<(Entity Factor, bool Underneath)>(); + foreach (var (factor, underneath) in FactorsOfTheIntegrand(expr)) + { + if (!factor.ContainsNode(x)) + { + constant = underneath ? constant / factor : constant * factor; + continue; + } + var (@base, exponent) = factor is Powf(var powerBase, var power) ? (powerBase, power) : (factor, (Entity)Number.Integer.One); + if (exponent.ContainsNode(x)) + return null; + if (underneath) + exponent = (-exponent).InnerSimplified; + var index = powers.FindIndex(pair => pair.Base == @base); + if (index >= 0) + powers[index] = (@base, (powers[index].Exponent + exponent).InnerSimplified); + else if (exponent is Number.Integer) + polynomial.Add((factor, underneath)); + else + powers.Add((@base, exponent)); + } + if (powers.Count != 2) + return null; + // A whole power of a base that turned out to be one of the two goes into its exponent. + for (var i = polynomial.Count - 1; i >= 0; i--) + { + var (factor, underneath) = polynomial[i]; + var (@base, exponent) = factor is Powf(var powerBase, Number.Integer power) ? (powerBase, (Entity)power) : (factor, (Entity)Number.Integer.One); + var index = powers.FindIndex(pair => pair.Base == @base); + if (index < 0) + { + if (underneath) + return null; + continue; + } + powers[index] = (@base, (powers[index].Exponent + (underneath ? -exponent : exponent)).InnerSimplified); + polynomial.RemoveAt(i); + } + var ((first, a), (second, b)) = (powers[0], powers[1]); + if (a is Number.Integer || b is Number.Integer) + return null; + if ((a + b).Simplify() is not Number.Integer { EInteger: var sum } || sum.CompareTo(EInteger.FromInt32(-2)) > 0 + || sum.CompareTo(EInteger.FromInt32(-MaximumLinearPowerSum)) < 0) + return null; + var k = -sum.ToInt32Checked(); + if (!TreeAnalyzer.TryGetPolyLinear(first, x, out var b1, out var a1) || !TreeAnalyzer.TryGetPolyLinear(second, x, out var b2, out var a2)) + return null; + var d = (a1 * b2 - a2 * b1).InnerSimplified; + if (d.Evaled is Number.Complex { IsZero: true }) + return null; + Entity product = Number.Integer.One; + foreach (var (factor, _) in polynomial) + product *= factor; + if (!TreeAnalyzer.TryGetPolynomial(product, x, out var coefficients) + || coefficients.Keys.Any(power => power.Sign < 0 || power.CompareTo(EInteger.FromInt32(k - 2)) > 0)) + return null; + var t = Variable.CreateUnique(expr, "t_quotient"); + // Whole powers written out, so that what is read as a polynomial in t is one. + static Entity Raised(Entity @base, int power) => power switch + { + 0 => Number.Integer.One, + 1 => @base, + _ => MathS.Pow(@base, power), + }; + Entity inT = Number.Integer.Zero; + foreach (var pair in coefficients) + { + var j = pair.Key.ToInt32Checked(); + inT += pair.Value * Raised(a1 - a2 * t, j) * Raised(b2 * t - b1, k - 2 - j); + } + if (!TreeAnalyzer.TryGetPolynomial(inT.InnerSimplified, t, out var ofT)) + return null; + Entity polynomialInX = Number.Integer.Zero; + foreach (var pair in ofT) + { + var (i, s) = (pair.Key.ToInt32Checked(), pair.Value); + if (i < 0 || i > k - 2) + return null; + polynomialInX += s / (a + (i + 1)) * Raised(first, i + 1) * Raised(second, k - 1 - i); + } + var answer = (-constant * MathS.Pow(d, 1 - k)).InnerSimplified * MathS.Pow(first, a) * MathS.Pow(second, b) * polynomialInX.InnerSimplified; + return answer.Nodes.Any(node => node == MathS.NaN) ? null : answer; + } + + /// + /// The most negative sum of the two exponents + /// takes: the polynomial + /// it writes has k - 1 terms, each a product of powers summing to k. + /// + private const int MaximumLinearPowerSum = 12; + /// /// The constant k with equal to k times /// , read off one monomial with both written as polynomials diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs b/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs index 223a8a250..9fd6f149c 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs @@ -762,6 +762,9 @@ private static Entity Normalized(Entity expr, Entity.Variable x) => // of the powers raised by one -- `e^x x^2 ln(x)^2 (3 + (3 + x) ln(x))` -- read off // the sum, for symbolic exponents too. Before the split, which loses it. if ((answer = IndefiniteIntegralSolver.SolveAsTheDerivativeOfAProductOfPowers(expr, x)) is { }) return answer; + // Powers of two linears whose exponents sum to a whole number below -2, beside a polynomial of + // low degree, by t = L1/L2: a power of t beside a polynomial. The rule above reads the sum -2. + if ((answer = IndefiniteIntegralSolver.SolveTwoLinearPowersWhoseExponentsSumToAWholeNumber(expr, x)) is { }) return answer; // A rational function of x and of a tower of exponentials and logarithms over it, // by the Risch-Norman ansatz: a rational function of the same plus logarithms of // the denominator's factors, the coefficients unknown and solved for. Before the diff --git a/Sources/Tests/UnitTests/Calculus/TwoLinearPowersIntegralTest.cs b/Sources/Tests/UnitTests/Calculus/TwoLinearPowersIntegralTest.cs new file mode 100644 index 000000000..987c520a1 --- /dev/null +++ b/Sources/Tests/UnitTests/Calculus/TwoLinearPowersIntegralTest.cs @@ -0,0 +1,50 @@ +// +// 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 +{ + /// + /// Powers of two linears whose exponents sum to a whole number below -2, beside a + /// polynomial of degree at most two less than its negative: under t = L1/L2 a power of + /// t beside a polynomial, so the antiderivative is the two powers times a polynomial. + /// Rubi's (a + b x)^m (c + d x)^n files. + /// #718 + /// + [Trait("Area", "Calculus")] + public sealed class TwoLinearPowersIntegralTest + { + [Theory] + [InlineData("(a + b*x)^m*(c + d*x)^(-3 - m)")] + [InlineData("x^(n - 4)/(a + b*x)^n")] + [InlineData("(a + b*x)^m*(c + d*x)^(-5 - m)*(p + q*x)^3")] + [InlineData("(a + b*x)*(c + d*x)^(n - 4)/(p + q*x)^n")] + [InlineData("(a + b*x)^m*(c + d*x)^(-4 - m)*(p + q*x)*(g + h*x)")] + [InlineData("(a + b*x)^(5/2)/(c + d*x)^(11/2)")] + public void UnderTheQuotientOfTheLinears(string integrand) + { + var integral = integrand.ToEntity().Integrate("x"); + Assert.DoesNotContain("integral(", integral.Stringize()); + Entity Pinned(Entity e) => e.Substitute("a", 1.3).Substitute("b", 0.9).Substitute("c", 0.7) + .Substitute("d", 1.1).Substitute("g", 0.8).Substitute("h", 0.6).Substitute("m", 0.37) + .Substitute("n", 1.21).Substitute("p", 0.4).Substitute("q", 1.3); + var derivative = Pinned(integral.Substitute("C", 0)).Differentiate("x"); + var original = Pinned(integrand.ToEntity()); + foreach (var at in new[] { -2.4, -1.5, 0.3, 0.8, 1.5, 2.4 }) + { + var want = original.Substitute("x", at).EvalNumerical(); + var got = derivative.Substitute("x", at).EvalNumerical(); + Assert.True(Math.Abs((double)(got - want).RealPart) + Math.Abs((double)(got - want).ImaginaryPart) + < 1e-9 * Math.Max(1, Math.Abs((double)want.RealPart) + Math.Abs((double)want.ImaginaryPart)), + $"d/dx of the antiderivative of {integrand} is {got} at x = {at}, where the integrand is {want}"); + } + } + } +}