diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index 66bfaa7fd..441e6c0df 100644 --- a/BREAKING-CHANGES.md +++ b/BREAKING-CHANGES.md @@ -281,6 +281,19 @@ divided by Smith's method now, which scales by the larger part and squares nothi |---|---|---| | `(1.98e119 + i)/(3.39e125 + i)` with the downcasting off | `0` | `5.8407…e-7` | +### A product of powers of a sine and a cosine beside a non-polynomial is linearised + +**Answers where there were none.** `x^3 (1 - c^2 x^2)^(5/2)/(a + b asin(c x))` was declined, and its square ran out of +time (Rubi's 5.1.4). Under `x = sin(u)/c` they become `sin(u)^p cos(u)^q` beside a power of `a + b u`. A product of +powers of a sine and a cosine beside such a factor is now written as a sum of sines and cosines of multiples of the +angle, each of which the integrator answers +([#718](https://github.com/asc-community/AngouriMath/issues/718)). + +| Input | Was (2.5.0) | Now | +|---|---|---| +| `"x^3*(1-c^2*x^2)^(5/2)/(a+b*asin(c*x))".ToEntity().Integrate("x")` | `integral(...)` | 508 characters | +| `"x^3*(1-c^2*x^2)^(5/2)/(a+b*asin(c*x))^2".ToEntity().Integrate("x")` | `integral(...)` | 1,470 characters | + ### A piecewise is defined where the case it takes is **Domain conditions that were false.** The domain condition of a piecewise asked every case's predicate diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs index 5e97ce0ba..7ef7cd080 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs @@ -26347,6 +26347,104 @@ when two.Equals(EInteger.FromInt32(2)) && e1 == MathS.e && e2 == MathS.e return answer.Nodes.Any(node => node == MathS.NaN) ? null : answer; } + /// + /// A product of whole powers of a sine and a cosine of one linear times something else in x + /// that is no polynomial, sin(u)^3 cos(u)^6/(a + b u), linearised: the product is a sum + /// of cos(k u) and sin(k u) by the product-to-sum formulas, and each term beside + /// the rest its own question. + /// + /// + /// sin(3 u)/(a + b u) is answered in sines and cosine integrals, and the product of + /// powers was not: Rubi's 5.1.4 x^3 (1 - c^2 x^2)^(5/2)/(a + b asin(c x)), which under + /// u = asin(c x) is that product, was declined. Last but one, so that everything + /// answering the product as written goes first. + /// https://github.com/asc-community/AngouriMath/issues/718 + /// + internal static Entity? SolveAProductOfPowersOfASineAndACosineByLinearising(Entity expr, Entity.Variable x) + { + if (!Integration.AnsweringTheQuestionAsked) + return null; + static IEnumerable<(Entity Factor, bool Below)> Factors(Entity e, bool below) => e switch + { + Mulf(var l, var r) => Factors(l, below).Concat(Factors(r, below)), + Divf(var l, var r) => Factors(l, below).Concat(Factors(r, !below)), + _ => new[] { (e, below) } + }; + Entity? argument = null; + int sines = 0, cosines = 0; + Entity rest = Number.Integer.One; + foreach (var (factor, below) in Factors(expr, false)) + { + var (function, power) = factor is Powf(var b0, Number.Integer { EInteger: var p0 }) && b0 is Sinf or Cosf && p0.Sign > 0 && p0.CanFitInInt32() + ? (b0, p0.ToInt32Unchecked()) : (factor, 1); + if (!below && function is Sinf or Cosf && function.ContainsNode(x)) + { + var y0 = function is Sinf(var ys) ? ys : ((Cosf)function).Argument; + if (argument is not null && argument != y0) + return null; + argument = y0; + if (function is Sinf) sines += power; else cosines += power; + } + else + rest = below ? rest / factor : rest * factor; + } + if (argument is null || sines + cosines < 2 || sines + cosines > 24 + || !TreeAnalyzer.TryGetPolyLinear(argument, x, out _, out _) + || !rest.ContainsNode(x) || rest.Nodes.Any(node => node is Sinf or Cosf or Tanf or Cotanf or Secantf or Cosecantf) + || TreeAnalyzer.TryGetPolynomial(rest, x, out var asPolynomial) && asPolynomial.Keys.All(key => key.Sign >= 0)) + return null; + // cos(k y) and sin(k y) coefficients, multiplied by sin y and cos y one at a time. + var cosine = new Dictionary { [0] = ERational.One }; + var sine = new Dictionary(); + static void Add(Dictionary into, int k, ERational value) + => into[k] = (into.TryGetValue(k, out var old) ? old + value : value).ToLowestTerms(); + var half = ERational.Create(1, 2); + for (var step = 0; step < sines + cosines; step++) + { + var bySine = step < sines; + var (nextCosine, nextSine) = (new Dictionary(), new Dictionary()); + foreach (var (k, value) in cosine.Select(pair => (pair.Key, pair.Value))) + if (bySine) + { // cos(k y) sin y = (sin((k + 1) y) - sin((k - 1) y))/2 + Add(nextSine, k + 1, value * half); + if (k - 1 > 0) Add(nextSine, k - 1, -value * half); + else if (k - 1 < 0) Add(nextSine, 1 - k, value * half); + } + else + { // cos(k y) cos y = (cos((k + 1) y) + cos((k - 1) y))/2 + Add(nextCosine, k + 1, value * half); + Add(nextCosine, System.Math.Abs(k - 1), value * half); + } + foreach (var (k, value) in sine.Select(pair => (pair.Key, pair.Value))) + if (bySine) + { // sin(k y) sin y = (cos((k - 1) y) - cos((k + 1) y))/2 + Add(nextCosine, System.Math.Abs(k - 1), value * half); + Add(nextCosine, k + 1, -value * half); + } + else + { // sin(k y) cos y = (sin((k + 1) y) + sin((k - 1) y))/2 + Add(nextSine, k + 1, value * half); + if (k - 1 > 0) Add(nextSine, k - 1, value * half); + else if (k - 1 < 0) Add(nextSine, 1 - k, -value * half); + } + (cosine, sine) = (nextCosine, nextSine); + } + Entity answer = Number.Integer.Zero; + foreach (var (terms, isSine) in new[] { (cosine, false), (sine, true) }) + foreach (var (k, value) in terms.OrderBy(term => term.Key).Select(pair => (pair.Key, pair.Value))) + { + if (value.IsZero) + continue; + var multiple = k == 1 ? argument : Number.Integer.Create(k) * argument; + Entity piece = k == 0 ? rest : (isSine ? MathS.Sin(multiple) : MathS.Cos(multiple)) * rest; + if (Integration.ComputeAsAQuestionOfItsOwn(piece, x, integrateByParts: true) is not { } integrated) + return null; + answer += Number.Rational.Create(value.Numerator, value.Denominator) * integrated; + } + answer = answer.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 fd7ca2fdb..f716bb36f 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs @@ -1364,6 +1364,9 @@ private static Entity WithItsPiecewiseSumsCombined(Entity expr) if ((answer = IndefiniteIntegralSolver.SolvePowersOfOneLinearAsOnePower(expr, x, integrateByParts)) is { }) return answer; // And beside a polynomial too, where everything before has declined: `x ln(e (f (a + b x)^p (c + d x)^q)^r)`. if ((answer = IndefiniteIntegralSolver.SolveARationalFunctionTimesALogarithmOfPowersOfLinears(expr, x, integrateByParts, besideAPolynomialToo: true)) is { }) return answer; + // A product of powers of a sine and a cosine beside something that is no polynomial, + // linearised, where everything answering it as written has declined: `sin(u)^3 cos(u)^6/(a + b u)`. + if ((answer = IndefiniteIntegralSolver.SolveAProductOfPowersOfASineAndACosineByLinearising(expr, x)) is { }) return answer; // Last of the rewrites, because it is the only one that fires on an integrand nothing // is wrong with -- it clears a parameter rather than a shape -- so everything that // can answer the problem as written gets to try first. diff --git a/Sources/Tests/UnitTests/Calculus/AProductOfPowersOfASineAndACosineLinearisedIntegralTest.cs b/Sources/Tests/UnitTests/Calculus/AProductOfPowersOfASineAndACosineLinearisedIntegralTest.cs new file mode 100644 index 000000000..582df104c --- /dev/null +++ b/Sources/Tests/UnitTests/Calculus/AProductOfPowersOfASineAndACosineLinearisedIntegralTest.cs @@ -0,0 +1,57 @@ +// +// 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 product of powers of a sine and a cosine beside something that is no polynomial, linearised: + /// Rubi's 5.1.4 x^3 (1 - c^2 x^2)^(5/2)/(a + b asin(c x)), which was declined. + /// #718 + /// + [Trait("Area", "Calculus")] + public sealed class AProductOfPowersOfASineAndACosineLinearisedIntegralTest + { + 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("x^3*(1-c^2*x^2)^(5/2)/(a+b*asin(c*x))")] + [InlineData("x^3*(1-c^2*x^2)^(5/2)/(a+b*asin(c*x))^2")] + [InlineData("sin(3+2*x)^3*cos(3+2*x)^6/(a+b*x)")] + public void IsAnsweredAndDifferentiatesBack(string integrand) + { + // |c x| below 1 at every sampled x, where the inverse sine is real. + DifferentiatesBack(integrand, ("a", 1.3), ("b", 0.7), ("c", 0.3)); + DifferentiatesBack(integrand, ("a", -2), ("b", 1.5), ("c", -0.2)); + } + } +}