En aras de la exhaustividad, demostraré la afirmación de un modo más elemental. Empecemos por demostrar que $f^2$ es Riemann -integrable:
Sabemos que $f:[a,b]\to \mathbb{R}$ es Riemann -integrable por lo que se deduce que $|f|$ es Riemann -integrable también. Sostiene que $|f|^2=f^2$ por lo que si $\int\limits_a^b f(x)^2dx $ existe entonces tenemos la igualdad $\int\limits_a^b |f(x)|^2dx =\int\limits_a^b f(x)^2dx$ .
En $|f|$ es Riemann -integrable (y por tanto acotada) se deduce que para una $\epsilon>0$ existe una partición $P:=\{t_0,t_1\cdots, t_n\}$ de $[a,b]$ tal que:
$$ \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]\}. $$
Desde $|f([a,b])|\geq 0$ sabemos que $\sup(|f|^2)=\sup(|f|)^2$ y $\inf(|f|^2)=\inf(|f|)^2$ respectivamente.
Utilizamos estos resultados para dar un límite superior de la diferencia del Darboux -suma de $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]\}. $$
Esto demuestra que $f^2$ es Riemann -integrable. Aplicando las reglas ordinarias de Riemann -integrables (en relación con la adición de funciones y la multiplicación de factores) y la sugerencia de que $fg=\frac{1}{4}((f+g)^2-(f-g)^2)$ se ve inmediatamente que $fg$ es Riemann -integrable.