Symbolic powers of a ± a sin beside a power of the cosine are integrated through the sine

BREAKING-CHANGES.md

@@ -886,6 +886,26 @@ either now; and the secant quotient was simplified to `cos/(b^(5/2) d)` through
 `(1/cos)^(3/2) = cos^(-3/2)`, which is wrong for every negative cosine and cancelled against the
 same wrong step on the denominator. On the 1774-problem independent suites this is 1707 to 1705, the two `asin` rows, with the sign error gone.
 
+### Symbolic powers of `a ± a sin` beside a power of the cosine are integrated through the sine
+
+**Answers where there were none.** `(a + a sin(h + f x))^m sqrt(c - c sin(h + f x))` and its kin,
+Rubi's `(a + b sin)^m (c + d sin)^n` files with `a^2 = b^2`, `c^2 = d^2` and a symbolic exponent,
+were declined or past the budget. `(1 + sin(y))(1 - sin(y))` is `cos(y)^2`, so under `u = sin(y)`
+each such power, a power of `g cos(y)` beside them and `dy` itself are powers of `1 + u` and
+`1 - u`, up to a factor constant on each interval, and the rest is asked in `u`; the cosine's the
+same way. The answer is the antiderivative in `u` times that factor, the powers as written over
+the form they were rewritten to. Where the question in `u` is `(1 + u)^A (1 - u)^B` with `A + B`
+a whole number below `-2` it is declined still
+([#718](https://github.com/asc-community/AngouriMath/issues/718)).
+
+| Input | Was (2.5.0) | Now |
+|---|---|---|
+| `"(a + a*sin(h + f*x))^m*sqrt(c - c*sin(h + f*x))".ToEntity().Integrate("x")` | `integral(...)` | the integrand times `(1 + sin(h + f x))/(f (m + 1/2) cos(h + f x))`, written unreduced |
+| `"cos(h + f*x)^2*(a + a*sin(h + f*x))^m/sqrt(c - c*sin(h + f*x))".ToEntity().Integrate("x")` | `integral(...)` | the integrand times `(1 + sin(h + f x))/(f (m + 3/2) cos(h + f x))`, written unreduced |
+| `"(g*cos(h + f*x))^(1 - 2*m)*(a + a*sin(h + f*x))^m*(c - c*sin(h + f*x))^(m - 1)".ToEntity().Integrate("x")` | `integral(...)` | the integrand times `-(1 - sin(h + f x)) ln(1 - sin(h + f x))/(f cos(h + f x))` |
+| `"(a + a*sin(h + f*x))^m*(c - c*sin(h + f*x))^(-1 - m)".ToEntity().Integrate("x")` | `integral(...)` | the integrand times `(1 + sin(h + f x))(1 - sin(h + f x))/(f (1 + 2 m) cos(h + f x))`, written unreduced |
+| `"(a + a*cos(h + f*x))^m*sqrt(c - c*cos(h + f*x))".ToEntity().Integrate("x")` | `integral(...)` | the integrand times `-(1 + cos(h + f x))/(f (m + 1/2) sin(h + f x))`, written unreduced |
+
 ### A power of x comes out of a fractional power of a sum whose every term has it
 
 `(a x^j + b x^n)^p` with a fractional `p`, the power of x common to every term inside the power,

Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs

@@ -3600,6 +3600,159 @@ Entity Step(ERational power)
             return answer.Nodes.Any(node => node == MathS.NaN) ? null : answer;
         }
 
+        /// <summary>
+        /// Powers of <c>a ± a sin(y)</c> with a symbolic exponent, beside a power of
+        /// <c>g cos(y)</c> and anything else in the sine alone, under <c>u = sin(y)</c>:
+        /// <c>(1 + sin(y))(1 - sin(y))</c> is <c>cos(y)^2</c>, so each factor is a power of
+        /// <c>1 + u</c> or of <c>1 - u</c> up to a factor constant on each interval where it is
+        /// defined, and <c>dy = du/cos(y)</c> is one more pair of half powers. The cosine's the
+        /// same way, by <c>u = cos(y)</c> beside a power of <c>g sin(y)</c>.
+        /// </summary>
+        /// <remarks>
+        /// <para>
+        /// Rubi's <c>(a + b sin)^m (c + d sin)^n</c> files with <c>a^2 = b^2</c> and
+        /// <c>c^2 = d^2</c> hold about a hundred and twenty problems with a symbolic exponent,
+        /// <c>(a + a sin(e + f x))^m sqrt(c - c sin(e + f x))</c> the plainest, and none was
+        /// answered: the half angle at which <c>1 ± sin(y)</c> is a square, the rule before this
+        /// one, wants numeric powers, and the substitution search spent the budget on them. What
+        /// is asked in <c>u</c> is <c>(1 + u)^A (1 - u)^B R(u)</c>, which is elementary where one
+        /// of <c>A</c>, <c>B</c> and <c>A + B</c> is a whole number and the rest is a
+        /// polynomial, and is declined where it is not.
+        /// </para>
+        /// <para>
+        /// The answer is <c>K G(sin(y))/F</c>, for <c>G</c> the antiderivative in <c>u</c>,
+        /// <c>F</c> the slope of <c>y</c> and <c>K</c> the powers as written over the form they
+        /// were rewritten to and over <c>cos(y)</c>: a quotient whose logarithmic derivative is
+        /// zero, so the answer holds whatever <c>a</c>, <c>g</c> and the exponents are --
+        /// <c>(a (1 + u))^m</c> is not <c>a^m (1 + u)^m</c> for every <c>a</c>, and the quotient
+        /// of the two is constant on each interval either way. At the question asked only, and
+        /// only where an exponent is symbolic: numeric half powers are the half angle's.
+        /// https://github.com/asc-community/AngouriMath/issues/718
+        /// </para>
+        /// </remarks>
+        internal static Entity? SolveSymbolicPowersOfOnePlusMinusASineThroughTheSine(Entity expr, Entity.Variable x, bool integrateByParts)
+        {
+            if (!Integration.AnsweringTheQuestionAsked)
+                return null;
+            Entity? argument = null;
+            foreach (var node in expr.Nodes)
+            {
+                if (TrigonometricArgument(node) is not { } thisArgument || !thisArgument.ContainsNode(x))
+                    continue;
+                if (argument is null)
+                    argument = thisArgument;
+                else if (argument != thisArgument)
+                    return null;
+            }
+            if (argument is null || !TreeAnalyzer.TryGetPolyLinear(argument, x, out var rate, out _)
+                || rate.ContainsNode(x) || TreeAnalyzer.IsZero(rate))
+                return null;
+            var sine = MathS.Sin(argument);
+            var cosine = MathS.Cos(argument);
+            // Each factor as a base and an exponent, a power of a power read as one power: the
+            // quotient of the two spellings is constant on each interval, and K carries it.
+            var factors = new List<(Entity Factor, Entity Base, Entity Exponent)>();
+            foreach (var factor in Mulf.LinearChildren(expr))
+            {
+                Entity @base = factor, exponent = Number.Integer.One;
+                while (@base is Powf(var inner, var power))
+                {
+                    @base = inner;
+                    exponent = power * exponent;
+                }
+                if (exponent.ContainsNode(x))
+                    return null;
+                factors.Add((factor, @base, exponent.InnerSimplified));
+            }
+            // Which function the sums are of: all of one kind.
+            bool? ofTheSine = null;
+            foreach (var (_, @base, _) in factors)
+                if (@base.ContainsNode(x) && ReadAsOnePlusMinusAFunction(@base, sine, cosine) is var (_, _, isSine))
+                {
+                    if (ofTheSine is { } kind && kind != isSine)
+                        return null;
+                    ofTheSine = isSine;
+                }
+            if (ofTheSine is not { } sineKind)
+                return null;
+            var function = sineKind ? sine : cosine;
+            var complement = sineKind ? cosine : sine;
+            var u = Variable.CreateUnique(expr, "u_one_plus_minus");
+            // The exponents of 1 + f, of 1 - f and of the complement, the factors that carry
+            // them as they stand, and the rest in u.
+            Entity plus = Number.Integer.Zero, minus = Number.Integer.Zero, ofTheComplement = Number.Integer.Zero;
+            Entity asItIs = Number.Integer.One, rest = Number.Integer.One;
+            var symbolic = false;
+            foreach (var (factor, @base, exponent) in factors)
+            {
+                if (!factor.ContainsNode(x))
+                {
+                    rest *= factor;
+                    continue;
+                }
+                if (ReadAsOnePlusMinusAFunction(@base, sine, cosine) is var (_, isPlus, _))
+                {
+                    if (isPlus)
+                        plus += exponent;
+                    else
+                        minus += exponent;
+                    asItIs *= factor;
+                    symbolic |= exponent is not Number;
+                    continue;
+                }
+                if (AsAPowerOfTheComplement(@base) is { } sign)
+                {
+                    ofTheComplement += sign * exponent;
+                    asItIs *= factor;
+                    symbolic |= exponent is not Number;
+                    continue;
+                }
+                var inU = factor.Substitute(function, u);
+                if (inU.ContainsNode(x))
+                    return null;
+                rest *= inU;
+            }
+            if (!symbolic)
+                return null;
+            // dy = du/complement takes one more half power of each.
+            var half = Number.Rational.Create(1, 2);
+            var a = (plus + ofTheComplement * half - half).Simplify();
+            var b = (minus + ofTheComplement * half - half).Simplify();
+            // A zero exponent leaves its factor out: `(1 - u)^0` is one provided `u` is not one, and
+            // the condition would stand between the question and the rules that answer it.
+            Entity PowerOf(Entity @base, Entity exponent) => exponent == Number.Integer.Zero ? Number.Integer.One : MathS.Pow(@base, exponent);
+            var integrand = (PowerOf(1 + u, a) * PowerOf(1 - u, b) * rest).InnerSimplified;
+            if (Integration.ComputeAsAQuestionOfItsOwn(integrand, u, integrateByParts) is not { } result)
+                return null;
+            // d sin(y) = cos(y) dy and d cos(y) = -sin(y) dy.
+            Entity constant = asItIs * PowerOf(1 + function, (-a).InnerSimplified) * PowerOf(1 - function, (-b).InnerSimplified) / (rate * complement);
+            if (!sineKind)
+                constant = -constant;
+            var answer = constant * result.Substitute(u, function);
+            return answer.Nodes.Any(node => node == MathS.NaN) ? null : answer;
+
+            // 1 for a constant times the complement, -1 for a constant times its reciprocal
+            // function, null for anything else.
+            int? AsAPowerOfTheComplement(Entity @base)
+            {
+                int? found = null;
+                foreach (var part in Mulf.LinearChildren(@base))
+                {
+                    if (!part.ContainsNode(x))
+                        continue;
+                    if (found is not null)
+                        return null;
+                    if (part == complement)
+                        found = 1;
+                    else if (sineKind ? part is Secantf(var s) && s == argument : part is Cosecantf(var c) && c == argument)
+                        found = -1;
+                    else
+                        return null;
+                }
+                return found;
+            }
+        }
+
         /// <summary>
         /// Half-odd powers of <c>a ± a sec(y)</c>, beside powers of <c>d sec(y)</c> or
         /// <c>d cos(y)</c> and anything rational in the sine and cosine of <c>y</c>, by the

Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs

@@ -848,6 +848,10 @@ private static Entity Normalized(Entity expr, Entity.Variable x) =>
             // cosine, by the half angle at which they are squares: `1 + sin(y)` is `2 sin(u)^2`. Before
             // the substitution search, which spent twenty seconds on the radicals of the sine.
             if ((answer = IndefiniteIntegralSolver.SolveByTheHalfAngleWhereOnePlusASineIsASquare(expr, x, integrateByParts)) is { }) return answer;
+            // And symbolic powers of `a ± a sin(y)` beside a power of `g cos(y)`, by u = sin(y), where
+            // `(1 + u)(1 - u)` is the cosine's square: the half angle wants numeric powers, and the
+            // substitution search spent the budget on these.
+            if ((answer = IndefiniteIntegralSolver.SolveSymbolicPowersOfOnePlusMinusASineThroughTheSine(expr, x, integrateByParts)) is { }) return answer;
             // And a half-odd power of a +- a sec(y), which is that square over cos(y): by the half-angle
             // tangent, in which the whole is rational beside one root.
             if ((answer = IndefiniteIntegralSolver.SolveByTheHalfAngleTangentBesideAHalfOddPowerOfOnePlusASecant(expr, x, integrateByParts)) is { }) return answer;

Sources/Tests/UnitTests/Calculus/PowersOfOnePlusAndMinusASineIntegralTest.cs

@@ -0,0 +1,50 @@
+//
+// 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>
+    /// Powers of <c>a ± a sin(y)</c> with a symbolic exponent, beside a power of <c>g cos(y)</c> and
+    /// a function of the sine, under <c>u = sin(y)</c>, where <c>(1 + u)(1 - u)</c> is the cosine's
+    /// square; and the cosine's the same way. Rubi's <c>(a + b sin)^m (c + d sin)^n</c> files.
+    /// <a href="https://github.com/asc-community/AngouriMath/issues/718">#718</a>
+    /// </summary>
+    [Trait("Area", "Calculus")]
+    public sealed class PowersOfOnePlusAndMinusASineIntegralTest
+    {
+        [Theory]
+        [InlineData("(a + a*sin(h + f*x))^m*sqrt(c - c*sin(h + f*x))")]
+        [InlineData("(a + a*sin(h + f*x))^m*(c - c*sin(h + f*x))^(5/2)*(p + q*sin(h + f*x)^2)")]
+        [InlineData("cos(h + f*x)^2*(a + a*sin(h + f*x))^m/sqrt(c - c*sin(h + f*x))")]
+        [InlineData("(g*cos(h + f*x))^(1 - 2*m)*(a + a*sin(h + f*x))^m*(c - c*sin(h + f*x))^(m - 1)")]
+        [InlineData("(a + a*sin(h + f*x))^m*(c - c*sin(h + f*x))^(-1 - m)")]
+        [InlineData("(a - a*sin(h + f*x))^m*(c + c*sin(h + f*x))^n*(b*(m - n) + b*(1 + m + n)*sin(h + f*x))")]
+        [InlineData("(a + a*cos(h + f*x))^m*sqrt(c - c*cos(h + f*x))")]
+        public void ThroughTheSine(string integrand)
+        {
+            var integral = integrand.ToEntity().Integrate("x");
+            Assert.DoesNotContain("integral(", integral.Stringize());
+            Entity Pinned(Entity e) => e.Substitute("a", 1.3).Substitute("b", 0.9).Substitute("c", 0.7)
+                .Substitute("h", 0.4).Substitute("f", 1.3).Substitute("g", 0.8).Substitute("m", 0.37)
+                .Substitute("n", 1.21).Substitute("p", 1.1).Substitute("q", 0.6);
+            var derivative = Pinned(integral.Substitute("C", 0)).Differentiate("x");
+            var original = Pinned(integrand.ToEntity());
+            foreach (var at in new[] { -2.1, -0.9, 0.3, 0.7, 1.6, 2.2 })
+            {
+                var want = original.Substitute("x", at).EvalNumerical();
+                var got = derivative.Substitute("x", at).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) + Math.Abs((double)want.ImaginaryPart)),
+                    $"d/dx of the antiderivative of {integrand} is {got} at x = {at}, where the integrand is {want}");
+            }
+        }
+    }
+}