diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index 71c1b2b42..77824a4f1 100644 --- a/BREAKING-CHANGES.md +++ b/BREAKING-CHANGES.md @@ -596,6 +596,20 @@ unevaluated, is answered by this route as well. | `"(a*x+2*b*n*ln(c*x^n))/(a*x^2+b*x*ln(c*x^n)^2)".ToEntity().Integrate("x")` | `integral(...)` | 37 characters | | `"tanh(a+2*ln(x))^2/x^3".ToEntity().Integrate("x")` | `integral(...)` | 7,796 characters | +### A quotient over a sum of two powers of even quadratics is split over the roots of its binomial + +**Answers where there were none.** `(1 + x^2)^2/(a x^6 + b (1 + x^2)^3)` was left unevaluated. The rule that splits +`w^m/(a + b w^n)` over the `n` roots of the binomial, for an even quadratic `w`, took a constant beside the power +only. A second power of an even quadratic `v` is read as well. `w^m v^(n - 1 - m)/(a v^n + b w^n)` is the same +sum in `t = w/v`, and each root's `1/(w - r v)` is again an even quadratic +([#718](https://github.com/asc-community/AngouriMath/issues/718)). + +| Input | Was (2.5.0) | Now | +|---|---|---| +| `"(1+x^2)^2/(a*x^6+b*(1+x^2)^3)".ToEntity().Integrate("x")` | `integral(...)` | 878 characters | +| `"x^2*(1+x^2)/(a*x^6+b*(1+x^2)^3)".ToEntity().Integrate("x")` | `integral(...)` | 896 characters | +| `"(1+x^2)/(a*x^4+b*(1+x^2)^2)".ToEntity().Integrate("x")` | `integral(...)` | 289 characters | + ### A symbolic cubic with a repeated root is written over its roots **Answers where there were none.** `1/(x^3 - 3b x + 2b^(3/2))` was left unevaluated: its coefficients are symbolic, so diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs index e41ab1075..941345a65 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs @@ -20066,6 +20066,14 @@ Piecewise piecewise when piecewise.Cases.Any() => piecewise.Cases.First().Expres /// 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. + /// + /// Two powers of even quadratics, w^m v^(n - 1 - m)/(a v^n + b w^n), are the same sum in + /// t = w/v: the integrand is t^m/(v (a + b t^n)), and 1/(v (t - r_j)) is + /// 1/(w - r_j v), again an even quadratic P + Q x^2, whose integral + /// arctan(sqrt(Q) x/sqrt(P))/(sqrt(P) sqrt(Q)) differentiates back for any P and + /// Q off zero. Rubi's 1.3.1 (1 + x^2)^2/(a x^6 + b (1 + x^2)^3), a sum of two cubes, + /// was declined. + /// /// https://github.com/asc-community/AngouriMath/issues/718 /// internal static Entity? SolveAPowerOfAnEvenQuadraticOverABinomialInIt(Entity expr, Entity.Variable x, bool integrateByParts) @@ -20081,20 +20089,19 @@ Piecewise piecewise when piecewise.Cases.Any() => piecewise.Cases.First().Expres { if (!factor.ContainsNode(x)) constant = under ? constant / factor : constant * factor; - else if (!under && above is null) - above = factor; - else if (under && below is null) + else if (!under) + above = above is null ? factor : above * factor; + else if (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; + // a + b w^n: one term in x, a whole power of an even quadratic, beside terms free of it; or + // a v^n + b w^n, two such powers, a power of x^2 among them written as one of x. + var powers = new List<(Entity Scale, Entity Base, int Degree)>(); Entity free = Number.Integer.Zero; - var degree = 0; foreach (var term in Sumf.LinearChildren(below)) { if (!term.ContainsNode(x)) @@ -20102,37 +20109,78 @@ Piecewise piecewise when piecewise.Cases.Any() => piecewise.Cases.First().Expres free = free + term; continue; } - if (inner is not null) + if (powers.Count == 2) return null; + Entity termScale = Number.Integer.One; + Entity? termBase = null; + var termDegree = 0; 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) + termScale = termScale * factor; + else if (termBase is null && factor is Powf(var w, Number.Integer { EInteger: var n }) && n.CompareTo(EInteger.FromInt32(2)) >= 0 + && n.CompareTo(EInteger.FromInt32(16)) <= 0) { - inner = w; - degree = n.ToInt32Unchecked(); + termBase = w; + termDegree = n.ToInt32Unchecked(); } else return null; } + if (termBase is null) + return null; + powers.Add((termScale, termBase, termDegree)); + } + if (powers.Count == 2) + { + if (!VanishesIdentically(free)) + return null; + // x^(2 n) beside v^n is (x^2)^n. + for (var k = 0; k < 2; k++) + if (powers[k].Base == x && powers[k].Degree % 2 == 0 && powers[1 - k].Base != x) + powers[k] = (powers[k].Scale, MathS.Pow(x, Number.Integer.Create(2)), powers[k].Degree / 2); } - if (inner is null || VanishesIdentically(free) || Functions.PartialFractions.HoldsTheImaginaryUnit(expr) + if (powers.Count == 0 || powers.Any(power => power.Degree > 8) || powers.Count == 2 && powers[0].Degree != powers[1].Degree) + return null; + var (scale, inner, degree) = powers[0]; + Entity? other = powers.Count == 2 ? powers[1].Base : null; + if (other is not null) + free = powers[1].Scale; + Entity otherP = Number.Integer.One, otherQ = Number.Integer.Zero; + if (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. + if (other is not null) + { + if (!TreeAnalyzer.TryGetPolyQuadratic(other, x, out var readQ, out var otherLinear, out var readP) + || !VanishesIdentically(otherLinear) || VanishesIdentically(readQ)) + return null; + (otherP, otherQ) = (readP, readQ); + } + // The numerator a constant multiple of w^m v^(n - 1 - m), m below n, which is w^m alone + // for a binomial: then t^m/(V (a + b t^n)) for t = w/v, split over the roots in t, and + // 1/(V (t - r)) is 1/(w - r v). var power = -1; Entity multiple = Number.Integer.One; - if (above is null) + if (above is null && other is null) power = 0; - else - for (var m = 1; m < degree && power < 0; m++) + else if (above is not null) + for (var m = other is null ? 1 : 0; 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)) + if (other is not null) + { + // No zeroth powers: InnerSimplified writes w^0 as 1 provided w is not 0. + var rest = degree - 1 - m; + wm = m == 0 ? MathS.Pow(other, Number.Integer.Create(rest)) + : rest == 0 ? wm : wm * MathS.Pow(other, Number.Integer.Create(rest)); + } + var top = other is null ? 2 * m : 2 * (degree - 1); + // Expanded for two powers: the reader takes a sum of monomials, and v^k is not one. + var (aboveAsRead, wmAsRead) = other is null ? (above, wm) : (above.Expand(), wm.Expand()); + if (!TreeAnalyzer.TryGetPolynomial(aboveAsRead, x, out var aboveRead) || !TreeAnalyzer.TryGetPolynomial(wmAsRead, x, out var wmRead) + || !aboveRead.TryGetValue(EInteger.FromInt32(top), out var lead) || !wmRead.TryGetValue(EInteger.FromInt32(top), out var wLead)) continue; var ratio = (lead / wLead).InnerSimplified; if (!VanishesIdentically((above - ratio * wm).Expand())) @@ -20144,16 +20192,20 @@ Piecewise piecewise when piecewise.Cases.Any() => piecewise.Cases.First().Expres return null; var a = free.InnerSimplified; var b = scale.InnerSimplified; + Entity? each = null; 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 + if (other is null) { - Piecewise piecewise when piecewise.Cases.Any() => piecewise.Cases.First().Expression, - Providedf provided => provided.Expression, - _ => node, - }); + if (Integration.ComputeAsAQuestionOfItsOwn(1 / (namedP + q * MathS.Pow(x, Number.Integer.Create(2))), x, integrateByParts) is not { } asked) + return null; + // See the remarks: the arctangent arm for every root. + each = asked.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++) @@ -20162,7 +20214,20 @@ Piecewise piecewise when piecewise.Cases.Any() => piecewise.Cases.First().Expres // 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); + Entity over; + if (each is not null) + over = each.Substitute(namedP, p - root); + else + { + // 1/(P + Q x^2) for w - r v: arctan(sqrt(Q) x/sqrt(P))/(sqrt(P) sqrt(Q)), whose + // derivative is 1/(P + Q x^2) for any P, Q off zero, the two roots squaring back. + var atP = Functions.PartialFractions.InLowestTermsOverTheSymbols(p - root * otherP); + var atQ = Functions.PartialFractions.InLowestTermsOverTheSymbols(q - root * otherQ); + if (VanishesIdentically(atP) || VanishesIdentically(atQ)) + return null; + over = MathS.Arctan(MathS.Sqrt(atQ) * x / MathS.Sqrt(atP)) / (MathS.Sqrt(atP) * MathS.Sqrt(atQ)); + } + sum = sum + MathS.Pow(root, Number.Integer.Create(power + 1)) * over; } var answer = (-constant * multiple / (degree * a) * sum).InnerSimplified; return Functions.PartialFractions.DerivativeHoldsAtSampledPoints(answer, expr, x) ? answer : null; diff --git a/Sources/Tests/UnitTests/Calculus/ASumOfTwoPowersOfEvenQuadraticsIntegralTest.cs b/Sources/Tests/UnitTests/Calculus/ASumOfTwoPowersOfEvenQuadraticsIntegralTest.cs new file mode 100644 index 000000000..1d65d6c32 --- /dev/null +++ b/Sources/Tests/UnitTests/Calculus/ASumOfTwoPowersOfEvenQuadraticsIntegralTest.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 quotient over a sum of two whole powers of even quadratics, a v^n + b w^n, split over the + /// roots of a + b t^n in t = w/v: Rubi's 1.3.1 (1 + x^2)^2/(a x^6 + b (1 + x^2)^3), + /// which was left unevaluated. + /// #718 + /// + [Trait("Area", "Calculus")] + public sealed class ASumOfTwoPowersOfEvenQuadraticsIntegralTest + { + 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)^2/(a*x^6+b*(1+x^2)^3)")] + [InlineData("x^2*(1+x^2)/(a*x^6+b*(1+x^2)^3)")] + [InlineData("x^4/(a*x^6+b*(1+x^2)^3)")] + [InlineData("(2+3*x^2)^2/(a*x^6+b*(2+3*x^2)^3)")] + [InlineData("(1+x^2)/(a*x^4+b*(1+x^2)^2)")] + public void IsAnsweredAndDifferentiatesBack(string integrand) + { + DifferentiatesBack(integrand, ("a", 2), ("b", 3)); + DifferentiatesBack(integrand, ("a", -5), ("b", 0.5)); + } + } +}