Integrate the reciprocal of a binomial in a high power of a sine, cosine or hyperbolic one over its roots

BREAKING-CHANGES.md

@@ -396,6 +396,26 @@ in Rubi's 1.1.3.4 come out shorter as well, and those that took two seconds take
 | `"1/(x^8*(a + b*x^3)*(c + d*x^3))".ToEntity().Integrate("x")` | `integral(...)` | 1,217 characters |
 | `"1/(x^6*(a + b*x^4)*(c + d*x^4))".ToEntity().Integrate("x")` | `integral(...)` | 795 characters |
 
+### The reciprocal of a binomial in a high power of a sine, cosine or hyperbolic one is split over its roots
+
+**Answers where there were none.** `1/(a + b sin(x)^5)` was declined, and so were twenty more of Rubi's 4.1.7, 4.2.7,
+6.1.7 and 6.2.7. With `s` the function, `n` from 5 to 12, `q` 2 for an even `n` and 1 for an odd one, and `r_k` the
+`n/q`-th roots of `-a/b`, `1/(a + b s^n)` is `(q/(n a))` times the sum of `r_k/(r_k - s^q)`, and each term is an
+arctangent, or an inverse hyperbolic tangent, of the tangent of the argument or of its half
+([#718](https://github.com/asc-community/AngouriMath/issues/718)).
+
+| Input | Was (2.5.0) | Now |
+|---|---|---|
+| `"1/(a+b*sin(x)^5)".ToEntity().Integrate("x")` | `integral(...)` | 2,110 characters |
+| `"1/(a-b*sin(x)^6)".ToEntity().Integrate("x")` | `integral(...)` | 648 characters |
+| `"1/(a+b*sin(x)^8)".ToEntity().Integrate("x")` | `integral(...)` | 569 characters |
+| `"1/(a+b*cos(x)^5)".ToEntity().Integrate("x")` | `integral(...)` | 2,362 characters |
+| `"1/(a-b*cos(x)^8)".ToEntity().Integrate("x")` | `integral(...)` | 841 characters |
+| `"1/(a+b*sinh(x)^5)".ToEntity().Integrate("x")` | `integral(...)` | 2,678 characters |
+| `"1/(a+b*sinh(x)^6)".ToEntity().Integrate("x")` | `integral(...)` | 1,139 characters |
+| `"1/(a+b*cosh(x)^5)".ToEntity().Integrate("x")` | `integral(...)` | 3,607 characters |
+| `"1/(a-b*cosh(x)^8)".ToEntity().Integrate("x")` | `integral(...)` | 1,441 characters |
+
 ### A whole power of an inverse function of `c + d x^2` is integrated by its recurrence
 
 **Answers where there were none.** `(a + b acos(1 + d x^2))^4` was declined, and so were twenty-five more of Rubi's

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

@@ -25940,6 +25940,98 @@ static bool IsASumOfPowersOfX(Entity sum, Entity.Variable x)
             return answer.Nodes.Any(node => node == MathS.NaN) ? null : answer;
         }
 
+        /// <summary>
+        /// The reciprocal of <c>a + b s^n</c>, <c>s</c> a sine or a cosine of a linear and a symbol in
+        /// <c>a</c> or <c>b</c>, over the roots of the binomial: with <c>q = 2</c> for an even <c>n</c>
+        /// and <c>1</c> for an odd one, <c>m = n/q</c> and <c>r_k</c> the <c>m</c>-th roots of
+        /// <c>-a/b</c>, <c>1/(a + b s^n) = (1/(m a)) sum r_k/(r_k - s^q)</c>, and each term is an
+        /// arctangent of <c>tan(y)</c> or of <c>tan(y/2)</c>.
+        /// </summary>
+        /// <remarks>
+        /// Each term is written <c>F(K u)/K</c> with only <c>K^2</c> fixed, so its derivative is the
+        /// same whichever root <c>K</c> is: <c>(atan(K tan(y))/K)' = 1/(cos(y)^2 + K^2 sin(y)^2)</c>.
+        /// A root that makes <c>K^2</c> zero or infinite -- <c>1/(1 + sin(y)^5)</c> has <c>-1</c> -- is
+        /// another shape, and declines. <c>n</c> up to <c>4</c>, and numbers for <c>a</c> and
+        /// <c>b</c>, are answered by the rules before it. Rubi's 4.1.7 and 4.2.7
+        /// <c>1/(a + b sin(x)^5)</c>, <c>1/(a - b cos(x)^6)</c> and the like were declined.
+        /// https://github.com/asc-community/AngouriMath/issues/718
+        /// </remarks>
+        internal static Entity? SolveAReciprocalOfABinomialInAPowerOfASineOverItsRoots(Entity expr, Entity.Variable x)
+        {
+            if (!TryReadAsQuotient(expr, out var numerator, out var denominator) || numerator.ContainsNode(x))
+                return null;
+            // The hyperbolic sine and cosine are parsed as `(e^y -+ e^(-y))/2`, and read back from it.
+            static (Entity Argument, bool Sine, bool Hyperbolic)? Read(Entity node) => node switch
+            {
+                Sinf(var y) => (y, true, false),
+                Cosf(var y) => (y, false, false),
+                Divf(Sumf(Powf(var e1, var y), Powf(var e2, var minusY)), Number.Integer { EInteger: var two })
+                    when two.Equals(EInteger.FromInt32(2)) && e1 == MathS.e && e2 == MathS.e
+                    && Functions.PartialFractions.IsZeroAsAValue((y + minusY).InnerSimplified) => (y, false, true),
+                Divf(Minusf(Powf(var e1, var y), Powf(var e2, var minusY)), Number.Integer { EInteger: var two })
+                    when two.Equals(EInteger.FromInt32(2)) && e1 == MathS.e && e2 == MathS.e
+                    && Functions.PartialFractions.IsZeroAsAValue((y + minusY).InnerSimplified) => (y, true, true),
+                _ => null
+            };
+            var functions = denominator.Nodes.Where(node => node.ContainsNode(x) && Read(node) is not null).Distinct().ToList();
+            if (functions.Count != 1)
+                return null;
+            var function = functions[0];
+            var (argument, sine, hyperbolic) = Read(function)!.Value;
+            if (!TreeAnalyzer.TryGetPolyLinear(argument, x, out var slope, out _) || Functions.PartialFractions.IsZeroAsAValue(slope))
+                return null;
+            var s = Variable.CreateUnique(expr, "s_bin");
+            var inS = denominator.Substitute(function, s);
+            if (inS.ContainsNode(x) || !TreeAnalyzer.TryGetPolynomial(inS, s, out var read) || read.Count != 2
+                || !read.ContainsKey(EInteger.Zero) || read.Keys.Max() is not { } top || !top.CanFitInInt32())
+                return null;
+            var n = top.ToInt32Unchecked();
+            // Past what the rules before it answer: they write `1/(a + b cosh(x)^3)` through `e^x` in
+            // 655 characters, where this writes 1491.
+            if (n < 5 || n > 12)
+                return null;
+            var a = read[EInteger.Zero];
+            var b = read[top];
+            if (!(a + b).Vars.Any() || Functions.PartialFractions.IsZeroAsAValue(a) || Functions.PartialFractions.IsZeroAsAValue(b))
+                return null;
+            var q = n % 2 == 0 ? 2 : 1;
+            var m = n / q;
+            var half = q == 2 ? argument : argument / 2;
+            var u = hyperbolic ? MathS.Hyperbolic.Tanh(half) : MathS.Tan(half);
+            var radius = MathS.Pow(-a / b, Number.Rational.Create(1, m));
+            Entity sum = Number.Integer.Zero;
+            for (var k = 0; k < m; k++)
+            {
+                // The root of unity as a cosine and a sine: `(-1)^(2/3)` is 1 here, the real cube
+                // root of 1, where it means the principal `e^(2 pi i/3)` to Rubi.
+                var turn = 2 * MathS.pi * k / m;
+                var r = (k == 0 ? radius : radius * (MathS.Cos(turn) + MathS.i * MathS.Sin(turn))).InnerSimplified;
+                // K^2, the factor in front, and the argument of the arctangent over K.
+                var (square, front, shift) = (q, sine, hyperbolic) switch
+                {
+                    (2, true, false) => ((r - 1) / r, (Entity)Number.Integer.One, (Entity)Number.Integer.Zero),
+                    (2, false, _) => (r / (r - 1), r / (r - 1), Number.Integer.Zero),
+                    (1, true, false) => ((MathS.Sqr(r) - 1) / MathS.Sqr(r), 2, -1 / r),
+                    (2, true, true) => ((r + 1) / r, Number.Integer.One, Number.Integer.Zero),
+                    (1, true, true) => (1 + 1 / MathS.Sqr(r), 2, 1 / r),
+                    _ => ((r - 1) / (r + 1), 2 * r / (r + 1), Number.Integer.Zero),
+                };
+                square = square.InnerSimplified;
+                if (Functions.PartialFractions.IsZeroAsAValue(square) || Functions.PartialFractions.IsZeroAsAValue((1 / square).InnerSimplified))
+                    return null;
+                var root = MathS.Sqrt(square);
+                // An even power under tan(y), where K multiplies u; an odd one under tan(y/2), where it divides.
+                // The hyperbolic ones are the same with the inverse hyperbolic tangent.
+                Entity Inverse(Entity z) => hyperbolic ? MathS.Hyperbolic.Artanh(z) : MathS.Arctan(z);
+                sum += q == 2 ? front * Inverse(root * u) / root : front * Inverse((u + shift) / root) / root;
+            }
+            var answer = (numerator / (m * a * slope) * sum).InnerSimplified;
+            // Refused where it differs, not where it cannot be told: the intervals do not settle
+            // the arctangents of complex arguments at every sampled point, and `1/(a + b sin(x)^5)`,
+            // right at every point checked exactly, was declined for too few comparisons.
+            return Functions.PartialFractions.DerivativeDiffersWhereReal(answer, expr, x) ? 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

@@ -968,6 +968,8 @@ private static Entity WithItsPiecewiseSumsCombined(Entity expr)
             // among the coefficients, split over the quadratic's two roots: under the
             // half-angle it is a quartic in t with a symbol in every coefficient.
             if ((answer = IndefiniteIntegralSolver.SolveARationalFunctionOfASineOverASymbolicQuadratic(expr, x, integrateByParts)) is { }) return answer;
+            // And over a binomial in a higher power of it, over the binomial's roots.
+            if ((answer = IndefiniteIntegralSolver.SolveAReciprocalOfABinomialInAPowerOfASineOverItsRoots(expr, x)) is { }) return answer;
             // Fractional powers of `a ± a sin(y)` beside anything rational in the sine and
             // 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.

Sources/Tests/UnitTests/Calculus/AReciprocalOfABinomialInAPowerOfASineIntegralTest.cs

@@ -0,0 +1,65 @@
+//
+// 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>
+    /// The reciprocal of <c>a + b s^n</c>, <c>s</c> a sine or a cosine and <c>n</c> past 4, or a
+    /// hyperbolic sine or cosine, over the roots of the binomial: Rubi's 4.1.7, 4.2.7, 6.1.7 and 6.2.7. Pinned with <c>-a/b</c> of either sign, so the
+    /// real root among them is positive or negative.
+    /// <a href="https://github.com/asc-community/AngouriMath/issues/718">#718</a>
+    /// </summary>
+    [Trait("Area", "Calculus")]
+    public sealed class AReciprocalOfABinomialInAPowerOfASineIntegralTest
+    {
+        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("1/(a+b*sin(x)^5)")]
+        [InlineData("1/(a-b*sin(x)^6)")]
+        [InlineData("1/(a+b*sin(x)^8)")]
+        [InlineData("1/(a-b*cos(x)^5)")]
+        [InlineData("1/(a+b*cos(x)^6)")]
+        [InlineData("3/(a+b*sin(2*x+1)^7)")]
+        [InlineData("1/(a+b*cosh(x)^7)")]
+        [InlineData("1/(a-b*cosh(x)^5)")]
+        [InlineData("1/(a+b*sinh(x)^6)")]
+        [InlineData("1/(a+b*sinh(2*x+1)^5)")]
+        [InlineData("1/(a-b*cosh(x)^8)")]
+        public void IsAnsweredAndDifferentiatesBack(string integrand)
+        {
+            DifferentiatesBack(integrand, ("a", 0.6), ("b", 1.25));
+            DifferentiatesBack(integrand, ("a", -0.75), ("b", 0.4));
+        }
+    }
+}