I'm not familiar with the algorithm, but it looks like a good addition to eval_sum_symbolic method. Currently, if the sequence term is a sum it
splits it in two and calls itself for each part.
If the input was a polynomial, after some number of recursive calls (apparently, equal to the number of terms), it gets reduced to a monomial, for which the sum is
evaluated with Faulhaber formula.
So, implementing a direct approach that works with the entire polynomial would improve the performance. I imagine it would start with
and appear near the beginning of eval_sum_symbolic.