[solution]
**Step 1: Establish that $h := g - f$ is zero except on a null set.**
Define $h: [a, b] \to \mathbb{R}$ by $h(x) := g(x) - f(x)$. By hypothesis, $h(x) = 0$ for all $x \in [a, b] \setminus E$, where $E$ has Lebesgue measure zero. The function $h$ is bounded because both $f$ and $g$ are bounded; let $M := \sup_{[a,b]} |h|$.
**Step 2: Show $h \in \mathcal{R}([a, b])$ with integral zero using the Darboux criterion.**
Let $\varepsilon > 0$. Since $E$ has Lebesgue measure zero, there exists a countable collection of open intervals $\{(a_k, b_k)\}_{k=1}^\infty$ with $E \subseteq \bigcup_{k=1}^\infty (a_k, b_k)$ and $\sum_{k=1}^\infty (b_k - a_k) < \varepsilon / (2M)$. By compactness of $[a, b]$, finitely many of these intervals cover $E \cap [a, b]$, say $(a_1, b_1), \ldots, (a_N, b_N)$ with $\sum_{k=1}^N (b_k - a_k) < \varepsilon / (2M)$.
**Step 3: Construct a partition subordinate to the cover.**
Let $P$ be any partition of $[a, b]$ whose points include all endpoints $a_1, b_1, \ldots, a_N, b_N$ (after intersecting with $[a, b]$). Separate the subintervals of $P$ into two groups: those $I_j$ entirely contained in some $(a_k, b_k)$, and those $I_j$ that are disjoint from $E$. On the latter, $h = 0$ identically, so $M_j - m_j = 0$. On the former, $M_j - m_j \le 2M$. The total length of the former subintervals is at most $\sum_{k=1}^N (b_k - a_k) < \varepsilon / (2M)$. Therefore
\begin{align*}
U(h, P) - L(h, P) &= \sum_{\text{near } E} (M_j - m_j) \Delta x_j + \sum_{\text{away from } E} 0 \cdot \Delta x_j \le 2M \cdot \frac{\varepsilon}{2M} = \varepsilon.
\end{align*}
By the Darboux criterion, $h \in \mathcal{R}([a, b])$. Since $0 \le U(h, P)$ and $L(h, P) \le 0$ can be arranged (the function $h$ takes value $0$ on most of $[a, b]$), passing to the limit gives $\int_a^b h\, dx = 0$.
**Step 4: Conclude by linearity.**
Since $g = f + h$ and both $f$ and $h$ are in $\mathcal{R}([a, b])$, linearity of the Riemann integral gives
\begin{align*}
\int_a^b g(x)\, dx = \int_a^b f(x)\, dx + \int_a^b h(x)\, dx = \int_a^b f(x)\, dx + 0 = \int_a^b f(x)\, dx.
\end{align*}
This result shows that the Riemann integral "does not see" modifications on null sets — a property shared with, and ultimately formalised more completely by, the Lebesgue integral.
[/solution]