Gauss Quadrature Proof. Let $x_0<x_1< ... < x_n$ be the roots of an n+1 degree orthonormal polynomial $\phi_{n+1}$ with respect to the inner product:
$$\langle g,h \rangle = \int_a^bw(x)g(x)h(x)dx$$
and $$p_n= \sum_{j =  0}^nf(x_j)L_{n,j} \in \mathbb{P_n}$$
be the lagrange interpolating polynomial for the given data.
Then
$$ \int_a^bw(x)f(x) dx \approx \sum_{j = 0}^n \int^b_a w(x)L_{n,j}(x)f(x_j) $$
Proof:(The bits in red I have added to try and understand it)
Let $f \in \mathbb{P_{2n+1}}$ such that $f(x) = q(x) \phi_{n+1} + r(x). Then;$
$$\int^b_aw(x)f(x)dx = \color{red}{\int^b_aw(x)q(x)\phi_{n+1}dx + \int^b_aw(x)r(x)dx} = \color{green}{\sum^n_{j=0}\int^b_a w(x)L_{n,j}(x)r(x_j)}$$ 
The first red integral I know is zero, since $\langle\phi_{n+1},q_n\rangle = 0; \forall q_n \in \mathbb{P_n} $ the second integral is said to be "computed exactly". $\color{red}{\text{why is that so?}}$
 A: Antwon(theswan),
As you say, let $f \in \mathbb{P_{2n+1}}$; then $f$ factors by long division into $f(x) = q(x) \phi_{n+1} + r(x)$. But then $f(x)$ is of degree $2n+1$, and $\phi_{n+1}$ is of degree $n+1$ and $q(x),r(x)$ are of degree at most $n$.
Now we take the $x_j$ to be the zeroes of $\phi_{n+1}$, so you have $f(x_j) = r(x_j)$ for $j = 0, \dots, n$.
I think you made it at least this far since you understand the factoring of the integral and understand the orthogonality principle behind the first term of the factorization being $0$.
Now the second integral is where the "clever" factoring comes into play.
We have:
$\int^b_a w(x)f(x)dx =^* \int^b_a w(x)r(x) =^{**} \color{red}{\sum^n_{j=0} P_j r(x_j)} = \sum^n_{j=0} \int^b_a w(x) L_{n,j}(x)r(x_j)$
where the first equivalence ($*$) comes from the orthogonality (which you recognize) $\langle\phi_{n+1},q_n\rangle = 0$ and the exactness of ($**$), which pertains to your question, comes from the Newton-Cotes formulas, which tell us that ($**$) is exact for any set of roots ($x_j$) for polynomials of degree at most $n$. And since we have chosen our roots wisely, i.e. the ($x_j$) of $\phi_{n+1}$, we can use what we deduced earlier, that $f(x_j) = r(x_j)$, which leads to the result.
I think your hang-up was due to the fact that the second integral in red is of degree at most $n$. For this reason, Newton-Cotes can be applied to draw exactness. This concept should jive with the idea of Gauss 'improving' the exactness of quadrature rules for polynomials up to $2n+1$ from the Newton-Cotes formulas, and this proof that you cite is using the fact that we already know exactness holds for one of these polynomials up to degree $n$.
NOTE: I highlighted the formula in red because I am having difficulty deciphering your notation. Explicitly, $\int^b_a w(x)f(x) dx \approx \sum^n_{j=0} P_j f(x_j)$ will be exact, if $P_j = \int^b_a w(x) p_j(x) dx$ for $j = 0, \dots, n$ and the $p_j$ are the Lagrange interpolating polynomials.
