For sake of completeness I will prove the statement in a more elementary way. Let's start by showing that $f^2$ is Riemann-integrable:
We know that $f:[a,b]\to \mathbb{R}$ is Riemann-integrable so it follows that $|f|$ is Riemann-integrable as well. It holds that $|f|^2=f^2$, so if $\int\limits_a^b f(x)^2dx $ exists then we have the equality $\int\limits_a^b |f(x)|^2dx =\int\limits_a^b f(x)^2dx$.
As $|f|$ is Riemann-integrable (and hence bounded) it follows that for an arbitrary $\epsilon>0$ there exists a partitition $P:=\{t_0,t_1\cdots, t_n\}$ of $[a,b]$ such that:
$$
\sum\limits_{i=1}^{n}(M_i-m_i)(t_i-t_{i-1})<\frac{\epsilon}{2\sup(|f|)},\\
\text{where } M_i:=\sup\{|f(x)|\mid x\in[t_{i-1},t_i]\}\\
\text{and } m_i:=\inf\{|f(x)|\mid x\in[t_{i-1},t_i]\}.
$$
Since $|f([a,b])|\geq 0$, we know that $\sup(|f|^2)=\sup(|f|)^2$ and $\inf(|f|^2)=\inf(|f|)^2$, respectively.
We use these results to give an upper bound of the difference of the Darboux-sums of $f^2$:
$$
\sum\limits_{i=1}^{n}(M'_i-m'_i)(t_i-t_{i-1})=\sum\limits_{i=1}^{n}(M_i^2-m_i^2)(t_i-t_{i-1})=\sum\limits_{i=1}^{n}(M_i-m_i)(M_i+m_i)(t_i-t_{i-1})\\ \leq 2\sup(|f|)\sum\limits_{i=1}^{n}(M_i-m_i)(t_i-t_{i-1})<2\sup(|f|)\frac{\epsilon}{2\sup(|f|)}=\epsilon\\
\text{where } M'_i:=\sup\{f(x)^2\mid x\in[t_{i-1},t_i]\}\\
\text{and } m'_i:=\inf\{f(x)^2\mid x\in[t_{i-1},t_i]\}.
$$
This shows that $f^2$ is Riemann-integrable. Applying the ordinary rules of Riemann-integrable functions (regarding addition of functions and multiplying factors) and the hint that $fg=\frac{1}{4}((f+g)^2-(f-g)^2)$ one immediately sees that $fg$ is Riemann-integrable.