Integrate a product of powers of a sine and a cosine beside a non-polynomial by linearising it

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

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;
         }
 
+        /// <summary>
+        /// A product of whole powers of a sine and a cosine of one linear times something else in x
+        /// that is no polynomial, <c>sin(u)^3 cos(u)^6/(a + b u)</c>, linearised: the product is a sum
+        /// of <c>cos(k u)</c> and <c>sin(k u)</c> by the product-to-sum formulas, and each term beside
+        /// the rest its own question.
+        /// </summary>
+        /// <remarks>
+        /// <c>sin(3 u)/(a + b u)</c> is answered in sines and cosine integrals, and the product of
+        /// powers was not: Rubi's 5.1.4 <c>x^3 (1 - c^2 x^2)^(5/2)/(a + b asin(c x))</c>, which under
+        /// <c>u = asin(c x)</c> 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
+        /// </remarks>
+        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<int, ERational> { [0] = ERational.One };
+            var sine = new Dictionary<int, ERational>();
+            static void Add(Dictionary<int, ERational> 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<int, ERational>(), new Dictionary<int, ERational>());
+                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;
+        }
+
         /// <summary>
         /// A rational function of a sine, or of a cosine, over a quadratic in it with a symbol
         /// among the coefficients: <c>sin(x)/(a + b sin(x) + c sin(x)^2)</c>, and the same with

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.

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
+{
+    /// <summary>
+    /// A product of powers of a sine and a cosine beside something that is no polynomial, linearised:
+    /// Rubi's 5.1.4 <c>x^3 (1 - c^2 x^2)^(5/2)/(a + b asin(c x))</c>, which was declined.
+    /// <a href="https://github.com/asc-community/AngouriMath/issues/718">#718</a>
+    /// </summary>
+    [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));
+        }
+    }
+}