diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index eb039db20..4e194e312 100644 --- a/BREAKING-CHANGES.md +++ b/BREAKING-CHANGES.md @@ -3052,6 +3052,25 @@ polynomial exponent in `u` now, and the question asked in `u`. Rubi's 2.3 | `"F^(a + b/(c + d*x))/(c + d*x)".ToEntity().Integrate("x")` | `integral(...)` | an exponential integral | | `"F^(a + b/(c + d*x)^3)*(c + d*x)^2".ToEntity().Integrate("x")` | `integral(...)` | an exponential integral of the cube, and the exponential | +### An exponential of the reciprocal of a linear over another linear, and its elementary moments, are integrated under `u = 1/(c + d x)` first + +**Answers where there were none, and sooner.** `F^(e + f (a + b x)/(c + d x))/(g + h x)^4` was declined by 2.5.0, +and the unreleased master answered it in ten seconds and 12 GB, the third power in six; `F^(a + b/(c + d x)^3)/(c + d x)^13`, +`^16` and `^19` ran past ten seconds there, Rubi's 2.3. The substitution `u = 1/(c + d x)` that answers each in +hundredths of a second was made only after every other rule had declined. It is made first now where another linear +stands below the bar, and where the integral is elementary: `F^(a + b/(c + d x)^k)/(c + d x)^n` with `n - 1` a positive +multiple of `k >= 2`, a polynomial times an exponential under `w = u^k`. A power of `c + d x` below the bar arrives +written monic, `c/d + x` beside the `c + d x` of the exponent, and is read as a constant multiple of it by value. Two +answers 2.5.0 gave are respelled, `F^(a + b/(c + d x)^2)/(c + d x)^3` and `F^(a + b/(c + d x)^3)/(c + d x)^4`, from 58 +characters to 74 ([#233](https://github.com/asc-community/AngouriMath/issues/233)). + +| Input | Was (2.5.0) | Now | +|---|---|---| +| `"F^(m + f*(a + b*x)/(c + d*x))/(g + h*x)^4".ToEntity().Integrate("x")` | `integral(...)` | 1,800 characters | +| `"F^(m + f*(a + b*x)/(c + d*x))/(g + h*x)^3".ToEntity().Integrate("x")` | `integral(...)` | 969 characters | +| `"F^(a + b/(c + d*x)^3)/(c + d*x)^13".ToEntity().Integrate("x")` | `integral(...)` | 262 characters | +| `"F^(a + b/(c + d*x)^2)/(c + d*x)^3".ToEntity().Integrate("x")` | 58 characters | 74 characters | + ### A power of a monomial in a symbolic power of `x` is integrated as the power of `x` it is **Answers where there were none.** `(c x^n)^b` is written `c^b x^(n b)`, and the exponent `n b` was diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs index 13fcdeb3d..bd20612e7 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs @@ -26055,6 +26055,113 @@ private static bool TryReadAPowerTimesARadicalQuadratic(Entity expr, Entity.Vari return (factor * answerInU.Substitute(u, 1 / linear)).InnerSimplified; } + /// + /// An exponential of a constant plus a multiple of the reciprocal of a linear over a whole + /// power of another linear, F^(a + b/(c + d x))/(g + h x)^n, under + /// 's substitution u = 1/(c + d x), + /// asked before the rules that do not read it: in u it is + /// -d^(n - 1) u^(n - 2) F^(a + b u)/((d g - c h) u + h)^n, the exponential integral of a + /// shifted argument, which and + /// answer. The exponent may be written + /// e + f (a + b x)/(c + d x). Rubi's 2.3. + /// https://github.com/asc-community/AngouriMath/issues/233 + /// + /// + /// In its place at the end of the dispatch the substitution was reached only after every rule + /// above it had declined, and F^(e + f (a + b x)/(c + d x))/(g + h x)^n took 2 to 10 + /// seconds and 0.4 to 12 GB for n = 1 to 4, past the budget for 3 and 4, for the answer + /// it writes in a twentieth of a second. Only for an exponent linear in u and another + /// linear below the bar whose slope in u, d g - c h, is not zero as a value. A + /// multiple of c + d x is not another linear, though a power of c + d x alone + /// arrives as one, written monic: its elementary cases are the ansatz's, which writes them + /// shorter, and 's. + /// Another linear above the bar only, and a Gaussian in u, keep the rule's place. + /// + internal static Entity? SolveAnExponentialInTheReciprocalOfALinearOverAnotherLinear(Entity expr, Entity.Variable x, bool integrateByParts) + { + if (ReadAnExponentialInTheReciprocalOfALinear(expr, x) is not { ExponentInU: var exponentInU, U: var u, Beside: (_, var below, _) } + || !below.ContainsNode(u) + || !TreeAnalyzer.TryGetPolyLinear(exponentInU, u, out var rate, out _) || rate is null || VanishesIdentically(rate)) + return null; + // Each linear below the bar a linear in u, its slope not zero as a value. + foreach (var factor in Mulf.LinearChildren(below)) + { + var linearInU = factor is Powf(var raised, Number.Integer) ? raised : factor; + if (!TreeAnalyzer.TryGetPolyLinear(linearInU, u, out var slopeInU, out _) || slopeInU is null || VanishesIdentically(slopeInU)) + return null; + } + return SolveAnExponentialInTheReciprocalOfALinear(expr, x, integrateByParts); + } + + /// + /// An exponential of a constant plus a multiple of a power of the reciprocal of a linear, + /// over a power of the linear, where the integral is elementary: + /// F^(a + b/(c + d x)^k)/(c + d x)^n with k >= 2 and n - 1 a positive + /// multiple of k. Under 's + /// u = 1/(c + d x) it is -(1/d) u^(n - 2) F^(a + b u^k), and under w = u^k a + /// polynomial in w times F^(a + b w), which the question in u is answered by. + /// Rubi's 2.3. https://github.com/asc-community/AngouriMath/issues/233 + /// + /// + /// In the substitution's place at the end of the dispatch these were reached after every rule + /// above it had declined: F^(a + b/(c + d x)^2)/(c + d x)^13 took 6 seconds and 6 GB, + /// F^(a + b/(c + d x)^3)/(c + d x)^13 and higher ran past the budget, and the questions + /// in u take a hundredth of a second. The power below the bar arrives written monic, + /// c/d + x beside the c + d x of the exponent, and is read as a constant multiple + /// of the linear by value. Only where the integral is elementary, and only for a linear other + /// than the variable: asked here for every power, the exponential integral and the error + /// function came out longer or slower than where they are answered now, a linear exponent in + /// u, which the ansatz answers in the base, three times longer, and the moments of + /// f^(a + b/x^2), answered in a hundredth of a second already, respelled. + /// + internal static Entity? SolveAnElementaryMomentOfAnExponentialInAPowerOfTheReciprocal(Entity expr, Entity.Variable x, bool integrateByParts) + { + if (ReadAnExponentialInTheReciprocalOfALinear(expr, x) is not var (constant, @base, exponentInU, u, linear, slope, power, beside) + || linear == x || power.Sign >= 0 + || !TreeAnalyzer.TryGetPolynomial(exponentInU, u, out var monomials) || monomials.Count > 2) + return null; + // A + B u^k: one power past the constant, its coefficient not zero as a value. + var raised = monomials.Where(pair => pair.Key.Sign != 0).ToList(); + if (raised.Count != 1 || monomials.Keys.Any(degree => degree.Sign < 0) || VanishesIdentically(raised[0].Value)) + return null; + var k = raised[0].Key; + // u^(n - 2) du with n - 1 a positive multiple of k is a whole power of w = u^k times dw. + var belowMinusOne = power.Negate().Subtract(EInteger.One); + var powerInU = belowMinusOne.Subtract(EInteger.One); + if (k.CompareTo(EInteger.FromInt32(2)) < 0 || k.CompareTo(EInteger.FromInt32(8)) > 0 + || belowMinusOne.Sign <= 0 || !belowMinusOne.Remainder(k).IsZero || powerInU.CompareTo(EInteger.FromInt32(32)) > 0) + return null; + var factor = -constant / slope; + // Another linear only as a multiple of the linear, which is how the power below the bar + // arrives, written monic: `c/d + x` beside the `c + d x` of the exponent is + // `(d^0 c - c) u + 1` in u, a constant whose slope is zero as a value and not as written. + // Its value at u = 0, and d to the power the reader kept apart. Bare, since the `d^0` in + // it is 1 provided d is not 0, which the substitution already assumes. + if (beside is var (above, below, slopePower)) + { + if (!IsAConstantInU(above) || !IsAConstantInU(below)) + return null; + var ofTheSlope = slopePower - EInteger.One; + var multiple = Functions.PartialFractions.Bare((above.Substitute(u, Number.Integer.Zero) / below.Substitute(u, Number.Integer.Zero)).InnerSimplified); + factor = ofTheSlope.IsZero ? -constant + : ofTheSlope.Sign > 0 ? -constant * MathS.Pow(slope, Number.Integer.Create(ofTheSlope)) + : -constant / MathS.Pow(slope, Number.Integer.Create(-ofTheSlope)); + if (multiple != Number.Integer.One) + factor *= multiple; + } + var exponential = MathS.Pow(@base, exponentInU); + var question = powerInU.IsZero ? exponential : MathS.Pow(u, Number.Integer.Create(powerInU)) * exponential; + if (Integration.ComputeAsAQuestionOfItsOwn(question, u, integrateByParts) is not { } answerInU + || answerInU.Nodes.Any(node => node is Integralf || node == MathS.NaN)) + return null; + return (factor * answerInU.Substitute(u, 1 / linear)).InnerSimplified; + + // A product of whole powers of linears in u, each of slope zero as a value. + bool IsAConstantInU(Entity product) => !product.ContainsNode(u) || Mulf.LinearChildren(product).All(part => + !part.ContainsNode(u) || TreeAnalyzer.TryGetPolyLinear(part is Powf(var raisedLinear, Number.Integer) ? raisedLinear : part, u, out var slopeInU, out _) + && (slopeInU is null || VanishesIdentically(slopeInU))); + } + /// /// read as a constant times F^E times a whole power of a /// linear L, with every x in E inside a reciprocal power of L: diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs b/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs index ac3141bc9..4cc23dd9f 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs @@ -841,11 +841,17 @@ private static Entity WithItsPiecewiseSumsCombined(Entity expr) if ((answer = IndefiniteIntegralSolver.SolveByWritingAPowerOfAnExponentialAsAMultipleOfOne(expr, x, integrateByParts)) is { }) return answer; // An exponential of a quadratic in 1/L beside a power of L, onto the Gaussian under u = 1/L. if ((answer = IndefiniteIntegralSolver.SolveAGaussianInAReciprocal(expr, x)) is { }) return answer; + // And of the reciprocal of a linear over another linear, the same substitution's, ahead of + // its place below: `F^(e + f (a + b x)/(c + d x))/(g + h x)^4` ran past the budget there. + if ((answer = IndefiniteIntegralSolver.SolveAnExponentialInTheReciprocalOfALinearOverAnotherLinear(expr, x, integrateByParts)) is { }) return answer; // And a trigonometric function of the reciprocal of a linear beside a power of it, under u = 1/L. if ((answer = IndefiniteIntegralSolver.SolveAFunctionOfTheReciprocalOfALinear(expr, x, integrateByParts)) is { }) return answer; // And one of `a + b/x` beside any rational function, under `x = 1/u`: // `sin(c + d/x)/(e + f x)` is `-sin(c + d u)/(u (e u + f)) du`. if ((answer = IndefiniteIntegralSolver.SolveByTheReciprocalOfAFunctionOfTheReciprocal(expr, x, integrateByParts)) is { }) return answer; + // And a power of the linear alone beside an exponential of a power of its reciprocal, where + // that is elementary: `F^(a + b/(c + d x)^3)/(c + d x)^13` ran past the budget below. + if ((answer = IndefiniteIntegralSolver.SolveAnElementaryMomentOfAnExponentialInAPowerOfTheReciprocal(expr, x, integrateByParts)) is { }) return answer; // An exponential of a multiple of a logarithm is a power of the argument, which is // how every inverse hyperbolic function under an exponential arrives. if ((answer = IndefiniteIntegralSolver.SolveByFoldingAnExponentialOfALogarithm(expr, x, integrateByParts)) is { }) return answer; diff --git a/Sources/Tests/UnitTests/Calculus/AnExponentialInTheReciprocalOfALinearOverAnotherLinearIntegralTest.cs b/Sources/Tests/UnitTests/Calculus/AnExponentialInTheReciprocalOfALinearOverAnotherLinearIntegralTest.cs new file mode 100644 index 000000000..c94bda036 --- /dev/null +++ b/Sources/Tests/UnitTests/Calculus/AnExponentialInTheReciprocalOfALinearOverAnotherLinearIntegralTest.cs @@ -0,0 +1,94 @@ +// +// 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 +{ + /// + /// An exponential of a constant plus a multiple of 1/L, L = c + d x, over a whole + /// power of another linear M = g + h x, and of a power of 1/L over a power of + /// L where that is elementary. Under u = 1/L the first is an exponential of a + /// linear over a power of a linear in u, the exponential integral of a shifted argument, + /// and the second u^(n - 2) F^(a + b u^k) with n - 1 a multiple of k. The + /// substitution was reached only after every other rule had declined: + /// F^(e + f (a + b x)/(c + d x))/(g + h x)^n took 4 to 12 GB for n = 2 to 4, and + /// F^(a + b/(c + d x)^3)/(c + d x)^13 ran past the budget. Rubi's 2.3. + /// #233 + /// + [Trait("Area", "Calculus")] + public sealed class AnExponentialInTheReciprocalOfALinearOverAnotherLinearIntegralTest + { + [Theory] + [InlineData("F^(m + f*(a + b*x)/(c + d*x))/(g + h*x)", 1.1)] + [InlineData("F^(m + f*(a + b*x)/(c + d*x))/(g + h*x)^2", 1.1)] + [InlineData("F^(m + f*(a + b*x)/(c + d*x))/(g + h*x)^3", 1.1)] + [InlineData("F^(m + f*(a + b*x)/(c + d*x))/(g + h*x)^4", 1.1)] + [InlineData("F^(m + f*(a + b*x)/(c + d*x))/(g + h*x)^3", -0.9)] + [InlineData("F^(a + b/(c + d*x))/(k + f*x)^3", 1.1)] + [InlineData("e^(k/(c + d*x))/(a + b*x)^3", -0.9)] + [InlineData("f^(c/(a + b*x))/x^3", 1.1)] + [InlineData("F^(a + b/(c + d*x))*(p + q*x)/(g + h*x)^2", 1.1)] + [InlineData("F^(a + b/(c + d*x))/((g + h*x)*(p + q*x))", 1.1)] + [InlineData("e^(1/(1 + x))/(2 + x)^2", 1.1)] + [InlineData("F^(a + b/(c + d*x)^3)/(c + d*x)^13", 1.1)] + [InlineData("F^(a + b/(c + d*x)^3)/(c + d*x)^4", 1.1)] + [InlineData("F^(a + b/(c + d*x)^2)/(c + d*x)^7", 1.1)] + [InlineData("F^(a + b/(c + d*x)^2)/(c + d*x)^3", 1.1)] + [InlineData("e^(k/(1 + 2*x)^2)/(1 + 2*x)^5", -0.9)] + public void DifferentiatesBack(string integrand, double k) + { + var integral = integrand.ToEntity().Integrate("x"); + var text = integral.Stringize(); + Assert.DoesNotContain("integral(", text); + Assert.True(text.Length < 5000, $"{text.Length} characters of answer for {integrand}"); + Entity Pinned(Entity e) => e.Substitute("a", 1.3).Substitute("b", 0.7).Substitute("c", 0.4).Substitute("d", 0.3) + .Substitute("f", 2.3).Substitute("g", 0.6).Substitute("h", 1.9).Substitute("m", 0.5).Substitute("F", 1.7) + .Substitute("p", 0.9).Substitute("q", 0.45).Substitute("k", k); + var derivative = Pinned(integral.Substitute("C", 0)).Differentiate("x"); + var original = Pinned(integrand.ToEntity()); + var compared = 0; + foreach (var at in new[] { -0.9, -0.4, 0.3, 0.8, 1.3, 2.9 }) + { + var want = original.Substitute("x", at).EvalNumerical(); + var got = derivative.Substitute("x", at).EvalNumerical(); + if (want.IsNaN) + continue; + compared++; + 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}"); + } + Assert.True(compared >= 5, $"only {compared} points could be compared for {integrand}"); + } + + /// + /// Answered by the substitution before the rules that decline it: in its old place the + /// declines cost 3.9, 7.9 and 11.7 GB for the second, third and fourth power, and 6.2 and + /// 5.6 GB for the moments. Spelled in + /// symbols no other test here uses, since an integral already answered on the thread is + /// read from the cache and allocates nothing. + /// + [Theory] + [InlineData("G^(n + s*(r + t*x)/(u + v*x))/(w + y*x)^2")] + [InlineData("G^(n + s*(r + t*x)/(u + v*x))/(w + y*x)^3")] + [InlineData("G^(n + s*(r + t*x)/(u + v*x))/(w + y*x)^4")] + [InlineData("G^(n + s/(u + v*x)^2)/(u + v*x)^13")] + [InlineData("G^(n + s/(u + v*x)^3)/(u + v*x)^10")] + public void IsAnsweredBeforeTheSearch(string integrand) + { + var expression = integrand.ToEntity(); + var before = GC.GetAllocatedBytesForCurrentThread(); + var text = expression.Integrate("x").Stringize(); + var allocated = GC.GetAllocatedBytesForCurrentThread() - before; + Assert.DoesNotContain("integral(", text); + Assert.True(allocated < 1_000_000_000L, $"integrating {integrand} allocated {allocated} bytes"); + } + } +}