From c267decf8a29ed76cb60462c6798a24f5a8d15d2 Mon Sep 17 00:00:00 2001 From: Rafael Vuijk Date: Wed, 21 Oct 2026 05:43:21 +0000 Subject: [PATCH] Leave a half-angle question whose logarithms want a sum over roots to the second pass 1/(1 + sin(x)^5) is asked one level down as (1 + sin^5)^(-1), and the half-angle substitution asks 2 (t^2 + 1)^4/((t^2 + 1)^5 + 32 t^5), whose denominator is (t + 1)^2 times an irreducible octic: its logarithms want residues of degree eight, which only the second pass writes, as a sum over roots. The first pass declined it after partial fractions had tried every split and spelling for 28 seconds, and the second answered in half a second. In the first pass the half-angle substitution now declines at once, setting ASumOverRootsWouldAnswer, where the question in t has rational coefficients, Hermite's reduction, the resultant and its factors show a residue field of degree past two, and no product of the denominator's factors is a polynomial in (t + s)^k with k >= 2 -- the binomials, biquadratics and blocks the closed-form rules read. The residue fields come from a helper LogarithmicPart now shares. Rubi's 4.1.7:326 goes from 29.4 s and 56 GB to 2.3 s and 0.75 GB, and 329 and sin(x)/(1 + sin(x)^5) likewise, with the same answers, character for character; 394 solved rows of 4.1.7, 86 of 4.2.7 and others answer as before. Part of #233. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_012sonx8iAspMiwRwokT1Ura --- .../Integration/IndefiniteIntegralSolver.cs | 14 ++ .../Continuous/Integration/RothsteinTrager.cs | 148 +++++++++++++++--- ...verRootsLeftToTheSecondPassIntegralTest.cs | 66 ++++++++ 3 files changed, 204 insertions(+), 24 deletions(-) create mode 100644 Sources/Tests/UnitTests/Calculus/ASumOverRootsLeftToTheSecondPassIntegralTest.cs diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs index 5c1f510a0..77dfaeeae 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs @@ -39246,6 +39246,20 @@ static bool NoRemainder(Entity rest) if (integrand.Vars.Any(symbol => symbol != t) && WithTheLeadingCoefficientDividedOut(integrand, t) is { } normalized) integrand = normalized; + // In the first pass, a rational function of t with rational coefficients whose + // logarithms want a residue of degree above two, by no factor a closed form reads, is + // the second pass's: see RothsteinTrager.WantsASumOverRootsAlone. Asked anyway, every + // split and spelling of it was searched before the decline, and Rubi's + // `1/(1 + sin(x)^5)` spent thirty seconds there where the second pass answers it in + // half of one. + if (!Integration.SumsOverRootsAllowed && !integrand.Vars.Any(symbol => symbol != t) + && TryReadAsQuotient(integrand, out var aboveInT, out var belowInT) + && RothsteinTrager.WantsASumOverRootsAlone(aboveInT, belowInT, t)) + { + Integration.ASumOverRootsWouldAnswer = true; + return null; + } + return Integration.ComputeIndefiniteIntegral(integrand, t, integrateByParts) is { } result ? result.Substitute(t, MathS.Tan(argument / 2)) : null; diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/RothsteinTrager.cs b/Sources/AngouriMath/Functions/Continuous/Integration/RothsteinTrager.cs index 9e68cec5a..177cd8ebd 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/RothsteinTrager.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/RothsteinTrager.cs @@ -182,31 +182,8 @@ private static (Entity Rational, RationalPolynomial Above, RationalPolynomial Be private static Entity? LogarithmicPart(RationalPolynomial above, RationalPolynomial below, Variable x, bool sumsOverRoots = true) { var derivative = Derivative(below); - // The resultant in t of D and A - t D' has degree at most that of D in t, and is - // found by evaluating it at that many rational t and interpolating: each value is - // one numeric Sylvester determinant. As a determinant with polynomial entries - // eliminated fraction-free it was thirteen seconds for a denominator of degree six - // with four-digit coefficients, Timofeev's `(1 + 2x)/((4 + 4x + 3x^2) sqrt(x^2 + 6x - 1))` - // under Euler's substitution, asked several times over. - if (ResultantInT(below, above, derivative) is not { } resultantPolynomial || resultantPolynomial.IsZero) - return null; - if (!sumsOverRoots && HasAFactorPastQuadraticModuloAPrime(resultantPolynomial)) + if (ResidueFields(above, below, derivative, declineOnAFactorPastQuadratic: !sumsOverRoots) is not { } factors) return null; - if (SquareFree(resultantPolynomial) is not { } residueParts) - return null; - - // Every residue field first, and the residues only then. A residue in a field of degree - // above two is written as a sum over roots, below, and where there is one, every - // residue that is not rational goes into that sum with it: the conjugate pair of a - // quadratic factor beside a quartic one was twelve seconds of gcds over the extension, - // and a sum over roots needs no arithmetic in any extension at all. - var factors = new List<(IReadOnlyList Factors, int Multiplicity)>(); - foreach (var part in residueParts) - { - if (PolynomialFactorization.FactorPrimitive(ToInteger(part.Factor).PrimitivePart()) is not { } irreducibles) - return null; - factors.Add((irreducibles, part.Multiplicity)); - } var overRoots = factors.Any(part => part.Factors.Any(irreducible => irreducible.Factor.Degree > 2)); if (overRoots && !sumsOverRoots) return null; @@ -291,6 +268,129 @@ private static bool HasAFactorPastQuadraticModuloAPrime(RationalPolynomial resul return false; } + /// + /// The irreducible factors over the rationals of each squarefree part of the resultant whose + /// roots are the residues of over , with the + /// part's multiplicity: every residue field, before any residue is computed in one. + /// where the resultant or a factorization declines. + /// + private static List<(IReadOnlyList Factors, int Multiplicity)>? ResidueFields( + RationalPolynomial above, RationalPolynomial below, RationalPolynomial derivative, bool declineOnAFactorPastQuadratic = false) + { + // The resultant in t of D and A - t D' has degree at most that of D in t, and is + // found by evaluating it at that many rational t and interpolating: each value is + // one numeric Sylvester determinant. As a determinant with polynomial entries + // eliminated fraction-free it was thirteen seconds for a denominator of degree six + // with four-digit coefficients, Timofeev's `(1 + 2x)/((4 + 4x + 3x^2) sqrt(x^2 + 6x - 1))` + // under Euler's substitution, asked several times over. + if (ResultantInT(below, above, derivative) is not { } resultantPolynomial || resultantPolynomial.IsZero) + return null; + // Seen modulo a prime before the squarefree decomposition costs anything, where a field + // past the second degree is a decline: see HasAFactorPastQuadraticModuloAPrime. + if (declineOnAFactorPastQuadratic && HasAFactorPastQuadraticModuloAPrime(resultantPolynomial)) + return null; + if (SquareFree(resultantPolynomial) is not { } residueParts) + return null; + + // Every residue field first, and the residues only then. A residue in a field of degree + // above two is written as a sum over roots, below, and where there is one, every + // residue that is not rational goes into that sum with it: the conjugate pair of a + // quadratic factor beside a quartic one was twelve seconds of gcds over the extension, + // and a sum over roots needs no arithmetic in any extension at all. + var factors = new List<(IReadOnlyList Factors, int Multiplicity)>(); + foreach (var part in residueParts) + { + if (PolynomialFactorization.FactorPrimitive(ToInteger(part.Factor).PrimitivePart()) is not { } irreducibles) + return null; + factors.Add((irreducibles, part.Multiplicity)); + } + return factors; + } + + /// + /// Whether the logarithms of the integral of over + /// , a proper fraction of polynomials with rational + /// coefficients that reads, have a residue in a field of degree + /// above two, so that writes them only as a sum over roots, and + /// no product of the denominator's irreducible factors past the second degree is a + /// polynomial in (x + s)^k with k at least two. + /// + /// + /// The residues of a rational function are its own, whatever spelling or split it is asked + /// in, and every antiderivative has them as the coefficients of its logarithms. The rules + /// that write such a residue in closed form read a block of the denominator that is a + /// binomial, a biquadratic or a polynomial in a power of x, at its roots of unity or + /// in nested square roots, after the shift that takes away its second term where there is + /// one; each is a polynomial in (x + s)^k, and no product of the denominator's + /// irreducible factors may be one. The integral is then nothing the first pass writes, + /// which keeps sums over roots for the second. + /// + internal static bool WantsASumOverRootsAlone(Entity numerator, Entity denominator, Variable x) + { + if (!PolynomialFactoring.TryGetRationalCoefficients(numerator, x, 1, 0, MaxDegree, out var aboveCoefficients) + || !PolynomialFactoring.TryGetRationalCoefficients(denominator, x, 1, 1, MaxDegree, out var belowCoefficients)) + return false; + var above = RationalPolynomial.Create(aboveCoefficients); + var below = RationalPolynomial.Create(belowCoefficients); + if (above.IsZero || above.Degree >= below.Degree) + return false; + var scale = ERational.One.Divide(below.Leading); + above = above.ScaleBy(scale); + below = below.ScaleBy(scale); + if (HermiteReduce(above, below, x) is not var (_, logAbove, logBelow) || logAbove.IsZero || logBelow.Degree > MaxDegree + || ResidueFields(logAbove, logBelow, Derivative(logBelow)) is not { } fields + || !fields.Any(part => part.Factors.Any(irreducible => irreducible.Factor.Degree > 2))) + return false; + // Every block the splits could hand a closed form: each product of distinct irreducible + // factors, the whole denominator's among them. `x^5 + 1` is `(x + 1)` times a quartic + // that is no binomial, and is a binomial whole. + if (PolynomialFactorization.FactorPrimitive(ToInteger(below).PrimitivePart()) is not { } irreducibles + || irreducibles.Count > MostFactorsCombined) + return false; + for (var subset = 1; subset < 1 << irreducibles.Count; subset++) + { + IntegerPolynomial? block = IntegerPolynomial.One; + for (var i = 0; i < irreducibles.Count && block is not null; i++) + if ((subset & (1 << i)) != 0) + block = block.Multiply(irreducibles[i].Factor); + if (block is null || block.Degree > 2 && InAPowerOfAShift(block)) + return false; + } + return true; + } + + /// The most irreducible factors whose products reads. + private const int MostFactorsCombined = 8; + + /// + /// Whether , of degree n, is a polynomial in (x + s)^k + /// for some k of two or more, with s = a_(n-1)/(n a_n): written in + /// y = x + s, the powers of y it has, past the constant, share a factor. + /// + private static bool InAPowerOfAShift(IntegerPolynomial poly) + { + var n = poly.Degree; + var shift = ERational.Create(poly[n - 1], poly[n].Multiply(EInteger.FromInt32(n))).ToLowestTerms(); + // Horner in y - s: g = g (y - s) + a_i, from the leading coefficient down. + var g = new ERational[n + 1]; + for (var i = 0; i <= n; i++) + g[i] = ERational.Zero; + g[0] = ERational.FromEInteger(poly[n]); + for (var i = n - 1; i >= 0; i--) + { + for (var j = n - i; j >= 1; j--) + g[j] = g[j - 1].Subtract(g[j].Multiply(shift)); + g[0] = ERational.FromEInteger(poly[i]).Subtract(g[0].Multiply(shift)); + } + var common = 0; + for (var j = 1; j <= n; j++) + if (!g[j].IsZero) + common = Gcd(common, j); + return common >= 2; + + static int Gcd(int a, int b) => b == 0 ? a : Gcd(b, a % b); + } + /// /// The logarithms at the poles whose residues are roots of , /// as a sum over those poles: sum(A(r)/D'(r) ln(x - r), r in { r : E(r) = 0 }), diff --git a/Sources/Tests/UnitTests/Calculus/ASumOverRootsLeftToTheSecondPassIntegralTest.cs b/Sources/Tests/UnitTests/Calculus/ASumOverRootsLeftToTheSecondPassIntegralTest.cs new file mode 100644 index 000000000..f483139fa --- /dev/null +++ b/Sources/Tests/UnitTests/Calculus/ASumOverRootsLeftToTheSecondPassIntegralTest.cs @@ -0,0 +1,66 @@ +// +// 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 +{ + /// + /// The half-angle tangent makes of 1/(1 + sin(x)^5) a rational function whose logarithms + /// want a residue of degree eight, which only the second pass writes, as a sum over roots; the + /// first pass searched every split of it for thirty seconds before declining. Rubi's 4.1.7. + /// #233 + /// + [Trait("Area", "Calculus")] + public sealed class ASumOverRootsLeftToTheSecondPassIntegralTest + { + private static string AnsweredAndCheckedBack(string integrand) + { + var integral = integrand.ToEntity().Integrate("x"); + var text = integral.Stringize(); + Assert.DoesNotContain("integral(", text); + // b > c, so that the root is real. + Entity Pinned(Entity e) => e.Substitute("b", 2.7).Substitute("c", 0.4); + var derivative = Pinned(integral.Substitute("C", 0)).Differentiate("x"); + var original = Pinned(integrand.ToEntity()); + var compared = 0; + foreach (var at in new[] { -1.9, -0.7, 0.3, 0.8, 1.3, 2.1 }) + { + var want = original.Substitute("x", at).EvalNumerical(); + var got = derivative.Substitute("x", at).EvalNumerical(); + if (want.IsNaN) + continue; + compared++; + Assert.True(Math.Abs((double)(got - want).RealPart) + Math.Abs((double)(got - want).ImaginaryPart) + < 1e-9 * Math.Max(1, Math.Abs((double)want.RealPart) + Math.Abs((double)want.ImaginaryPart)), + $"d/dx of the antiderivative of {integrand} is {got} at x = {at}, where the integrand is {want}"); + } + Assert.True(compared >= 5, $"only {compared} points could be compared for {integrand}"); + return text; + } + + // A residue of degree above two, by a factor no closed form reads: the second pass's sum over + // roots, as before. + [Theory] + [InlineData("1/(1 + sin(x)^5)")] + [InlineData("1/(1 - sin(x)^5)")] + [InlineData("sin(x)/(1 + sin(x)^5)")] + [InlineData("1/(1 + sin(x)^3)")] + [InlineData("1/(1 + sin(2*x)^5)")] + public void LeavesTheSumOverRootsToTheSecondPass(string integrand) + => Assert.Contains("sum(", AnsweredAndCheckedBack(integrand)); + + // In the half-angle tangent `1/(1 + cos(x)^3)` is `(1 + t^2)^2/(1 + 3 t^4)`, over a binomial + // the first pass writes at its roots: it is still answered in closed form. + [Theory] + [InlineData("1/(1 + cos(x)^3)")] + public void KeepsTheClosedFormOfAnEvenBlock(string integrand) + => Assert.DoesNotContain("sum(", AnsweredAndCheckedBack(integrand)); + } +}