Solution (source code)

= Solution

For fixed <interpolation nodes> $c_0,\ldots,c_{s-1}$, interpolate any $p\in\mathbb P_{s-1}[x]$ as
$$
p(x)=\sum_{k=0}^{s-1}p(c_k)\ell_k(x),
\qquad
\ell_k(x)=\prod_{j\ne k}\frac{x-c_j}{c_k-c_j}.
$$
Integrating shows that the unique weights making the <quadrature rule> exact through degree $s-1$ are
$$
\boxed{w_k=\int_a^b\ell_k(x)\,dx
=\int_a^b\prod_{j\ne k}\frac{x-c_j}{c_k-c_j}\,dx}.
$$

Now let the nodes be the zeros of the degree-$s$ <orthogonal polynomial> $q_s$. Given any $p\in\mathbb P_{2s-1}[x]$, <polynomial division> gives
$$
p=q_sh+r,
\qquad \deg h\le s-1,quad \deg r\le s-1.
$$
Orthogonality gives $\int_a^bq_sh=0$, while $q_s(c_k)=0$ makes the quadrature sum for $q_sh$ vanish. The rule is already exact on $r$, so it is exact on $p$. This is the defining precision property of <Gaussian quadrature>.