Consider the $$$n$$$ degree polynomial $$$P(x) = \sum_{i=0}^n p_ix^i$$$. I call the polynomial $$$Q(x)$$$ it's "prefix sum polynomial" if it satisfies the following condition:

Essentially, $$$Q(x)$$$ is a prefix sum on the values of $$$P(x)$$$. So,

Now, all we need is a formula for the prefix sum of $$$y^i$$$. This is given by Faulhaber's Formula:

where $$$B_j$$$ is the $$$j$$$th Bernoulli number. The formula only applies for positive values of $$$i$$$, so we have to handle $$$[x^0]P(x)$$$ separately. Moreover, we can tell from this formula that the degree of $$$Q(x)$$$ will be $$$n+1$$$.

Just using the formula above is enough to find $$$Q(x)$$$ in $$$O(n^2)$$$, but let us go further. From what we know so far, the coefficients of $$$Q(x)$$$ would be:

This looks suspiciously like a convolution. Let us define $$$A(x)$$$ and $$$B(x)$$$ as:

Then,

We can include the contribution of $$$[x^0]P(x)$$$ after this by adding it to $$$[x^0]Q(x)$$$ and $$$[x^1]Q(x)$$$. This way, we can calculate $$$Q(x)$$$ in $$$O(n \log n)$$$ using FFT.

Finally, we need to tackle the problem of calculating the Bernoulli numbers. We could use the recursive formula on the Wikipedia page, but that's $$$O(n^2)$$$. Instead, we can use the exponential generating function of the Bernoulli numbers:

Now, we can solve the entire problem in $$$O(n \log{n})$$$.

I'd like to thank neal since I only found this technique by looking at his submission for ARC 104 E. I later found out that this same idea is used in the problem SERSUM by chemthan, but I decided to make the blog anyway to introduce it to more people.