The direct calculation of the derivatives isn't that hard. But one can also quickly see the values of the Riemann tensor for a sphere – and similar simple shapes – by using the definition of the Riemann tensor via the parallel transport of vectors.
$$\delta V^\alpha = R_{\alpha\beta\gamma\delta}V^\beta d\Sigma^{\gamma\delta}$$
Around a point of the sphere $S^d$, the transport around an area given by $d\Sigma^{\gamma\delta}$ for fixed values of the indices (locally orthonormal basis) allows you to see that all of this is happening on an $S^2$ only. The other dimensions are unaffected. That's why you get $R_{\alpha\beta\gamma\delta}$ equal to $(1/a^2)$ times $g_{\alpha\gamma}g_{\beta\delta}-g_{\alpha\delta}g_{\beta\gamma}$. Effectively, the antisymmetrized pair of indices $\alpha\beta$ has to be the same as the pair $\gamma\delta$. I didn't assume anything special about the point; all points on a sphere are equally good by a symmetry. So the Ansatz for the Riemann tensor has to hold everywhere.
Note that it's multiplied by $1/a^2$, the inverse squared radius of the sphere. Many of your formulae omit it; moreover, you are using a confusing symbol $R$ for the radius which looks like the Ricci scalar – a different thing.
In $d=2$, the Riemann tensor only has one independent component and the formula for the Riemann tensor in terms of the metric tensor above actually holds for any surface if $1/a^2$ is replaced by $R/2$. Note that the two-sphere has (Ricci scalar) $R=2/a^2$. Also, the Ricci tensor is $R_{ij}=Rg_{ij}/2$ in $d=2$ so that the vacuum Einstein equations are obeyed identically.
For $d=3$, the Riemann tensor has 3 independent components, just like the Ricci tensor, so the Riemann tensor may be written in terms of the Ricci tensor. That's not true for higher dimensions, either.
There is a relatively fast approach to computing the Riemann tensor, Ricci tensor and Ricci scalar given a metric tensor known as the Cartan method or method of moving frames. Given a line element,
$$ds^2 = g_{\mu\nu}dx^\mu dx^\nu$$
you pick an orthonormal basis $e^a = e^a_\mu dx^\mu$ such that $ds^2 = \eta_{ab}e^a e^b$. The first Cartan structure equation,
$$de^a + \omega^a_b \wedge e^b = 0$$
allows one to solve for the spin connection components $\omega^a_b$ from which one can compute the Ricci tensor in the orthonormal basis:
$$R^a_b = d\omega^a_b + \omega^a_c \wedge \omega^c_b.$$
The entire process simply requires exterior differentiation of the basis and spin connection. The Riemann components may be deduced from the relation,
$$R^a_b = R^a_{bcd} \, e^c \wedge e^d$$
possibly with a factor of $\frac12$ depending on your conventions. To convert back to the coordinate basis, one must simply contract with the basis back:
$$R^\mu_{\nu \lambda \kappa} = (e^{-1})^\mu_a \, R^a_{bcd}\, e^b_\nu \, e^c_\lambda \, e^d_\kappa.$$
For an explicit calculation see my previous answers here, here and here. The gravitational physics lectures at pirsa.org also provide explicit examples. As for using computer algebra systems, if all you're looking to do is compute curvature tensors, Hartle's textbook for Mathematica is your best option or the GREAT package. If you'd like to do more advanced stuff like perturbation theory, then xAct is required.
Best Answer
$g^{\alpha\beta}$ is symmetric in $\alpha$ and $\beta$, while $R_{\alpha\beta\gamma\mu}$ is anti-symmetric in $\alpha$ and $\beta$, so the contraction $g^{\alpha\beta}R_{\alpha\beta\gamma\mu}$ is necessarily $0$, and cannot be $R_{\gamma\mu}$.
Moreover, it is not correct to say, that if the contraction of $2$ tensors with another tensor (here the metric tensor) are equals, then the $2$ tensors are equal. For instance, if you take a tensor $T_1$, and a tensor $T_2$, with $T_2-T_1$ anti-symmetric in lower indices $\alpha, \beta$ and contract the tensors $T_1$ and $T_2$ with a tensor symmetric in upper indices $\alpha, \beta$, you will get the same result.