diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index a2e79d456..5b1e84346 100644 --- a/BREAKING-CHANGES.md +++ b/BREAKING-CHANGES.md @@ -447,6 +447,17 @@ 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 an even quadratic over a binomial in it is split over the binomial's roots + +**Answers where there were none.** `(1 - x^2)/(a + b (1 - x^2)^4)` was declined. `w^m/(a + b w^n)` for `w = p + q x^2` and +`m < n` is now split over the `n` roots `r` of `a + b w^n` into `1/(p - r + q x^2)`, each in closed form. The roots +past the first are complex, so the answer carries `i` ([#718](https://github.com/asc-community/AngouriMath/issues/718)). + +| Input | Was (2.5.0) | Now | +|---|---|---| +| `"(1-x^2)/(a+b*(1-x^2)^4)".ToEntity().Integrate("x")` | `integral(...)` | 570 characters | +| `"(1-x^2)/(a+b*(-1+x^2)^4)".ToEntity().Integrate("x")` | `integral(...)` | 539 characters | + ### A power of x over a sum of two cubes is split over the cube roots of unity **Answers where there were none.** `x^k/((3 a + b x^2)^3 + 27 a^2 c x^3)` was declined: the integrator's resultant needs diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs index f78fa02c3..dec1de4e6 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs @@ -19677,6 +19677,121 @@ Piecewise piecewise when piecewise.Cases.Any() => piecewise.Cases.First().Expres return Functions.PartialFractions.DerivativeHoldsAtSampledPoints(answer, expr, x) ? answer : null; } + /// + /// A power of an even quadratic w = p + q x^2 over a binomial in it, w^m/(a + b w^n) + /// with 0 <= m < n, split over the n roots r_j of a + b w^n: + /// w^m/(a + b w^n) = -(1/(n a)) sum_j r_j^(m + 1)/(w - r_j), and each + /// 1/(p - r_j + q x^2) is integrated in closed form. + /// + /// + /// The residue of w^m/(a + b w^n) at r_j is r_j^m/(n b r_j^(n - 1)), and + /// b r_j^n = -a makes it -r_j^(m + 1)/(n a). The roots are + /// (-a/b)^(1/n) times the n-th roots of unity. The one integral 1/(P + q x^2) is asked with + /// P named, and splits on the sign of P q, which means nothing for a complex + /// P: its first arm, the arctangent, differentiates back for any P q != 0, and + /// the answer is checked below. Rubi's 1.3.1 (1 - x^2)/(a + b (1 - x^2)^4) was declined. + /// https://github.com/asc-community/AngouriMath/issues/718 + /// + internal static Entity? SolveAPowerOfAnEvenQuadraticOverABinomialInIt(Entity expr, Entity.Variable x, bool integrateByParts) + { + if (!Integration.AnsweringTheQuestionAsked || !expr.Vars.Any(symbol => symbol != x)) + return null; + var (numerator, denominator) = Functions.SingleQuotient.Of(expr); + Entity constant = Number.Integer.One; + Entity? above = null; + Entity? below = null; + foreach (var (side, under) in new[] { (numerator, false), (denominator, true) }) + foreach (var factor in Mulf.LinearChildren(side)) + { + if (!factor.ContainsNode(x)) + constant = under ? constant / factor : constant * factor; + else if (!under && above is null) + above = factor; + else if (under && below is null) + below = factor; + else + return null; + } + if (below is null) + return null; + // a + b w^n: one term in x, a whole power of an even quadratic, beside terms free of it. + Entity? inner = null; + Entity scale = Number.Integer.One; + Entity free = Number.Integer.Zero; + var degree = 0; + foreach (var term in Sumf.LinearChildren(below)) + { + if (!term.ContainsNode(x)) + { + free = free + term; + continue; + } + if (inner is not null) + return null; + foreach (var factor in Mulf.LinearChildren(term)) + { + if (!factor.ContainsNode(x)) + scale = scale * factor; + else if (inner is null && factor is Powf(var w, Number.Integer { EInteger: var n }) && n.CompareTo(EInteger.FromInt32(2)) >= 0 + && n.CompareTo(EInteger.FromInt32(8)) <= 0) + { + inner = w; + degree = n.ToInt32Unchecked(); + } + else + return null; + } + } + if (inner is null || VanishesIdentically(free) || Functions.PartialFractions.HoldsTheImaginaryUnit(expr) + || !TreeAnalyzer.TryGetPolyQuadratic(inner, x, out var q, out var linearTerm, out var p) + || !VanishesIdentically(linearTerm) || VanishesIdentically(q)) + return null; + // The numerator a constant multiple of w^m, m below n. + var power = -1; + Entity multiple = Number.Integer.One; + if (above is null) + power = 0; + else + for (var m = 1; m < degree && power < 0; m++) + { + var wm = MathS.Pow(inner, Number.Integer.Create(m)); + if (!TreeAnalyzer.TryGetPolynomial(above, x, out var aboveRead) || !TreeAnalyzer.TryGetPolynomial(wm, x, out var wmRead) + || !aboveRead.TryGetValue(EInteger.FromInt32(2 * m), out var lead) || !wmRead.TryGetValue(EInteger.FromInt32(2 * m), out var wLead)) + continue; + var ratio = (lead / wLead).InnerSimplified; + if (!VanishesIdentically((above - ratio * wm).Expand())) + continue; + power = m; + multiple = ratio; + } + if (power < 0) + return null; + var a = free.InnerSimplified; + var b = scale.InnerSimplified; + var namedP = Variable.CreateUnique(expr, "P_b"); + if (Integration.ComputeAsAQuestionOfItsOwn(1 / (namedP + q * MathS.Pow(x, Number.Integer.Create(2))), x, integrateByParts) is not { } each) + return null; + // See the remarks: the arctangent arm for every root. + each = each.Replace(node => node switch + { + Piecewise piecewise when piecewise.Cases.Any() => piecewise.Cases.First().Expression, + Providedf provided => provided.Expression, + _ => node, + }); + var radius = MathS.Pow(Functions.PartialFractions.InLowestTermsOverTheSymbols(-a / b), Number.Rational.Create(1, degree)); + Entity sum = Number.Integer.Zero; + for (var j = 0; j < degree; j++) + { + // cos + i sin rather than (-1)^(2 j/n): an odd denominator's power of -1 is read as + // the real root, and (-1)^(2/3) is 1 (#204). + var turn = 2 * MathS.pi * Number.Rational.Create(j, degree); + Entity root = j == 0 ? radius : radius * (MathS.Cos(turn) + MathS.i * MathS.Sin(turn)).InnerSimplified; + sum = sum + MathS.Pow(root, Number.Integer.Create(power + 1)) * each.Substitute(namedP, p - root); + } + var answer = (-constant * multiple / (degree * a) * sum).InnerSimplified; + 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: diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs b/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs index 27552674c..0b3cbfba1 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs @@ -681,6 +681,8 @@ private static Entity WithItsPiecewiseSumsCombined(Entity expr) // A power of x over a sextic that is a sum of two cubes, split over the cube roots of // unity into three quadratics before the partial fractions try the symbolic sextic. if ((answer = IndefiniteIntegralSolver.SolveAPowerOfXOverASumOfTwoCubes(expr, x, integrateByParts)) is { }) return answer; + // A power of an even quadratic over a binomial in it, split over the binomial's roots. + if ((answer = IndefiniteIntegralSolver.SolveAPowerOfAnEvenQuadraticOverABinomialInIt(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/APowerOfAnEvenQuadraticOverABinomialInItIntegralTest.cs b/Sources/Tests/UnitTests/Calculus/APowerOfAnEvenQuadraticOverABinomialInItIntegralTest.cs new file mode 100644 index 000000000..fc1050c21 --- /dev/null +++ b/Sources/Tests/UnitTests/Calculus/APowerOfAnEvenQuadraticOverABinomialInItIntegralTest.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 power of an even quadratic w = p + q x^2 over a binomial in it, + /// w^m/(a + b w^n), split over the n roots of a + b w^n into + /// 1/(p - r + q x^2). Rubi's 1.3.1 (1 - x^2)/(a + b (1 - x^2)^4), and an odd + /// n, whose roots of unity are written as cos + i sin. + /// #718 + /// + [Trait("Area", "Calculus")] + public sealed class APowerOfAnEvenQuadraticOverABinomialInItIntegralTest + { + private static void DifferentiatesBack(string integrand, params (string, double)[] pins) + { + var integral = integrand.ToEntity().Integrate("x"); + Assert.DoesNotContain("integral(", integral.Stringize()); + + 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("(1-x^2)/(a+b*(1-x^2)^4)")] + [InlineData("(1-x^2)/(a+b*(-1+x^2)^4)")] + [InlineData("1/(a+b*(2+3*x^2)^3)")] + [InlineData("(2+3*x^2)^2/(a+b*(2+3*x^2)^3)")] + public void IsAnsweredAndDifferentiatesBack(string integrand) + { + DifferentiatesBack(integrand, ("a", 2), ("b", 3)); + DifferentiatesBack(integrand, ("a", -1.5), ("b", 2)); + DifferentiatesBack(integrand, ("a", 5), ("b", -0.5)); + } + } +}