[proofplan]
We use the Darboux criterion for [Riemann integrability](/page/Riemann%20Integrability). Since $f$ and $g$ are Riemann integrable, they are bounded and admit partitions whose upper-lower Darboux gaps are arbitrarily small. The key estimate is that the oscillation of $fg$ on each subinterval is bounded by a boundedness constant for $g$ times the oscillation of $f$, plus a boundedness constant for $f$ times the oscillation of $g$. Summing this estimate over a sufficiently fine common refinement gives an arbitrarily small Darboux gap for $fg$.
[/proofplan]
[step:Choose boundedness constants and Darboux partitions for $f$ and $g$]
If $a=b$, then $[a,b]=\{a\}$ and every [bounded function](/page/Bounded%20Function) on $[a,b]$ has zero upper-lower Darboux gap for the partition $\{a\}$. Hence $fg$ is Riemann integrable on $[a,b]$ by the Darboux criterion. For the rest of the proof, assume $a<b$.
Because $f$ and $g$ are Riemann integrable on $[a,b]$, both functions are bounded. Define constants $A,B \in [0,\infty)$ by
\begin{align*}
A := \sup_{x \in [a,b]} |f(x)|
\end{align*}
and
\begin{align*}
B := \sup_{x \in [a,b]} |g(x)|.
\end{align*}
For a bounded function $h: [a,b] \to \mathbb{R}$ and a partition $Q=\{q_0,q_1,\dots,q_m\}$ with $a=q_0<q_1<\cdots<q_m=b$, define the upper and lower Darboux sums by
\begin{align*}
U(h,Q) := \sum_{j=1}^{m}\left(\sup_{x \in [q_{j-1},q_j]} h(x)\right)(q_j-q_{j-1})
\end{align*}
and
\begin{align*}
L(h,Q) := \sum_{j=1}^{m}\left(\inf_{x \in [q_{j-1},q_j]} h(x)\right)(q_j-q_{j-1}).
\end{align*}
Let $\varepsilon > 0$. By the Darboux criterion for Riemann integrability, there are partitions $P_f$ and $P_g$ of $[a,b]$ such that
\begin{align*}
U(f,P_f)-L(f,P_f) < \frac{\varepsilon}{2(B+1)}
\end{align*}
and
\begin{align*}
U(g,P_g)-L(g,P_g) < \frac{\varepsilon}{2(A+1)}.
\end{align*}
Let $P$ be the common refinement obtained by taking the union of the partition points of $P_f$ and $P_g$. Refinement decreases upper Darboux sums and increases lower Darboux sums, so
\begin{align*}
U(f,P)-L(f,P) \leq U(f,P_f)-L(f,P_f) < \frac{\varepsilon}{2(B+1)}
\end{align*}
and
\begin{align*}
U(g,P)-L(g,P) \leq U(g,P_g)-L(g,P_g) < \frac{\varepsilon}{2(A+1)}.
\end{align*}
[/step]
[step:Control the oscillation of the product on each subinterval]
Write the partition $P$ as
\begin{align*}
P = \{x_0,x_1,\dots,x_n\}
\end{align*}
with
\begin{align*}
a=x_0 < x_1 < \cdots < x_n=b.
\end{align*}
For each $i \in \{1,\dots,n\}$, define the closed subinterval
\begin{align*}
I_i := [x_{i-1},x_i].
\end{align*}
For a bounded function $h: [a,b] \to \mathbb{R}$, define its oscillation on $I_i$ by
\begin{align*}
\operatorname{osc}_{I_i}(h) := \sup_{x \in I_i} h(x)-\inf_{x \in I_i} h(x).
\end{align*}
For any $x,y \in I_i$, we have
\begin{align*}
|(fg)(x)-(fg)(y)| = |f(x)g(x)-f(y)g(y)|.
\end{align*}
Adding and subtracting $f(y)g(x)$ and applying the triangle inequality gives
\begin{align*}
|f(x)g(x)-f(y)g(y)| \leq |g(x)|\,|f(x)-f(y)|+|f(y)|\,|g(x)-g(y)|.
\end{align*}
Since $|f(y)| \leq A$ and $|g(x)| \leq B$, and since $x,y \in I_i$, this implies
\begin{align*}
|(fg)(x)-(fg)(y)| \leq B\,\operatorname{osc}_{I_i}(f)+A\,\operatorname{osc}_{I_i}(g).
\end{align*}
For every bounded function $h: I_i \to \mathbb{R}$ on the nonempty interval $I_i$, the oscillation satisfies
\begin{align*}
\operatorname{osc}_{I_i}(h)=\sup_{x,y \in I_i}|h(x)-h(y)|.
\end{align*}
Applying this identity to $h=fg$ and taking the supremum over all $x,y \in I_i$ gives
\begin{align*}
\operatorname{osc}_{I_i}(fg) \leq B\,\operatorname{osc}_{I_i}(f)+A\,\operatorname{osc}_{I_i}(g).
\end{align*}
[guided]
The product $fg$ can vary on $I_i$ for two reasons: $f$ can vary while $g$ stays bounded, and $g$ can vary while $f$ stays bounded. The estimate separates these two effects.
Fix $i \in \{1,\dots,n\}$ and take arbitrary points $x,y \in I_i$. We compare $(fg)(x)$ and $(fg)(y)$ by inserting the intermediate term $f(y)g(x)$:
\begin{align*}
f(x)g(x)-f(y)g(y) = g(x)(f(x)-f(y))+f(y)(g(x)-g(y)).
\end{align*}
Taking absolute values and using the triangle inequality, we obtain
\begin{align*}
|f(x)g(x)-f(y)g(y)| \leq |g(x)|\,|f(x)-f(y)|+|f(y)|\,|g(x)-g(y)|.
\end{align*}
The constants $A$ and $B$ were chosen so that $|f(z)| \leq A$ and $|g(z)| \leq B$ for every $z \in [a,b]$. Therefore
\begin{align*}
|f(x)g(x)-f(y)g(y)| \leq B\,|f(x)-f(y)|+A\,|g(x)-g(y)|.
\end{align*}
Because $x,y \in I_i$, the differences $|f(x)-f(y)|$ and $|g(x)-g(y)|$ are bounded by the corresponding oscillations on $I_i$. Hence
\begin{align*}
|f(x)g(x)-f(y)g(y)| \leq B\,\operatorname{osc}_{I_i}(f)+A\,\operatorname{osc}_{I_i}(g).
\end{align*}
Since this holds for every pair $x,y \in I_i$, taking the supremum over all such pairs gives
\begin{align*}
\operatorname{osc}_{I_i}(fg) \leq B\,\operatorname{osc}_{I_i}(f)+A\,\operatorname{osc}_{I_i}(g).
\end{align*}
This is the central estimate: boundedness controls the size of the factors, while integrability controls the total oscillation of each factor across the partition.
[/guided]
[/step]
[step:Sum the oscillation estimate to make the Darboux gap of $fg$ small]
For each $i \in \{1,\dots,n\}$, set
\begin{align*}
\Delta x_i := x_i-x_{i-1}.
\end{align*}
Using the definition of upper and lower Darboux sums in terms of oscillation, we have
\begin{align*}
U(fg,P)-L(fg,P)=\sum_{i=1}^{n} \operatorname{osc}_{I_i}(fg)\,\Delta x_i.
\end{align*}
Applying the oscillation estimate from the previous step on each $I_i$ gives
\begin{align*}
U(fg,P)-L(fg,P) \leq \sum_{i=1}^{n}\bigl(B\,\operatorname{osc}_{I_i}(f)+A\,\operatorname{osc}_{I_i}(g)\bigr)\Delta x_i.
\end{align*}
Distributing the finite sum,
\begin{align*}
U(fg,P)-L(fg,P) \leq B\bigl(U(f,P)-L(f,P)\bigr)+A\bigl(U(g,P)-L(g,P)\bigr).
\end{align*}
By the choice of $P$,
\begin{align*}
U(fg,P)-L(fg,P) < B\frac{\varepsilon}{2(B+1)}+A\frac{\varepsilon}{2(A+1)}.
\end{align*}
Since $B/(B+1) \leq 1$ and $A/(A+1) \leq 1$, it follows that
\begin{align*}
U(fg,P)-L(fg,P) < \varepsilon.
\end{align*}
[/step]
[step:Apply the Darboux criterion to conclude Riemann integrability]
We have shown that for every $\varepsilon > 0$ there exists a partition $P$ of $[a,b]$ such that
\begin{align*}
U(fg,P)-L(fg,P) < \varepsilon.
\end{align*}
The function $fg: [a,b] \to \mathbb{R}$ is bounded because $|f(x)g(x)| \leq AB$ for every $x \in [a,b]$. Therefore the Darboux criterion implies that $fg$ is Riemann integrable on $[a,b]$.
[/step]
