diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index 91cf2dca3..05ed4545d 100644 --- a/BREAKING-CHANGES.md +++ b/BREAKING-CHANGES.md @@ -253,6 +253,30 @@ then cancels before the division of the improper fraction | `"(a + b*x)^6/(a*c + (b*c + a*d)*x + b*d*x^2)^2".ToEntity().Integrate("x")` | past a minute | 3,592 characters | | `"csc(e + f*x)^3/(a + b*tan(e + f*x)^2)^2".ToEntity().Integrate("x")` | `integral(...)` | 7,992 characters | +### What a reduction leaves over its block is integrated in named coefficients + +**Answers where there were none, and sooner.** `(c + d x^3 + g x^6 + f x^9)/(x^4 (a + b x^3)^3)` took 42 s on the +unreleased master, and `x^5`, `x^7` and `x^8` in place of `x^4` took 35 to 50 s, Rubi's 1.1.3.8. Four polynomials +over `(a + b x^2 + c x^4)^2` and `^3` in 1.2.2.5 ran past two minutes, and +`(d + f x + g x^2 + h x^3 + k x^4)/(a + b x^2 + c x^4)^3` took 31 s. Two reductions -- the split at a power of x beside +a block, and the one that takes a polynomial over a power of a trinomial in `x^n` down to its first power -- +integrated what was left over the block with the coefficients they had computed, such as +`(a^2 g - 3 a b d + 6 b^2 c)/a^2`, and the substitution search spent up to a minute simplifying quotients in them +before another rule answered. Each such coefficient is a fresh symbol while the rest is integrated now, and is put +back in the answer. The answers over a power of the trinomial are large, 22,000 to 98,000 characters, and right. +One answer grew: `1/(x^7 (a + b x^3)(c + d x^3))` took 52 s on the unreleased master for 14,565 characters, and +takes 6 s for 265,499 now ([#718](https://github.com/asc-community/AngouriMath/issues/718)). + +| Input | Was (2.5.0) | Now | +|---|---|---| +| `"(c + d*x^3 + g*x^6 + f*x^9)/(x^4*(a + b*x^3)^3)".ToEntity().Integrate("x")` | `integral(...)` | 547 characters | +| `"(c + d*x^3 + g*x^6 + f*x^9)/(x^8*(a + b*x^3)^3)".ToEntity().Integrate("x")` | `integral(...)` | 1,859 characters | +| `"1/(x^7*(a + b*x^3)*(c + d*x^3))".ToEntity().Integrate("x")` | `integral(...)` | 265,499 characters | +| `"(d + f*x + g*x^2 + h*x^3 + k*x^4)/(a + b*x^2 + c*x^4)^3".ToEntity().Integrate("x")` | `integral(...)` | 31,531 characters | +| `"(d + f*x + g*x^2 + h*x^3 + k*x^4 + m*x^5 + n*x^6 + p*x^7)/(a + b*x^2 + c*x^4)^2".ToEntity().Integrate("x")` | `integral(...)` | 21,719 characters | +| `"(d + f*x + g*x^2 + h*x^3 + k*x^4 + m*x^5 + n*x^8 + p*x^11)/(a + b*x^2 + c*x^4)^3".ToEntity().Integrate("x")` | `integral(...)` | 77,389 characters | +| `"csc(c + d*x)^2/(a - b*sin(c + d*x)^4)^3".ToEntity().Integrate("x")` | `integral(...)` | 157,787 characters | + ### One linear twice is not two linear powers **Answers that had no value.** `(a c + b c x)^(-3 - 2 p) (f + g x) (a^2 + 2 a b x + b^2 x^2)^p` was answered diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs index bb53654c7..881634610 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs @@ -1361,7 +1361,7 @@ Entity Coefficient(Dictionary polynomial, int power) } // R = (N - P B)/x^k, whose coefficient of x^(m - k) is that of x^m in N - P B. - Entity remainder = Number.Integer.Zero; + var remainder = new Dictionary(); for (var m = k; m <= System.Math.Max(degreeAbove, degreeBelow + k - 1); m++) { var coefficient = Coefficient(above, m); @@ -1369,7 +1369,7 @@ Entity Coefficient(Dictionary polynomial, int power) coefficient -= series[j] * Coefficient(below, m - j); coefficient = Functions.PartialFractions.InLowestTermsOverTheSymbols(coefficient); if (!VanishesIdentically(coefficient)) - remainder += m == k ? coefficient : coefficient * MathS.Pow(x, m - k); + remainder[m - k] = coefficient; } // The part over x^k term by term, and the rest over the block as it is written. @@ -1379,7 +1379,7 @@ Entity Coefficient(Dictionary polynomial, int power) overX += j == k - 1 ? series[j] * IntegralPatterns.AntiderivativeLog(x) : series[j] * MathS.Pow(x, j - k + 1) / (j - k + 1); - if (remainder == Number.Integer.Zero) + if (remainder.Count == 0) return overX; return OverTheRestOfTheBlock(remainder, block, x, integrateByParts) is { } overTheBlock ? overX + overTheBlock @@ -1398,10 +1398,69 @@ Entity Coefficient(Dictionary polynomial, int power) /// https://github.com/asc-community/AngouriMath/issues/718 /// private static Entity? OverTheRestOfTheBlock(Entity remainder, Entity block, Entity.Variable x, bool integrateByParts) + => AtAnotherLinearOfTheBlock(remainder, block, x, integrateByParts) + ?? Integration.ComputeIndefiniteIntegral(remainder / block, x, integrateByParts); + + /// + /// for a rest + /// given by its coefficients, integrated a level down in named coefficients. + /// + private static Entity? OverTheRestOfTheBlock(Dictionary rest, Entity block, Entity.Variable x, bool integrateByParts) + { + Entity remainder = Number.Integer.Zero; + foreach (var pair in rest.OrderBy(pair => pair.Key)) + remainder += pair.Key == 0 ? pair.Value : pair.Value * (pair.Key == 1 ? x : MathS.Pow(x, pair.Key)); + return AtAnotherLinearOfTheBlock(remainder, block, x, integrateByParts) + ?? IntegratedInNamedCoefficients(rest, block, x, integrateByParts); + } + + /// + /// R/B split at a linear the block still holds, where that answer differentiates back. + /// + private static Entity? AtAnotherLinearOfTheBlock(Entity remainder, Entity block, Entity.Variable x, bool integrateByParts) => IntegrateOverAPowerOfALinearBesideABlock(remainder, block, x, integrateByParts) is { } atAnotherRoot && Functions.PartialFractions.DerivativeHoldsAtSampledPoints(atAnotherRoot, remainder / block, x) ? atAnotherRoot - : Integration.ComputeIndefiniteIntegral(remainder / block, x, integrateByParts); + : null; + + /// + /// The integral of R/B, for a polynomial R in x given by its coefficients, with each coefficient of + /// R that is more than a symbol or a number a fresh symbol while it is integrated, + /// and put back in the answer. + /// + /// + /// What a reduction leaves over its block has coefficients the reduction computed, such as + /// (a^2 g - 3 a b d + 6 b^2 c)/a^2, and the substitution search simplifies each + /// quotient it tries in them. Over (a/b + x^3)^3 that took 40 s to decline where three + /// symbols take it 45 ms, and over (a + b x^2 + c x^4) past a minute. + /// https://github.com/asc-community/AngouriMath/issues/718 + /// + private static Entity? IntegratedInNamedCoefficients(Dictionary rest, Entity block, Entity.Variable x, bool integrateByParts) + { + var named = new List<(Variable Name, Entity Value)>(); + Entity everything = block; + foreach (var value in rest.Values) + everything += value; + Entity inNames = Number.Integer.Zero; + foreach (var pair in rest.OrderBy(pair => pair.Key)) + { + var power = pair.Key; + var coefficient = pair.Value; + if (coefficient is not (Variable or Number)) + { + var name = Variable.CreateUnique(everything, "k_r"); + everything += name; + named.Add((name, coefficient)); + coefficient = name; + } + inNames += power == 0 ? coefficient : coefficient * (power == 1 ? x : MathS.Pow(x, power)); + } + if (Integration.ComputeIndefiniteIntegral(inNames / block, x, integrateByParts) is not { } answer) + return null; + foreach (var (name, value) in named) + answer = answer.Substitute(name, value); + return answer; + } /// /// N/(L^k B) for a linear L = g + h x other than x, beside a block B @@ -1897,13 +1956,13 @@ Entity InU(Dictionary polynomial) Entity answer = Number.Integer.Zero; foreach (var pair in overPowers.OrderByDescending(pair => pair.Key)) answer = answer + pair.Value / (pair.Key == 1 ? trinomial : MathS.Pow(trinomial, pair.Key)); - Entity rest = Number.Integer.Zero; - foreach (var pair in atThisPower.OrderBy(pair => pair.Key)) + var rest = new Dictionary(); + foreach (var pair in atThisPower) if (!VanishesIdentically(pair.Value)) - rest = rest + pair.Value * (pair.Key == 0 ? Number.Integer.One : pair.Key == 1 ? x : MathS.Pow(x, pair.Key)); - if (rest != Number.Integer.Zero) + rest[pair.Key] = pair.Value; + if (rest.Count > 0) { - if (Integration.ComputeIndefiniteIntegral(rest / trinomial, x, integrateByParts) is not { } overTheTrinomial) + if (IntegratedInNamedCoefficients(rest, trinomial, x, integrateByParts) is not { } overTheTrinomial) return null; answer = answer + overTheTrinomial; } diff --git a/Sources/Tests/UnitTests/Calculus/TheRestOverABlockInNamedCoefficientsIntegralTest.cs b/Sources/Tests/UnitTests/Calculus/TheRestOverABlockInNamedCoefficientsIntegralTest.cs new file mode 100644 index 000000000..93136e882 --- /dev/null +++ b/Sources/Tests/UnitTests/Calculus/TheRestOverABlockInNamedCoefficientsIntegralTest.cs @@ -0,0 +1,45 @@ +// +// 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 +{ + /// + /// What a reduction leaves over its block is integrated with its compound coefficients named: + /// with them as written, the substitution search spent up to a minute declining it, and these + /// did not finish. Rubi's 1.1.3.8 and 1.2.2.5. + /// #718 + /// + [Trait("Area", "Calculus")] + public sealed class TheRestOverABlockInNamedCoefficientsIntegralTest + { + [Theory] + [InlineData("(c + d*x^3 + g*x^6 + f*x^9)/(x^4*(a + b*x^3)^3)")] + [InlineData("(c + d*x^3 + g*x^6 + f*x^9)/(x^8*(a + b*x^3)^3)")] + [InlineData("(d*x + g*x^3)/(a + b*x^2 + c*x^4)^3")] + public void DifferentiatesBack(string integrand) + { + var integral = integrand.ToEntity().Integrate("x"); + Assert.DoesNotContain("integral(", integral.Stringize()); + Entity Pinned(Entity e) => e.Substitute("a", 2).Substitute("b", 3).Substitute("c", 5) + .Substitute("d", -1).Substitute("g", 7).Substitute("f", 1); + var derivative = Pinned(integral.Substitute("C", 0)).Differentiate("x"); + var original = Pinned(integrand.ToEntity()); + foreach (var at in new[] { -1.7, -0.6, 0.45, 1.3 }) + { + 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)), + $"d/dx of the antiderivative of {integrand} is {got} at x = {at}, where the integrand is {want}"); + } + } + } +}