Add GHZ circuit function

doc/source/api-reference/qibo.rst

@@ -401,6 +401,11 @@ with the last qubit as the control qubit and the first qubit as a target qubit.
 .. autofunction:: qibo.models.encodings.entangling_layer
 
 
+Greenberger-Horne-Zeilinger (GHZ) state
+"""""""""""""""""""""""""""""""""""""""
+
+.. autofunction:: qibo.models.encodings.ghz_state
+
 .. _error-mitigation:
 
 Error Mitigation

src/qibo/models/encodings.py

@@ -376,6 +376,35 @@ def entangling_layer(
     return circuit
 
 
+def ghz_state(nqubits: int, **kwargs):
+    """Generates an :math:`n`-qubit Greenberger-Horne-Zeilinger (GHZ) state that takes the form
+
+    .. math::
+        \\ket{\\text{GHZ}} = \\frac{\\ket{0}^{\\otimes n} + \\ket{1}^{\\otimes n}}{\\sqrt{2}}
+
+    where :math:`n` is the number of qubits.
+
+    Args:
+        nqubits (int): number of qubits :math:`n >= 2`.
+        kwargs (dict, optional): additional arguments used to initialize a Circuit object.
+            For details, see the documentation of :class:`qibo.models.circuit.Circuit`.
+
+    Returns:
+        :class:`qibo.models.circuit.Circuit`: Circuit that prepares the GHZ state.
+    """
+    if nqubits < 2:
+        raise_error(
+            ValueError,
+            f"nqubits given as {nqubits}. nqubits needs to be >= 2.",
+        )
+
+    circuit = Circuit(nqubits, **kwargs)
+    circuit.add(gates.H(0))
+    circuit.add(gates.CNOT(qubit, qubit + 1) for qubit in range(nqubits - 1))
+
+    return circuit
+
+
 def _generate_rbs_pairs(nqubits: int, architecture: str, **kwargs):
     """Generating list of indexes representing the RBS connections
 

tests/test_models_encodings.py

@@ -11,6 +11,7 @@
 from qibo.models.encodings import (
     comp_basis_encoder,
     entangling_layer,
+    ghz_state,
     phase_encoder,
     unary_encoder,
     unary_encoder_random_gaussian,
@@ -279,3 +280,24 @@ def test_circuit_kwargs(density_matrix):
 
     test = unary_encoder_random_gaussian(4, density_matrix=density_matrix)
     assert test.density_matrix is density_matrix
+
+
+@pytest.mark.parametrize("density_matrix", [False, True])
+@pytest.mark.parametrize("nqubits", [1, 2, 3, 4])
+def test_ghz_circuit(backend, nqubits, density_matrix):
+    if nqubits < 2:
+        with pytest.raises(ValueError):
+            GHZ_circ = ghz_state(nqubits, density_matrix=density_matrix)
+    else:
+        target = np.zeros(2**nqubits, dtype=complex)
+        target[0] = 1 / np.sqrt(2)
+        target[2**nqubits - 1] = 1 / np.sqrt(2)
+        target = backend.cast(target, dtype=target.dtype)
+
+        GHZ_circ = ghz_state(nqubits, density_matrix=density_matrix)
+        state = backend.execute_circuit(GHZ_circ).state()
+
+        if density_matrix:
+            target = backend.np.outer(target, backend.np.conj(target.T))
+
+        backend.assert_allclose(state, target)