For $d\in N^+, d>1$, how to calculate $\sum_{i_1=1}^n\sum_{i_2=i_1+1}^n...\sum_{i_d=i_{d-1}+1}^nx_{i_1}...x_{i_d}$ in O(n) time? For example, when $d=2$, we have $\sum_{i=1}^n\sum_{j=i+1}^nx_ix_j=\frac{1}{2} ((\sum x_i)^2-\sum x_i^2)$ and time complexity is O(n).
 A: The method that Fimpellizieri alluded to by finding the coefficients of a polynomial can be done in $O(n^2)$ time. A procedure could be as follows:


*

*Define the polynomial $p(x) = (x - x_1)(x-x_2)\cdots(x-x_n)$. The problem then becomes to find the coefficient of $x^{n-d}$ in $p(x)$. We choose to calculate $p(1)$, $p(\omega)$, $p(\omega^2), \cdots,$ $p(\omega^{n-1})$, a procedure that takes $O(n^2)$ time. Here $\omega = e^{2\pi i/n}$.

*Taking $p(x) = a_{n-1} x^{n-1} + a_{n-2}x^{n-2} + \cdots + a_1 x + a_0$, we now have obtained the matrix equation 
$$\begin{bmatrix}
1 & 1 & 1 & \cdots & 1 \\
1 & \omega & \omega^2 & \cdots & \omega^{n-1} \\
1 & \omega^2 & \omega^4 & \cdots & \omega^{2(n-1)} \\
\vdots & \vdots & \vdots & \ddots & \vdots \\
1 & \omega^{n-1} & \omega^{2(n - 1)} & \cdots & \omega^{(n-1)(n-1)}
\end{bmatrix} \begin{bmatrix} a_0 \\ a_1 \\ a_2 \\ \vdots \\ a_{n - 1} \end{bmatrix} = \begin{bmatrix} p(1) \\ p(\omega) \\ p(\omega^2) \\ \vdots \\ p(\omega^{n - 1}) \end{bmatrix}$$  But you might have noticed that the first matrix looks familiar - in fact, it is the discrete fourier transform matrix $\mathcal{F}$! It has the property that $\frac{1}{\sqrt{n}} \mathcal{F}$ is unitary, meaning that $\mathcal{F}^{-1} = \frac{1}{n} \mathcal{F}^* = \frac{1}{n} \bar {\mathcal{F}}$ as $\mathcal{F}$ is also symmetric. This implies that $\vec{a} = \frac{1}{n} \bar{\mathcal{F}} \vec{p}$, where $\vec{a}$ is the vector of $a_i$'s and $p$ is the vector of $p(\omega^i)$'s. Since we wish to pick out only the coefficient of $x^{n-d}$ (this is equal to the quantity you desire by Vieta's formulas), we only need to multiply the $(n-d)$th row of $\frac{1}{n}\bar{\mathcal{F}}$ and $\vec{p}$, which yields
$$a_{n-d} = \frac{1}{n} \sum_{k = 0}^{n-1} \bar{\omega}^k p(\omega^k)$$
This step can be accomplished in $O(n)$ time and recovers $a_{n-d}$, the quantity you desired.
I know that this does not meet your time constraints, but I would like to highlight a couple strengths of this method. First, if the values of $p(\omega)$, $p(\omega^2), \cdots,$ $p(\omega^{n-1})$ where given, then we could bypass step $2$ and this algorithm would be linear in $n$. Second, if we wanted to obtain all coefficients in the polynomial (i.e. calculate the quantity you desire for every $d < n$), this would be just as efficient a method as using an $O(n)$ algorithm to calculate every single quantity. In fact, if you wanted to obtain all coefficients in the polynomial, step 2 can be done in $O(n \text{log} n)$ time, by utilizing the Fast Fourier Transform.
