Integrate a power of a quotient of linears beside whole powers of them

BREAKING-CHANGES.md

@@ -154,6 +154,18 @@ integrand is then `x^(2n + 1) (c + d x)^n (2 c + 3 d x)`, the derivative of `x^(
 | `"x^(2*n)*(c+d*x)^n*(2*c*x+3*d*x^2)".ToEntity().Integrate("x")` | `integral(...)` | 61 characters |
 | `"x^(n-1)*(c+d*x^2)^n*(c*x+3*d*x^3)".ToEntity().Integrate("x")` | `integral(...)` | 65 characters |
 
+### A power of a quotient of linears beside whole powers of them is integrated under the quotient
+
+**Answers where there were none.** `e^(n acoth(a x))/(c - a^2 c x^2)^2` was declined, and so were its `^3` and `^4`, Rubi's
+7.4.2. The parser folds `e^(n acoth(a x))` into `((a x + 1)/(a x - 1))^(n/2)`, and under `t = (a x + 1)/(a x - 1)` -- the
+Moebius substitution, `p + q x = t (r + s x)` and `dx = (q r - p s)/(q - s t)^2 dt` -- a whole power of the quadratic
+is a polynomial in `t` times a power of it ([#233](https://github.com/asc-community/AngouriMath/issues/233)).
+
+| Input | Was (2.5.0) | Now |
+|---|---|---|
+| `"e^(n*acoth(a*x))/(c-a^2*c*x^2)^2".ToEntity().Integrate("x")` | `integral(...)` | three powers of `(a x + 1)/(a x - 1)`, `n/2 - 1` to `n/2 + 1` |
+| `"e^(n*acoth(a*x))/(c-a^2*c*x^2)^4".ToEntity().Integrate("x")` | `integral(...)` | seven powers of it, `n/2 - 3` to `n/2 + 3` |
+
 ### The last digits of an evaluated sine, cosine, tangent, arcsine or arctangent
 
 **Silent, and in the hundredth digit.** `EvalNumerical` at the default hundred digits gives the

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

@@ -709,6 +709,109 @@ private static bool IsADerivativeOfAFactorBelowTheBar(Entity expr, Entity.Variab
             return (constant * (MathS.Pow(x, raised) * logarithm - 2 * slope * integral) / raised).InnerSimplified;
         }
 
+        /// <summary>
+        /// A power of a quotient of two linears beside whole powers of the two, under the Moebius
+        /// substitution <c>t = (p + q x)/(r + s x)</c>: then <c>p + q x = t (r + s x)</c>,
+        /// <c>r + s x = D/(q - s t)</c> and <c>dx = D/(q - s t)^2 dt</c> with <c>D = q r - p s</c>, so
+        /// <c>t^alpha (p + q x)^j (r + s x)^k dx</c> is <c>D^(j + k + 1) t^(alpha + j) (q - s t)^(-(j + k + 2)) dt</c>,
+        /// a polynomial in <c>t</c> times a power of it wherever <c>j + k &lt;= -2</c>. A whole power of
+        /// a quadratic that is a multiple of the two linears' product counts as a power of each.
+        /// <c>e^(n acoth(a x))</c> arrives folded as <c>((a x + 1)/(a x - 1))^(n/2)</c>, and over
+        /// <c>(c - a^2 c x^2)^m</c> it is Rubi's 7.4.2.
+        /// https://github.com/asc-community/AngouriMath/issues/233
+        /// </summary>
+        /// <remarks>
+        /// Only for a power <c>alpha</c> that is not whole: a whole one makes a rational function,
+        /// which the partial fractions answer. <c>t^gamma</c> is differentiated to
+        /// <c>gamma t^(gamma - 1) t'</c> on the principal branch wherever <c>t</c> is not 0, the
+        /// negative half-line included, so the answer written back in the quotient holds where
+        /// the integrand is defined. A half-odd power of the quadratic is left alone: it is not
+        /// the product of the linears' powers on every side of their roots.
+        /// </remarks>
+        internal static Entity? SolveByTheMoebiusSubstitutionOfAQuotientOfLinears(Entity expr, Entity.Variable x)
+        {
+            Powf? power = null;
+            Entity constant = Number.Integer.One;
+            var others = new List<(Entity Base, EInteger Power)>();
+            foreach (var (factor, underneath) in FactorsOfTheIntegrand(expr))
+            {
+                if (!factor.ContainsNode(x))
+                {
+                    constant = underneath ? constant / factor : constant * factor;
+                    continue;
+                }
+                if (!underneath && power is null && factor is Powf(Divf(var above, var below), var alpha)
+                    && !alpha.ContainsNode(x) && alpha.Evaled is not Number.Integer
+                    && above.ContainsNode(x) && below.ContainsNode(x))
+                {
+                    power = (Powf)factor;
+                    continue;
+                }
+                var (@base, whole) = factor is Powf(var b, Number.Integer n) && n.EInteger.CanFitInInt32() ? (b, n.EInteger) : (factor, EInteger.One);
+                others.Add((@base, underneath ? -whole : whole));
+            }
+            if (power is null || others.Count == 0 || power.Base is not Divf(var top, var bottom)
+                || !TreeAnalyzer.TryGetPolyLinear(top, x, out var q, out var p) || TreeAnalyzer.IsZero(q)
+                || !TreeAnalyzer.TryGetPolyLinear(bottom, x, out var s, out var r) || TreeAnalyzer.IsZero(s))
+                return null;
+            var determinant = (q * r - p * s).InnerSimplified;
+            if (VanishesIdentically(determinant))
+                return null;
+            // Each factor a multiple of the top, of the bottom, or of their product.
+            var j = EInteger.Zero;
+            var k = EInteger.Zero;
+            foreach (var (@base, n) in others)
+            {
+                Entity multiple;
+                if (TreeAnalyzer.TryGetPolyLinear(@base, x, out var slope, out var intercept))
+                {
+                    if (VanishesIdentically(slope * p - intercept * q))
+                    {
+                        multiple = (slope / q).InnerSimplified;
+                        j += n;
+                    }
+                    else if (VanishesIdentically(slope * r - intercept * s))
+                    {
+                        multiple = (slope / s).InnerSimplified;
+                        k += n;
+                    }
+                    else
+                        return null;
+                }
+                else if (TreeAnalyzer.TryGetPolyQuadratic(@base, x, out var a2, out var b2, out var c2) && !TreeAnalyzer.IsZero(a2))
+                {
+                    multiple = (a2 / (q * s)).InnerSimplified;
+                    if (!VanishesIdentically(b2 - multiple * (p * s + q * r)) || !VanishesIdentically(c2 - multiple * p * r))
+                        return null;
+                    j += n;
+                    k += n;
+                }
+                else
+                    return null;
+                constant *= n.Equals(EInteger.One) ? multiple : MathS.Pow(multiple, Number.Integer.Create(n));
+            }
+            var degree = (j + k + EInteger.FromInt32(2)).Negate();
+            if (degree.Sign < 0 || degree.CompareTo(EInteger.FromInt32(24)) > 0 || !j.CanFitInInt32())
+                return null;
+            // D^(j + k + 1) t^(alpha + j) (q - s t)^degree, term by term in t.
+            var m = degree.ToInt32Checked();
+            Entity sum = Number.Integer.Zero;
+            Entity binomial = Number.Integer.One;
+            for (var i = 0; i <= m; i++)
+            {
+                var offset = j + EInteger.FromInt32(i + 1);
+                var raised = (offset.IsZero ? power.Exponent : power.Exponent + Number.Integer.Create(offset)).InnerSimplified;
+                if (raised.Evaled is Number.Complex { IsZero: true })
+                    return null;
+                // Free of x and small, so simplified: the multiples are otherwise left as
+                // `(1/(a (-a)))^(-2) (a - (-a))^(-3)`.
+                var coefficient = (binomial * MathS.Pow(q, Number.Integer.Create(m - i)) * MathS.Pow(-s, Number.Integer.Create(i)) / raised).Simplify();
+                sum += coefficient * MathS.Pow(power.Base, raised);
+                binomial = binomial * (m - i) / (i + 1);
+            }
+            return ((constant * MathS.Pow(determinant, Number.Integer.Create(j + k + EInteger.One))).Simplify() * sum).InnerSimplified;
+        }
+
         /// <summary>
         /// A quotient by the square of a base, by parts against the base's reciprocal, where the
         /// base's derivative makes the rest exact: <c>N/(M v^2)</c> is <c>g v'/v^2</c> with

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

@@ -784,6 +784,9 @@ private static Entity WithItsPiecewiseSumsCombined(Entity expr)
             if ((answer = IndefiniteIntegralSolver.SolveAPowerOverAPowerOfALogarithm(expr, x)) is { }) return answer;
             // And x^m artanh(A x/sqrt(d + A^2 x^2)), by parts, the logarithm's derivative read as the root it is.
             if ((answer = IndefiniteIntegralSolver.SolveByPartsAgainstAnArtanhWhoseDerivativeIsARoot(expr, x, integrateByParts)) is { }) return answer;
+            // A power of a quotient of linears beside whole powers of them, which is how e^(n acoth(a x))
+            // over a power of c - a^2 c x^2 arrives folded, under t = the quotient.
+            if ((answer = IndefiniteIntegralSolver.SolveByTheMoebiusSubstitutionOfAQuotientOfLinears(expr, x)) is { }) return answer;
             // A power of an exponential with a positive base is the exponential of the product,
             // exactly, and only that spelling is one the exponential rules read.
             if ((answer = IndefiniteIntegralSolver.SolveByFlatteningAPowerOfAnExponential(expr, x, integrateByParts)) is { }) return answer;

Sources/Tests/UnitTests/Calculus/APowerOfAQuotientOfLinearsBesideThemIntegralTest.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 a quotient of two linears beside whole powers of them, under
+    /// <c>t = (p + q x)/(r + s x)</c>: <c>e^(n acoth(a x))/(c - a^2 c x^2)^m</c> arrives as
+    /// <c>((a x + 1)/(a x - 1))^(n/2)</c> over the power of the quadratic, and is a polynomial in
+    /// <c>t</c> times a power of it. For <c>m</c> past 1 it was declined. Rubi's 7.4.2.
+    /// <a href="https://github.com/asc-community/AngouriMath/issues/233">#233</a>
+    /// </summary>
+    [Trait("Area", "Calculus")]
+    public sealed class APowerOfAQuotientOfLinearsBesideThemIntegralTest
+    {
+        [Theory]
+        [InlineData("e^(n*acoth(a*x))/(c - a^2*c*x^2)^2")]
+        [InlineData("e^(n*acoth(a*x))/(c - a^2*c*x^2)^3")]
+        [InlineData("e^(n*atanh(a*x))/(c - a^2*c*x^2)^2")]
+        [InlineData("((a*x + 1)/(a*x - 1))^(n/2)/((a*x + 1)^2*(a*x - 1)^2)")]
+        [InlineData("((x + 2)/(3*x - 1))^(1/3)/((x + 2)*(3*x - 1)^2)")]
+        public void UnderTheQuotient(string integrand)
+        {
+            var integral = integrand.ToEntity().Integrate("x");
+            var text = integral.Stringize();
+            Assert.DoesNotContain("integral(", text);
+            Entity Pinned(Entity e) => e.Substitute("a", 1.3).Substitute("c", 0.7).Substitute("n", 0.9);
+            var derivative = Pinned(integral.Substitute("C", 0)).Differentiate("x");
+            var original = Pinned(integrand.ToEntity());
+            var compared = 0;
+            // On both sides of the roots, -0.77 and 0.77 for a = 1.3, and of 1/3 and -2.
+            foreach (var at in new[] { -2.9, -1.4, -0.4, 0.2, 0.6, 1.3, 2.9 })
+            {
+                var want = original.Substitute("x", at).EvalNumerical();
+                var got = derivative.Substitute("x", at).EvalNumerical();
+                if (want.IsNaN)
+                    continue;
+                compared++;
+                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}");
+            }
+            Assert.True(compared >= 5, $"only {compared} points could be compared for {integrand}");
+        }
+    }
+}