diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index 4093449e2..422854f4c 100644 --- a/BREAKING-CHANGES.md +++ b/BREAKING-CHANGES.md @@ -614,6 +614,22 @@ polynomial in the symbols is written as the two quadratics first | `"cot(e + f*x)^3*sqrt(a + b*sec(e + f*x))".ToEntity().Integrate("x")` | `integral(...)`; past a minute on the unreleased master | 36,760 characters, a case for each sign of the roots | | `"sec(e + f*x)^2/(a*sin(e + f*x) + b*tan(e + f*x))^2".ToEntity().Integrate("x")` | `integral(...)`; past a minute on the unreleased master | 68,854 characters | +### A term with a power of `x` is integrated by parts with the rest of the sum in its remainder + +**Answers where there were none.** `sqrt(csc(x)) (x cos(x) - 4 sec(x) tan(x))`, Rubi's 4.7.7, was declined: neither +of its two terms is elementary alone. Where a sum's terms are not integrated one by one, a term `x^k f` whose `f` +has `x` only inside transcendental functions is now integrated by parts, and what parts leaves is asked together +with the other terms. Here that is `-2 sqrt(sin(x)) - 4 sqrt(sin(x))/cos(x)^2`, a power of the sine over a power of +the cosine, which is integrated in closed form +([#718](https://github.com/asc-community/AngouriMath/issues/718)). + +| Input | Was (2.5.0) | Now | +|---|---|---| +| `"sqrt(csc(x))*(x*cos(x) - 4*sec(x)*tan(x))".ToEntity().Integrate("x")` | `integral(...)` | `x * -(1 / sin(x)) ^ (-1/2) / (-1/2) + (-4) * sin(x) ^ (3/2) / cos(x)` | +| `"x*sin(x)/sqrt(cos(x)) + 4*sqrt(cos(x))/sin(x)^2".ToEntity().Integrate("x")` | `integral(...)` | `x * -sqrt(cos(x)) / (1/2) + (-4) / sin(x) * cos(x) ^ (3/2)` | +| `"x^2*cos(x)*e^sin(x) + 2*x*e^sin(x)".ToEntity().Integrate("x")` | `integral(...)` | `x ^ 2 * e ^ sin(x)` | +| `"(c + d*x)*erfc(a + b*x)^2".ToEntity().Integrate("x")` | 187 characters, reading `erfc` as a symbol; `integral(...)` on the unreleased master | 1,196 characters | + ### An odd power of the tangent beside a function of the secant is integrated in the secant **Answers where there were none.** `cot(c + d x)^3/(a + b sec(c + d x))^(3/2)` ran past a minute. In `u = sec(y)`, diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs index 0a1498898..2aedd1ee0 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs @@ -99,6 +99,8 @@ internal static class IndefiniteIntegralSolver // corpus for nothing. if (!Integration.AnsweringTheQuestionAsked) return null; + if (integrateByParts && ByPartsAgainstThePowersOfX(splitted, x) is { } byPartsOverTheSum) + return byPartsOverTheSum; // What is polynomial in a term is what has no root of the variable in it: a // logarithm or an exponential of x is an indeterminate as far as gathering and // cancelling go. Bronstein's @@ -143,6 +145,276 @@ internal static class IndefiniteIntegralSolver }); } + /// + /// A sum whose terms are not elementary apart, integrated by parts against the power + /// of in each term that has one, with what parts leaves and the + /// other terms asked together as one integrand: x^k f + g is + /// x^k F + ∫(g - k x^(k - 1) F) for F the antiderivative of f. + /// + /// + /// + /// sqrt(csc(x)) (x cos(x) - 4 sec(x) tan(x)), Rubi's 4.7.7, is + /// (2 x sqrt(sin(x)) - 4 sin(x)^(3/2)/cos(x))'. Split into its two terms, neither is + /// elementary: parts on the first leaves -2 sqrt(sin(x)), and that and the second + /// term are elementary only together, as sqrt(sin(x)) (-2 - 4/cos(x)^2). So the + /// remainder is gathered into one quotient with the factors its terms share taken out, + /// which is the spelling the rules read; and where it is zero, as for + /// x cos(x) e^sin(x) + e^sin(x), the parts are the answer. + /// + /// + /// Only where some term with a power of beside something else in + /// has an elementary F, and only at the question asked: the + /// remainder is a new question, and the powers of of the terms + /// integrated are lower in it, so asking it again descends. + /// + /// + private static Entity? ByPartsAgainstThePowersOfX(IReadOnlyList terms, Entity.Variable x) + { + // Once: the cofactors and the remainder are asked at the top, where this runs again, + // and `(c + d x) erfc(a + b x)^2` went from a 13-second decline to a hang that way. + if (integratingASumByParts) + { + Integration.DeclinedForItsScope(); + return null; + } + integratingASumByParts = true; + try + { + return ByPartsAgainstThePowersOfXOnce(terms, x); + } + finally + { + integratingASumByParts = false; + } + } + + /// + /// Whether a term of has a factor that is a polynomial + /// in of positive degree. + /// + private static bool HasAPolynomialFactorIn(Entity antiderivative, Entity.Variable x) + => Sumf.LinearChildren(antiderivative).Any(term => Mulf.LinearChildren(term).Any(factor => + factor.ContainsNode(x) && TreeAnalyzer.TryGetPolynomial(factor, x, out var powers) + && powers.Keys.All(power => power.Sign >= 0) && powers.Values.All(coefficient => !coefficient.ContainsNode(x)))); + + /// + /// Whether every occurrence of in is inside + /// the argument of a function that is not algebraic: cos(x)/sqrt(sin(x)), not + /// sqrt(x^2 + 3) or atanh(x)/(1 - x^2). + /// + private static bool OnlyInsideTranscendentalFunctions(Entity expr, Entity.Variable x) + => !expr.ContainsNode(x) + || expr is not (Sumf or Minusf or Mulf or Divf or Powf or Entity.Variable or Number) + || expr is not Entity.Variable && expr.DirectChildren.All(child => OnlyInsideTranscendentalFunctions(child, x)); + + /// Set while asks its questions. + [System.ThreadStatic] private static bool integratingASumByParts; + + private static Entity? ByPartsAgainstThePowersOfXOnce(IReadOnlyList terms, Entity.Variable x) + { + Entity parts = Number.Integer.Zero; + var remainder = new List(); + var any = false; + foreach (var term in terms) + { + var power = EInteger.Zero; + Entity rest = Number.Integer.One; + foreach (var factor in Mulf.LinearChildren(term)) + if (factor == x) + power = power.Add(EInteger.One); + else if (factor is Powf(var @base, Number.Integer { EInteger.Sign: > 0 } exponent) && @base == x) + power = power.Add(exponent.EInteger); + else + rest *= factor; + // Only a cofactor with x in it nowhere but inside transcendental functions. A power + // of x beside an algebraic one is what the algebraic reductions answer, and asking + // each such cofactor and the remainder took `(x^2 + 3)^2/(x + sqrt(x^2 + 3))` from a + // half-second decline past two minutes; and beside a rational function of x as + // well, `x^3 atanh(c x^2)/(1 - c^2 x^4)` in what parts leaves of + // `(a + b atanh(c x^2))^2/x^5`, the cofactor is a slow decline that doubled its time. + if (power.IsZero || !rest.ContainsNode(x) || !OnlyInsideTranscendentalFunctions(rest, x)) + { + remainder.Add(term); + continue; + } + // A term whose cofactor is not elementary goes into the remainder as it is: in + // `x^2 cos(x) e^sin(x) + 2 x e^sin(x)` the second is what parts on the first leaves. + // And one with a polynomial in x beside something in its antiderivative: what parts + // leaves then has the power it took away, and asking it does not descend -- + // `x erfc(a + b x)^2`, whose cofactor's antiderivative has `x erfc(a + b x)^2` in it. + if (Integration.ComputeAsAQuestionOfItsOwn(rest, x, integrateByParts: true) is not { } antiderivative + || HasAPolynomialFactorIn(antiderivative, x)) + { + remainder.Add(term); + continue; + } + var monomial = power.Equals(EInteger.One) ? (Entity)x : MathS.Pow(x, Number.Integer.Create(power)); + parts += monomial * antiderivative; + // No `x^0` for a first power: it would bring its `provided not x = 0` along. + var lower = power.Equals(EInteger.One) ? Number.Integer.One : MathS.Pow(x, Number.Integer.Create(power.Subtract(EInteger.One))); + remainder.Add(-Number.Integer.Create(power) * lower * antiderivative); + any = true; + } + if (!any) + return null; + var left = remainder.Aggregate((l, r) => l + r).InnerSimplified; + if (TreeAnalyzer.IsZero(left)) + return parts; + var gathered = OverTheLeastPowersOfTheSineAndTheCosine(left, x, out var closed); + if (closed is not null) + return parts + closed; + if (gathered is null) + { + var (above, below) = Functions.SingleQuotient.Of(left); + gathered = WithTheFactorsEveryTermSharesOut(above.InnerSimplified, x) / below; + } + return Integration.ComputeAsAQuestionOfItsOwn(gathered, x, integrateByParts: true) is { } integrated + ? parts + integrated + : null; + } + + /// + /// with each term read as a number times powers of + /// sin(u) and cos(u) for one argument u, and written as the least + /// power of each times the sum of what is left; where a term + /// is not such a product. + /// + /// + /// What parts leaves is spelled the way the antiderivative came back, and the other + /// terms the way they were written: (1/sin(x))^(-1/2) beside + /// sqrt(csc(x)) sec(x) tan(x) share nothing as written, and are + /// sqrt(sin(x)) and sqrt(sin(x))/cos(x)^2 as powers. Reading + /// (1/f)^p as f^(-p) holds wherever f is positive or p is + /// whole; where neither, the integrand is not real. + /// + /// And where the sum is the derivative of one such + /// product, which the rules otherwise answer in the half angle's tangent at four times + /// the size: (K sin(u)^(a + 1) cos(u)^(b + 1))' is + /// K u' sin(u)^a cos(u)^b ((a + b + 2) cos(u)^2 - (b + 1)), and + /// K u' sin(u)^a cos(u)^b ((a + 1) - (a + b + 2) sin(u)^2) the other way, so + /// sqrt(sin(x)) (-2 - 4/cos(x)^2) is (-4 sin(x)^(3/2)/cos(x))'. + /// + /// + private static Entity? OverTheLeastPowersOfTheSineAndTheCosine(Entity sum, Entity.Variable x, out Entity? antiderivative) + { + antiderivative = null; + Entity? argument = null; + var read = new List<(Entity Coefficient, Number.Rational Sine, Number.Rational Cosine)>(); + foreach (var term in TreeAnalyzer.GatherLinearChildrenOverSumAndExpand(sum, e => e.ContainsNode(x)) ?? new List { sum }) + if (Read(term) is { } powers) + read.Add(powers); + else + return null; + if (argument is null || read.Count < 2) + return null; + var leastSine = read.Select(r => r.Sine).Aggregate((l, r) => r < l ? r : l); + var leastCosine = read.Select(r => r.Cosine).Aggregate((l, r) => r < l ? r : l); + Entity Power(Entity @base, Number.Rational power) + => power == Number.Integer.Zero ? Number.Integer.One : power == Number.Integer.One ? @base : MathS.Pow(@base, power); + var slope = argument.Differentiate(x).InnerSimplified; + var two = Number.Integer.Create(2); + if (!slope.ContainsNode(x)) + { + Entity Total(System.Func<(Entity Coefficient, Number.Rational Sine, Number.Rational Cosine), bool> which) + => read.Where(which).Select(r => r.Coefficient).Aggregate((Entity)Number.Integer.Zero, (l, r) => l + r).InnerSimplified; + var (a, b) = (leastSine, leastCosine); + Entity? k = null; + if (read.All(r => r.Sine == a && (r.Cosine == b || r.Cosine == b + two)) && b + Number.Integer.One != Number.Integer.Zero) + { + var (squared, alone) = (Total(r => r.Cosine != b), Total(r => r.Cosine == b)); + var candidate = (-alone / ((b + Number.Integer.One) * slope)).InnerSimplified; + if (TreeAnalyzer.IsZero((squared - candidate * slope * (a + b + two)).InnerSimplified)) + k = candidate; + } + else if (read.All(r => r.Cosine == b && (r.Sine == a || r.Sine == a + two)) && a + Number.Integer.One != Number.Integer.Zero) + { + var (squared, alone) = (Total(r => r.Sine != a), Total(r => r.Sine == a)); + var candidate = (alone / ((a + Number.Integer.One) * slope)).InnerSimplified; + if (TreeAnalyzer.IsZero((squared + candidate * slope * (a + b + two)).InnerSimplified)) + k = candidate; + } + if (k is not null) + { + var (sinePower, cosinePower) = (a + Number.Integer.One, b + Number.Integer.One); + Entity product = k; + if (sinePower < Number.Integer.Zero) product /= Power(MathS.Sin(argument), -sinePower); else product *= Power(MathS.Sin(argument), sinePower); + if (cosinePower < Number.Integer.Zero) product /= Power(MathS.Cos(argument), -cosinePower); else product *= Power(MathS.Cos(argument), cosinePower); + antiderivative = product; + } + } + var left = read.Select(r => r.Coefficient * Power(MathS.Sin(argument), r.Sine - leastSine) * Power(MathS.Cos(argument), r.Cosine - leastCosine)) + .Aggregate((l, r) => l + r); + Entity above = left, below = Number.Integer.One; + if (leastSine < Number.Integer.Zero) below *= Power(MathS.Sin(argument), -leastSine); else above = Power(MathS.Sin(argument), leastSine) * above; + if (leastCosine < Number.Integer.Zero) below *= Power(MathS.Cos(argument), -leastCosine); else above = Power(MathS.Cos(argument), leastCosine) * above; + return below == Number.Integer.One ? above : above / below; + + (Entity Coefficient, Number.Rational Sine, Number.Rational Cosine)? Read(Entity node) + { + if (!node.ContainsNode(x)) + return (node, Number.Integer.Zero, Number.Integer.Zero); + switch (node) + { + case Sinf(var u): return Base(u, Number.Integer.One, Number.Integer.Zero); + case Cosf(var u): return Base(u, Number.Integer.Zero, Number.Integer.One); + case Secantf(var u): return Base(u, Number.Integer.Zero, Number.Integer.MinusOne); + case Cosecantf(var u): return Base(u, Number.Integer.MinusOne, Number.Integer.Zero); + case Tanf(var u): return Base(u, Number.Integer.One, Number.Integer.MinusOne); + case Cotanf(var u): return Base(u, Number.Integer.MinusOne, Number.Integer.One); + case Mulf(var l, var r): + return (Read(l), Read(r)) is ({ } a, { } b) ? ((a.Coefficient * b.Coefficient).InnerSimplified, a.Sine + b.Sine, a.Cosine + b.Cosine) : null; + case Divf(var l, var r): + return (Read(l), Read(r)) is ({ } n, { } d) ? ((n.Coefficient / d.Coefficient).InnerSimplified, n.Sine - d.Sine, n.Cosine - d.Cosine) : null; + case Powf(var @base, Number.Rational power) when Read(@base) is { } inner: + // A whole power of anything read; a fractional one only of a single + // power of the sine or of the cosine with nothing beside it, so that + // no root of a product is split. + if (power is Number.Integer) + return (MathS.Pow(inner.Coefficient, power).InnerSimplified, inner.Sine * power, inner.Cosine * power); + if (inner.Coefficient == Number.Integer.One && (inner.Sine == Number.Integer.Zero || inner.Cosine == Number.Integer.Zero)) + return (Number.Integer.One, inner.Sine * power, inner.Cosine * power); + return null; + default: + return null; + } + } + + (Entity Coefficient, Number.Rational Sine, Number.Rational Cosine)? Base(Entity u, Number.Rational sine, Number.Rational cosine) + { + if (argument is not null && argument != u) + return null; + argument = u; + return (Number.Integer.One, sine, cosine); + } + } + + /// + /// as the factors with in them that every + /// term has, as written, times the sum of what is left; as written where there are none. + /// + private static Entity WithTheFactorsEveryTermSharesOut(Entity sum, Entity.Variable x) + { + var terms = Sumf.LinearChildren(sum).Select(term => Mulf.LinearChildren(term).ToList()).ToList(); + if (terms.Count < 2) + return sum; + var shared = new List(); + foreach (var factor in terms[0]) + if (factor.ContainsNode(x) && terms.Skip(1).All(other => other.Contains(factor))) + { + shared.Add(factor); + foreach (var other in terms.Skip(1)) + other.Remove(factor); + } + if (shared.Count == 0) + return sum; + var firstLeft = new List(terms[0]); + foreach (var factor in shared) + firstLeft.Remove(factor); + Entity Product(List factors) => factors.Count == 0 ? Number.Integer.One : factors.Aggregate((l, r) => l * r); + var left = new[] { Product(firstLeft) }.Concat(terms.Skip(1).Select(Product)).Aggregate((l, r) => l + r); + return Product(shared) * left; + } + /// The largest term a sum is split over as written, before the expanding gather. private const int LargestTermTakenAsWritten = 60; diff --git a/Sources/Tests/UnitTests/Calculus/ATermIntegratedByPartsTakesTheRestOfTheSumTest.cs b/Sources/Tests/UnitTests/Calculus/ATermIntegratedByPartsTakesTheRestOfTheSumTest.cs new file mode 100644 index 000000000..ce6e7350d --- /dev/null +++ b/Sources/Tests/UnitTests/Calculus/ATermIntegratedByPartsTakesTheRestOfTheSumTest.cs @@ -0,0 +1,58 @@ +// +// 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 sum whose terms are not elementary apart, integrated by parts against the power of + /// x in a term, with what parts leaves asked together with the other terms: + /// sqrt(csc(x)) (x cos(x) - 4 sec(x) tan(x)) is + /// (2 x sqrt(sin(x)) - 4 sin(x)^(3/2)/cos(x))'. Rubi's 4.7.7. + /// #718 + /// + [Trait("Area", "Calculus")] + public sealed class ATermIntegratedByPartsTakesTheRestOfTheSumTest + { + [Theory] + [InlineData("sqrt(csc(x))*(x*cos(x) - 4*sec(x)*tan(x))", 100, "0.2 0.6 0.9 2.3")] + [InlineData("x*cos(x)/sqrt(sin(x)) - 4*sqrt(sin(x))/cos(x)^2", 100, "0.2 0.6 0.9 2.3")] + [InlineData("x*sin(x)/sqrt(cos(x)) + 4*sqrt(cos(x))/sin(x)^2", 100, "0.2 0.6 0.9 1.3")] + [InlineData("sqrt(csc(2*x + 1))*(x*cos(2*x + 1) - 2*sec(2*x + 1)*tan(2*x + 1))", 150, "0.2 0.35 0.6 0.9")] + [InlineData("x*cos(x)*e^sin(x) + e^sin(x)", 50, "0.35 1.1 2.3 4.4")] + [InlineData("x^2*cos(x)*e^sin(x) + 2*x*e^sin(x)", 50, "0.35 1.1 2.3 4.4")] + public void DifferentiatesBack(string integrand, int atMost, string points) + { + var integral = integrand.ToEntity().Integrate("x"); + var text = integral.Stringize(); + Assert.DoesNotContain("integral(", text); + Assert.True(text.Length < atMost, $"{text.Length} characters of answer for {integrand}"); + var derivative = integral.Substitute("C", 0).Differentiate("x"); + var original = integrand.ToEntity(); + foreach (var at in points.Split(' ')) + { + var want = original.Substitute("x", at.ToEntity()).EvalNumerical(); + Assert.True(Math.Abs((double)want.ImaginaryPart) < 1e-12, $"{integrand} is not real at x = {at}"); + var got = derivative.Substitute("x", at.ToEntity()).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}"); + } + } + + /// + /// What parts leaves is -2 sqrt(sin(x)), and beside 3 sqrt(sin(x)) it is + /// sqrt(sin(x)), which is not elementary: declined, not answered. + /// + [Fact] + public void ARemainderThatIsNotElementaryIsDeclined() + => Assert.Contains("integral(", "x*cos(x)/sqrt(sin(x)) + 3*sqrt(sin(x))".ToEntity().Integrate("x").Stringize()); + } +}