def vectorized_interaction(s, alpha, kprop):
    im = complex(0, 1)
I = np.array([[1, 0, 0], [0, 1, 0], [0, 0, 1]])
k2 = kprop * kprop
sd = s[:, (None)] - s
Rv = np.sqrt((sd ** 2).sum(2))
nv = sd / Rv[:, :, (None)]
kRv = Rv * kprop
Av = 1.0 / (kRv * kRv) - im / kRv
nxnv = nv[:, :, :, (None)] * nv[:, :, (None), :]
P = (3 * Av[:, :, (None), (None)] - 1) * nxnv + (1 - Av[:, :, (None), (None)]
    ) * I
multv = -alpha * (k2 * np.exp(im * kRv)) / Rv
outv = P * multv[:, :, (None), (None)]
outv[np.eye(N, dtype=bool)] = I
return outv.transpose(0, 2, 1, 3).reshape(N * 3, -1)
Cython speedup isn't as large as expected