\mu(x)$, $x\in \mbox{supp}(\mu)$. Set a system of nodes ${\bf x}_n=(x_{1}, \ldots, x_{n})$ such that $\displaystyle P_n(x)=\prod_{j=1}^n(x-x_{j})=q_{\mu_n,n}$. This means that ${\bf x}_n$ is the system of $n$ nodes corresponding to the Gaussian quadrature rule for the measure $\mu_n$. Given an arbitrary polynomial $p\in \Pi_{m(n)}$, we have that $$\label{supergauss}
q_{m(n)}(x)p(x)-\sum_{j=1}^n q_{m(n)}(x_{j})p(x_{j})L_{j,n}(x)=q_{\mu_n,n}(x)\mathcal{P}_{n-1}(x),$$ where $\displaystyle L_{j,n}(x):=\prod_{\underset{k\not=j}{k=1}}^n \frac{x-x_{k}}{x_{j}-x_{k}},\quad j=1,\ldots,n,$ and $ \mathcal{P}_{n-1}$ is a certain polynomial with $ \displaystyle \deg \mathcal{P}_{n-1}=n- (m(n)-\deg p)-1\leq n-1.$

Observe that $\displaystyle \int p(x) \, \mbox{d}\,  \mu (x)= \int q_{m(n)}(x)p(x) \mbox{d} \mu_n(x).$ Hence from (\[supergauss\]) we obtain $$\int p(x) \, \mbox{d}\,  \mu (x)- \int \sum_{j=1}^n q_{m(n)}(x_{j})p(x_{j})L_{j,n}(x)\,  \mbox{d}\,  \mu_n (x)=\int q_{\mu_n,n}(x)\mathcal{P}_{n-1}(x)\, \mbox{d}\, \mu_n(x),$$ which vanishes because $q_{\mu_n,n}$ satisfies the orthogonality relations for $\mu_n$ as in (\[ortoargentino\]). We conclude then $$\int p(x)\, \mbox{d}\, \mu(x) =\sum_{j=1}^n p(x_{j,n})q_{m(n)}(x_{j,n}) \int L_{j,n}(x)\frac{\mbox{d}\,  \lambda_0 (x)}{q_{m(n)}(x)}=\sum_{j=1}^n w_{j,n} p(x_{j,n}).$$ This is an interpolatory integration rule with degree of exactness $m(n)$, where the weights can be defined via $$\label{lambda}
w_{j,n}=q_{m(n)}(x_{j})\int L_{j,n}(x)\, \mbox{d} \mu_n(x)=q_{m(n)}(x_{j}) \widetilde{w}_{j,n}, \quad j=1,\ldots,n.$$ The numbers $\widetilde{w}_{j,n}$, $j=1, \ldots,n$ are the weights corresponding to a Gaussian quadrature rule, which are all positive. Since $q_{m(n)}$ is also positive the weights $w_{j,n}>0$. According to Pólya’s condition a sequence of these rules of integration is convergent.

Let us consider $\displaystyle {\bf 