diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index 20f9a22be..15c014970 100644 --- a/BREAKING-CHANGES.md +++ b/BREAKING-CHANGES.md @@ -497,6 +497,24 @@ and 7.3.2. `x^m F(x^n)^p`, with `F` a sum holding an inverse tangent, inverse hy | `"(a+b*atanh(c*x^(3/2)))^2/x^4".ToEntity().Integrate("x")` | `integral(...)` | 330 characters | | `"(a+b*ln(c*(d+k/x^(2/3))^n))^2/x^5".ToEntity().Integrate("x")` | `integral(...)` | 2,034 characters | +### A polynomial times an artanh or arcoth and the logarithm of its quadratic is integrated in closed form + +**Answers where there were none.** `x^4 (a + b atanh(c x)) (d + e ln(1 - c^2 x^2))` was declined, and so were the +other powers of x up to `x^5` and their arcoth twins, eleven of Rubi's 7.3.4 and 7.4.1. With `y = c x`, `T` the +inverse function and `B = ln(1 - y^2)`, the derivatives `T' = 1/(1 - y^2)` and `B' = -2 y/(1 - y^2)` are rational, and +the integral of `y^j T B` is `p0 T B + p1 T + p2 B + q1 T^2 + q2 B^2 + r` for polynomials `p0`, `p1`, `p2` and `r` and +constants `q1` and `q2`, the constants chosen so that each division by `1 - y^2` is exact. Only the two derivatives +are read, so the answer holds for the artanh and the arcoth alike and for any constant factor inside the logarithm. +Only where the product `T B` is in the integrand; without it the rules after this one answer +([#233](https://github.com/asc-community/AngouriMath/issues/233)). + +| Input | Was (2.5.0) | Now | +|---|---|---| +| `"(a+b*atanh(c*x))*(d+k*ln(1-c^2*x^2))".ToEntity().Integrate("x")` | `integral(...)` | 354 characters | +| `"x^4*(a+b*atanh(c*x))*(d+k*ln(1-c^2*x^2))".ToEntity().Integrate("x")` | `integral(...)` | 721 characters | +| `"x^5*(a+b*acoth(c*x))*(d+k*ln(1-c^2*x^2))".ToEntity().Integrate("x")` | `integral(...)` | 717 characters | +| `"(a+b*acoth(c*x))*(d+k*ln(1-c^2*x^2))".ToEntity().Integrate("x")` | `integral(...)` | 354 characters | + ### A power of x times an inverse tangent and the logarithm of its quadratic is integrated by parts **Answers where there were none.** `x^4 (a + b atan(c x)) (d + e ln(1 + c^2 x^2))` ran out of time, and so did three more diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs index 91e57daa8..2785988ff 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs @@ -28924,6 +28924,216 @@ static void Add(Dictionary into, int k, ERational value) return answer.Nodes.Any(node => node == MathS.NaN) ? null : answer; } + /// + /// A polynomial in x times a linear combination of 1, T, B and T B, + /// T = artanh(c x) or arcoth(c x) and B = ln(K (1 - c^2 x^2)), in closed form. + /// With y = c x, T' = 1/(1 - y^2) and B' = -2 y/(1 - y^2), the integral of + /// y^j T B is p0 T B + p1 T + p2 B + q1 T^2 + q2 B^2 + r, where + /// p0 = (y^(j + 1) - [j odd])/(j + 1), q1 = p0(1), q2 = p0(1)/4, and + /// p1, p2 and r are the integrals of (2 y p0 - 2 q1)/(1 - y^2), + /// (4 y q2 - p0)/(1 - y^2) and (2 y p2 - p1)/(1 - y^2), polynomials. + /// + /// + /// Matching the coefficients of T B, T, B and 1 asks that the three + /// numerators vanish at y = 1 and y = -1: the constant of p0 makes the first + /// two do, and the constants of p1 and p2 the third. Only the derivatives of + /// T and B are read, which are rational and the same on every branch, so the + /// answer holds for either inverse function and every K, in the generic case + /// c != 0. The coefficients are rational in y, and y^k is c^k x^k. + /// The rules after it ran out of time on Rubi's 7.3.4 and 7.4.1 + /// x^m (a + b artanh(c x))(d + e ln(1 - c^2 x^2)) and the arcoth twins. + /// Read only where T B is there: without it the rules after it answer. + /// https://github.com/asc-community/AngouriMath/issues/233 + /// + internal static Entity? SolveAPolynomialTimesAnInverseHyperbolicTangentAndTheLogarithmOfItsQuadratic(Entity expr, Entity.Variable x) + { + var logarithms = expr.Nodes.Where(node => node is Logf(var @base, var of) && @base == MathS.e && of.ContainsNode(x)).Distinct().ToList(); + if (logarithms.Count != 2) + return null; + // ln(N/D), N a multiple of 1 + c x and D one of 1 - c x, is 2 T: artanh(c x) and + // arcoth(c x) arrive so. The other is B. + Entity? inverse = null, square = null, c = null; + foreach (var logarithm in logarithms) + { + var (above, below) = Functions.SingleQuotient.Of(((Logf)logarithm).Antilogarithm); + if (!below.ContainsNode(x)) + { + if (square is not null) + return null; + square = logarithm; + continue; + } + if (inverse is not null + || !TreeAnalyzer.TryGetPolyLinear(above, x, out var n1, out var n0) + || !TreeAnalyzer.TryGetPolyLinear(below, x, out var d1, out var d0) + || Functions.PartialFractions.IsZeroAsAValue(n0) || Functions.PartialFractions.IsZeroAsAValue(d0) + || Functions.PartialFractions.IsZeroAsAValue(n1)) + return null; + // A numeric c is left to the rules after it, which answer it for an even power of x: + // `atanh(2 x) ln(1 - 4 x^2)` was. + var slope = (n1 / n0).InnerSimplified; + if (!slope.Vars.Any() || !Functions.PartialFractions.IsZeroAsAValue((d1 / d0 + slope).InnerSimplified)) + return null; + (inverse, c) = (logarithm, slope); + } + if (inverse is null || square is null || c is null + || !TreeAnalyzer.TryGetPolyQuadratic(((Logf)square).Antilogarithm, x, out var g2, out var g1, out var g0) + || Functions.PartialFractions.IsZeroAsAValue(g0) || !Functions.PartialFractions.IsZeroAsAValue(g1) + || !Functions.PartialFractions.IsZeroAsAValue((g2 / g0 + MathS.Sqr(c)).InnerSimplified)) + return null; + // The polynomials beside 1, T, B and T B, kind = [T] + 2 [B]. + var t = Variable.CreateUnique(expr, "t_ahl"); + var s = Variable.CreateUnique(expr, "s_ahl"); + var written = expr.Substitute(inverse, 2 * t).Substitute(square, s); + const int mostDegree = 20; + var beside = new Entity?[4, mostDegree + 1]; + static int? Small(EInteger power, int most) => power.Sign >= 0 && power.CanFitInInt32() && power.ToInt32Unchecked() <= most ? power.ToInt32Unchecked() : null; + if (!TreeAnalyzer.TryGetPolynomial(written, t, out var inT)) + return null; + foreach (var byT in inT) + { + if (Small(byT.Key, 1) is not { } i || !TreeAnalyzer.TryGetPolynomial(byT.Value, s, out var inS)) + return null; + foreach (var byS in inS) + { + if (Small(byS.Key, 1) is not { } j || !TreeAnalyzer.TryGetPolynomial(byS.Value, x, out var inX)) + return null; + foreach (var byX in inX) + { + var coefficient = byX.Value; + if (Small(byX.Key, mostDegree) is not { } k || coefficient.ContainsNode(x) || coefficient.ContainsNode(t) || coefficient.ContainsNode(s)) + return null; + var kind = i + 2 * j; + beside[kind, k] = beside[kind, k] is { } before ? before + coefficient : coefficient; + } + } + } + for (var k = 0; k <= mostDegree; k++) + for (var kind = 0; kind < 4; kind++) + if (beside[kind, k] is { } coefficient && Functions.PartialFractions.IsZeroAsAValue(coefficient)) + beside[kind, k] = null; + if (!Enumerable.Range(0, mostDegree + 1).Any(k => beside[3, k] is not null)) + return null; + // y^k/c^(j + 1) of the integral of y^j, as c^(k - j - 1) x^k. + Entity TimesAPowerOfC(Entity factor, int exponent) => exponent switch + { + 0 => factor, + 1 => factor * c, + > 1 => factor * MathS.Pow(c, Number.Integer.Create(exponent)), + -1 => factor / c, + _ => factor / MathS.Pow(c, Number.Integer.Create(-exponent)) + }; + var collected = new Entity?[6, mostDegree + 2]; + for (var kind = 0; kind < 4; kind++) + for (var j = 0; j <= mostDegree; j++) + { + if (beside[kind, j] is not { } coefficient) + continue; + if (IntegralOfAPowerBesideAnInverseHyperbolicTangentAndItsLogarithm(j, kind) is not { } integral) + return null; + for (var atom = 0; atom < 6; atom++) + for (var k = 0; k < integral[atom].Length; k++) + { + if (integral[atom][k].IsZero) + continue; + if (k > mostDegree + 1) + return null; + var term = TimesAPowerOfC(Number.Rational.Create(integral[atom][k]) * coefficient, k - j - 1); + collected[atom, k] = collected[atom, k] is { } before ? before + term : term; + } + } + var inverseFunction = inverse / 2; + Entity[] atoms = { inverseFunction * square, inverseFunction, square, MathS.Sqr(inverseFunction), MathS.Sqr(square), Number.Integer.One }; + Entity answer = Number.Integer.Zero; + for (var atom = 0; atom < 6; atom++) + { + Entity polynomial = Number.Integer.Zero; + for (var k = 0; k <= mostDegree + 1; k++) + if (collected[atom, k] is { } coefficient && coefficient.InnerSimplified is var folded + && !Functions.PartialFractions.IsZeroAsAValue(folded)) + polynomial += folded * (k == 0 ? Number.Integer.One : k == 1 ? x : MathS.Pow(x, Number.Integer.Create(k))); + if (polynomial != Number.Integer.Zero) + answer += polynomial * atoms[atom]; + } + answer = answer.InnerSimplified; + if (answer.Nodes.Any(node => node == MathS.NaN) || !Functions.PartialFractions.DerivativeHoldsAtSampledPoints(answer, expr, x)) + return null; + return answer; + } + + /// + /// The integral of y^j times 1, T, B or T B for + /// 0 to 3, where T' = 1/(1 - y^2) and B' = -2 y/(1 - y^2), + /// as the rational polynomials in y beside T B, T, B, T^2, + /// B^2 and 1, the two squares' constants; null where a division by + /// 1 - y^2 leaves a remainder, which the constants are chosen to prevent. + /// + private static ERational[][]? IntegralOfAPowerBesideAnInverseHyperbolicTangentAndItsLogarithm(int j, int kind) + { + static ERational[] Monomial(int degree, ERational coefficient) + { + var p = Enumerable.Repeat(ERational.Zero, degree + 1).ToArray(); + p[degree] = coefficient; + return p; + } + static ERational[] Plus(ERational[] p, ERational[] q) + { + var sum = Enumerable.Repeat(ERational.Zero, System.Math.Max(p.Length, q.Length)).ToArray(); + for (var k = 0; k < p.Length; k++) + sum[k] = sum[k].Add(p[k]); + for (var k = 0; k < q.Length; k++) + sum[k] = sum[k].Add(q[k]); + return sum; + } + static ERational[] Times(ERational factor, ERational[] p) => p.Select(a => a.Multiply(factor)).ToArray(); + static ERational[] TimesY(ERational[] p) => new[] { ERational.Zero }.Concat(p).ToArray(); + static ERational At(ERational[] p, ERational y) + { + var value = ERational.Zero; + for (var k = p.Length - 1; k >= 0; k--) + value = value.Multiply(y).Add(p[k]); + return value; + } + static ERational[] Integral(ERational[] p) => + new[] { ERational.Zero }.Concat(p.Select((a, k) => a.Divide(ERational.FromInt32(k + 1)))).ToArray(); + // n/(1 - y^2), q_k = n_k + q_(k - 2), and null unless (1 - y^2) q is n. + static ERational[]? OverOneMinusTheSquare(ERational[] n) + { + var q = Enumerable.Repeat(ERational.Zero, System.Math.Max(n.Length - 2, 0)).ToArray(); + for (var k = 0; k < q.Length; k++) + q[k] = n[k].Add(k >= 2 ? q[k - 2] : ERational.Zero); + var back = Plus(q, Times(ERational.FromInt32(-1), TimesY(TimesY(q)))); + for (var k = 0; k < System.Math.Max(back.Length, n.Length); k++) + if ((k < back.Length ? back[k] : ERational.Zero).CompareTo(k < n.Length ? n[k] : ERational.Zero) != 0) + return null; + return q; + } + var one = ERational.One; + var minusOne = ERational.FromInt32(-1); + var two = ERational.FromInt32(2); + var four = ERational.FromInt32(4); + var rise = Monomial(j + 1, one.Divide(ERational.FromInt32(j + 1))); + var zero = new[] { ERational.Zero }; + // p0(1) = -p0(-1): for an odd j, y^(j + 1) is even, and p0 is (y^(j + 1) - 1)/(j + 1). + var p0 = kind != 3 ? zero : j % 2 == 1 ? Plus(rise, Monomial(0, minusOne.Divide(ERational.FromInt32(j + 1)))) : rise; + var q1 = At(p0, one); + var q2 = q1.Divide(four); + if (OverOneMinusTheSquare(Plus(Times(two, TimesY(p0)), Monomial(0, q1.Multiply(two).Negate()))) is not { } p1Prime + || OverOneMinusTheSquare(Plus(Monomial(1, q2.Multiply(four)), Times(minusOne, p0))) is not { } p2Prime) + return null; + var p1 = Plus(kind == 1 ? rise : zero, Integral(p1Prime)); + var p2 = Plus(kind == 2 ? rise : zero, Integral(p2Prime)); + // The constants that make 2 y p2 - p1 vanish at 1 and at -1. + var c1 = At(p2, one).Subtract(At(p2, minusOne)).Subtract(At(p1, one).Add(At(p1, minusOne)).Divide(two)); + var c2 = At(p1, one).Subtract(At(p1, minusOne)).Subtract(At(p2, one).Multiply(two)).Subtract(At(p2, minusOne).Multiply(two)).Divide(four); + p1 = Plus(p1, Monomial(0, c1)); + p2 = Plus(p2, Monomial(0, c2)); + if (OverOneMinusTheSquare(Plus(Times(two, TimesY(p2)), Times(minusOne, p1))) is not { } rPrime) + return null; + var r = Plus(kind == 0 ? rise : zero, Integral(rPrime)); + return new[] { p0, p1, p2, new[] { q1 }, new[] { q2 }, r }; + } + /// /// A function of an inverse function of a linear L = c + d x times powers of multiples of /// L, (c e + d e x)^3 (a + b arsinh(c + d x))^2, under u = L: each power of diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs b/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs index aac987b4a..59ed76c5a 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs @@ -695,6 +695,10 @@ private static Entity WithItsPiecewiseSumsCombined(Entity expr) // A power of x times an inverse tangent and the logarithm of its quadratic, by parts with the // two remainders that have no antiderivative apiece paired, before the rules that timed out. if ((answer = IndefiniteIntegralSolver.SolveAPowerOfXTimesAnArctangentAndTheLogarithmOfItsQuadratic(expr, x)) is { }) return answer; + // A polynomial times 1, an inverse hyperbolic tangent or cotangent of c x, the logarithm of + // 1 - c^2 x^2 and their product, in closed form, before the rules that timed out on + // `x^3 (a + b artanh(c x))(d + e ln(1 - c^2 x^2))`. + if ((answer = IndefiniteIntegralSolver.SolveAPolynomialTimesAnInverseHyperbolicTangentAndTheLogarithmOfItsQuadratic(expr, x)) is { }) return answer; // A rational function over the quadratic of an inverse tangent beside it, split, before the // rules that timed out on `x^4 atan(a x)/(c + a^2 c x^2)`. if ((answer = IndefiniteIntegralSolver.SolveARationalFunctionOverTheQuadraticOfAnArctangentBesideIt(expr, x)) is { }) return answer; diff --git a/Sources/Tests/UnitTests/Calculus/APolynomialTimesAnInverseHyperbolicTangentAndTheLogarithmOfItsQuadraticIntegralTest.cs b/Sources/Tests/UnitTests/Calculus/APolynomialTimesAnInverseHyperbolicTangentAndTheLogarithmOfItsQuadraticIntegralTest.cs new file mode 100644 index 000000000..3a17a9438 --- /dev/null +++ b/Sources/Tests/UnitTests/Calculus/APolynomialTimesAnInverseHyperbolicTangentAndTheLogarithmOfItsQuadraticIntegralTest.cs @@ -0,0 +1,74 @@ +// +// 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 +{ + /// + /// A polynomial times an inverse hyperbolic tangent or cotangent of c x and the logarithm of + /// 1 - c^2 x^2, in closed form: Rubi's 7.3.4 and 7.4.1 + /// x^m (a + b artanh(c x))(d + e ln(1 - c^2 x^2)), which timed out. + /// #233 + /// + [Trait("Area", "Calculus")] + public sealed class APolynomialTimesAnInverseHyperbolicTangentAndTheLogarithmOfItsQuadraticIntegralTest + { + private static void DifferentiatesBack(string integrand, params (string, double)[] pins) + { + var integral = integrand.ToEntity().Integrate("x"); + var text = integral.Stringize(); + Assert.DoesNotContain("integral(", text); + Assert.True(text.Length < 5000, $"{text.Length} characters of answer for {integrand}"); + + 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"); + var compared = 0; + foreach (var at in new[] { -2.3, -0.7, 0.4, 1.1, 1.9, 3.2 }) + { + var want = original.Substitute("x", at).EvalNumerical(); + if (want.IsNaN) + continue; + var got = derivative.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}"); + compared++; + } + Assert.True(compared >= 5, $"only {compared} points could be compared for {integrand}"); + } + + [Theory] + [InlineData("x^4*(a+b*atanh(c*x))*(d+k*ln(1-c^2*x^2))")] + [InlineData("x^3*(a+b*atanh(c*x))*(d+k*ln(1-c^2*x^2))")] + [InlineData("x^2*(a+b*atanh(c*x))*(d+k*ln(1-c^2*x^2))")] + [InlineData("x*(a+b*atanh(c*x))*(d+k*ln(1-c^2*x^2))")] + [InlineData("(a+b*atanh(c*x))*(d+k*ln(1-c^2*x^2))")] + [InlineData("x^5*(a+b*acoth(c*x))*(d+k*ln(1-c^2*x^2))")] + [InlineData("x^4*(a+b*acoth(c*x))*(d+k*ln(1-c^2*x^2))")] + [InlineData("x^3*(a+b*acoth(c*x))*(d+k*ln(1-c^2*x^2))")] + [InlineData("x^2*(a+b*acoth(c*x))*(d+k*ln(1-c^2*x^2))")] + [InlineData("x*(a+b*acoth(c*x))*(d+k*ln(1-c^2*x^2))")] + [InlineData("(x^2+1)*atanh(c*x)*ln(1-c^2*x^2)")] + [InlineData("x^7*atanh(c*x)*ln(1-c^2*x^2)")] + [InlineData("atanh(c*x)*ln(2-2*c^2*x^2)")] + public void IsAnsweredAndDifferentiatesBack(string integrand) + { + DifferentiatesBack(integrand, ("a", 1.3), ("b", 0.7), ("c", 0.3), ("d", 2), ("k", -0.4)); + DifferentiatesBack(integrand, ("a", -2), ("b", 1.5), ("c", -1.2), ("d", 0.5), ("k", 3)); + } + } +}