diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index f6dd61c5c..0e797a69d 100644 --- a/BREAKING-CHANGES.md +++ b/BREAKING-CHANGES.md @@ -336,6 +336,23 @@ quotient of two such linears the sum is `(b - d t^2)^2 + (c t^2 - a)^2`. Rubi's | `"(c + d*tan(x))^(3/2)/(a + b*tan(x))^3".ToEntity().Integrate("x")` | `integral(...)` | the same | | `"1/((a + b*tan(x))^(3/2)*(c + d*tan(x))^(3/2))".ToEntity().Integrate("x")` | `integral(...)`; past a minute on the unreleased master | in the root of the quotient of the two, in a second | +### A quotient in `x^2` over a power of a linear in `x^2` and a biquadratic is split in `x^2` + +**Shorter answers, sooner.** `sqrt(c + d tan(x)) (A + B tan(x) + C tan(x)^2)/(a + b tan(x))^3` was +declined in 2.5.0 and answered in 3.3 million characters since; it is answered in seven thousand. +Under `u = tan(x)` and `t = sqrt(c + d u)` it is a polynomial in `t^2` over +`(a d + b (t^2 - c))^3 ((t^2 - c)^2 + d^2)`, which was split in `t`, the sum of two squares over its +conjugates. Such a quotient, one linear in `x^2` to a power beside one biquadratic, is split in +`s = x^2 - r` now, `r` the linear's root: in powers of `s` it is the split over a power of `x^2` and a +biquadratic that the rule above makes, and the terms in `1/(x^2 - r)^j` go by their reduction to +`atan(x/sqrt(-r))/sqrt(-r)` ([#718](https://github.com/asc-community/AngouriMath/issues/718)). + +| Input | Was (2.5.0) | Now | +|---|---|---| +| `"sqrt(c + d*tan(x))*(A + B*tan(x) + C*tan(x)^2)/(a + b*tan(x))^3".ToEntity().Integrate("x")` | `integral(...)`; 3,299,915 characters after 24 seconds on the unreleased master | 6,912 characters | +| `"sqrt(c + d*tan(x))/(a + b*tan(x))^2".ToEntity().Integrate("x")` | `integral(...)`; 40,540 characters on the unreleased master | 1,465 characters | +| `"(c + d*tan(x))^(3/2)/(a + b*tan(x))^3".ToEntity().Integrate("x")` | `integral(...)`; 71,293 characters on the unreleased master | 2,786 characters | + ### A quotient in `x^2` over a power of `x` and a biquadratic is split in `x^2` **Shorter answers, sooner.** `cot(x)^(13/2) (a + b tan(x))^(5/2) (A + B tan(x))` was declined in 2.5.0 diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs index 49cf0c67e..729347128 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs @@ -2485,59 +2485,214 @@ over is Entity.Powf(var @base, var power) ? if (Functions.PartialFractions.IsZeroAsAValue(discriminant)) return null; - // The numerator in w = x^2, divided down by w^k (c w^2 + b w + a) to a remainder of a - // degree below k + 2. + // The numerator in w = x^2 over w^k (c w^2 + b w + a). var degree = above.Keys.Max() / 2; - var inW = new Entity[System.Math.Max(degree, k + 1) + 1]; - for (var j = 0; j < inW.Length; j++) + var inW = new Entity[degree + 1]; + for (var j = 0; j <= degree; j++) inW[j] = above.TryGetValue(2 * j, out var at) ? at : Number.Integer.Zero; + var (polynomial, pole, d, e) = SplitOverAPowerAndAQuadratic(inW, k, a, b, c); + Entity polynomialPart = Number.Integer.Zero; + for (var j = polynomial.Length - 1; j >= 0; j--) + if (polynomial[j] != Number.Integer.Zero) + polynomialPart += polynomial[j] * MathS.Pow(x, 2 * j); + + Entity total = Number.Integer.Zero; + if (polynomialPart != Number.Integer.Zero) + { + if (Integration.ComputeIndefiniteIntegral(polynomialPart, x, integrateByParts) is not { } whole) + return null; + total += whole; + } + // r_m x^(2m - 2k), whose exponent is odd and so never minus one. + for (var m = 0; m < k; m++) + { + if (pole[m] == Number.Integer.Zero) + continue; + var exponent = 2 * (m - k) + 1; + total += pole[m] * MathS.Pow(x, exponent) / exponent; + } + return ArctangentsOverABiquadratic(total, d, e, a, b, c, discriminant, x); + } + + /// + /// A polynomial in x^2 over a power of a linear in x^2 and a biquadratic, with a + /// symbol in it: N(x^2)/(L(x^2)^k Q(x^2)), L(w) = l (w - r), split in + /// s = x^2 - r. Written in powers of s, the quotient is the one + /// splits over w^k Q, and the + /// terms in s^(-j) go by the reduction + /// int dx/(x^2 + p)^j = x/(2 p (j - 1) (x^2 + p)^(j - 1)) + (2 j - 3)/(2 p (j - 1)) int dx/(x^2 + p)^(j - 1), + /// p = -r, down to atan(x/sqrt(p))/sqrt(p), which holds for every complex + /// p but zero, as the biquadratic's own terms do. + /// + /// + /// Under t = sqrt(c + d u) Rubi's sqrt(c + d u)/((a + b u)^n (1 + u^2)), which the + /// tangent substitution makes of sqrt(c + d tan(x))/(a + b tan(x))^n, is + /// t^2/((a d + b (t^2 - c))^n ((t^2 - c)^2 + d^2)) up to a constant. Split in t, the + /// sum of two squares went over its conjugates, and the answer to + /// sqrt(c + d tan(x)) (A + B tan(x) + C tan(x)^2)/(a + b tan(x))^3 ran to 3.5 million + /// characters. Only one linear in x^2, to a power, beside one biquadratic: with + /// r zero it is 's, and a root + /// r shared with the biquadratic declines. + /// https://github.com/asc-community/AngouriMath/issues/718 + /// + internal static Entity? SolveAnEvenQuotientOverAPowerOfALinearInTheSquareAndABiquadratic(Entity expr, Entity.Variable x, bool integrateByParts) + { + if (!TryReadAsQuotient(expr, out var numerator, out var denominator) || !denominator.Vars.Any(v => v != x)) + return null; + Entity constant = Number.Integer.One; + (Entity Zeroth, Entity First, int Power)? linear = null; + (Entity A, Entity B, Entity C)? quadratic = null; + foreach (var factor in Mulf.LinearChildren(denominator)) + { + var (@base, power) = factor is Powf(var raised, Number.Integer { EInteger: var n }) && n.Sign > 0 && n.CanFitInInt32() + ? (raised, n.ToInt32Unchecked()) : (factor, 1); + if (!@base.ContainsNode(x)) + { + constant *= factor; + continue; + } + if (!TreeAnalyzer.TryGetPolynomial(@base, x, out var read) || read.Count == 0 + || read.Values.Any(coefficient => coefficient.ContainsNode(x)) + || read.Keys.Any(p => p.Sign < 0 || !p.IsEven || !p.CanFitInInt32())) + return null; + Entity At(int p) => read.TryGetValue(EInteger.FromInt32(p), out var value) ? value : Number.Integer.Zero; + var top = read.Keys.Max()!.ToInt32Unchecked(); + if (top == 2 && linear is null) + linear = (At(0), At(2), power); + else if (top == 4 && quadratic is null && power == 1) + quadratic = (At(0), At(2), At(4)); + else + return null; + } + if (linear is not { } l || quadratic is not { } quartic) + return null; + if (!TreeAnalyzer.TryGetPolynomial(numerator, x, out var above) || above.Count == 0 + || above.Values.Any(coefficient => coefficient.ContainsNode(x)) + || above.Keys.Any(p => p.Sign < 0 || !p.IsEven || !p.CanFitInInt32())) + return null; + var (a, b, c) = quartic; + if (Functions.PartialFractions.IsZeroAsAValue(a) || Functions.PartialFractions.IsZeroAsAValue(c) + || Functions.PartialFractions.IsZeroAsAValue(l.First)) + return null; + var r = Functions.PartialFractions.InLowestTermsOverTheSymbols(-l.Zeroth / l.First); + if (Functions.PartialFractions.IsZeroAsAValue(r)) + return null; + var discriminant = Functions.PartialFractions.Bare((b * b - 4 * a * c).Simplify()); + if (Functions.PartialFractions.IsZeroAsAValue(discriminant)) + return null; + // The biquadratic in s: Q(r + s) = Q(r) + Q'(r) s + c s^2, and Q(r) is not zero, or the + // linear's root is one of the biquadratic's. + var atRoot = Functions.PartialFractions.InLowestTermsOverTheSymbols(a + b * r + c * r * r); + if (Functions.PartialFractions.IsZeroAsAValue(atRoot)) + return null; + var slopeAtRoot = Functions.PartialFractions.InLowestTermsOverTheSymbols(b + 2 * c * r); + + // N(r + s) in powers of s. + var degree = above.Keys.Max()!.ToInt32Unchecked() / 2; + var inW = new Entity[degree + 1]; + for (var j = 0; j <= degree; j++) + inW[j] = above.TryGetValue(EInteger.FromInt32(2 * j), out var at) ? at : Number.Integer.Zero; + var inS = new Entity[degree + 1]; + for (var m = 0; m <= degree; m++) + { + Entity sum = Number.Integer.Zero; + EInteger binomial = EInteger.One; + for (var j = m; j <= degree; j++) + { + if (j > m) + binomial = binomial * j / (j - m); + if (inW[j] != Number.Integer.Zero) + sum += inW[j] * Number.Integer.Create(binomial) * MathS.Pow(r, Number.Integer.Create(j - m)); + } + inS[m] = Functions.PartialFractions.InLowestTermsOverTheSymbols(sum); + } + var k = l.Power; + var (polynomial, pole, dInS, eInS) = SplitOverAPowerAndAQuadratic(inS, k, atRoot, slopeAtRoot, c); + var square = MathS.Pow(x, 2) - r; + + Entity total = Number.Integer.Zero; Entity polynomialPart = Number.Integer.Zero; + for (var j = 0; j < polynomial.Length; j++) + if (polynomial[j] != Number.Integer.Zero) + polynomialPart += polynomial[j] * MathS.Pow(square, j); + if (polynomialPart != Number.Integer.Zero) + { + if (Integration.ComputeIndefiniteIntegral(polynomialPart.Expand(), x, integrateByParts) is not { } whole) + return null; + total += whole; + } + // pole[m] s^(m - k) is a term in 1/(x^2 - r)^j with j = k - m, by the reduction. + var p = -r; + var root = MathS.Sqrt(p); + Entity reduced = MathS.Arctan(x / root) / root; + var byPower = new Entity[k + 1]; + byPower[1] = reduced; + for (var j = 2; j <= k; j++) + byPower[j] = x / (2 * p * (j - 1) * MathS.Pow(square, j - 1)) + Number.Integer.Create(2 * j - 3) / (2 * p * (j - 1)) * byPower[j - 1]; + for (var m = 0; m < k; m++) + if (pole[m] != Number.Integer.Zero) + total += pole[m] * byPower[k - m]; + // (d + e s)/Q in w = x^2 is (d - e r + e w)/Q(w). + var d = Functions.PartialFractions.InLowestTermsOverTheSymbols(dInS - eInS * r); + total = ArctangentsOverABiquadratic(total, d, eInS, a, b, c, discriminant, x); + return total / (constant * MathS.Pow(l.First, k)); + } + + /// + /// N(s)/(s^k (a + b s + c s^2)), the coefficients of N in , + /// as a polynomial in s, the terms in s^(m - k) for m below k, and + /// (d + e s)/(a + b s + c s^2): N divided down by s^k times the quadratic first, + /// then the remainder R over it expanded at s = 0 by + /// r_m = (R_m - b r_(m-1) - c r_(m-2))/a, and d + e s what + /// (R - Q sum r_m s^m)/s^k leaves, every lower coefficient cancelling by the recurrence. + /// + private static (Entity[] Polynomial, Entity[] Pole, Entity D, Entity E) SplitOverAPowerAndAQuadratic( + Entity[] above, int k, Entity a, Entity b, Entity c) + { + var degree = above.Length - 1; + var inS = new Entity[System.Math.Max(degree, k + 1) + 1]; + for (var j = 0; j < inS.Length; j++) + inS[j] = j <= degree ? above[j] : Number.Integer.Zero; + var polynomial = new Entity[System.Math.Max(degree - k - 1, 0)]; + for (var j = 0; j < polynomial.Length; j++) + polynomial[j] = Number.Integer.Zero; for (var j = degree; j >= k + 2; j--) { - var lead = Functions.PartialFractions.InLowestTermsOverTheSymbols(inW[j] / c); + var lead = Functions.PartialFractions.InLowestTermsOverTheSymbols(inS[j] / c); if (lead == Number.Integer.Zero) continue; - polynomialPart += lead * MathS.Pow(x, 2 * (j - k - 2)); - inW[j - 1] = Functions.PartialFractions.InLowestTermsOverTheSymbols(inW[j - 1] - lead * b); - inW[j - 2] = Functions.PartialFractions.InLowestTermsOverTheSymbols(inW[j - 2] - lead * a); + polynomial[j - k - 2] = lead; + inS[j - 1] = Functions.PartialFractions.InLowestTermsOverTheSymbols(inS[j - 1] - lead * b); + inS[j - 2] = Functions.PartialFractions.InLowestTermsOverTheSymbols(inS[j - 2] - lead * a); } - // The remainder R over w^k Q is sum r_m w^(m - k) over m < k, the expansion of R/Q at - // w = 0, plus (d + e w)/Q: r_m = (R_m - b r_(m-1) - c r_(m-2))/a, and d + e w is what - // (R - Q sum r_m w^m)/w^k leaves, every lower coefficient cancelling by the recurrence. var pole = new Entity[k]; for (var m = 0; m < k; m++) { - var term = inW[m]; + var term = inS[m]; if (m >= 1) term -= b * pole[m - 1]; if (m >= 2) term -= c * pole[m - 2]; pole[m] = Functions.PartialFractions.InLowestTermsOverTheSymbols(term / a); } - var dTerm = inW[k]; + if (k == 0) + return (polynomial, pole, inS[0], inS[1]); + var dTerm = inS[k]; if (k >= 1) dTerm -= b * pole[k - 1]; if (k >= 2) dTerm -= c * pole[k - 2]; - var eTerm = inW[k + 1]; + var eTerm = inS[k + 1]; if (k >= 1) eTerm -= c * pole[k - 1]; - var d = k == 0 ? inW[0] : Functions.PartialFractions.InLowestTermsOverTheSymbols(dTerm); - var e = k == 0 ? (degree >= 1 ? inW[1] : Number.Integer.Zero) : Functions.PartialFractions.InLowestTermsOverTheSymbols(eTerm); + return (polynomial, pole, Functions.PartialFractions.InLowestTermsOverTheSymbols(dTerm), + Functions.PartialFractions.InLowestTermsOverTheSymbols(eTerm)); + } + /// + /// plus the antiderivative of (d + e x^2)/(a + b x^2 + c x^4) by + /// the two roots in x^2, being b^2 - 4 a c and not zero. + /// + private static Entity ArctangentsOverABiquadratic(Entity total, Entity d, Entity e, Entity a, Entity b, Entity c, Entity discriminant, Entity.Variable x) + { var q = MathS.Sqrt(discriminant); var firstRoot = (-b + q) / (2 * c); var secondRoot = (-b - q) / (2 * c); - Entity total = Number.Integer.Zero; - if (polynomialPart != Number.Integer.Zero) - { - if (Integration.ComputeIndefiniteIntegral(polynomialPart, x, integrateByParts) is not { } whole) - return null; - total += whole; - } - // r_m x^(2m - 2k), whose exponent is odd and so never minus one. - for (var m = 0; m < k; m++) - { - if (pole[m] == Number.Integer.Zero) - continue; - var exponent = 2 * (m - k) + 1; - total += pole[m] * MathS.Pow(x, exponent) / exponent; - } // `1/(x^2 - r)` is `atan(x/s)/s` with `s = sqrt(-r)` for every complex `r` but zero: // `d/dx atan(x/s)/s` is `1/(s^2 + x^2)` whatever `s` is, and for a positive `r` the // arctangent of an imaginary argument is the hyperbolic one, `-atanh(x/sqrt(r))/sqrt(r)`, diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs b/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs index c87ab1d23..15073706b 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs @@ -1017,6 +1017,8 @@ private static Entity Normalized(Entity expr, Entity.Variable x) => // in x^2 written with the root of the discriminant: partial fractions read a // written factor, and `a + b x^2 + c x^4` is written as one. if ((answer = IndefiniteIntegralSolver.SolveAnEvenPolynomialOverASymbolicBiquadratic(expr, x, integrateByParts)) is { }) return answer; + // And over a power of a linear in x^2 beside the biquadratic, split in x^2 minus its root. + if ((answer = IndefiniteIntegralSolver.SolveAnEvenQuotientOverAPowerOfALinearInTheSquareAndABiquadratic(expr, x, integrateByParts)) is { }) return answer; if ((answer = IndefiniteIntegralSolver.SolveByPartialFractions(expr, x, integrateByParts)) is { }) return answer; // A whole negative power of a polynomial of several terms among the factors, // written below the bar and asked again: the gathering on the way in writes diff --git a/Sources/Tests/UnitTests/Calculus/PowerOfALinearInTheSquareBesideABiquadraticIntegralTest.cs b/Sources/Tests/UnitTests/Calculus/PowerOfALinearInTheSquareBesideABiquadraticIntegralTest.cs new file mode 100644 index 000000000..259c4e884 --- /dev/null +++ b/Sources/Tests/UnitTests/Calculus/PowerOfALinearInTheSquareBesideABiquadraticIntegralTest.cs @@ -0,0 +1,54 @@ +// +// 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 polynomial in x^2 over a power of a linear in x^2 and a biquadratic, split in + /// x^2 minus the linear's root, the powers of 1/(x^2 - r) by their reduction to an + /// arctangent. Under the tangent and t = sqrt(c + d u), Rubi's + /// sqrt(c + d tan(x))/(a + b tan(x))^n comes to one; split in t, the answer to the + /// third row ran to 3.5 million characters. Rubi's 4.3.2.1 and 4.3.4.2. + /// #718 + /// + /// + /// Real where tan(x), a + b tan(x) and c + d tan(x) are positive, which the + /// points are. The pinned answer is differentiated as it is. + /// + [Trait("Area", "Calculus")] + public sealed class PowerOfALinearInTheSquareBesideABiquadraticIntegralTest + { + [Theory] + [InlineData("x^2/((k + x^2)^2*(a + b*x^2 + c*x^4))")] + [InlineData("sqrt(c + d*x)/((a + b*x)^3*(1 + x^2))")] + [InlineData("sqrt(c + d*tan(x))*(A + B*tan(x) + M*tan(x)^2)/(a + b*tan(x))^3")] + [InlineData("(c + d*tan(x))^(3/2)/(a + b*tan(x))^3")] + public void SplitInTheSquareLessTheRoot(string integrand) + { + var integral = integrand.ToEntity().Integrate("x"); + var text = integral.Stringize(); + Assert.DoesNotContain("integral(", text); + Assert.True(text.Length < 10000, $"{text.Length} characters of answer for {integrand}"); + Entity Pinned(Entity e) => e.Substitute("a", 1.3).Substitute("b", 0.7).Substitute("c", 1.1).Substitute("d", 0.6) + .Substitute("k", 0.9).Substitute("A", 0.4).Substitute("B", 1.1).Substitute("M", 0.3); + var derivative = Pinned(integral.Substitute("C", 0)).Differentiate("x"); + var original = Pinned(integrand.ToEntity()); + foreach (var at in new[] { 0.2, 0.5, 0.9, 1.2 }) + { + 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}"); + } + } + } +}