diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index 5051710b2..ab2fe2684 100644 --- a/BREAKING-CHANGES.md +++ b/BREAKING-CHANGES.md @@ -559,6 +559,23 @@ of the quadratic, parts leaves the quadratic one power lower, which is answered | `"(b + 2*c*x)/((a + b*x + c*x^2)^2*sqrt(d + k*x))".ToEntity().Integrate("x")` | `integral(...)` | 1,075 characters | | `"(b + 2*c*x)*(d + k*x)^(3/2)/(a + b*x + c*x^2)^3".ToEntity().Integrate("x")` | `integral(...)` | 4,054 characters | +### A rational function over a power of a symbolic quadratic is reduced a power at a time + +**Answers where there were none.** `1/((d + k x) (a + b x + c x^2)^4)` ran past a minute, the partial fractions taking the +symbolic quadratic apart over its roots and its powers with them. The linear is taken off at its root now, and the power +of the quadratic is lowered one at a time by the recurrence for the integral of `1/Q^k`, down to the first power, which +is written in closed form. A quadratic that factors over linears is left to the partial fractions +([#718](https://github.com/asc-community/AngouriMath/issues/718)). + +| Input | Was (2.5.0) | Now | +|---|---|---| +| `"x^8/(a + b*x + c*x^2)^4".ToEntity().Integrate("x")` | `integral(...)` | 8,511 characters | +| `"1/((d + k*x)*(a + b*x + c*x^2)^4)".ToEntity().Integrate("x")` | `integral(...)` | 9,091 characters | +| `"1/((d + k*x)^2*(a + b*x + c*x^2)^5)".ToEntity().Integrate("x")` | `integral(...)` | 25,767 characters | +| `"1/((d + k*x)^2*(a + c*x^2)^4)".ToEntity().Integrate("x")` | `integral(...)` | 9,599 characters | +| `"(f + g*x)/((d + k*x)^2*(a + b*x + c*x^2)^3)".ToEntity().Integrate("x")` | `integral(...)` | 21,943 characters | +| `"(A + B*x + C*x^2)/((d + k*x)^3*(a + c*x^2)^3)".ToEntity().Integrate("x")` | `integral(...)` | 25,201 characters | + ### An odd power of a quadratic's derivative over a power of the quadratic is integrated in the quadratic **Shorter answers, and answers where there were none.** `(b + 2c x)^9/(a + b x + c x^2)^3` was answered through the diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs index 3a7a8a3b4..3f60a9126 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs @@ -19342,6 +19342,208 @@ Entity TakeOut(Entity node) return inTermsOfU.Substitute(u, Number.Integer.One / x); } + /// + /// A rational function over a whole power of a symbolic quadratic, beside at most one + /// linear to a whole power, with the power of the quadratic lowered one at a time rather + /// than the quadratic taken apart: the linear comes off first, as + /// N/(L^j Q^n) = N(p)/(Q(p)^n L^j) + M/(L^(j - 1) Q^n) at its root p, and + /// what is left over Q^k is S/Q^(k - 1) + (alpha Q' + beta)/Q^k, whose second + /// part is -alpha/((k - 1) Q^(k - 1)) and the recurrence + /// integral 1/Q^k = (2 c x + b)/((k - 1) D Q^(k - 1)) + 2 c (2 k - 3)/((k - 1) D) integral 1/Q^(k - 1) + /// with D = 4 a c - b^2, down to the first power, which is the integrator's. + /// + /// + /// Rubi's 1.2.1.2 1/((d + e x) (a + b x + c x^2)^4) and x^8/(a + b x + c x^2)^4 + /// ran past a minute: the partial fractions take a symbolic quadratic apart over its + /// roots, and its powers with it. D and Q(p) are named while the reduction + /// runs and written back after, so that each stays one atom. + /// https://github.com/asc-community/AngouriMath/issues/718 + /// + internal static Entity? SolveARationalFunctionOverAPowerOfASymbolicQuadratic(Entity expr, Entity.Variable x, bool integrateByParts) + { + // At the question asked only: below a substitution the routes it pre-empts answered + // `(a + b x + c x^2)^3/((1 - d x)^(3/2) (1 + d x)^(3/2))` in 2,917 characters, and + // this in 67,746. + if (!Integration.AnsweringTheQuestionAsked) + return null; + var (numerator, denominator) = Functions.SingleQuotient.Of(expr); + if (!TreeAnalyzer.TryGetPolynomial(numerator, x, out var top) || top.Keys.Any(k => k.Sign < 0 || !k.CanFitInInt32())) + return null; + Entity? quadratic = null; + var power = 0; + Entity? linear = null; + var linearPower = 0; + Entity constant = Number.Integer.One; + foreach (var factor in Mulf.LinearChildren(denominator)) + { + if (!factor.ContainsNode(x)) + { + constant = constant * factor; + continue; + } + var (@base, whole) = factor is Powf(var b0, Number.Integer { EInteger: var w }) && w.Sign > 0 && w.CanFitInInt32() + ? (b0, w.ToInt32Unchecked()) : (factor, 1); + if (TreeAnalyzer.TryGetPolyLinear(@base, x, out var slope, out _) && !VanishesIdentically(slope)) + { + if (linear is not null && linear != @base) + return null; + linear = @base; + linearPower += whole; + continue; + } + if (quadratic is null && whole >= 2 && TreeAnalyzer.TryGetPolyQuadratic(@base, x, out var leading, out _, out _) && !VanishesIdentically(leading)) + { + quadratic = @base; + power = whole; + continue; + } + return null; + } + if (quadratic is null || !quadratic.Vars.Any(symbol => symbol != x) + || !TreeAnalyzer.TryGetPolyQuadratic(quadratic, x, out var a2, out var a1, out var a0)) + return null; + var discriminant = (4 * a2 * a0 - a1 * a1).InnerSimplified; + if (VanishesIdentically(discriminant)) + return null; + // Real coefficients only, whose cases the table's 1/Q reads; and with no linear, not a + // binomial in x^2, which the binomial rules answer: `x^22/(a + b x^2)^10` was theirs in + // a second and ten levels of this ran past the budget. + if (Functions.PartialFractions.HoldsTheImaginaryUnit(quadratic) || linear is null && VanishesIdentically(a1)) + return null; + // A quadratic that factors over linears is the partial fractions' to take apart, and + // they do it cheaply: `(d + e x)^10/((d + e x) (a e + c d x))^4`, written out, went from + // two seconds past the budget here. + // It does where b^2 - 4 a c is the square of a polynomial in the symbols, which is + // what is asked: the factoring does not find `(d + e x) (a e + c d x)` written out. + var opposite = Functions.PartialFractions.Bare((a1 * a1 - 4 * a2 * a0).Expand().InnerSimplified); + var symbols = opposite.Vars.Distinct().OrderBy(v => v.Name, System.StringComparer.Ordinal).ToArray(); + var indices = new Dictionary(); + for (var i = 0; i < symbols.Length; i++) + indices[symbols[i]] = i; + if (symbols.Length <= MultivariatePolynomial.MaxVariables + && MultivariatePolynomial.TryParse(opposite, indices) is { } parsed && parsed.TrySquareRoot() is not null) + return null; + var named = Variable.CreateUnique(expr, "D_q"); + var atTheRoot = Variable.CreateUnique(expr + named, "Q_p"); + var theRoot = Variable.CreateUnique(expr + named + atTheRoot, "x_p"); + Entity d = named; + + static Entity[] Trim(Entity[] poly) + { + var length = poly.Length; + while (length > 0 && poly[length - 1] == Number.Integer.Zero) + length--; + return poly.Take(length).ToArray(); + } + static Entity[] Plus(Entity[] left, Entity[] right) + { + var sum = new Entity[System.Math.Max(left.Length, right.Length)]; + for (var i = 0; i < sum.Length; i++) + sum[i] = ((i < left.Length ? left[i] : Number.Integer.Zero) + (i < right.Length ? right[i] : Number.Integer.Zero)).InnerSimplified; + return sum; + } + static Entity[] Scaled(Entity[] poly, Entity by) => poly.Select(c => (c * by).InnerSimplified).ToArray(); + // Each coefficient in lowest terms over the symbols, so that the next step works on it + // as written small. + static Entity[] Reduced(Entity[] poly) => poly.Select(Functions.PartialFractions.InLowestTermsOverTheSymbols).ToArray(); + + var degree = top.Keys.Count == 0 ? 0 : top.Keys.Max()!.ToInt32Unchecked(); + var start = new Entity[degree + 1]; + for (var i = 0; i <= degree; i++) + start[i] = top.TryGetValue(EInteger.FromInt32(i), out var coefficient) ? coefficient : Number.Integer.Zero; + + // The integrand as a sum of polynomials over L^j Q^k, one per (j, k), merged as they + // arrive so that each stays of small degree. + var terms = new Dictionary<(int Linear, int Quadratic), Entity[]>(); + void Add(int j, int k, Entity[] poly) + => terms[(j, k)] = terms.TryGetValue((j, k), out var already) ? Plus(already, poly) : poly; + Add(linearPower, power, start); + Entity answer = Number.Integer.Zero; + Entity m = Number.Integer.One, p = theRoot; + if (linear is not null && !TreeAnalyzer.TryGetPolyLinear(linear, x, out m!, out _)) + return null; + // Q - Q(p) = (x - p) (c (x + p) + b). + var cofactor = new Entity[] { (a2 * p + a1).InnerSimplified, a2 }; + + while (terms.Count > 0) + { + var (j, k) = terms.Keys.OrderByDescending(key => key.Linear).ThenByDescending(key => key.Quadratic).First(); + var poly = Reduced(terms[(j, k)]); + terms.Remove((j, k)); + if (Trim(poly).Length == 0) + continue; + if (j > 0) + { + // N/(L^j Q^k) = N(p)/(Q(p)^k m^j (x - p)^j) + N1/(m L^(j-1) Q^k) + // - N(p) (c (x + p) + b) sum_i Q(p)^(k-1-i)/(Q(p)^k m L^(j-1) Q^(k-i)), + // with N1 = (N - N(p))/(x - p). + Entity atP = Number.Integer.Zero; + for (var i = poly.Length - 1; i >= 0; i--) + atP = (atP * p + poly[i]).InnerSimplified; + var c = Functions.PartialFractions.InLowestTermsOverTheSymbols(atP / MathS.Pow(atTheRoot, Number.Integer.Create(k)) / MathS.Pow(m, Number.Integer.Create(j))); + answer = answer + (j == 1 ? c * MathS.Ln(linear!) : c * MathS.Pow(linear!, Number.Integer.Create(1 - j)) * MathS.Pow(m, Number.Integer.Create(j - 1)) / (1 - j)); + var n1 = new Entity[System.Math.Max(poly.Length - 1, 0)]; + Entity carry = Number.Integer.Zero; + for (var i = poly.Length - 1; i >= 1; i--) + { + carry = (poly[i] + carry * p).InnerSimplified; + n1[i - 1] = carry; + } + Add(j - 1, k, Scaled(n1, 1 / m)); + for (var i = 0; i < k; i++) + Add(j - 1, k - i, Scaled(cofactor, (-atP * MathS.Pow(atTheRoot, Number.Integer.Create(-1 - i)) / m).InnerSimplified)); + continue; + } + if (k == 0) + { + for (var i = 0; i < poly.Length; i++) + answer = answer + poly[i] * MathS.Pow(x, Number.Integer.Create(i + 1)) / (i + 1); + continue; + } + // P/Q^k: P = S Q + t1 x + t0, and t1 x + t0 = alpha Q' + beta. + var remainder = poly.ToArray(); + var quotient = Enumerable.Repeat((Entity)Number.Integer.Zero, System.Math.Max(remainder.Length - 2, 0)).ToArray(); + for (var i = remainder.Length - 1; i >= 2; i--) + { + var factor = (remainder[i] / a2).InnerSimplified; + quotient[i - 2] = factor; + remainder[i] = Number.Integer.Zero; + remainder[i - 1] = (remainder[i - 1] - factor * a1).InnerSimplified; + remainder[i - 2] = (remainder[i - 2] - factor * a0).InnerSimplified; + } + var t1 = remainder.Length > 1 ? remainder[1] : Number.Integer.Zero; + var t0 = remainder.Length > 0 ? remainder[0] : Number.Integer.Zero; + var alpha = Functions.PartialFractions.InLowestTermsOverTheSymbols(t1 / (2 * a2)); + var beta = Functions.PartialFractions.InLowestTermsOverTheSymbols(t0 - alpha * a1); + Add(0, k - 1, quotient); + if (k == 1) + { + // Over the first power in closed form: alpha ln(Q), and beta times the table's + // integral of 1/Q. Asked of the integrator instead, a linear over Q with these + // coefficients took forty seconds. + if (IntegralPatterns.TryStandardIntegrals(1 / quadratic, x) is not { } overQ) + return null; + answer = answer + alpha * MathS.Ln(quadratic) + beta * overQ; + continue; + } + var lower = MathS.Pow(quadratic, Number.Integer.Create(k - 1)); + answer = answer - alpha / ((k - 1) * lower) + beta * (2 * a2 * x + a1) / ((k - 1) * d * lower); + Add(0, k - 1, new[] { (beta * 2 * a2 * (2 * k - 3) / ((k - 1) * d)).InnerSimplified }); + } + Entity withoutTheLinear = Number.Integer.Zero; + Entity scale = Number.Integer.One; + answer = ((withoutTheLinear + scale * answer) / constant).Substitute(named, discriminant); + if (linear is not null && TreeAnalyzer.TryGetPolyLinear(linear, x, out var mm, out var nn)) + { + var root = (-nn / mm).InnerSimplified; + // Q(p) as one quotient in lowest terms, (c d^2 - b d e + a e^2)/e^2, rather than + // Q written at -d/e wherever it stood. + answer = answer.Substitute(atTheRoot, Functions.PartialFractions.InLowestTermsOverTheSymbols(a2 * root * root + a1 * root + a0)).Substitute(theRoot, root); + } + answer = answer.InnerSimplified; + return Functions.PartialFractions.DerivativeHoldsAtSampledPoints(answer, expr, x) ? answer : null; + } + /// /// A power of a linear beside the derivative of a quadratic, over a whole power of that /// quadratic, by parts once: Q' L^m/Q^n with u = L^m and diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs b/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs index 49f6fbd56..a02c2360c 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs @@ -669,6 +669,9 @@ private static Entity WithItsPiecewiseSumsCombined(Entity expr) // A power of a linear beside the derivative of a quadratic over a whole power of it, by // parts once, before anything takes the symbolic quadratic apart. if ((answer = IndefiniteIntegralSolver.SolveAPowerOfALinearBesideTheDerivativeOfAQuadraticOverItsPower(expr, x, integrateByParts)) is { }) return answer; + // A rational function over a whole power of a symbolic quadratic, the power lowered one + // at a time, before the partial fractions take the quadratic apart over its roots. + if ((answer = IndefiniteIntegralSolver.SolveARationalFunctionOverAPowerOfASymbolicQuadratic(expr, x, integrateByParts)) is { }) return answer; // A root or a power of a product beside a power of one of its factors, the product // split: `(a + b x^2)^(5/2)/sqrt(a^2 - b^2 x^4)` is `K (a + b x^2)^2/sqrt(a - b x^2)`. // Early, since what reads the product whole is slower and longer: Rubi's 1.2.1.2:2218, diff --git a/Sources/Tests/UnitTests/Calculus/ARationalFunctionOverAPowerOfASymbolicQuadraticIntegralTest.cs b/Sources/Tests/UnitTests/Calculus/ARationalFunctionOverAPowerOfASymbolicQuadraticIntegralTest.cs new file mode 100644 index 000000000..6404579ea --- /dev/null +++ b/Sources/Tests/UnitTests/Calculus/ARationalFunctionOverAPowerOfASymbolicQuadraticIntegralTest.cs @@ -0,0 +1,60 @@ +// +// 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 rational function over a whole power of a symbolic quadratic, beside one linear: the + /// linear taken off at its root and the power of the quadratic lowered one at a time. Rubi's + /// 1.2.1.2 1/((d + e x) (a + b x + c x^2)^4), written here with k for e, + /// which is Euler's number. + /// #718 + /// + [Trait("Area", "Calculus")] + public sealed class ARationalFunctionOverAPowerOfASymbolicQuadraticIntegralTest + { + private static void DifferentiatesBack(string integrand, params (string, double)[] pins) + { + var integral = integrand.ToEntity().Integrate("x"); + Assert.DoesNotContain("integral(", integral.Stringize()); + + // Pinned before it is differentiated: the answer is a piecewise on the sign of the + // discriminant, read where it is a number. + var pinned = integral.Substitute("C", 0); + Entity original = integrand.ToEntity(); + foreach (var (name, value) in pins) + { + pinned = pinned.Substitute(name, value); + original = original.Substitute(name, value); + } + var derivative = pinned.Differentiate("x"); + foreach (var at in new[] { -2.3, -0.7, 0.4, 1.9 }) + { + var got = derivative.Substitute("x", at).EvalNumerical(); + var want = original.Substitute("x", at).EvalNumerical(); + var difference = Math.Abs((double)(got - want).RealPart) + Math.Abs((double)(got - want).ImaginaryPart); + var scale = Math.Max(1.0, Math.Abs((double)want.RealPart) + Math.Abs((double)want.ImaginaryPart)); + Assert.True(difference / scale < 1e-8, + $"d/dx of the antiderivative of {integrand} is {got} at x = {at}, where the integrand is {want}"); + } + } + + [Theory] + [InlineData("x^8/(a+b*x+c*x^2)^4")] + [InlineData("1/((d+k*x)*(a+b*x+c*x^2)^4)")] + [InlineData("1/((d+k*x)^2*(a+c*x^2)^4)")] + public void IsAnsweredAndDifferentiatesBack(string integrand) + { + DifferentiatesBack(integrand, ("a", 2), ("b", 1.5), ("c", 1.25), ("d", 1.25), ("k", 2)); + DifferentiatesBack(integrand, ("a", -1), ("b", 2), ("c", 2), ("d", -0.5), ("k", 1.5)); + } + } +}