def decomposed_diffusion_program(qubits: List[int]) -> Program:
    """
    Constructs the diffusion operator used in Grover's Algorithm, acted on both sides by an
    a Hadamard gate on each qubit. Note that this means that the matrix representation of this
    operator is diag(1, -1, ..., -1). In particular, this decomposes the diffusion operator, which
    is a :math:`2**{len(qubits)}\times2**{len(qubits)}` sparse matrix, into
     :math:`\mathcal{O}(len(qubits)**2) single and two qubit gates.

    See C. Lavor, L.R.U. Manssur, and R. Portugal (2003) `Grover's Algorithm: Quantum Database
    Search`_ for more information.

    .. _`Grover's Algorithm: Quantum Database Search`: https://arxiv.org/abs/quant-ph/0301079

    :param qubits: A list of ints corresponding to the qubits to operate on.
                   The operator operates on bistrings of the form
                   ``|qubits[0], ..., qubits[-1]>``.
    """
    program = Program()
    if len(qubits) == 1:
        program.inst(Z(qubits[0]))
    else:
        program.inst([X(q) for q in qubits])
        program.inst(H(qubits[-1]))
        program.inst(RZ(-np.pi, qubits[0]))
        program += (ControlledProgramBuilder()
                              .with_controls(qubits[:-1])
                              .with_target(qubits[-1])
                              .with_operation(X_GATE)
                              .with_gate_name(X_GATE_LABEL).build())
        program.inst(RZ(-np.pi, qubits[0]))
        program.inst(H(qubits[-1]))
        program.inst([X(q) for q in qubits])
    return program