code: def decomposed_diffusion_program(qubits: List[int]) -> Program:
      
      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
docs: 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]>``.