A power of a cosine and a sine plus their amplitude is integrated

BREAKING-CHANGES.md

@@ -620,6 +620,22 @@ are right on both sides of zero ([#718](https://github.com/asc-community/Angouri
 | `"x^2/(1 + x^4)^(3/4)".ToEntity().Integrate("x")` | `integral(...)` | logarithms and an arctangent of `(1 + x^4)^(1/4)/x` |
 | `"x^6*(3 + 4*x^4)^(1/4)".ToEntity().Integrate("x")` | `integral(...)` | the same in `(3 + 4 x^4)^(1/4)/x` |
 
+### A power of a cosine and a sine plus their amplitude is integrated
+
+**Answers where there were none.** `sqrt(5 + 4 cos(x) + 3 sin(x))` was declined, with the rest of
+Rubi's 4.7.7 powers of `a + b cos(y) + c sin(y)` with `a^2 = b^2 + c^2`, several after searches past
+the budget. `b cos(y) + c sin(y)` is `R cos(y - phi)` for `R = sqrt(b^2 + c^2)`, so the base is
+`a (1 ± cos(y - phi))`, a square of the half angle; a half-odd power of it, or a negative whole one,
+is integrated in closed form now, through `T = b sin(y) - c cos(y)`, with no `phi` in the answer
+([#718](https://github.com/asc-community/AngouriMath/issues/718)).
+
+| Input | Was (2.5.0) | Now |
+|---|---|---|
+| `"sqrt(5 + 4*cos(x) + 3*sin(x))".ToEntity().Integrate("x")` | `integral(...)` | `2 (4 sin(x) - 3 cos(x))/sqrt(5 + 4 cos(x) + 3 sin(x))` |
+| `"1/sqrt(5 + 4*cos(x) + 3*sin(x))".ToEntity().Integrate("x")` | `integral(...)` | `ln((1 + q)/(1 - q))/sqrt(10)` for `q = (4 sin(x) - 3 cos(x))/(sqrt(10) sqrt(5 + 4 cos(x) + 3 sin(x)))` |
+| `"(-5 + 4*cos(x) + 3*sin(x))^(3/2)".ToEntity().Integrate("x")` | `integral(...)` | `T sqrt(S)/(3/2) - (40/3) T/sqrt(S)` for `T = 4 sin(x) - 3 cos(x)` and `S` the base, imaginary as the integrand is |
+| `"1/(b*cos(g + f*x) + c*sin(g + f*x) - sqrt(b^2 + c^2))^3".ToEntity().Integrate("x")` | `integral(...)` | `b sin(g + f x) - c cos(g + f x)` times powers of the base |
+
 ### A rational function of the sine and cosine with symbols in it is integrated by the half angle
 
 `sin(x)^2/(a + b cos(x))` was left unevaluated while `sin(x)^2/(2 + 3 cos(x))` was answered. Under

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

@@ -3550,6 +3550,148 @@ Entity Step(ERational power)
             return Step(n).InnerSimplified;
         }
 
+        /// <summary>
+        /// <c>(a + b cos(y) + c sin(y))^n</c> for <c>a^2 = b^2 + c^2</c>, <c>y</c> linear in
+        /// <c>x</c> and <c>n</c> half-odd or a negative whole number, in closed form.
+        /// </summary>
+        /// <remarks>
+        /// <para>
+        /// <c>b cos(y) + c sin(y)</c> is <c>R cos(y - phi)</c> for <c>R = sqrt(b^2 + c^2)</c>, so
+        /// the base is <c>a (1 ± cos(y - phi))</c>, a square of the half angle as
+        /// <c>a + a sin(y)</c> is, and <see cref="SolveAHalfPowerOfOnePlusASine"/>'s identity
+        /// holds for it. Rubi's <c>sqrt(5 + 4 cos(x) + 3 sin(x))</c> and the rest of its 4.7.7
+        /// with <c>a^2 = b^2 + c^2</c> were declined or past the budget. The answer needs no
+        /// <c>phi</c>: with <c>S</c> the base and <c>T = b sin(y) - c cos(y)</c>,
+        /// <c>T' = S - a</c>, <c>S' = -T</c> and <c>T^2 = S (2a - S)</c>, the last of which is
+        /// <c>a^2 = b^2 + c^2</c>, and from them
+        /// </para>
+        /// <code>
+        ///     d/dy (T S^(n-1)) = n S^n - a (2n - 1) S^(n-1)
+        /// </code>
+        /// <para>
+        /// which steps <c>n</c> down to <c>1/2</c>, where the integral is <c>2T/sqrt(S)</c>, or
+        /// up to <c>-1/2</c> or <c>-1</c>, where it is
+        /// <c>(2/sqrt(2a)) atanh(T/(sqrt(2a) sqrt(S)))</c> and <c>T/(a S)</c>. Each base is
+        /// checked by differentiating it with those three identities alone.
+        /// </para>
+        /// <para>
+        /// For a real <c>a</c> of either sign: the derivations use only
+        /// <c>S^(3/2) = S sqrt(S)</c> and <c>(sqrt(2a) sqrt(S))^2 = 2a S</c>, which hold for
+        /// the principal powers, and the argument of the hyperbolic arctangent stays between
+        /// <c>-1</c> and <c>1</c>, one less its square being <c>S/(2a)</c>, which is not negative
+        /// since <c>S</c> has the sign of <c>a</c>. An antiderivative on every interval between
+        /// the zeros of <c>S</c>, where <c>T</c> changes sign. <c>a^2 = b^2 + c^2</c> is decided,
+        /// not assumed, and the cosine and the sine are both there: one alone is
+        /// <see cref="SolveAHalfPowerOfOnePlusASine"/>'s.
+        /// https://github.com/asc-community/AngouriMath/issues/718
+        /// </para>
+        /// </remarks>
+        internal static Entity? SolveAPowerOfACosineAndASinePlusTheirAmplitude(Entity expr, Entity.Variable x)
+        {
+            // A constant times one power of the base, the power above the bar or below it.
+            if (!TryReadAsQuotient(expr, out var above, out var below))
+                (above, below) = (expr, Number.Integer.One);
+            Entity constant = Number.Integer.One;
+            Entity? @base = null;
+            ERational? exponent = null;
+            foreach (var (side, sign) in new[] { (above, 1), (below, -1) })
+                foreach (var factor in Mulf.LinearChildren(side))
+                {
+                    if (!factor.ContainsNode(x))
+                    {
+                        constant = sign > 0 ? constant * factor : constant / factor;
+                        continue;
+                    }
+                    if (@base is not null)
+                        return null;
+                    (@base, exponent) = factor is Powf(var inner, var power) && power.Evaled is Number.Rational rational
+                        ? (inner, rational.ERational)
+                        : (factor, ERational.One);
+                    if (sign < 0)
+                        exponent = exponent.Negate();
+                }
+            if (@base is null || exponent is null)
+                return null;
+            var n = exponent;
+            // Half-odd, or a negative whole number, and of a modest size: a positive whole power
+            // is a polynomial in the sine and cosine.
+            var isHalfOdd = n.Denominator.Equals(EInteger.FromInt32(2));
+            if (!isHalfOdd && !(n.Denominator.Equals(EInteger.One) && n.Sign < 0)
+                || n.Abs().CompareTo(ERational.FromInt32(MaximumAmplitudePower)) > 0)
+                return null;
+            // a + b cos(y) + c sin(y), read off the sum.
+            Entity a = Number.Integer.Zero;
+            Entity? b = null, c = null, argument = null;
+            foreach (var term in Sumf.LinearChildren(@base))
+            {
+                if (!term.ContainsNode(x))
+                {
+                    a = a == Number.Integer.Zero ? term : a + term;
+                    continue;
+                }
+                Entity coefficient = Number.Integer.One;
+                Entity? function = null;
+                foreach (var factor in Mulf.LinearChildren(term))
+                    if (!factor.ContainsNode(x))
+                        coefficient = coefficient == Number.Integer.One ? factor : coefficient * factor;
+                    else if (function is null && factor is Sinf or Cosf)
+                        function = factor;
+                    else
+                        return null;
+                var inner = function?.DirectChildren.First();
+                if (inner is null || argument is not null && inner != argument)
+                    return null;
+                argument = inner;
+                if (function is Cosf)
+                {
+                    if (b is not null) return null;
+                    b = coefficient;
+                }
+                else
+                {
+                    if (c is not null) return null;
+                    c = coefficient;
+                }
+            }
+            if (b is null || c is null || argument is null || a == Number.Integer.Zero)
+                return null;
+            if (!TreeAnalyzer.TryGetPolyLinear(argument, x, out var rate, out _) || rate.ContainsNode(x)
+                || rate.Evaled is Number.Complex { IsZero: true })
+                return null;
+            // a^2 = b^2 + c^2, decided rather than assumed.
+            var excess = (MathS.Sqr(a) - MathS.Sqr(b) - MathS.Sqr(c)).InnerSimplified;
+            if (excess.Evaled is not Number.Complex { IsZero: true }
+                && !(excess.Vars.Any() && Functions.PartialFractions.Bare(excess.Simplify()).Evaled is Number.Complex { IsZero: true }))
+                return null;
+            if (a.Evaled is Number.Complex { IsZero: true })
+                return null;
+
+            var s = @base;
+            var t = b * MathS.Sin(argument) - c * MathS.Cos(argument);
+            var half = ERational.Create(1, 2);
+            Entity Integral(ERational power)
+            {
+                if (power.CompareTo(half) == 0)
+                    return 2 * t / MathS.Sqrt(s);
+                if (power.CompareTo(half.Negate()) == 0)
+                    return 2 / MathS.Sqrt(2 * a) * MathS.Hyperbolic.Artanh(t / (MathS.Sqrt(2 * a) * MathS.Sqrt(s)));
+                if (power.CompareTo(ERational.FromInt32(-1)) == 0)
+                    return t / (a * s);
+                var current = Number.Rational.Create(power);
+                if (power.Sign > 0)
+                    return t * MathS.Pow(s, Number.Rational.Create(power.Subtract(ERational.One))) / current
+                        + a * (2 * current - 1) / current * Integral(power.Subtract(ERational.One));
+                return ((current + 1) * Integral(power.Add(ERational.One)) - t * MathS.Pow(s, current)) / (a * (2 * current + 1));
+            }
+            return (constant * Integral(n) / rate).InnerSimplified;
+        }
+
+        /// <summary>
+        /// The largest power, either way, that <see cref="SolveAPowerOfACosineAndASinePlusTheirAmplitude"/>
+        /// steps through: each step is a term of the answer.
+        /// </summary>
+        private const int MaximumAmplitudePower = 8;
+
         /// <summary>
         /// Fractional powers of <c>a ± a sin(y)</c>, or of <c>a ± a cos(y)</c>, beside anything
         /// rational in the sine and cosine of <c>y</c>, by the half angle at which they are

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

@@ -858,6 +858,9 @@ 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 a power of `a + b cos(y) + c sin(y)` with `a^2 = b^2 + c^2`, the same square
+            // turned by a phase, in closed form: the substitution search spent the budget on it.
+            if ((answer = IndefiniteIntegralSolver.SolveAPowerOfACosineAndASinePlusTheirAmplitude(expr, x)) 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.

Sources/Tests/UnitTests/Calculus/CosineAndSinePlusTheirAmplitudeIntegralTest.cs

@@ -0,0 +1,54 @@
+//
+// 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 power of <c>a + b cos(y) + c sin(y)</c> with <c>a^2 = b^2 + c^2</c>, which is
+    /// <c>a (1 ± cos(y - phi))</c>, a square of the half angle: half-odd powers either way and
+    /// negative whole ones, in closed form through <c>T = b sin(y) - c cos(y)</c>. Rubi's 4.7.7.
+    /// <a href="https://github.com/asc-community/AngouriMath/issues/718">#718</a>
+    /// </summary>
+    /// <remarks>
+    /// Compared as complex numbers on both sides of zero and away from the zeros of the base:
+    /// with <c>a</c> negative the base is not positive anywhere, and the integrand is imaginary
+    /// on the whole line.
+    /// </remarks>
+    [Trait("Area", "Calculus")]
+    public sealed class CosineAndSinePlusTheirAmplitudeIntegralTest
+    {
+        [Theory]
+        [InlineData("sqrt(5 + 4*cos(x) + 3*sin(x))")]
+        [InlineData("(5 + 4*cos(x) + 3*sin(x))^(5/2)")]
+        [InlineData("1/sqrt(5 + 4*cos(x) + 3*sin(x))")]
+        [InlineData("1/(5 + 4*cos(x) + 3*sin(x))^(3/2)")]
+        [InlineData("(-5 + 4*cos(x) + 3*sin(x))^(3/2)")]
+        [InlineData("1/(b*cos(g + f*x) + c*sin(g + f*x) + sqrt(b^2 + c^2))^(5/2)")]
+        [InlineData("(b*cos(g + f*x) + c*sin(g + f*x) - sqrt(b^2 + c^2))^(3/2)")]
+        [InlineData("1/(b*cos(g + f*x) + c*sin(g + f*x) - sqrt(b^2 + c^2))^3")]
+        public void IsWrittenThroughTheDerivativeOfItsBase(string integrand)
+        {
+            var integral = integrand.ToEntity().Integrate("x");
+            Assert.DoesNotContain("integral(", integral.Stringize());
+            Entity Pinned(Entity e) => e.Substitute("b", 0.9).Substitute("c", 0.7).Substitute("g", 0.4).Substitute("f", 1.3);
+            var derivative = Pinned(integral.Substitute("C", 0)).Differentiate("x");
+            var original = Pinned(integrand.ToEntity());
+            foreach (var at in new[] { -1.9, -1.1, -0.4, 0.4, 1.1, 1.9 })
+            {
+                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}");
+            }
+        }
+    }
+}