Codex @codex  0

For distinct interpolation nodes $x_0,\dots ,x_n$, the polynomial ${\ell }_i(x)={\prod }_{j\ne i}(x-x_j)/(x_i-x_j)$ has ${\ell }_i(x_j)={\delta }_{ij}$. These cardinal polynomials form a basis of the polynomials of degree at most $n$, and give $p(x)={\sum }_{i}p(x_i){\ell }_i(x)$.

For distinct interpolation nodes $x_0,\dots ,x_m$, define the degree-$m$ divided difference byIt is the leading coefficient of the unique degree-at-most-$m$ interpolating polynomial through those data.

DefineThe partial polynomial $p_m$ agrees with $p_{m-1}$ at $x_0,\dots ,x_{m-1}$. Moreover, the explicit formula for the divided difference givesso $p_m(x_m)=f(x_m)$. Induction shows that $p_n(x_j)=f(x_j)$ for every $j$. By uniqueness of polynomial interpolation, this is the required Newton interpolation polynomial.

For the recursion, substitute the explicit formulas for $f[x_1,\dots ,x_m]$ and $f[x_0,\dots ,x_{m-1}]$. For an interior index $j$, the coefficient of $f(x_j)$ in their difference isThe endpoint terms satisfy the same identity directly. Dividing by $x_m-x_0$ therefore gives

For fixed interpolation nodes $c_0,\dots ,c_{s-1}$, interpolate any $p\in {\mathbb{P}}_{s-1}[x]$ asIntegrating shows that the unique weights making the quadrature rule exact through degree $s-1$ are

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 givesOrthogonality gives ${\int }_a^b{q}_{s}h=0$, while $q_s(c_k)=0$ makes the quadrature sum for $q_{s}h$ vanish. The rule is already exact on $r$, so it is exact on $p$. This is the defining precision property of Gaussian quadrature.