As far as I know, it is only the fully covariant form $R_{ijkl}$ or the fully contravariant form $R^{ijkl}$ that have the nice symmetries reducing the number of independent components to 21. (And then the algebraic Bianchi identity reduces the number to 20.)
These tensors are antisymmetric in the first two indices, antisymmetric in the last two indices, and symmetric when the first two are swapped with the last two.
This means that to get the independent components, you can first restrict $ij$ and $kl$ to be the six pairs 01, 02, 03, 12, 13, and 23. You don't have to consider 00, 11, 22, and 33 because these components are zero by antisymmetry. And you don't have to consider, say, 10, because you can swap it to get 01, with a minus sign, again by antisymmetry.
Next, by the symmetry of swapping the first two with the last two, you can restrict $ij$ to be "less than or equal to" $kl$, giving the following 21 components:
$$\begin{matrix}
0101 && 0102 && 0103 && 0112 && 0113 && 0123 \\
&& 0202 && 0203 && 0212 && 0213 && 0223 \\
&& && 0303 && 0312 && 0313 && 0323 \\
&& && && 1212 && 1213 && 1223 \\
&& && && && 1313 && 1323 \\
&& && && && && 2323
\end{matrix}$$
The algebraic Bianchi identity in four dimensions is a single relation between the three components in which all the indices are different:
$$R_{0123}-R_{0213}+R_{0312}=0$$
so you can consider any one of these three to be redundant and not calculate it "the hard way".
So I think you need to give up on the idea that you are only interested in ${R^i}_{jkl}$. You should instead calculate the 20 components $R_{ijkl}$ above, use the Bianchi identity to get the 21st, and then use symmetry and antisymmetry to get the remaining 236. For example,
$$R_{2103}=-R_{0312}$$
because starting with 2103 you can reverse the first two to get 1203, introducing a minus sign because of the antisymmetry, and then swap the first two with the last two to get 0312, which is in the list above.
So, to calculate a component of $R_{ijkl}$, do you have to first compute four components of ${R^i}_{jkl}$ using the standard formula and then contract them with the metric tensor? No! There is a formula (from Wikpedia's List of formulas in Riemannian geometry) which you can use to directly calculate fully-covariant components:
$$R_{iklm}=\frac{1}{2}\left(\frac{\partial^2 g_{im}}{\partial x^k \partial x^l}+\frac{\partial^2 g_{kl}}{\partial x^i \partial x^m}-\frac{\partial^2 g_{il}}{\partial x^k \partial x^m}-\frac{\partial^2 g_{km}}{\partial x^i \partial x^l}\right)+g_{np}({\Gamma^n}_{kl}{\Gamma^p}_{im}-{\Gamma^n}_{km}{\Gamma^p}_{il})$$
So a good strategy would be to
Calculate first derivatives of the metric components, remembering that the metric is symmetric, and cache them.
Use the cached first derivatives to calculate the Christoffel symbols, remembering that they are symmetric on the lower indices, and cache them.
Use the cached first derivatives to calculate second derivatives of the metric, remembering that they are symmetric in the partials, and cache them.
Calculate the 20 independent fully-covariant components of the Riemann tensor using the formula above.
Calculate everything else from these.
If you can show that the Einstein tensor of the Kerr metric vanishes, that will be a good test case.