The previous code was entirely scalar and quite messy. By rewriting the
code in terms of Eigen3 vectors it becomes more expressive and clear. It
also becomes much faster as we can now leverage SIMD implicitly through
Eigen's internal optimisation engine.

Note: This commit may introduce very slightly different numerical
results. This is not a bug, rather it is a consequence of Eigen's
internal optimisation heuristics and the way IEEE 754 floating-point
mathematics works.