I just finished a course on numerical mathematics, and have become quite interested in interpolation and how it ties into numerical integration. What suprised me while studiyng quadrature was the fact that Simpson's rule calculates integrals of cubic polynomials exactly, even though it is based on quadratic interpolation. Seeing how the structure and symmetry of equidistant points could lead to higher precision was kind of mind-blowing, and made me wonder about an "optimal" node arrangement for numerical integration.
Let's say we interpolate a $d^{\text{th}}$ degree polynomial $f \in \mathbb{P}_d$ at $n$ points $x_1, \dots x_n$ by an $(n-1)^{\text{th}}$ degree polynomial $p_{n-1} \in \mathbb{P}_{n-1}$. Then the interpolation error $f - p_{n-1}$ satisfies $$ f - p_{n-1} \in \mathbb{P}_d \quad \text{and} \quad (f-p_{n-1})(x_i) = 0, \quad i = 1, \dots, n. $$ That is, the error is a $d^{\text{th}}$ degree polynomial with a least $n$ roots, and can thus be written as $$ (f-p_{n-1})(x) = p_{d-n}(x) \prod_{i=1}^n (x - x_i), $$ where $p_{d-n}$ is some polynomial of degree $d-n$. Now, for the integral of the error to be zero independent of $p_{d-n}$, I suppose the only way is to choose the node polynomial to be orthogonal to $p_{d-n}$, which can only be achieved in general if the nodal polynomial is a (possibly scaled) Legendre polynomial of degree strictly greater than $d-n$. Thus, we need $d-n < n$, which implies $d \leq 2n-1$, and $ \prod_{i=1}^n (x - x_i) = cL_n(x)$ (scaled Legendre polynomial of degree $n$).
Now my question is: is this sufficient to deduce that we should pick $x_1, \dots, x_n$ as the $n$ roots of $L_n$? And more generally: is $2n-1$ really the maximum precision given $n$ points?