docs: def symorth(S):
      "Symmetric orthogonalization"
      E,U = np.linalg.eigh(S)
      n = len(E)
      Shalf = np.identity(n,'d')
      for i in range(n):
          Shalf[i,i] /= np.sqrt(E[i])
      return simx(Shalf,U,True)