Efficiently evaluate a triple nested summation How do you efficiently evaluate the following nested sum, perhaps as a product of matrices and/or vectors:
$$
\sum_{i=1}^{p}\sum_{j=1}^{p}\sum_{k=1}^{p}\int\alpha_{i}\alpha_{j}\alpha_{k}\;\pi_{i}(s)\pi_{j}(s)\pi_{k}(s) \; d s
$$
where $\pi_j(s)$ represents a $j$-th order orthogonal polynomial and $\alpha_j$ represents its coefficients. 
I'm looking for an expansion along the lines of:
$$
\sum_{i=1}^{p}\int\alpha_{i}\pi_{i}(s)ds=
\mathbf{a^{T}}\mathbf{P}\mathbf{w}
$$
where for an $N-point$ quadrature rule, this becomes:
$$
\mathbf{a}=\left[\begin{array}{c} \alpha_{1}\\
\vdots\\
\alpha_{p}
\end{array}\right],\;\mathbf{w}=\left[\begin{array}{c} w_{1}\\
\vdots\\
w_{N}
\end{array}\right],\;\mathbf{P}_{ij}=\pi_{i}(s_{j}),
$$
where $i=1,\ldots, p$ and $j=1,\ldots, N$. 
Here the vector $\mathbf{w}$ are quadrature weights with corresponding points $s_1, \cdots, s_N$. Then the vector-matrix-vector product above is an evaluation of a bilinear form.
 A: Since the indices are independent, this reduces to $(\sum\limits_{i=1}^p a_i)(\sum\limits_{j=1}^p b_j)(\sum\limits_{k=1}^p c_k)$
I suppose one way of expressing it would be $A^TV B^TV C^TV$, where V is the column matrix with $v_i = 1 \,\, \forall i$
A: I think the problem got a bit over-simplified in its editing.  Instead of using $\alpha_i$ for essentially all the polynomial coefficients, I think we should treat the case of coefficients $a_i, b_i, c_i$ for $i=1,\ldots,p$.
In any case this "more general case" can be used for the special case.
Given that we are willing to have the integration done by quadrature, there is a significant simplification to the complexity over just evaluating the triple summation.
Let $f(s) = \sum_{i=1}^p a_i \pi_i(s)$, $g(s) = \sum_{i=1}^p b_i \pi_i(s)$, and $h(s) = \sum_{i=1}^p c_i \pi_i(s)$.  Then the revised expression is the integral of the triple product:
$$ \int f(s) g(s) h(s) ds = \sum_{i=1}^p \sum_{j=1}^p \sum_{k=1}^p \int a_i b_j c_k \; \pi_i(s) \pi_j(s) \pi_k(s) ds $$
A quadrature rule for the left hand side can take advantage of recursive rules for evaluating the "orthogonal polynomials" $\pi_i(s)$ at successive orders $i$ for some fixed quadrature node $s$.
It is suggested in the revised Question that an $N$-point rule will be used to evaluate/approximate the integral.  We will therefore have operational complexity of $N$ times the complexity of evaluating integrand $f(s) g(s) h(s)$, plus some linear complexity $O(N)$ for multiplying the evaluations by weights and summing up the resulting terms.
The evaluation of $f(s) g(s) h(s)$ comes down to evaluating each factor, and these share/amortize the evaluations of underlying orthogonal polynomials $\pi_i(s)$.  If the orthogonal polynomials are "from the Legendre family" (as suggested in the OP's Comment on the Question), then for each quadrature point $s$ we can loop through the orders $m=1,\ldots,p$ as follows:
$$ \pi_0(s) = 1 \; , \; \pi_1(s) = s $$
$$ (m+1) \pi_{m+1}(s) = (2m+1) s \;\pi_m(s) - m\;\pi_{m-1}(s) $$
On this loop we can piggyback the accumulation of sums $f(s), g(s), h(s)$ with the prescribed coefficients $a_i, b_i, c_i$ respectively.  Therefore the complexity of evaluating the integrand at one quadrature point is $O(p)$, which subsumes the constant overhead of finally multiplying the three accumulated factors, $f(s)g(s)h(s)$.
The overall complexity of this quadrature is therefore $O(pN)$, whereas a naive nested looping scheme will cost $O(p^3N)$.
