code: def ghz_state(qubits: Union[int, Qubits]) -> State:
      
      N, qubits = qubits_count_tuple(qubits)
      ket = np.zeros(shape=[2] * N)
      ket[(0, ) * N] = 1 / sqrt(2)
      ket[(1, ) * N] = 1 / sqrt(2)
      return State(ket, qubits)
docs: Return a GHZ state on N qubits