From 9871064e6662fdc856105598d34c630b440cdd28 Mon Sep 17 00:00:00 2001 From: Rafael Vuijk Date: Tue, 20 Oct 2026 21:36:55 +0000 Subject: [PATCH] Integrate a polynomial times an artanh and the logarithm of its quadratic in closed form x^4 (a + b atanh(c x)) (d + e ln(1 - c^2 x^2)) and its kin up to x^5, with their arcoth twins, were declined. 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. 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, 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. The rule applies only where the product T B is in the integrand. Rubi pocket (352 rows of 7.3.4 and 7.4.1): 338 -> 347 solved. The moves re-run alone gain 10 rows (7.3.4:649-653, 7.4.1:361-363, 367, 368). Two rows master answers in 23-24 s timed out alone; at a 60 s budget both builds answer them with the same allocation to 0.01%. Part of #233. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_012sonx8iAspMiwRwokT1Ura --- BREAKING-CHANGES.md | 18 ++ .../Integration/IndefiniteIntegralSolver.cs | 210 ++++++++++++++++++ .../Integration/Integration.Definition.cs | 4 + ...dTheLogarithmOfItsQuadraticIntegralTest.cs | 74 ++++++ 4 files changed, 306 insertions(+) create mode 100644 Sources/Tests/UnitTests/Calculus/APolynomialTimesAnInverseHyperbolicTangentAndTheLogarithmOfItsQuadraticIntegralTest.cs 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)); + } + } +}