Comment by papa2fire

Following your suggestion, I’ve implemented the PSA algorithm (the optimized one) in C++. It might be tricky to understand without an explanation, especially since (1) it’s generalized for any dimensionality, and (2) I manually implemented the basic transformation ("oscillating shift") instead of using NumPy’s reshape+roll+flatten, which I believe makes it easier to grasp. But you can verify that it returns the same results as the naive algorithm. For the quiniela case, it shouldn’t take more than 2-3 seconds to compute the sum of neighbors for all points. https://github.com/petopello/PSA/blob/main/PSA.cpp