Nice. I am not surprised that the approach presented here was discussed in Julia circles already some 5 years ago. So it was about half as fast, if I interpret the thread correctly.
I'd be very happy if this can of worms could be reopened, but am currently not active enough in julia dev to champion it.
Also somebody would need figure out something very clever for AVX2 / NEON. (AVX512 has the required instructions)
Also I can't imagine the mess with GPU -- if rand statistics differ widely between CPU and GPU that's a no-go, and I don't know what works well on which GPUs.