Integrate half an odd power of an inverse hyperbolic cosine of c + d x^2 in error functions

BREAKING-CHANGES.md

@@ -396,6 +396,20 @@ 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 |
 
+### Half an odd power of an inverse hyperbolic cosine of `c + d x^2` is integrated in error functions
+
+**Answers where there were none.** `(a + b acosh(1 + d x^2))^(1/2)` was declined, and so were thirteen more of Rubi's
+7.2.5, from the power `-7/2` to `5/2`. With `P = a + b u`, the integrals of `e^(+-u/2) P^(-1/2)` are error functions of
+`sqrt(P)/sqrt(2 b)`; the power `1/2` follows by parts and the rest by the recurrence for the whole powers, run up and
+down ([#718](https://github.com/asc-community/AngouriMath/issues/718)).
+
+| Input | Was (2.5.0) | Now |
+|---|---|---|
+| `"(a+b*acosh(1+d*x^2))^(1/2)".ToEntity().Integrate("x")` | `integral(...)` | 625 characters |
+| `"(a+b*acosh(-1+d*x^2))^(5/2)".ToEntity().Integrate("x")` | `integral(...)` | 868 characters |
+| `"1/(a+b*acosh(1+d*x^2))^(1/2)".ToEntity().Integrate("x")` | `integral(...)` | 418 characters |
+| `"1/(a+b*acosh(-1+d*x^2))^(7/2)".ToEntity().Integrate("x")` | `integral(...)` | 1,098 characters |
+
 ### A difference of two cubes is split over the cube roots of unity where its answer cannot be sampled
 
 **Answers where there were none.** `1/(a - b cosh(x)^3)` was declined, which Rubi's 6.2.7 has. Under `u = e^x` it is

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

@@ -28963,10 +28963,14 @@ static bool Through(Entity node, Entity.Variable x)
         /// </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.IsZero || power.EInteger.Abs().CompareTo(EInteger.FromInt32(12)) > 0)
-                return null;
-            var n = power.EInteger.ToInt32Unchecked();
+            var (@base, power) = expr is Powf(var b0, Number.Rational p0) ? (b0, p0) : (expr, (Number.Rational)Number.Integer.One);
+            // A whole power, or half an odd one.
+            var twice = (power * 2).InnerSimplified;
+            if (twice is not Number.Integer { EInteger: var twiceE } || twiceE.IsZero || twiceE.Abs().CompareTo(EInteger.FromInt32(24)) > 0)
+                return null;
+            var doubled = twiceE.ToInt32Unchecked();
+            var half = doubled % 2 != 0;
+            var n = doubled / 2;
             // 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)
@@ -29006,6 +29010,8 @@ static bool Through(Entity node, Entity.Variable x)
             var h = root / (d * x);
             var e = (Entity)Number.Integer.Create(sign);
             var a_bu = a + b * inverse;
+            if (half)
+                return HalfPower();
             if (n < 0)
                 return NegativePower();
             // I(k) by the recurrence, from I(0) = x and I(1) = x (a + b u) + 2 e b H.
@@ -29025,6 +29031,53 @@ static bool Through(Entity node, Entity.Variable x)
                 return null;
             return integrals[n].InnerSimplified;
 
+            // Half an odd power of the inverse hyperbolic cosine, through u as below: with
+            // P = a + b u, the integrals of e^(+-u/2) P^(-1/2) are the error functions of
+            // z = sqrt(P)/sqrt(2 b), (2 sqrt(pi)/sqrt(2 b)) e^(-+a/(2 b)) erfi(z) and erf(z), z written
+            // so that z sqrt(2 b) is the root of P itself. I(1/2) is by parts from them, and the
+            // recurrence runs up and down from I(-1/2) and I(1/2). The circular functions would want
+            // the Fresnel integrals, which have no node.
+            Entity? HalfPower()
+            {
+                if (inverse is not Logf || sign != -1 || squareAtZero != 1)
+                    return null;
+                var plus = Functions.PartialFractions.IsZeroAsAValue((c - 1).InnerSimplified);
+                var theta = inverse / 2;
+                Entity Exp(Entity z) => MathS.Pow(MathS.e, z);
+                Entity Sinh(Entity z) => (Exp(z) - Exp(-z)) / 2;
+                Entity Cosh(Entity z) => (Exp(z) + Exp(-z)) / 2;
+                // g = sinh, h = g' = cosh for c = 1; g = cosh, h = sinh for c = -1.
+                var kappa = x / (plus ? Sinh(theta) : Cosh(theta));
+                var z = MathS.Sqrt(a_bu) / MathS.Sqrt(2 * b);
+                var front = 2 * MathS.Sqrt(MathS.pi) / MathS.Sqrt(2 * b);
+                var upward = front * Exp(-a / (2 * b)) * MathS.Erfi(z);
+                var downward = front * Exp(a / (2 * b)) * MathS.Erf(z);
+                // The integral over u of cosh(theta) P^(-1/2), or of sinh(theta) P^(-1/2).
+                Entity OverRoot(bool ofCosh) => ofCosh ? (upward + downward) / 2 : (upward - downward) / 2;
+                var root = MathS.Sqrt(a_bu);
+                var halves = new Dictionary<int, Entity>
+                {
+                    [-1] = kappa / 2 * OverRoot(plus),
+                    // By parts, with the antiderivative of h(theta) in u being 2 h~(theta), h~ the other one.
+                    [1] = kappa / 2 * (2 * (plus ? Sinh(theta) : Cosh(theta)) * root - b * OverRoot(!plus)),
+                };
+                Entity Step(int t) => (Entity)Number.Rational.Create(t, 2);
+                for (var t = 3; t <= doubled; t += 2)
+                {
+                    var m = Step(t);
+                    halves[t] = x * MathS.Pow(a_bu, m) + 2 * e * m * b * h * MathS.Pow(a_bu, m - 1)
+                        - 4 * e * MathS.Sqr(b) * m * (m - 1) * halves[t - 4];
+                }
+                for (var t = -3; t >= doubled; t -= 2)
+                {
+                    var m = Step(t + 4);
+                    halves[t] = (x * MathS.Pow(a_bu, m) + 2 * e * m * b * h * MathS.Pow(a_bu, m - 1) - halves[t + 4])
+                        / (4 * e * MathS.Sqr(b) * m * (m - 1));
+                }
+                var answer = halves[doubled].InnerSimplified;
+                return Functions.PartialFractions.DerivativeDiffersWhereReal(answer, expr, x) ? null : answer;
+            }
+
             // A negative power, through u: w - c is a constant times g(theta)^2 with theta = u/2 + beta,
             // so x = kappa g(theta) with kappa = x/g(theta) constant on each side of 0, and
             // dx = (kappa/2) g'(theta) du. Then I(-1) and I(-2) are sine and cosine integrals of

Sources/Tests/UnitTests/Calculus/APowerOfAnInverseFunctionOfAQuadraticIntegralTest.cs

@@ -76,6 +76,20 @@ public void ANegativePowerIsAnsweredAndDifferentiatesBack(string integrand, int
             DifferentiatesBack(integrand, ("a", -0.75), ("b", -0.7), ("d", 0.15 * sign));
         }
 
+        // Half an odd power of the inverse hyperbolic cosine, in the error functions of
+        // sqrt(a + b u)/sqrt(2 b).
+        [Theory]
+        [InlineData("(a+b*acosh(1+d*x^2))^(1/2)")]
+        [InlineData("(a+b*acosh(1+d*x^2))^(5/2)")]
+        [InlineData("1/(a+b*acosh(1+d*x^2))^(3/2)")]
+        [InlineData("(a+b*acosh(-1+d*x^2))^(3/2)")]
+        [InlineData("1/(a+b*acosh(-1+d*x^2))^(7/2)")]
+        public void HalfAnOddPowerIsAnsweredAndDifferentiatesBack(string integrand)
+        {
+            DifferentiatesBack(integrand, ("a", 2), ("b", 1.25), ("d", 0.3));
+            DifferentiatesBack(integrand, ("a", 1.5), ("b", 0.7), ("d", 1.5));
+        }
+
         [Fact]
         public void WithNumbersIsAnswered()
             => DifferentiatesBack("asin(1-2*x^2)^2");