Add GHZ circuit function

doc/source/api-reference/qibo.rst

@@ -129,6 +129,13 @@ and the second will act to ``(1, 2)`` corresponding to
 
     [Y(1), Z(2), CNOT(1, 2), H(1), H(2)]
 
+Greenberger-Horne-Zeilinger (GHZ) state
+"""""""""""""""""""""""""""""""""""""""
+
+.. autoclass:: qibo.models.encodings.ghz_state
+    :members:
+    :member-order: bysource
+
 .. _applicationspecific:
 
 Quantum Fourier Transform (QFT)

src/qibo/models/encodings.py

@@ -470,3 +470,31 @@ def _parametrized_two_qubit_gate(gate, q0, q1, params=None):
         return gate(q0, q1, *params)
 
     return gate(q0, q1)
+
+
+def ghz_state(nqubits: int, **kwargs):
+    """Generates an :math:`n`-qubits Greenberger-Horne-Zeilinger (GHZ) state that takes the form
+
+    .. math::
+        |GHZ\\rangle_n = \\frac{|0\\rangle ^{\\otimes n} + |1\\rangle ^{\\otimes n}}{\\sqrt{2}}
+
+    where :math:`n` is the number of qubits.
+
+    Args:
+        nqubits (int): number of qubits, nqubits >= 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`: GHZ state
+    """
+    if nqubits < 2:
+        raise_error(
+            ValueError,
+            f"nqubits given as {nqubits}. nqubits needs to be >= 2.",
+        )
+    else:
+        circuit = Circuit(nqubits)
+        circuit.add(gates.H(0))
+        circuit.add(gates.CNOT(qubit, qubit + 1) for qubit in range(nqubits - 1))
+    return circuit

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,25 @@ 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))
+            state = backend.np.outer(state, backend.np.conj(state.T))
+
+        backend.assert_allclose(state, target)