Let $x_0,\cdots,x_k$ be roots of the Legendre polynomial $L_{k+1}(x)$. Show that for any $y\in(-1,1)$, the error of the Gauss-Legendre quadrature rule using $x_0,\cdots,x_k$ for approximating $\int_{-1}^{1}f(x)\,dx$ equals
$\int_{-1}^1 f[x_0,\cdots,x_k,y,x](x-x_0)\cdots(x-x_k)(x-y)\,dx$
Attempt
I presume the error term for $f(x) \approx \sum_{k=0}^{n}c_k(x-x_k)$ is $f[x_0,\dots,x_k,x](x-x_0)\dots(x-x_k)$, but the expected outcome differs by some terms. I am not sure why.
We need to use a trick here and reorder the terms. Notice that
$$\begin{aligned} f[x_0,\dots,x_k,y,x] &= f[y,x_0,x_1,\dots,x_k,x] \\ &= \frac{f[x_0,\dots,x_k,x]-f[x_0,\dots,x_k,y]}{x-y} \end{aligned}$$
We can write the function $f(x)$ as
$$f(x)=q(x)L_{k+1}(x)+r(x)$$
where $r(x)$ has degree at most $k$ and $L_{k+1}(x)$ has degree $k+1$. The degree of precision of the Gauss-Legendre quadrature rule is $2k+1$ that interpolates $k+1$ points, $x_0,\dots,x_k$. If $q(x)$ has degree at most $k$, then by orthonormality of the quadrature rule, $\int_{-1}^1q(x)L_{k+1}(x)\, dx=0$ and hence this term becomes our error term. We see that this is precisely the case when
$$\scriptsize \int_{-1}^1 f[x_0,\dots,x_k,y,x](x-x_0)\dots(x-x_k)(x-y)\, dx = \int_{-1}^1(f[x_0,\dots,x_k,x]-f[x_0,\dots,x_k,y])\cdot(x-x_0)\dots(x-x_k)\, dx $$
where $(f[x_0,\dots,x_k,x]-f[x_0,\dots,x_k,y])$ is a constant.