diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index e6e6df82f..ad8d396dc 100644 --- a/BREAKING-CHANGES.md +++ b/BREAKING-CHANGES.md @@ -1185,6 +1185,24 @@ arctangent it is ([#718](https://github.com/asc-community/AngouriMath/issues/718 | `"(A + B*x + C*x^2)*sqrt(c + d*x)/((a + b*x)^3*sqrt(e + f*x))".ToEntity().Integrate("x")` | `integral(...)` | the roots over powers of `a + b x`, a logarithm, and an arctangent | | `"x^2*sqrt(a + 2*x)*sqrt(c - 3*x)".ToEntity().Integrate("x")` | `integral(...)` | a cubic times the two roots, and an arctangent of `sqrt(6)/2 sqrt(a + 2 x)/sqrt(c - 3 x)` | +### A negative power of a linear beside a root of another is integrated by the recurrence on the power + +**Answers where there were none, and faster ones, written differently.** `1/((a + b x)^4 sqrt(c + d x))` +was declined by 2.5.0, and the unreleased master answered it in fourteen seconds, by the substitution +for the root, over a power of a quadratic with symbols in it; with a cubic above it, past a minute. +A whole power `m <= -2` of one linear beside a power that is not whole of another, times a +polynomial, is integrated now by the recurrence +`int L1^m L2^n = L1^(m + 1) L2^(n + 1)/((m + 1) D) - (m + n + 2) d/((m + 1) D) int L1^(m + 1) L2^n`, +`D = b c - a d`, down to `int L2^n/L1`, the polynomial written in powers of `L1` first. The +answers are the same functions as master's where master gave one, written as powers of the two +linears ([#718](https://github.com/asc-community/AngouriMath/issues/718)). + +| Input | Was (2.5.0) | Now | +|---|---|---| +| `"1/((a + b*x)^4*sqrt(c + d*x))".ToEntity().Integrate("x")` | `integral(...)` | powers of `a/b + x` beside `sqrt(c + d x)`, and the arctangent or logarithm `1/((a + b x) sqrt(c + d x))` integrates to | +| `"(p + q*x + r*x^2 + s*x^3)/((a + b*x)^5*sqrt(c + d*x))".ToEntity().Integrate("x")` | `integral(...)` | the same, term by term | +| `"x/((a + b*x)^4*sqrt(c + d*x))".ToEntity().Integrate("x")` | `integral(...)` | the same | + ### A polynomial over a power of a linear beside the root of a quadratic is integrated `P(x) Q^(m/2)/(g + h x)^k`, a polynomial over a power of a linear beside a half-odd power of a diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs index c367f55ed..0de0501a1 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs @@ -9424,6 +9424,154 @@ private static Entity WithWholePowersOfProductsDistributed(Entity integrand) return integrand; } + /// + /// A whole power m <= -2 of a linear beside a power that is not whole of another, + /// times a polynomial: 1/((a + b x)^4 sqrt(c + d x)), + /// (A + B x + C x^2)/((a + b x)^3 (c + d x)^(3/2)). The polynomial is written in powers + /// of the first linear, and each negative power is raised to -1 by the recurrence + /// int L1^m L2^n = L1^(m + 1) L2^(n + 1)/((m + 1) D) - (m + n + 2) d/((m + 1) D) int L1^(m + 1) L2^n, + /// L1 = a + b x, L2 = c + d x, D = b c - a d; what is left are the + /// questions int L2^n/L1 and int L1^k L2^n for k >= 0, which the + /// integrator answers in a step. + /// + /// + /// + /// The recurrence is the derivative of L1^(m + 1) L2^(n + 1), which is + /// L1^m L2^n ((m + 1) D + (m + n + 2) d L1), read back. Under v = sqrt(c + d x) + /// the same integrand is a rational function over a power of a quadratic with symbols in + /// it, and the reduction of that took four seconds at the square and twenty at the fourth + /// power; and the partial fractions of a rational function beside a root are terms of + /// exactly this shape, so it is their cost too. + /// + /// + /// A closed reduction to smaller questions of another shape, so at any depth. + /// https://github.com/asc-community/AngouriMath/issues/718 + /// + /// + internal static Entity? SolveANegativePowerOfALinearBesideAPowerOfAnotherByTheRecurrence(Entity expr, Entity.Variable x, bool integrateByParts) + { + Entity constant = Number.Integer.One; + var linears = new List<(Entity Base, Entity Exponent)>(); + Entity polynomial = Number.Integer.One; + // The factors with their exponents, a whole power of a product read through: the question + // comes back as `((a/b + x)^4 sqrt(c + d x))^(-1)` once a constant is out of the bar. + var factors = new List<(Entity Base, Entity Exponent)>(); + void Gather(Entity node, int times) + { + switch (node) + { + case Mulf(var left, var right): + Gather(left, times); + Gather(right, times); + break; + case Divf(var numerator, var denominator): + Gather(numerator, times); + Gather(denominator, -times); + break; + case Powf(Mulf or Divf or Powf, Number.Integer whole) when whole.EInteger.CanFitInInt32() && System.Math.Abs(whole.EInteger.ToInt32Checked()) <= MaximumRecurrenceDepth: + Gather(((Powf)node).Base, times * whole.EInteger.ToInt32Checked()); + break; + case Powf(var powerBase, var power): + factors.Add((powerBase, times == 1 ? power : (power * times).InnerSimplified)); + break; + default: + factors.Add((node, Number.Integer.Create(times))); + break; + } + } + Gather(expr, 1); + foreach (var (@base, exponent) in factors) + { + if (!@base.ContainsNode(x)) + { + constant *= MathS.Pow(@base, exponent); + continue; + } + if (exponent.ContainsNode(x)) + return null; + if (TreeAnalyzer.TryGetPolyLinear(@base, x, out var slope, out _) && !TreeAnalyzer.IsZero(slope)) + { + var index = linears.FindIndex(pair => pair.Base == @base); + if (index >= 0) + linears[index] = (@base, (linears[index].Exponent + exponent).InnerSimplified); + else + linears.Add((@base, exponent)); + continue; + } + if (exponent is not Number.Integer { EInteger.Sign: > 0 } whole || !whole.EInteger.CanFitInInt32()) + return null; + polynomial *= whole == Number.Integer.One ? @base : MathS.Pow(@base, whole); + } + // The whole power below the bar, the power that is not whole, and the linears that are + // whole powers above it, which belong to the polynomial. + (Entity Base, int Power)? below = null; + (Entity Base, Entity Exponent)? beside = null; + foreach (var (@base, exponent) in linears) + { + if (exponent is Number.Integer whole && whole.EInteger.CanFitInInt32()) + { + var power = whole.EInteger.ToInt32Checked(); + if (power >= 0) + polynomial *= power == 1 ? @base : MathS.Pow(@base, power); + else if (power <= -2 && power >= -MaximumRecurrenceDepth && below is null) + below = (@base, power); + else + return null; + } + else if (beside is null && exponent is not Number.Integer) + beside = (@base, exponent); + else + return null; + } + if (below is not var (first, m) || beside is not var (second, n)) + return null; + if (!TreeAnalyzer.TryGetPolyLinear(first, x, out var b, out var a) || !TreeAnalyzer.TryGetPolyLinear(second, x, out var d, out var c)) + return null; + var determinant = (b * c - a * d).InnerSimplified; + if (determinant.Evaled is Number.Complex { IsZero: true }) + return null; + // The polynomial in powers of the first linear: x = (t - a)/b. + var t = Variable.CreateUnique(expr, "t_linear"); + if (!TreeAnalyzer.TryGetPolynomial(polynomial.Substitute(x, (t - a) / b), t, out var inT) + || inT.Keys.Any(power => power.Sign < 0 || !power.CanFitInInt32())) + return null; + // int L1^k L2^n for k from -1 down to m, each from the one above it. + var lowest = m + inT.Keys.Min(power => power.ToInt32Checked()); + var raised = new Dictionary(); + if (lowest <= -1) + { + if (Integration.ComputeIndefiniteIntegral(MathS.Pow(second, n) / first, x, integrateByParts) is not { } overTheFirst + || overTheFirst.Nodes.Any(node => node is Integralf)) + return null; + raised[-1] = overTheFirst; + for (var k = -2; k >= lowest; k--) + raised[k] = MathS.Pow(first, k + 1) * MathS.Pow(second, n + 1) / ((k + 1) * determinant) + - (k + n + 2) * d / ((k + 1) * determinant) * raised[k + 1]; + } + Entity answer = Number.Integer.Zero; + foreach (var pair in inT) + { + var (k, coefficient) = (m + pair.Key.ToInt32Checked(), pair.Value); + if (!raised.TryGetValue(k, out var integral)) + { + var question = k == 0 ? MathS.Pow(second, n) : MathS.Pow(first, k) * MathS.Pow(second, n); + if (Integration.ComputeIndefiniteIntegral(question, x, integrateByParts) is not { } answered + || answered.Nodes.Any(node => node is Integralf)) + return null; + integral = answered; + } + answer += coefficient.InnerSimplified * integral; + } + return constant * answer; + } + + /// + /// The most negative whole power + /// takes down + /// the recurrence: each step is one term more in the answer. + /// + private const int MaximumRecurrenceDepth = 12; + /// /// An integrand holding a fractional power of something linear in the variable, /// turned into a rational function by u^q = a*x + b. diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs b/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs index 223a8a250..fe96dd383 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs @@ -758,6 +758,10 @@ private static Entity Normalized(Entity expr, Entity.Variable x) => // degree-8 polynomial by expanding, both right and only one worth reading. Three // tests pinned exactly that and caught this when the split was put in front of the // substitution wholesale. The expanding call stays where it was, below. + // A whole power below -1 of a linear beside a power that is not whole of another, by the + // recurrence that raises the whole power to -1: closed, where the substitution for the + // root makes a rational function over a power of a symbolic quadratic. + if ((answer = IndefiniteIntegralSolver.SolveANegativePowerOfALinearBesideAPowerOfAnotherByTheRecurrence(expr, x, integrateByParts)) is { }) return answer; // A product of powers times a sum that is the derivative of the product with some // 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. diff --git a/Sources/Tests/UnitTests/Calculus/LinearPowerRecurrenceIntegralTest.cs b/Sources/Tests/UnitTests/Calculus/LinearPowerRecurrenceIntegralTest.cs new file mode 100644 index 000000000..0f3b70e78 --- /dev/null +++ b/Sources/Tests/UnitTests/Calculus/LinearPowerRecurrenceIntegralTest.cs @@ -0,0 +1,48 @@ +// +// 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 whole power below -1 of a linear beside a power that is not whole of another, times + /// a polynomial, by the recurrence that raises the whole power to -1. Rubi's + /// P(x) (a + b x)^m (c + d x)^n files. + /// #718 + /// + [Trait("Area", "Calculus")] + public sealed class LinearPowerRecurrenceIntegralTest + { + [Theory] + [InlineData("1/((a + b*x)^4*sqrt(c + d*x))")] + [InlineData("(p + q*x + r*x^2 + s*x^3)/((a + b*x)^5*sqrt(c + d*x))")] + [InlineData("(c + d*x)^(3/2)/(a + b*x)^3")] + [InlineData("1/(x^3*(c + d*x)^(5/2))")] + [InlineData("x/((a + b*x)^4*sqrt(c + d*x))")] + public void ByTheRecurrence(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("p", 0.4).Substitute("q", 1.3).Substitute("r", 0.8) + .Substitute("s", 0.6); + 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}"); + } + } + } +}