My manual implementation was originally intended to be used as a reference implementation to see, how the running estimates are updated and how the batchnorm layer works in general.
The performance difference is expected, as I used native methods without specific (and faster) kernels.
The difference in after a couple of iterations is most likely created due to the accumulation of rounding errors created by the limited floating point precision.
It’s also an indication for this, since you see a reduction in these errors using float64.
What’s your use case, that you would like to use this implementation in a real model?