Solution (source code)

= Solution
{c}

Set $W=\operatorname{diag}(w_0,\ldots,w_n)$. For column indices $j,k\in\{1,\ldots,n+1\}$, the product $P_{j-1}P_{k-1}$ has <polynomial degree> at most $2n$, so exactness of <Gauss-Legendre quadrature> and <orthonormality> give
$$
(A^TWA)_{jk}=\sum_{i=0}^nw_iP_{j-1}(x_i)P_{k-1}(x_i)
=\int_{-1}^1P_{j-1}(x)P_{k-1}(x)\,dx=\delta_{jk}.
$$
Thus $Q=W^{1/2}A$ is a square <matrix> with $Q^TQ=I$, and therefore is an <orthogonal matrix>: $QQ^T=I$ as well. It follows that
$$
\boxed{AA^T=W^{-1},\qquad
(AA^T)_{ij}=\begin{cases}
w_{i-1}^{-1},&i=j,\\
0,&i\ne j.
\end{cases}}
$$
Equivalently, the rows of the <generalized Vandermonde matrix> are <orthogonal> and
$$
\sum_{k=0}^nP_k(x_i)^2=\frac1{w_i}.
$$
The normalization matters: here $P_{n+1}=\sqrt{(2n+3)/2}\,L_{n+1}$, where $L_{n+1}$ is the conventional <Legendre polynomial> used above. Thus a formula containing only the nodes and the question's <orthonormal> <polynomials> is
$$
\boxed{(AA^T)_{ij}=\delta_{ij}
\frac{(1-x_{i-1}^2)[P_{n+1}'(x_{i-1})]^2}{2n+3}.}
$$
This identity is the <discrete orthogonality of Gaussian quadrature>.