A rational integral over the roots of its denominator

BREAKING-CHANGES.md

@@ -1699,6 +1699,31 @@ where it is worked out, is not affected.
 | `"sum(1/(k - x), k, 1, n)".ToEntity().DomainCondition` | `not k - x = 0` | `forall k in { k : k in ZZ and 1 <= k and k <= n } : not k - x = 0` |
 | `"integral(1/(t - x), t, 0, 1)".ToEntity().DomainCondition` | `not t - x = 0` | `forall t in [0; 1] : not t - x = 0` |
 
+### A rational function whose residues are of degree three or more is integrated over its poles
+
+**Wider.** The logarithmic part of a rational integral is a logarithm at each pole, weighted by the
+residue there. Where the residues are rational or roots of a quadratic they were written out, as
+logarithms and arctangents. Where they lie in a field of degree three or more, the rational
+integrator declined, and `1/(x^3 + x + 1)` was left unevaluated. It is written now as a sum over the
+poles, `sum(A(r)/D'(r) ln(x - r), r in { r : E(r) = 0 })`, with the node above: `E` is the factor
+of the denominator whose poles have those residues, the gcd of `D` and `D'^n R(A/D')` over the
+rationals. The integrand of [#1285](https://github.com/asc-community/AngouriMath/issues/1285),
+`sqrt(x)/(1 + x + x^4)`, is `2u^2/(1 + u^2 + u^8)` under `u = sqrt(x)`, and is answered through it.
+
+**Only where nothing else answers.** The integral is asked again with sums over roots allowed once
+the whole of it has come back unevaluated, so an integrand another rule writes in closed form keeps
+that form. Asked with the other rules, the sum answered each of Jeffrey's three terms of
+`(-1 + 4 cos(x) + 5 cos(x)^2)/(-1 - 4 cos(x) - 3 cos(x)^2 + 4 cos(x)^3)` with a sum over the roots of
+a sextic, where the whole is one arctangent. A denominator with a symbol among its coefficients is
+still declined. Both columns measured on a build, `v2.5.0` against this change.
+
+| | Was (2.5.0) | Is |
+|---|---|---|
+| `"1/(x^3 + x + 1)".Integrate("x")` | `integral(1 / (x ^ 3 + x + 1), x)` | `sum(1 / (3 * r ^ 2 + 1) * ln(x - r), r in { r : r ^ 3 + r + 1 = 0 }) + C` |
+| `"sqrt(x)/(1 + x + x^4)".Integrate("x")` | `integral(sqrt(x) / (1 + x + x ^ 4), x)` | `sum(2 * r ^ 2 / (8 * r ^ 7 + 2 * r) * ln(x ^ (1/2) - r), r in { r : r ^ 8 + r ^ 2 + 1 = 0 }) + C` |
+| `"1/(1 - x^4 + x^8)".Integrate("x")` | `integral(1 / (1 - x ^ 4 + x ^ 8), x)` | `sum(1 / (8 * r ^ 7 + (-4) * r ^ 3) * ln(x - r), r in { r : r ^ 8 + -r ^ 4 + 1 = 0 }) + C` |
+| `"1/(x^4 + a*x + 1)".Integrate("x")` | `integral(1 / (x ^ 4 + a * x + 1), x)` | the same |
+
 ### `(a + b asech(c x))/(d + e x)^2` is integrated, and a root written apart with `|x|` no longer needs a parity
 
 Rubi's 7.5.1 with a symbolic linear below the bar ran for ten minutes without an answer. By

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

@@ -446,14 +446,54 @@ private static Entity Normalized(Entity expr, Entity.Variable x) =>
                 descentTruncated = true;
                 return null;
             }
-            if (descentDepth == 0)
+            if (descentDepth != 0)
+                return ComputeIndefiniteIntegralGuarded(expr, x, integrateByParts);
+            descentTruncated = false;
+            inProgress?.Clear();
+            ASumOverRootsWouldAnswer = false;
+            var answer = ComputeIndefiniteIntegralGuarded(expr, x, integrateByParts);
+            if (answer is not null || !ASumOverRootsWouldAnswer || SumsOverRootsAllowed)
+                return answer;
+            // Once more with sums over roots, and with a memo of its own, since the first pass
+            // remembered every one of those rational functions as declined. Last, and only where
+            // nothing else answered: asked with the rest, a sum over roots answered each of
+            // Jeffrey's three terms, where the whole is one arctangent the rational integrator
+            // writes when it is asked the whole.
+            // https://github.com/asc-community/AngouriMath/issues/1285
+            var (memo, stamp) = (answered, answeredUnder);
+            SumsOverRootsAllowed = true;
+            answered = null;
+            descentTruncated = false;
+            inProgress?.Clear();
+            try
             {
-                descentTruncated = false;
-                inProgress?.Clear();
+                answer = ComputeIndefiniteIntegralGuarded(expr, x, integrateByParts);
             }
-            return ComputeIndefiniteIntegralGuarded(expr, x, integrateByParts);
+            finally
+            {
+                SumsOverRootsAllowed = false;
+                (answered, answeredUnder) = (memo, stamp);
+            }
+            // In place of the decline the first pass remembered, which a question asked again
+            // would otherwise be answered with, without the rule ever being reached to say that a
+            // second pass could answer it.
+            if (answer is not null && memo is not null && stamp is not null && SettingsState.StillHolds(stamp))
+                memo[(Normalized(expr, x), x, integrateByParts, true)] = answer;
+            return answer;
         }
 
+        /// <summary>
+        /// Whether the rational integrator may answer with a sum over the roots of a factor of the
+        /// denominator: only in the second pass of a question nothing else answered.
+        /// </summary>
+        [System.ThreadStatic] internal static bool SumsOverRootsAllowed;
+
+        /// <summary>
+        /// Set where the rational integrator declined only because the residues lie in a field of
+        /// degree above two, which a sum over roots would have written.
+        /// </summary>
+        [System.ThreadStatic] internal static bool ASumOverRootsWouldAnswer;
+
         /// <summary>
         /// <see cref="ComputeIndefiniteIntegral"/> past its entry check: the cycle guard, the
         /// depth, and the memo.

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

@@ -39,8 +39,10 @@ namespace AngouriMath.Functions.Algebra
     /// <c>a ± w</c> with <c>w^2</c> rational, and the gcd is computed in <c>Q(w)</c>: for
     /// <c>w</c> real the pair is two logarithms with a root in them, and for <c>w = i b</c> it
     /// is <c>a ln(P^2 + b^2 Q^2) + b LogToAtan(P, b Q)</c>, Rioboo's continuous arctangent
-    /// form. A factor of higher degree would put the residues in a field this does not do
-    /// arithmetic in, and is declined.
+    /// form. A factor of higher degree puts the residues in a field this does no arithmetic in,
+    /// and those logarithms are written as a sum over the roots of the factor of the denominator
+    /// they belong to, <c>sum(A(r)/D'(r) ln(x - r), r in { r : E(r) = 0 })</c>
+    /// (https://github.com/asc-community/AngouriMath/issues/1285).
     /// </para>
     /// <para>
     /// Bronstein, <i>Symbolic Integration I</i>, §2.2 (HermiteReduce), §2.5
@@ -64,8 +66,8 @@ internal static class RothsteinTrager
         /// <summary>
         /// The integral of <paramref name="numerator"/> over <paramref name="denominator"/>,
         /// both polynomials in <paramref name="x"/> with rational coefficients and the fraction
-        /// proper; <see langword="null"/> where they are not, or where a residue lies in a field
-        /// of degree above two.
+        /// proper; <see langword="null"/> where they are not, or where the answer does not
+        /// differentiate back to the integrand at the sampled points.
         /// </summary>
         internal static Entity? Integrate(Entity numerator, Entity denominator, Variable x)
         {
@@ -148,23 +150,54 @@ private static (Entity Rational, RationalPolynomial Above, RationalPolynomial Be
             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 declines the whole, and the conjugate pair of a quadratic factor
-            // beside it was twelve seconds of gcds over the extension before the quartic factor
-            // was reached and declined it.
-            var factored = new List<(IReadOnlyList<SquareFreeDecomposition.SquareFreePart> Factors, int Multiplicity)>();
+            // 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<SquareFreeDecomposition.SquareFreePart> Factors, int Multiplicity)>();
             foreach (var part in residueParts)
             {
-                var factors = PolynomialFactorization.FactorPrimitive(ToInteger(part.Factor).PrimitivePart());
-                if (factors is null || factors.Any(irreducible => irreducible.Factor.Degree > 2))
+                if (PolynomialFactorization.FactorPrimitive(ToInteger(part.Factor).PrimitivePart()) is not { } irreducibles)
                     return null;
-                factored.Add((factors, part.Multiplicity));
+                factors.Add((irreducibles, part.Multiplicity));
+            }
+            var overRoots = factors.Any(part => part.Factors.Any(irreducible => irreducible.Factor.Degree > 2));
+            // Only once nothing else has answered the question: see
+            // Integration.SumsOverRootsAllowed. Until then this declines, as it always did.
+            if (overRoots && !Integration.SumsOverRootsAllowed)
+            {
+                Integration.ASumOverRootsWouldAnswer = true;
+                return null;
+            }
+            var factored = new List<(IReadOnlyList<SquareFreeDecomposition.SquareFreePart> Factors, int Multiplicity)>();
+            var summed = RationalPolynomial.One;
+            var summedRoots = 0;
+            foreach (var (irreducibles, multiplicity) in factors)
+            {
+                if (!overRoots)
+                {
+                    factored.Add((irreducibles, multiplicity));
+                    continue;
+                }
+                factored.Add((irreducibles.Where(irreducible => irreducible.Factor.Degree == 1).ToList(), multiplicity));
+                foreach (var irreducible in irreducibles.Where(irreducible => irreducible.Factor.Degree > 1))
+                {
+                    summed = summed.Multiply(RationalPolynomial.FromInteger(irreducible.Factor));
+                    summedRoots += irreducible.Factor.Degree * multiplicity;
+                }
             }
             Entity total = Integer.Create(0);
-            foreach (var (factors, multiplicity) in factored)
+            if (overRoots)
             {
-                var part = (Multiplicity: multiplicity, Factors: factors);
-                foreach (var irreducible in factors)
+                if (OverTheRoots(summed, summedRoots, above, below, derivative, x) is not { } sum)
+                    return null;
+                total = sum;
+            }
+            foreach (var (kept, multiplicity) in factored)
+            {
+                var part = (Multiplicity: multiplicity, Factors: kept);
+                foreach (var irreducible in kept)
                 {
                     var f = irreducible.Factor;
                     Entity? term = f.Degree switch
@@ -181,6 +214,42 @@ private static (Entity Rational, RationalPolynomial Above, RationalPolynomial Be
             return total;
         }
 
+        /// <summary>
+        /// The logarithms at the poles whose residues are roots of <paramref name="residues"/>,
+        /// as a sum over those poles: <c>sum(A(r)/D'(r) ln(x - r), r in { r : E(r) = 0 })</c>,
+        /// where <c>E</c> is the factor of <c>D</c> they are the roots of. <see langword="null"/>
+        /// where <c>E</c> does not have the <paramref name="roots"/> the residues' multiplicities
+        /// count.
+        /// </summary>
+        /// <remarks>
+        /// The residue at a simple pole <c>b</c> of <c>A/D</c> is <c>A(b)/D'(b)</c>, a root of the
+        /// resultant, and <c>D'(b)</c> is not zero there. So the poles whose residues are roots of
+        /// <c>R</c> are the common roots of <c>D</c> and the polynomial <c>D'^n R(A/D')</c>, and
+        /// <c>E</c> is their gcd over the rationals: no arithmetic in the field the residues lie in
+        /// is needed. The sum is the one Rothstein–Trager's theorem writes over the roots of the
+        /// resultant, with its terms taken pole by pole. <c>E</c> has a factor of degree above two,
+        /// since a pole in a field of degree at most two has its residue there too, so the sum is
+        /// left standing rather than written out in radicals, unless each such factor is a
+        /// binomial, whose roots are written.
+        /// https://github.com/asc-community/AngouriMath/issues/1285
+        /// </remarks>
+        private static Entity? OverTheRoots(RationalPolynomial residues, int roots,
+            RationalPolynomial above, RationalPolynomial below, RationalPolynomial derivative, Variable x)
+        {
+            // D'^n R(A/D') = sum over k of r_k A^k D'^(n - k).
+            var n = residues.Degree;
+            var cleared = RationalPolynomial.Zero;
+            for (var k = 0; k <= n; k++)
+                cleared = cleared.Add(above.Pow(k).Multiply(derivative.Pow(n - k)).ScaleBy(residues[k]));
+            var poles = RationalPolynomial.FromInteger(IntegerPolynomial.Gcd(ToInteger(below), ToInteger(cleared)).PrimitivePart());
+            if (poles.Degree != roots)
+                return null;
+            // The polynomials are in x alone, so any other name is free.
+            var r = MathS.Var(x.Name == "r" ? "w" : "r");
+            var summand = above.ToEntity(r) / derivative.ToEntity(r) * MathS.Ln(x - r);
+            return MathS.Sum(summand, r, new Set.ConditionalSet(r, poles.ToEntity(r).Equalizes(Integer.Zero)));
+        }
+
         /// <summary><c>c ln(gcd(D, A - c D'))</c> for a rational residue <paramref name="c"/>.</summary>
         private static Entity? RationalResidue(
             ERational c, int degree, RationalPolynomial above, RationalPolynomial below, RationalPolynomial derivative, Variable x)