Do variations obey the product rule? I have been trying to derive the Einstein equation from the Einstein-Hilbert action
$$ S[g_{\mu \nu}] = \frac{1}{16 \pi} \int_M \text{d}^4x \sqrt{-g}R $$
The standard derivation states that the variation $\delta S =0$ when we vary the metric components. In this derivation, we use the fact that $\delta (\sqrt{-g}R)= \delta \sqrt{-g} R + \sqrt{-g}\delta R$. To me this seems quite obvious but I thought I would try and prove this.
My understanding is that if we have an action
$$ S[q] = \int \text{d}t L(q,\dot{q},t)$$
we vary it as 
$$ \delta S[q] = \int \text{d}t \delta L(q,\dot{q},t)$$
so 
$$ \delta L= \frac{\partial L}{\partial q} \delta q + \frac{\partial L}{\partial \dot{q}} \delta \dot{q} \\
= \left( \frac{\partial L}{\partial q} - \frac{\text{d}}{\text{d}t} \frac{\partial L}{\partial \dot{q}} \right) \delta q + \frac{\text{d}}{\text{d}t} \left( \frac{\partial L}{\partial \dot{q}} \delta q \right)$$
As we integrate over $\delta L$, we ignore the last term as this produces a boundary term which we can take to vanish, therefore I say
$$ \delta L = \left( \frac{\partial L}{\partial q} - \frac{\text{d}}{\text{d}t} \frac{\partial L}{\partial \dot{q}} \right) \delta q $$
Okay, now let's say my Lagrangian is a product: $L(q,\dot{q},t) = f(q,\dot{q},t)g(q,\dot{q},t)$. Plugging this into the above formula for the variation, I have
$$ \frac{\partial L}{\partial q} = \frac{\partial f}{\partial q} g + \frac{\partial g }{\partial q} f $$
$$ \frac{\partial L}{\partial \dot{q}} = \frac{\partial f}{\partial \dot{q}}g + \frac{\partial g}{\partial \dot{q}} f$$
$$ \frac{\text{d}}{\text{d}t} \frac{\partial L}{\partial \dot{q}} = \left( \frac{\text{d}}{\text{d}t} \frac{\partial f}{\partial \dot{q}} \right) g + \left( \frac{\text{d}}{\text{d}t}\frac{\partial g}{\partial \dot{q}} \right) f + \frac{\partial f}{\partial \dot{q}} \frac{\text{d}g}{\text{d}t} + \frac{\partial q}{\partial \dot{q}} \frac{\text{d}f}{\text{d}t}$$
so the variation of this Lagrangian is 
$$ \delta L = \left( \frac{\partial f}{\partial q} - \frac{\text{d}}{\text{d}t} \frac{\partial f}{\partial \dot{q}} \right)g \delta q + \left( \frac{\partial g}{\partial q} - \frac{\text{d}}{\text{d}t} \frac{\partial g}{\partial \dot{q}} \right)f \delta q -  \frac{\partial f}{\partial \dot{q}} \frac{\text{d}g}{\text{d}t}\delta q - \frac{\partial q}{\partial \dot{q}} \frac{\text{d}f}{\text{d}t}\delta q$$
or 
$$ \delta L = g \delta f + f \delta g - \frac{\partial f}{\partial \dot{q}} \frac{\text{d}g}{\text{d}t}\delta q - \frac{\partial g}{\partial \dot{q}} \frac{\text{d}f}{\text{d}t}\delta q $$
Now I can't seem to get rid of those horrible extra terms - I can't see how they would produce a boundary term when integrated. Maybe my understanding was incorrect and variations do not obey the product rule? Many standard resources suggest that varying the Einstein-Hilbert action obeys the product rule... is this just an exception?

My question:
How can I show that variations obey $\delta (fg) = f \delta g + g \delta f$
 A: It is true that
$$ \delta L =  f \delta g  + g \delta f,$$
but the expressions for variations are:
\begin{align} \delta f&=\frac{\partial f}{\partial q} \delta q+ \frac{\partial f}{\partial \dot q}\delta \dot q = \left[ \frac{\partial f}{\partial q} - \frac{d}{dt}\left(\frac{\partial f}{\partial \dot q} \right) \right] \delta q \ + \ \frac{d}{dt}\left(  \frac{\partial f}{\partial \dot q} \delta q\right)\\ \delta g&=\frac{\partial g}{\partial q} \delta q+ \frac{\partial g}{\partial \dot q}\delta \dot q = \left[ \frac{\partial g}{\partial q} - \frac{d}{dt}\left(\frac{\partial g}{\partial \dot q} \right) \right] \delta q \ + \ \frac{d}{dt}\left(  \frac{\partial g}{\partial \dot q} \delta q\right)\\\delta L&=\frac{\partial L}{\partial q} \delta q+ \frac{\partial L}{\partial \dot q}\delta \dot q = \left[ \frac{\partial f}{\partial q}g+\frac{\partial g}{\partial q}f - \frac{d}{dt}\left(\frac{\partial f}{\partial \dot q}g + \frac{\partial g}{\partial \dot q}f \right) \right] \delta q \ + \ \frac{d}{dt}\left(  \frac{\partial f}{\partial \dot q}g \delta q + \frac{\partial g}{\partial \dot q}f \delta q\right)\end{align}
In each case, the variation is equal to the Euler-Lagrange equation times $\delta q$, plus a total derivative. The total derivatives are important!
