diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index 5f7c251b9..6e0870f08 100644 --- a/BREAKING-CHANGES.md +++ b/BREAKING-CHANGES.md @@ -101,6 +101,25 @@ with the slope in another place, `F^(a + b/(c + d x)^2)/(c + d x)^2` in 15 ms wh | `"(k*g+f*g*x)^(p-1)/ln(d*(k+f*x)^p)".ToEntity().Integrate("x")` | `integral(...)` | 150 characters, `Ei(ln(d (k + f x)^p))` | | `"F^(a+b/(c+d*x)^2)/(c+d*x)^2".ToEntity().Integrate("x")` | `integral(...)` | `-F^a sqrt(pi) erf(sqrt(-b ln F)/(c + d x))/(2 d sqrt(-b ln F))` | +### An even quotient over a power of a biquadratic beside another factor is split in `x^2` + +**Answers where there were none.** `1/((x^2 - 1)(c + d x^2 + f x^4)^2)` ran past the budget, and so did Rubi's 6.1.7 +`csch(c + d x)/(a - b sinh(c + d x)^4)^2`, which is that shape under `w = cosh(c + d x)`, and 1.2.2.3's +`1/((d + e x^2)(a + b x^2 + c x^4)^2)` took eleven seconds and 10.8 GB. Split in x, the series at the +roots of `x^2 - 1` came to coefficients of a hundred thousand nodes. An even rational function with a symbol in it, over a +power of a biquadratic of three terms beside another factor in `x^2`, is now split in `w = x^2`, and each term is asked in +x. While the split is solved, each factor is made monic and its compound coefficients are named. Every other factor +stands once, and not where the numerator vanishes at a linear's root: under the tangent `sinh(u)^4/(a - b sinh(u)^4)^2` +shares a power of `1 - w` with its numerator, and is left to the rules that answered it. Checked by its derivative +([#718](https://github.com/asc-community/AngouriMath/issues/718)). + +| Input | Was (2.5.0) | Now | +|---|---|---| +| `"1/((x^2-1)*(c+d*x^2+f*x^4)^2)".ToEntity().Integrate("x")` | `integral(...)` | 2,164 characters | +| `"x^2/((x^2+1)*(c+d*x^2+f*x^4)^2)".ToEntity().Integrate("x")` | `integral(...)` | 2,062 characters, where before this change master gave 13,714 in 3.8 s | +| `"csch(c+d*x)/(a-b*sinh(c+d*x)^4)^2".ToEntity().Integrate("x")` | `integral(...)` | 2,511 characters | +| `"1/((d+k*x^2)*(a+b*x^2+c*x^4)^2)".ToEntity().Integrate("x")` | `integral(...)` | 9,968 characters, in 0.9 s, where before this change master gave 18,154 in 11 s and 10.8 GB | + ### A polynomial over a symbolic quartic that is two quadratics in a hidden variable is split over them **Answers where there were none.** `1/(a + b x + c x^2 + b x^3 + a x^4)` was declined, and so were Rubi's 1.3.1 diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs index 50e63673b..4277f4aa4 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs @@ -4631,6 +4631,117 @@ over is Entity.Powf(var @base, var power) ? return total / (constant * MathS.Pow(l.First, k)); } + /// + /// An even rational function with a symbol in it, over a power past the first of a biquadratic of + /// three terms beside another factor in x^2, split over its written factors in w = x^2 and each + /// term asked in x. 1/((x^2 - 1)(c + d x^2 + f x^4)^2) is + /// A/(w - 1) + (B + C w)/Q(w) + (D + E w)/Q(w)^2, and every term is a shape the rules + /// answer at once; split in x, the series at the roots of x^2 - 1 ran to coefficients + /// of a hundred thousand nodes, with powers like c^40 f^58, and past the budget. That + /// is Rubi's 6.1.7 csch(u)/(a - b sinh(u)^4)^2 under w = cosh(u). Checked by + /// its derivative. + /// https://github.com/asc-community/AngouriMath/issues/718 + /// + internal static Entity? SolveAnEvenQuotientOverAPowerOfABiquadraticBesideAnotherFactorInTheSquare(Entity expr, Entity.Variable x, bool integrateByParts) + { + if (!TryReadAsQuotient(expr, out var numerator, out var denominator) || !(numerator + denominator).Vars.Any(v => v != x)) + return null; + var w = Variable.CreateUnique(expr, "w"); + // Each coefficient that is more than a symbol or a number a fresh symbol while the split is + // solved and the terms integrated, and put back after: in `a - b (x^2 - 1)^2` they are + // `a - b`, `2 b` and `-b`, and the split took 3 s and 1.4 GB where with three symbols it takes 0.1. + var named = new List<(Variable Name, Entity Value)>(); + Entity everything = expr + w; + Entity Named(Entity coefficient) + { + if (coefficient is Variable or Number) + return coefficient; + foreach (var (already, value) in named) + if (value == coefficient) + return already; + var name = Variable.CreateUnique(everything, "k_q"); + everything += name; + named.Add((name, coefficient)); + return name; + } + // A polynomial in x with only even powers, written in w, its degree in w, how many terms it + // has and its leading coefficient; monic, where asked, with that coefficient out in front and + // the others over it in lowest terms, so that the names come back as `(b - a)/b`. + (Entity InW, int Degree, int Terms, Entity Leading)? InTheSquare(Entity polynomial, bool monic) + { + if (!TreeAnalyzer.TryGetPolynomial(polynomial, x, out var read) || read.Count == 0) + return null; + foreach (var pair in read) + if (pair.Key.Sign < 0 || !pair.Key.IsEven || !pair.Key.CanFitInInt32() || pair.Value.ContainsNode(x)) + return null; + var top = read.Keys.Max()!; + var leading = monic ? read[top] : Number.Integer.One; + Entity inW = Number.Integer.Zero; + foreach (var (power, coefficient) in read.OrderBy(pair => pair.Key).Select(pair => (pair.Key, pair.Value))) + { + var half = power.ToInt32Unchecked() / 2; + // Folded first, so that a number the reading left as `3 (-2)` is not a name. + var inNames = power == top && monic ? Number.Integer.One + : Named(monic ? Functions.PartialFractions.InLowestTermsOverTheSymbols(coefficient / leading) : coefficient.InnerSimplified); + inW += half == 0 ? inNames : inNames * (half == 1 ? w : MathS.Pow(w, half)); + } + return (inW, top.ToInt32Unchecked() / 2, read.Count, leading); + } + Entity denominatorInW = Number.Integer.One; + Entity leadingOfTheDenominator = Number.Integer.One; + var linearsInW = new List(); + var factorsInX = 0; + var aRepeatedBiquadratic = false; + foreach (var factor in Mulf.LinearChildren(denominator)) + { + if (!factor.ContainsNode(x)) + { + denominatorInW *= factor; + continue; + } + var (@base, power) = factor is Powf(var raised, Number.Integer { EInteger: var n }) && n.Sign > 0 && n.CanFitInInt32() + ? (raised, n.ToInt32Unchecked()) : (factor, 1); + if (InTheSquare(@base, monic: true) is not { } read || read.Degree is not (1 or 2)) + return null; + factorsInX++; + // With the middle term: over a power of `a + b x^4` the split in x answers at once. + var repeated = read.Degree == 2 && read.Terms == 3 && power >= 2; + aRepeatedBiquadratic |= repeated; + // Every other factor once: under the tangent Rubi's `sinh(u)^4/(a - b sinh(u)^4)^2` is + // `w^2 (w - 1)^4/((1 - w)^3 Q^2)`, which shares the power of `1 - w` with its numerator, + // and split as written it was 1.46 million nodes after twenty-two seconds. + if (!repeated && power != 1) + return null; + if (read.Degree == 1) + linearsInW.Add(read.InW); + denominatorInW *= power == 1 ? read.InW : MathS.Pow(read.InW, power); + leadingOfTheDenominator *= power == 1 ? read.Leading : MathS.Pow(read.Leading, power); + } + if (factorsInX < 2 || !aRepeatedBiquadratic || InTheSquare(numerator, monic: false) is not { } above) + return null; + // Nor where the numerator vanishes at the root of one of the linears, which it then shares. + foreach (var linear in linearsInW) + if (TreeAnalyzer.TryGetPolyLinear(linear, w, out var slope, out var intercept) + && Functions.PartialFractions.IsZeroAsAValue(above.InW.Substitute(w, -intercept / slope))) + return null; + if (!Functions.PartialFractions.TrySplitOverWrittenFactors(above.InW, denominatorInW, w, out var decomposition, aRepeatedQuadraticIsRead: true) + || decomposition.Complexity > 20_000) + return null; + // The split is a sum over a constant: each term of the sum asked in x, as a question of its own. + var (terms, overAll) = decomposition is Divf(var sum, var by) && !by.ContainsNode(w) ? (sum, by) : (decomposition, Number.Integer.One); + Entity total = Number.Integer.Zero; + foreach (var term in Sumf.LinearChildren(terms)) + { + if (Integration.ComputeAsAQuestionOfItsOwn(term.Substitute(w, MathS.Pow(x, 2)), x, integrateByParts) is not { } part) + return null; + total += part; + } + total /= overAll * leadingOfTheDenominator; + foreach (var (name, value) in named) + total = total.Substitute(name, value); + return Functions.PartialFractions.DerivativeHoldsAtSampledPoints(total, expr, x) ? total : null; + } + /// /// 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 diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs b/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs index b2acbfc9d..df5b1726b 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs @@ -1423,6 +1423,8 @@ internal static bool ASumOverRootsWouldAnswer 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; + // And over a power of a biquadratic beside another factor in x^2, split in x^2 itself. + if ((answer = IndefiniteIntegralSolver.SolveAnEvenQuotientOverAPowerOfABiquadraticBesideAnotherFactorInTheSquare(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/AnEvenQuotientOverAPowerOfABiquadraticIntegralTest.cs b/Sources/Tests/UnitTests/Calculus/AnEvenQuotientOverAPowerOfABiquadraticIntegralTest.cs new file mode 100644 index 000000000..c79aaf126 --- /dev/null +++ b/Sources/Tests/UnitTests/Calculus/AnEvenQuotientOverAPowerOfABiquadraticIntegralTest.cs @@ -0,0 +1,55 @@ +// +// 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 System.Threading; +using System.Threading.Tasks; +using AngouriMath.Extensions; +using Xunit; + +namespace AngouriMath.Tests.Calculus +{ + /// + /// An even rational function over a power of a biquadratic of three terms beside another factor + /// in x^2, split in w = x^2: 1/((x^2 - 1)(c + d x^2 + f x^4)^2) ran past the budget, + /// the series at the roots of x^2 - 1 coming to coefficients of a hundred thousand nodes, + /// and so did Rubi's 6.1.7 csch(u)/(a - b sinh(u)^4)^2, which is it under w = cosh(u). + /// #718 + /// + [Trait("Area", "Calculus")] + public sealed class AnEvenQuotientOverAPowerOfABiquadraticIntegralTest + { + [Theory] + [InlineData("1/((x^2 - 1)*(c + d*x^2 + f*x^4)^2)")] + [InlineData("x^2/((x^2 + 1)*(c + d*x^2 + f*x^4)^2)")] + // Rubi 6.1.7:285. + [InlineData("csch(c + d*x)/(a - b*sinh(c + d*x)^4)^2")] + public void IsAnsweredWithinTheGuard(string integrand) + { + Entity? integral = null; + // A thread of its own, so that the guard is not spent waiting behind other tests' work. + var task = Task.Factory.StartNew(() => integral = integrand.ToEntity().Integrate("x"), + CancellationToken.None, TaskCreationOptions.LongRunning, TaskScheduler.Default); + Assert.True(task.Wait(LimitTermination.Guard), $"{integrand} did not finish within {LimitTermination.Guard}"); + var text = integral!.Stringize(); + Assert.DoesNotContain("integral(", text); + Entity Pinned(Entity e) => e.Substitute("a", 1.3).Substitute("b", 0.7).Substitute("c", 0.4) + .Substitute("d", 0.9).Substitute("f", 0.8); + var derivative = Pinned(integral.Substitute("C", 0)).Differentiate("x"); + var original = Pinned(integrand.ToEntity()); + foreach (var at in new[] { -2.2, -0.9, 0.3, 1.4, 2.1 }) + { + var want = original.Substitute("x", at).EvalNumerical(); + var got = derivative.Substitute("x", at).EvalNumerical(); + Assert.False(got.IsNaN, $"d/dx of the antiderivative of {integrand} has no value at x = {at}, where the integrand is {want}"); + Assert.True(Math.Abs((double)(got - want).RealPart) + Math.Abs((double)(got - want).ImaginaryPart) + < 1e-9 * Math.Max(1, Math.Abs((double)want.RealPart) + Math.Abs((double)want.ImaginaryPart)), + $"d/dx of the antiderivative of {integrand} is {got} at x = {at}, where the integrand is {want}"); + } + } + } +}