Integrate a power of an inverse function of c + d x^2 by its recurrence

BREAKING-CHANGES.md

@@ -396,6 +396,25 @@ 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 |
 
+### 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
+5.1.5, 5.2.5, 7.1.5 and 7.2.5. For `u` the inverse sine or cosine of `c + d x^2` with `c^2 = 1`, or its inverse
+hyperbolic cosine, `integral (a + b u)^n` is `x (a + b u)^n + 2 e n b H (a + b u)^(n - 1)` less
+`4 e b^2 n (n - 1) integral (a + b u)^(n - 2)`, where `H u' = 2` and `e` is `1` for the circular functions and `-1` for
+the hyperbolic one. It closes at the first two powers
+([#718](https://github.com/asc-community/AngouriMath/issues/718)).
+
+| Input | Was (2.5.0) | Now |
+|---|---|---|
+| `"(a+b*acos(1+d*x^2))^4".ToEntity().Integrate("x")` | `integral(...)` | 301 characters |
+| `"(a+b*acos(-1+d*x^2))^2".ToEntity().Integrate("x")` | `integral(...)` | 151 characters |
+| `"(a+b*asin(1+d*x^2))^3".ToEntity().Integrate("x")` | `integral(...)` | 241 characters |
+| `"(a+b*acosh(-1+d*x^2))^2".ToEntity().Integrate("x")` | `integral(...)` | 214 characters |
+| `"(a+b*acosh(1+d*x^2))^3".ToEntity().Integrate("x")` | `integral(...)` | 334 characters |
+| `"acos(1+d*x^2)".ToEntity().Integrate("x")` | `integral(...)` | 76 characters |
+| `"asin(-1+d*x^2)".ToEntity().Integrate("x")` | `integral(...)` | 77 characters |
+
 ### A biquadratic with surds in its coefficients is split by its roots in the square
 
 **Answers where there were none.** `(sqrt(2) - x^2)/(1 - sqrt(2) x^2 + x^4)` was declined, and so was Rubi's 1.2.2.3

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

@@ -28844,6 +28844,85 @@ static bool Through(Entity node, Entity.Variable x)
             return Integration.ComputeAsTheSameQuestion(simplified, x, integrateByParts);
         }
 
+        /// <summary>
+        /// A whole power of <c>a + b u</c>, where <c>u</c> is the inverse sine or cosine of
+        /// <c>w = c + d x^2</c> with <c>c^2 = 1</c>, the inverse hyperbolic cosine of one, or the
+        /// inverse hyperbolic sine of one with <c>c^2 = -1</c>, by the recurrence
+        /// <c>integral (a + b u)^n = x (a + b u)^n + 2 e n b H (a + b u)^(n - 1) - 4 e b^2 n (n - 1) integral (a + b u)^(n - 2)</c>,
+        /// with <c>H = 1/(d x u'(w))</c>, <c>e = 1</c> for the circular functions and <c>-1</c> for the
+        /// hyperbolic ones. It closes at <c>n = 0</c> and <c>n = 1</c>, so nothing is asked below.
+        /// </summary>
+        /// <remarks>
+        /// Each <c>u</c> has <c>u'' = e w u'^3</c>, and with <c>c^2 = 1</c> (or <c>-1</c>) the radicand of
+        /// <c>u'</c> is <c>d x^2 (2 c + d x^2)</c> up to sign, so <c>H u' = 2</c> and
+        /// <c>(x (a + b u)^n + 2 e n b H (a + b u)^(n - 1))'</c> is <c>(a + b u)^n</c> plus a multiple of
+        /// <c>(a + b u)^(n - 2)</c>. Rubi's 5.1.5, 5.2.5 and 7.2.5 <c>(a + b acos(1 + d x^2))^4</c>
+        /// and the like were declined.
+        /// https://github.com/asc-community/AngouriMath/issues/718
+        /// </remarks>
+        internal static Entity? SolveAPowerOfAnInverseFunctionOfAQuadraticByTheRecurrence(Entity expr, Entity.Variable x)
+        {
+            var (@base, power) = expr is Powf(var b0, Number.Integer p0) ? (b0, p0) : (expr, Number.Integer.One);
+            if (power.EInteger.Sign <= 0 || power.EInteger.CompareTo(EInteger.FromInt32(12)) > 0)
+                return null;
+            var n = power.EInteger.ToInt32Unchecked();
+            // The inverse function, its argument w, the sign e, the radical H is written with, and
+            // whether w^2 is 1 or -1 at x = 0.
+            (Entity Node, Entity W, int Sign, Entity Root, int Square)? Read(Entity node)
+            {
+                switch (node)
+                {
+                    case Arcsinf(var w):
+                        return (node, w, 1, MathS.Sqrt(1 - MathS.Sqr(w)), 1);
+                    case Arccosf(var w):
+                        return (node, w, 1, -MathS.Sqrt(1 - MathS.Sqr(w)), 1);
+                    // The parser writes acosh(w) and asinh(w) as ln(w + sqrt(w^2 -+ 1)).
+                    case Logf(var logBase, Sumf(var w, Powf(var radicand, Number.Rational half)))
+                        when logBase == MathS.e && half == Number.Rational.Create(1, 2) && w.ContainsNode(x):
+                        foreach (var square in new[] { 1, -1 })
+                            if (VanishesIdentically((radicand - (MathS.Sqr(w) - square)).Expand()))
+                                return (node, w, -1, MathS.Sqrt(radicand), square);
+                        return null;
+                    default:
+                        return null;
+                }
+            }
+            var candidates = @base.Nodes.Select(Read).Where(read => read is not null).Select(read => read!.Value).ToList();
+            if (candidates.Count == 0 || candidates.Any(read => read.Node != candidates[0].Node))
+                return null;
+            var (inverse, argument, sign, root, squareAtZero) = candidates[0];
+            var t = Variable.CreateUnique(expr, "u_inv");
+            var inT = @base.Substitute(inverse, t);
+            if (inT.ContainsNode(x) || !TreeAnalyzer.TryGetPolyLinear(inT, t, out var b, out var a)
+                || Functions.PartialFractions.IsZeroAsAValue(b))
+                return null;
+            if (!TreeAnalyzer.TryGetPolynomial(argument, x, out var read) || read.Keys.Any(key => !key.IsZero && !key.Equals(EInteger.FromInt32(2)))
+                || !read.TryGetValue(EInteger.FromInt32(2), out var d) || d.ContainsNode(x) || Functions.PartialFractions.IsZeroAsAValue(d))
+                return null;
+            var c = read.TryGetValue(EInteger.Zero, out var c0) ? c0 : Number.Integer.Zero;
+            if (c.ContainsNode(x) || !Functions.PartialFractions.IsZeroAsAValue((MathS.Sqr(c) - squareAtZero).Expand().InnerSimplified))
+                return null;
+            var h = root / (d * x);
+            var e = (Entity)Number.Integer.Create(sign);
+            var a_bu = a + b * inverse;
+            // I(k) by the recurrence, from I(0) = x and I(1) = x (a + b u) + 2 e b H.
+            var integrals = new Entity[n + 1];
+            for (var k = 0; k <= n; k++)
+            {
+                Entity term = k == 0 ? x : x * MathS.Pow(a_bu, k) + 2 * e * k * b * h * (k == 1 ? Number.Integer.One : MathS.Pow(a_bu, k - 1));
+                if (k >= 2)
+                    term -= 4 * e * MathS.Sqr(b) * k * (k - 1) * integrals[k - 2];
+                integrals[k] = term;
+            }
+            // The recurrence is exact where H u' = 2 with the derivative the library takes, the one
+            // step that depends on a branch; checked here rather than at sampled points, which for
+            // `c = 1` and a positive `d` put `w` past 1, on the cut of the inverse sine and cosine,
+            // and declined `asin(1 + d x^2)^2`, right for the `d < 0` where it is real.
+            if (!Functions.PartialFractions.IsZeroAsAValue((h * inverse.Differentiate(x) - 2).InnerSimplified))
+                return null;
+            return integrals[n].InnerSimplified;
+        }
+
         /// <summary>
         /// A power of a constant multiple of a radicand that a logarithm holds under a root,
         /// written over that radicand: <c>acosh(c x)</c> is <c>ln(c x + sqrt(c^2 x^2 - 1))</c>,

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

@@ -664,6 +664,9 @@ private static Entity WithItsPiecewiseSumsCombined(Entity expr)
             // A sum free of the variable that comes to a number multiplied out, as that number,
             // before anything reads a coefficient or an exponent in it: `2 + 2 a - 2 (1 + a)`.
             if ((answer = IndefiniteIntegralSolver.SolveByFoldingConstantSumsToTheirNumbers(expr, x, integrateByParts)) is { }) return answer;
+            // A whole power of a linear in the inverse sine, cosine or hyperbolic cosine of
+            // `c + d x^2` with `c^2 = 1`, by its recurrence, which closes: by parts it went round.
+            if ((answer = IndefiniteIntegralSolver.SolveAPowerOfAnInverseFunctionOfAQuadraticByTheRecurrence(expr, x)) is { }) return answer;
             // A product of powers of quotients of linears with one linear in two of them, written
             // over each linear once with the constant that takes in front. First, since it answers
             // the question asked only, and the polynomial term below writes `1/u` as `u^(-1)` and

Sources/Tests/UnitTests/Calculus/APowerOfAnInverseFunctionOfAQuadraticIntegralTest.cs

@@ -0,0 +1,68 @@
+//
+// 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 whole power of <c>a + b u</c> for <c>u</c> the inverse sine, cosine or hyperbolic cosine of
+    /// <c>c + d x^2</c> with <c>c^2 = 1</c>, or the inverse hyperbolic sine of one with <c>c^2 = -1</c>,
+    /// by the recurrence that closes at the first two powers: Rubi's 5.1.5, 5.2.5 and 7.2.5. Pinned
+    /// with <c>d</c> of either sign, on the real domain of each and off it.
+    /// <a href="https://github.com/asc-community/AngouriMath/issues/718">#718</a>
+    /// </summary>
+    [Trait("Area", "Calculus")]
+    public sealed class APowerOfAnInverseFunctionOfAQuadraticIntegralTest
+    {
+        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}");
+            }
+        }
+
+        // d with the sign that keeps the argument real at every sampled x, and |d| below 2/2.3^2
+        // so that the inverse sine and cosine stay inside [-1, 1]; off it the exact checks hold too.
+        [Theory]
+        [InlineData("(a+b*acos(1+d*x^2))^4", -1)]
+        [InlineData("(a+b*acos(-1+d*x^2))^2", 1)]
+        [InlineData("(a+b*asin(1+d*x^2))^3", -1)]
+        [InlineData("(a+b*asin(-1+d*x^2))", 1)]
+        [InlineData("(a+b*acosh(-1+d*x^2))^2", -1)]
+        [InlineData("(a+b*acosh(1+d*x^2))^3", 1)]
+        [InlineData("(a+b*asinh(i+d*x^2))^2", 1)]
+        public void IsAnsweredAndDifferentiatesBack(string integrand, int sign)
+        {
+            DifferentiatesBack(integrand, ("a", 0.6), ("b", 1.25), ("d", 0.3 * sign));
+            DifferentiatesBack(integrand, ("a", -0.75), ("b", -0.7), ("d", 0.15 * sign));
+        }
+
+        [Fact]
+        public void WithNumbersIsAnswered()
+            => DifferentiatesBack("asin(1-2*x^2)^2");
+    }
+}