code: def _P(m):
      
      P = np.zeros((m**2,m**2), dtype=np.int64)
      for i in range(1, m**2 + 1):
          j = 1 + m*((i - 1) % m) + (i - 1)//m
          P[i-1, j-1] = 1
      return P
docs: Returns m^2 x m^2 permutation matrix that swaps rows i and j where
      j = 1 + m((i - 1) mod m) + (i - 1) div m, for i = 1 .. m^2