To avoid blackening the place, I will not use bold symbols -but the answer will be carried out in matrix notationform. Uppercase will denote matrices, lower case will denoteVectors are column vectors, a prime will denote the transpose.
Let a linear regression model
$$y = X_1b_1 + X_2b_2 + u_A \qquad [A]$$
The normal equations for the OLS estimator are
$$\begin{align} \left(X_1'X_1\right)b_1+\left(X_1'X_2\right)b_2=& X_1'y \qquad [1]\\ \\ \left(X_2'X_1\right)b_1+\left(X_2'X_2\right)b_2=& X_2'y \qquad [2]\\ \end{align}$$ Solving $[2]$ for $b_2$ we have $$[2]\rightarrow b_2= \left(X_2'X_2\right)^{-1}X_2'y-\left(X_2'X_2\right)^{-1}\left(X_2'X_1\right)b_1$$
Inserting this into $[1]$ we obtain $$\left(X_1'X_1\right)b_1+\left(X_1'X_2\right)\left(X_2'X_2\right)^{-1}X_2'y-\left(X_1'X_2\right)\left(X_2'X_2\right)^{-1}\left(X_2'X_1\right)b_1= X_1'y $$
Collecting terms w.r.t $b_1$ and $y$, $$X_1'\left[I-X_2\left(X_2'X_2\right)^{-1}X_2'\right]X_1b_1= X_1'\left[I-X_2\left(X_2'X_2\right)^{-1}X_2'\right]y$$ $$\Rightarrow X_1'M_2X_1b_1 = X_1'M_2y \qquad [3]$$
where $M_2$ is the "annihilator" or "residual maker"matrix related to $X_2$, namely the matrix that produces the residuals when a variable is regressed on $X_2$, by pre-multiplying this variable. This matrix is symmetric and idempotent, $M_2=M_2',\; M_2= M_2M_2$. So we can write
$$(M_2X_1)'(M_2X_1)b_1 = (M_2X_1)'y$$ $$\Rightarrow R_{1\sim2}'R_{1\sim2}b_1=R_{1\sim2}'y \Rightarrow \hat b_1 = \left(R_{1\sim2}'R_{1\sim2}\right)^{-1}R_{1\sim2}'y\qquad [4]$$
where $R_{1\sim2}$ denotes the residual vector from regressing $X_1$ on $X_2$.
This last formula is exactly the OLS formula from the regression model $$y= R_{1\sim2}d_1+u_B \qquad [B]$$
So eq. $[4]$ tells us that the coefficient estimate for $X_1$ that we will obtain in a multiple regression setting, will be exactly the same with what we will obtain if we regress the dependent variable on the residuals from the regression of $X_1$ on $X_2$.
Now consider the second case, regressing the residuals on the residuals. This is the model
$$R_{y\sim2} = R_{1\sim2}c_1+u_C \Rightarrow (M_2y)= (M_2X_1)c_1 +u_C \qquad [C]$$
The OLS estimator of $c$ is $$\hat c_1 = \left[(M_2X_1)'(M_2X_1)\right]^{-1}(M_2X_1)'(M_2y) \qquad [5]$$
By the properties of $M_2$ we have $$(M_2X_1)'(M_2y) = X_1'M_2'M_2y=X_1'M_2M_2y=X_1'M_2y=X_1'M_2'y=(M_2X_1)'y$$ Noting that $M_2X_1 = R_{1\sim2}$ eq. $[5]$ becomes
$$\hat c_1= \left(R_{1\sim2}'R_{1\sim2}\right)^{-1}R_{1\sim2}'y \qquad [6]$$
which is identical to eq. $[4]$, and so $ \hat c_1 = \hat d_1 =\hat b_1$. In other words the three models give mathematically identical results.
Let's now consider the issue of the estimator variance. Models $[B]$ and $[C]$ have the same regressor matrix so the question is what happens with the estimated error variances, $\sigma^2_B$ and $\sigma^2_C$. We will denote $M(r)_{1\sim2}$ the annihilator matrix of the regressor $R_{1\sim2}$. It has analogous properties as $M_2$ For model $[B]$ we have
$$u'_Bu_B = \left(M(r)_{1\sim2}y\right)'\left(M(r)_{1\sim2}y\right) = y'M(r)_{1\sim2}y \qquad [7]$$
while for model $[C]$ we have
$$u'_Cu_C = \left(M(r)_{1\sim2}(M_2y)\right)'\left(M(r)_{1\sim2}(M_2y)\right) = y'M_2M(r)_{1\sim2}M_2y \qquad [8]$$
Are the RHS of eq. $[7]$ and $[8]$ equal? I 'll leave that to the reader.