diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index 767f6bc37..927c9681c 100644 --- a/BREAKING-CHANGES.md +++ b/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 diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs index 12b4ccf64..22d434d73 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs @@ -28963,10 +28963,14 @@ static bool Through(Entity node, Entity.Variable x) /// 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 + { + [-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 diff --git a/Sources/Tests/UnitTests/Calculus/APowerOfAnInverseFunctionOfAQuadraticIntegralTest.cs b/Sources/Tests/UnitTests/Calculus/APowerOfAnInverseFunctionOfAQuadraticIntegralTest.cs index 0eb90e141..30284e6a1 100644 --- a/Sources/Tests/UnitTests/Calculus/APowerOfAnInverseFunctionOfAQuadraticIntegralTest.cs +++ b/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");