code: def productdiag(tensor: BKTensor) -> BKTensor:
        # DOCME: Explain
      N = rank(tensor)
      tensor = reshape(tensor, [2**(N//2), 2**(N//2)])
      tensor = np.diag(tensor)
      tensor = reshape(tensor, [2]*(N//2))
      return tensor
docs: Returns the matrix diagonal of the product tensor