A sum over the roots of a polynomial, so that a rational function with an irreducible denominator can be integrated

[Rafael-SOWNet]

The integral

integral sqrt(x)/(1 + x + x^4) dx

is elementary, and this library declines it. Under u = sqrt(x) it is 2u^2/(1 + u^2 + u^8), a rational function whose denominator is irreducible over the rationals, is not a biquadratic and is not a binomial — so none of the three splits the rational integrator has can take it apart, and it comes back unevaluated. Wolfram Alpha gives

sum_({ω: ω^8 + ω^2 + 1 = 0}) (ω log(sqrt(x) - ω))/(4 ω^6 + 1) + constant

which is the general form: a sum over the roots of a polynomial, with a logarithm at each root. It was the pinned boundary of PowerSubstitutionIntegralTest.WhereTheChainStops in #1284, and the review there is right that a boundary this library cannot represent an answer for should be an issue rather than a test (#1284 (comment)).

What the answer is

For P/Q with deg P < deg Q and Q squarefree, the antiderivative is

sum over roots r of Q of  (P(r)/Q'(r)) * ln(x - r)

exactly — the residue at each simple pole times the logarithm there. That is what the Rothstein–Trager theorem states; the Lazard–Rioboo–Trager refinement writes the same sum over the roots of the resultant R(t) = res_x(P - t Q', Q) rather than of Q, which groups the roots that share a residue and keeps the algebraic extension as small as it can be. Every computer algebra system that answers this integral answers it in that form. IntegrateAPolynomialOverABinomial (#1284) is the special case where Q = a x^n + b, whose roots are known in closed form and pair off into real logarithms and arctangents; this issue is the general case, where they are not.

What is missing is a node, not a rule

The rule is two lines once the answer can be written down. What cannot be written down today is a sum over the roots of a polynomial when those roots have no closed form. Three shapes to consider:

1. A RootSum(polynomial, variable, expression) node — the way Maxima's rootsum, Maple's RootOf/sum, and SymPy's RootSum(Poly, Lambda) do it. Evaluates numerically by finding the roots; simplifies when the polynomial factors; differentiates termwise under the sum, so the differentiate-back check still applies. Prints as the sum above.
2. A RootOf(polynomial, index) node and an ordinary Sumf over indices — finer, since a single root is then a value in its own right (a quintic's roots could be spoken of), but the index has to mean something stable, and ordering complex roots is exactly the mess every system with RootOf has to document.
3. Decline, and record the reason in the answer — which is what happens now, without the reason.

The first is the smallest change that makes the integral answerable, and it is also what a sum over roots is. It touches the printers, Evaled/EvalNumerical, Differentiate, Substitute and the parser, which is why it is an issue and not a PR.

What it would answer

Rubi's test suite has these by the dozen: every P/Q whose Q is irreducible of degree three or more and not a binomial. In the 463-problem sample run here, 1/(x^3 + x + 1) is the smallest, and x/(1 + x + x^5), (x^2 + 1)/(x^4 + x + 1) and the like sit behind it. With the node, the rational integrator is complete for squarefree denominators; a repeated factor is then the Hermite reduction on top of it, which is a rule again.

Part of #718.

[Happypig375]

need to consider if representing it as "sum on a set" (compare: max and min over a set) would make the expression easier to interact with compared to a dedicated rootsum node, which isn't how we write it out in math.

[Rafael-SOWNet]

A sum on a set reads better, and it is also the smaller change: max(expr, var in over) already takes a Set as its range, and a sum(expr, var in over) beside it would be the same shape with the same evaluation contract -- enumerate the set where it is finite and known, leave the node standing otherwise. A root sum is then not a node of its own at all:

sum(ω ln(sqrt(x) - ω)/(4ω^6 + 1), ω in { ω : ω^8 + ω^2 + 1 = 0 })

with the range a ConditionalSet on a polynomial equation, which exists and prints like that today. What falls out of that shape without being written for it:

• Simplification is the set's. Where the polynomial factors over the rationals the conditional set splits into a union, and a sum over a union of disjoint sets is a sum of sums; where a factor is linear or quadratic the roots are closed forms and that part of the sum is written out. The rational integrator's three splits become the special cases where the whole set is closed-form.
• Numeric evaluation is root finding, which a polynomial with numeric coefficients always has, so EvalNumerical never has to decline on the sum itself.
• Differentiation is termwise, since the range does not depend on the outer variable, and the differentiate-back check the integrator's tests use applies unchanged.
• No index and no ordering: the sum is commutative, so the question of which root is the first -- the thing every RootOf(p, i) design has to document -- does not arise. That is the argument against option 2 in the issue and it applies to a dedicated rootsum node too, which would have to say what order its terms come in when printed.

Two things the shape has to decide: a sum over a set that is not finite is not a value and should stay unevaluated rather than answer anything; and whether the existing sum(expr, var, from, to) becomes var in {from..to} or stays as it is (I would leave it -- the integer range is what people write, and the set form is the generalisation, not a replacement).

So: sum on a set, and the rule in the rational integrator is two lines once that node exists. I have not started on it -- this is still a proposal -- but that is the shape I would build if it is accepted.

[Happypig375]

Your analysis looks good to me. If we are to add integer ranges, we also need to think about whether {from..to} is how it's written by people.

[Rafael-SOWNet]

On the integer range: I would not add a {from..to} literal with this. sum(expr, k, from, to) is already how people write a sum over integers -- it is Σ with its bounds, and Maxima, Maple and SymPy all spell it that way -- so the set form is the generalisation beside it and not a replacement, and nothing in the proposal needs a range literal. {a..b} is how Haskell and Pascal write one; on paper it is {1, 2, ..., n} or ℤ ∩ [1, n], and the second of those the parser can already say, ZZ intersect [1; n], once a set that shape enumerates. If a literal is wanted later, {1, ..., n} with an ellipsis token is the one that reads as mathematics, and it can be its own small issue.

If this is the shape to build, would you put Accepted on it? I will take it from there -- the node, the parser, the printers, evaluation by root finding, termwise differentiation, and the two-line rule in the rational integrator, with the differentiate-back check on the Rubi problems it unlocks.