diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index ab2fe2684..1bbb3b9a1 100644 --- a/BREAKING-CHANGES.md +++ b/BREAKING-CHANGES.md @@ -447,6 +447,23 @@ case is most of its length ([#718](https://github.com/asc-community/AngouriMath/ | `"x^4/((a + b*x^2 + c*x^4)*sqrt(d + k*x^2))".ToEntity().Integrate("x")` | `integral(...)` | 28,551 characters | | `"1/(x^2*(a + b*x^2 + c*x^4)*sqrt(d + k*x^2))".ToEntity().Integrate("x")` | `integral(...)` | 10,194 characters | +### A power of a linear over a power of a quadratic is reduced a power at a time + +**Answers where there were none.** `sqrt(d + k x)/(a + b x + c x^2)^3` took thirty seconds for 26,437 characters on the +unreleased master, and `(d + k x)^(7/2)/(a + b x + c x^2)^2` ran past a minute. A power of a linear that is not whole, +times a linear, over a whole power of a quadratic is now reduced one power of the quadratic at a time: what is left each +time is the same shape with the quadratic's power one lower, and over the first power the integrator answers it +([#718](https://github.com/asc-community/AngouriMath/issues/718)). + +| Input | Was (2.5.0) | Now | +|---|---|---| +| `"(d + k*x)^(7/2)/(a + b*x + c*x^2)^2".ToEntity().Integrate("x")` | `integral(...)` | 4,691 characters | +| `"1/((d + k*x)^(3/2)*(a + b*x + c*x^2)^2)".ToEntity().Integrate("x")` | `integral(...)` | 2,789 characters | +| `"(d + k*x)^(7/2)/(a + b*x + c*x^2)^3".ToEntity().Integrate("x")` | `integral(...)` | 7,149 characters | +| `"sqrt(d + k*x)/(a + b*x + c*x^2)^3".ToEntity().Integrate("x")` | `integral(...)` | 4,819 characters | +| `"(A + B*x)/(x^(3/2)*(a + b*x + c*x^2)^2)".ToEntity().Integrate("x")` | `integral(...)` | 1,169 characters | +| `"(b + 2*c*x)*(d + k*x)^(3/2)/(a + b*x + c*x^2)^3".ToEntity().Integrate("x")` | `integral(...)` | 2,505 characters | + ### A power of the cosine beside a secant binomial is integrated in the cosine **Answers where there were none.** `cos(c + d x)^2/(a + b sec(c + d x))` was declined, and so were the cosine's powers diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs index 3f60a9126..ae9182ffe 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs @@ -19544,6 +19544,137 @@ void Add(int j, int k, Entity[] poly) return Functions.PartialFractions.DerivativeHoldsAtSampledPoints(answer, expr, x) ? answer : null; } + /// + /// A power of a linear that is not whole, times a linear, over a whole power of a quadratic, + /// with the power of the quadratic lowered one at a time: + /// integral L^m N Q^p = (alpha + beta x) L^(m + 1) Q^(p + 1) + integral L^m W Q^(p + 1) + /// for p <= -2, where alpha and beta make + /// N - (alpha + beta x) (p + 1) L Q' a multiple of Q and + /// W = (N - (alpha + beta x) (p + 1) L Q')/Q - beta L - (alpha + beta x) (m + 1) e is + /// linear again; at p = -1 the integrator answers what is left. + /// + /// + /// The derivative of the first term is L^m Q^p times + /// beta L Q + (alpha + beta x) (m + 1) e Q + (alpha + beta x) L (p + 1) Q', which with + /// L^m Q^(p + 1) W added is L^m Q^p N. (alpha + beta x) L Q' modulo + /// Q is (alpha + beta x) (r1 x + r0) with r1 = 2 c d - b e and + /// r0 = b d - 2 a e, for L = d + e x and Q = a + b x + c x^2, and + /// matching it to N is two equations in two unknowns. Rubi's 1.2.1.2 + /// (d + e x)^(7/2)/(a + b x + c x^2)^2 ran past a minute and + /// sqrt(d + e x)/(a + b x + c x^2)^3 took thirty seconds for 26,437 characters. + /// https://github.com/asc-community/AngouriMath/issues/718 + /// + internal static Entity? SolveAPowerOfALinearOverAPowerOfAQuadratic(Entity expr, Entity.Variable x, bool integrateByParts) + { + if (!Integration.AnsweringTheQuestionAsked) + return null; + var (numerator, denominator) = Functions.SingleQuotient.Of(expr); + Entity? quadratic = null; + var power = 0; + Entity? linear = null; + Entity exponent = Number.Integer.Zero; + Entity? other = null; + Entity constant = Number.Integer.One; + foreach (var (side, below) in new[] { (numerator, false), (denominator, true) }) + foreach (var factor in Mulf.LinearChildren(side)) + { + if (!factor.ContainsNode(x)) + { + constant = below ? constant / factor : constant * factor; + continue; + } + var (@base, raised) = factor is Powf(var b0, var raisedTo) && !raisedTo.ContainsNode(x) ? (b0, raisedTo) : (factor, (Entity)Number.Integer.One); + if (below && quadratic is null && raised is Number.Integer { EInteger: var whole } && whole.CompareTo(EInteger.FromInt32(2)) >= 0 + && whole.CanFitInInt32() && TreeAnalyzer.TryGetPolyQuadratic(@base, x, out var leading, out _, out _) && !VanishesIdentically(leading)) + { + quadratic = @base; + power = whole.ToInt32Unchecked(); + continue; + } + if (!TreeAnalyzer.TryGetPolyLinear(@base, x, out var slope, out _) || VanishesIdentically(slope)) + return null; + // The linear to a power that is not whole, either side; one more linear above, once. + if (raised is not Number.Integer && (linear is null || linear == @base)) + { + linear = @base; + exponent = below ? exponent - raised : exponent + raised; + continue; + } + if (!below && other is null && raised == Number.Integer.One) + { + other = @base; + continue; + } + return null; + } + exponent = exponent.InnerSimplified; + if (quadratic is null || linear is null || exponent is Number.Integer + || !TreeAnalyzer.TryGetPolyQuadratic(quadratic, x, out var c, out var b, out var a) + || !TreeAnalyzer.TryGetPolyLinear(linear, x, out var e, out var d) + || VanishesIdentically(b * b - 4 * a * c)) + return null; + // Not a positive power of x itself, which u = sqrt(x) answers sooner -- the + // reduction took `x^(9/2)/(a + b x + c x^2)^3` from 2.6 seconds to 8.3; not a quadratic + // over the imaginary unit; and not one that factors over linears, b^2 - 4 a c the + // square of a polynomial, which the partial fractions take apart cheaply. + if (linear == x && exponent is Number.Rational { IsPositive: true } || Functions.PartialFractions.HoldsTheImaginaryUnit(quadratic)) + return null; + var opposite = Functions.PartialFractions.Bare((b * b - 4 * a * c).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; + Entity nb = Number.Integer.Zero, na = constant; + if (other is not null) + { + if (!TreeAnalyzer.TryGetPolyLinear(other, x, out var ob, out var oa)) + return null; + nb = (constant * ob).InnerSimplified; + na = (constant * oa).InnerSimplified; + } + var r1 = (2 * c * d - b * e).InnerSimplified; + var r0 = (b * d - 2 * a * e).InnerSimplified; + Entity answer = Number.Integer.Zero; + for (var p = -power; p <= -2; p++) + { + // (p + 1) [r1 alpha + (r0 - r1 b/c) beta] = B, (p + 1) [r0 alpha - (r1 a/c) beta] = A. + var m11 = r1; + var m12 = r0 - r1 * b / c; + var m21 = r0; + var m22 = -r1 * a / c; + var determinant = Functions.PartialFractions.InLowestTermsOverTheSymbols((p + 1) * (m11 * m22 - m12 * m21)); + if (VanishesIdentically(determinant)) + return null; + var alpha = Functions.PartialFractions.InLowestTermsOverTheSymbols((nb * m22 - m12 * na) / determinant); + var beta = Functions.PartialFractions.InLowestTermsOverTheSymbols((m11 * na - m21 * nb) / determinant); + // N - (alpha + beta x) (p + 1) L Q', a cubic, over Q: its quotient is linear, read off + // the two highest coefficients, the remainder being zero by the choice of alpha, beta. + var c3 = -(p + 1) * beta * e * 2 * c; + var c2 = -(p + 1) * (beta * (d * 2 * c + e * b) + alpha * e * 2 * c); + var q1 = (c3 / c).InnerSimplified; + var q0 = ((c2 - q1 * b) / c).InnerSimplified; + // W = quotient - beta L - (alpha + beta x)(m + 1) e. + var wb = Functions.PartialFractions.InLowestTermsOverTheSymbols(q1 - beta * e - beta * (exponent + 1) * e); + var wa = Functions.PartialFractions.InLowestTermsOverTheSymbols(q0 - beta * d - alpha * (exponent + 1) * e); + answer = answer + (alpha + beta * x) * MathS.Pow(linear, (exponent + 1).InnerSimplified) * MathS.Pow(quadratic, Number.Integer.Create(p + 1)); + nb = wb; + na = wa; + } + // What is left asked with its two coefficients named: written out they are long, and the + // integrator took past a minute over them where with two symbols it takes a third of a + // second. + var namedA = Variable.CreateUnique(expr, "A_l"); + var namedB = Variable.CreateUnique(expr + namedA, "B_l"); + var left = (namedA + namedB * x) * MathS.Pow(linear, exponent) / quadratic; + if (Integration.ComputeAsAQuestionOfItsOwn(left, x, integrateByParts) is not { } rest) + return null; + answer = (answer + rest.Substitute(namedA, na).Substitute(namedB, nb)).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 a02c2360c..11c9fb5bd 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 power of a linear that is not whole, times a linear, over a whole power of a quadratic, + // the power of the quadratic lowered one at a time. + if ((answer = IndefiniteIntegralSolver.SolveAPowerOfALinearOverAPowerOfAQuadratic(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; diff --git a/Sources/Tests/UnitTests/Calculus/APowerOfALinearOverAPowerOfAQuadraticIntegralTest.cs b/Sources/Tests/UnitTests/Calculus/APowerOfALinearOverAPowerOfAQuadraticIntegralTest.cs new file mode 100644 index 000000000..1083b294d --- /dev/null +++ b/Sources/Tests/UnitTests/Calculus/APowerOfALinearOverAPowerOfAQuadraticIntegralTest.cs @@ -0,0 +1,59 @@ +// +// 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 half-odd power of a linear over a whole power of a quadratic, the power of the quadratic + /// lowered one at a time. Rubi's 1.2.1.2 sqrt(d + e x)/(a + b x + c x^2)^3, written here + /// with k for e, which is Euler's number. + /// #718 + /// + [Trait("Area", "Calculus")] + public sealed class APowerOfALinearOverAPowerOfAQuadraticIntegralTest + { + 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: what the integrator answers over the first power + // is piecewise on signs read where they are numbers. + 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("sqrt(d+k*x)/(a+b*x+c*x^2)^3")] + [InlineData("(d+k*x)^(7/2)/(a+b*x+c*x^2)^2")] + [InlineData("1/((d+k*x)^(3/2)*(a+b*x+c*x^2)^2)")] + 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)); + } + } +}