diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index 5d01b9202..e89c879dd 100644 --- a/BREAKING-CHANGES.md +++ b/BREAKING-CHANGES.md @@ -244,6 +244,21 @@ cache is now keyed on the downcasting setting too |---|---|---| | `MathS.FromString("1/3").Evaled` under downcasting off, after the same parse with it on | `1/3` | `0.3333…` (one hundred digits) | +### A polynomial in a sine over a power of a binomial in it is written in powers of the binomial + +**Answers where there were none.** `(A + B cos(c + d x) + C cos(c + d x)^2)/(a + b cos(c + d x))^4` ran out of time, and +so did twelve more of Rubi's 4.1.3.1, 4.2.3.1 and 4.2.4.2. `P(f)/(a + b f)^n`, for `f` a sine or cosine, is now +written in powers of `u = a + b f`, with coefficients `c_j = sum_i p_i C(i, j) (-a)^(i - j)/b^i`, and each power +`u^(j - n)` is integrated on its own +([#718](https://github.com/asc-community/AngouriMath/issues/718)). + +| Input | Was (2.5.0) | Now | +|---|---|---| +| `"(A+B*cos(c+d*x)+C*cos(c+d*x)^2)/(a+b*cos(c+d*x))^4".ToEntity().Integrate("x")` | `integral(...)` | 9,003 characters | +| `"(A+B*sin(e+f*x))*(c+d*sin(e+f*x))^3/(a+a*sin(e+f*x))^3".ToEntity().Integrate("x")` | `integral(...)` | 1,217 characters | +| `"(a+b*sin(e+f*x))^2*(A+B*sin(e+f*x))/(c+d*sin(e+f*x))^2".ToEntity().Integrate("x")` | `integral(...)` | 3,337 characters | +| `"cos(c+d*x)^4*(A+C*cos(c+d*x)^2)/(a+b*cos(c+d*x))^4".ToEntity().Integrate("x")` | `integral(...)` | 10,859 characters | + ### A difference of a negated term is printed as a sum **A printed form.** A solution and a difference were printed with the sign of a negative factor left inside: diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs index 079a84bee..b838c5e66 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs @@ -26205,6 +26205,87 @@ when two.Equals(EInteger.FromInt32(2)) && e1 == MathS.e && e2 == MathS.e return Functions.PartialFractions.DerivativeDiffersWhereReal(answer, expr, x) ? null : answer; } + /// + /// A polynomial in a sine or a cosine of a linear over a whole power of a binomial in it, + /// (A + B cos(y) + C cos(y)^2)/(a + b cos(y))^4, written in powers of the binomial: + /// with u = a + b f, f = (u - a)/b, and P(f)/u^n = sum c_j u^(j - n) with + /// c_j = sum_(i >= j) p_i C(i, j) (-a)^(i - j)/b^i, each power of the binomial its own + /// question. + /// + /// + /// The rules after it took the whole quotient through the half-angle and wrote 30 million + /// characters for (A + B cos(x) + C cos(x)^2)/(a + b cos(x))^4 in eighteen seconds, where + /// 1/(a + b cos(x))^4 alone is 44 thousand. Rubi's 4.1.4 and 4.2.4 + /// (a + b cos)^m (A + B cos + C cos^2) at m = -4 timed out. + /// https://github.com/asc-community/AngouriMath/issues/718 + /// + internal static Entity? SolveAPolynomialInASineOverAPowerOfABinomialInItInPowersOfTheBinomial(Entity expr, Entity.Variable x) + { + if (!Integration.AnsweringTheQuestionAsked) + return null; + var functions = expr.Nodes.Where(node => node is Sinf or Cosf && node.ContainsNode(x)).Distinct().ToList(); + if (functions.Count != 1) + return null; + var function = functions[0]; + var argument = function switch { Sinf(var y) => y, Cosf(var y) => y, _ => null! }; + if (!TreeAnalyzer.TryGetPolyLinear(argument, x, out var slope, out _) || Functions.PartialFractions.IsZeroAsAValue(slope)) + return null; + var s = Variable.CreateUnique(expr, "s_pow"); + var inS = expr.Substitute(function, s); + if (inS.ContainsNode(x)) + return null; + var (above, below) = Functions.SingleQuotient.Of(inS); + // Below the bar a factor free of s times one whole power of a binomial in s. + static IEnumerable Factors(Entity e) => e is Mulf(var l, var r) ? Factors(l).Concat(Factors(r)) : new[] { e }; + Entity constant = Number.Integer.One; + Entity? binomial = null; + var n = 0; + foreach (var factor in Factors(below)) + { + if (!factor.ContainsNode(s)) + { + constant *= factor; + continue; + } + if (binomial is not null) + return null; + (binomial, n) = factor is Powf(var b0, Number.Integer { EInteger: var power }) && power.CanFitInInt32() ? (b0, power.ToInt32Unchecked()) : (factor, 1); + } + if (binomial is null || n < 2 || n > 12 + || !TreeAnalyzer.TryGetPolyLinear(binomial, s, out var b, out var a) + || Functions.PartialFractions.IsZeroAsAValue(a) || Functions.PartialFractions.IsZeroAsAValue(b) + || !TreeAnalyzer.TryGetPolynomial(above, s, out var p) || p.Keys.Any(key => key.Sign < 0 || !key.CanFitInInt32())) + return null; + var degree = p.Keys.Max()!.ToInt32Unchecked(); + // The reciprocal of the power alone is what the rules after it answer, in a few milliseconds. + if (degree == 0 || degree > 12) + return null; + var u = a + b * function; + Entity answer = Number.Integer.Zero; + for (var j = 0; j <= degree; j++) + { + Entity c = Number.Integer.Zero; + var choose = EInteger.One; + for (var i = j; i <= degree; i++) + { + // C(i, j) from C(i - 1, j). + if (i > j) + choose = choose * i / (i - j); + if (p.TryGetValue(EInteger.FromInt32(i), out var pi)) + c += pi * Number.Integer.Create(choose) * (i == j ? Number.Integer.One : MathS.Pow(-a, Number.Integer.Create(i - j))) / (i == 0 ? Number.Integer.One : MathS.Pow(b, Number.Integer.Create(i))); + } + c = c.InnerSimplified; + if (Functions.PartialFractions.IsZeroAsAValue(c)) + continue; + var piece = j - n == 0 ? (Entity)Number.Integer.One : MathS.Pow(u, Number.Integer.Create(j - n)); + if (Integration.ComputeAsAQuestionOfItsOwn(piece, x, integrateByParts: true) is not { } integrated) + return null; + answer += c * integrated; + } + answer = (answer / constant).InnerSimplified; + return answer.Nodes.Any(node => node == MathS.NaN) ? null : answer; + } + /// /// A rational function of a sine, or of a cosine, over a quadratic in it with a symbol /// among the coefficients: sin(x)/(a + b sin(x) + c sin(x)^2), and the same with diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs b/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs index 382365c1e..9b1b2855b 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs @@ -680,6 +680,9 @@ private static Entity WithItsPiecewiseSumsCombined(Entity expr) // A power of a sine, cosine or hyperbolic one over a binomial in a power of it, over the // binomial's roots, before the rules that spent the budget on `sinh(y)^4/(a + b sinh(y)^3)`. if ((answer = IndefiniteIntegralSolver.SolveAPowerOverABinomialInAPowerOfASineOverItsRoots(expr, x)) is { }) return answer; + // A polynomial in a sine or a cosine over a whole power of a binomial in it, in powers of the + // binomial, before the half-angle wrote 30 million characters for `(A + B cos + C cos^2)/(a + b cos)^4`. + if ((answer = IndefiniteIntegralSolver.SolveAPolynomialInASineOverAPowerOfABinomialInItInPowersOfTheBinomial(expr, x)) is { }) return answer; // A rational function times the logarithm of a nested product of powers of linears, by parts, // its logarithm's derivative read off as written: `ln(e (f (a + b x)^p (c + d x)^q)^r)`. if ((answer = IndefiniteIntegralSolver.SolveARationalFunctionTimesALogarithmOfPowersOfLinears(expr, x, integrateByParts)) is { }) return answer; diff --git a/Sources/Tests/UnitTests/Calculus/APolynomialInASineOverAPowerOfABinomialInItIntegralTest.cs b/Sources/Tests/UnitTests/Calculus/APolynomialInASineOverAPowerOfABinomialInItIntegralTest.cs new file mode 100644 index 000000000..60fb55b54 --- /dev/null +++ b/Sources/Tests/UnitTests/Calculus/APolynomialInASineOverAPowerOfABinomialInItIntegralTest.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 polynomial in a sine or a cosine over a whole power of a binomial in it, written in powers of + /// the binomial: Rubi's 4.1.4 and 4.2.4 (a + b cos)^m (A + B cos + C cos^2) at m = -4, + /// which timed out. Pinned on the real domain. + /// #718 + /// + [Trait("Area", "Calculus")] + public sealed class APolynomialInASineOverAPowerOfABinomialInItIntegralTest + { + 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}"); + } + } + + // The parameters are p, q, r rather than Rubi's A, B, C: the helper takes C for the constant. + [Theory] + [InlineData("(p+q*cos(c+d*x)+r*cos(c+d*x)^2)/(a+b*cos(c+d*x))^4")] + [InlineData("cos(c+d*x)^3*(p+q*cos(c+d*x)+r*cos(c+d*x)^2)/(a+b*cos(c+d*x))^4")] + [InlineData("cos(c+d*x)^2*(p+r*cos(c+d*x)^2)/(a+b*cos(c+d*x))^4")] + [InlineData("(p+q*sin(c+d*x)+r*sin(c+d*x)^2)/(a+b*sin(c+d*x))^3")] + public void IsAnsweredAndDifferentiatesBack(string integrand) + { + DifferentiatesBack(integrand, ("a", 1.3), ("b", 0.7), ("c", 0.7), ("d", 1.07), ("p", 0.43), ("q", -0.45), ("r", 0.22)); + DifferentiatesBack(integrand, ("a", 2), ("b", -1.5), ("c", -0.3), ("d", 0.6), ("p", -1), ("q", 2), ("r", 0.5)); + } + } +}