def _recursive_builder(self, operation, gate_name, control_qubits, target_qubit):
        """Helper function used to define the controlled gate recursively. It uses the algorithm in
         the reference above. Namely it recursively constructs a controlled gate by applying a
         controlled square root of the gate, followed by a toffoli gate with len(control_qubits) - 1
         controls, applying the controlled adjoint of the square root, another toffoli with
         len(control_qubits) - 1 controls, and finally another controlled copy of the gate.

        :param numpy.ndarray operation: The matrix for the unitary to be controlled.
        :param String gate_name: The name for the gate being controlled.
        :param Sequence control_qubits: The qubits that are the controls.
        :param Qubit or Int target_qubit: The qubit that the gate should be applied to.
        :return: The intermediate Program being built.
        :rtype: Program
        """
        control_true = np.kron(ONE_PROJECTION, operation)
        control_false = np.kron(ZERO_PROJECTION, np.eye(2, 2))
        control_root_true = np.kron(ONE_PROJECTION, sqrtm(operation, disp=True))
        controlled_gate = control_true + control_false
        controlled_root_gate = control_root_true + control_false
        sqrt_name = self.format_gate_name(SQRT_PREFIX, gate_name)
        controlled_subprogram = pq.Program()
        control_gate = pq.Program()
        # For the base case, we check to see if there is just one control qubit.
        # If it is a CNOT we explicitly break the naming convention so as not to redefine the gate.
        if len(control_qubits) == 1:
            if gate_name == NOT_GATE_LABEL:
                control_name = CONTROL_PREFIX + gate_name
            else:
                control_name = self.format_gate_name(CONTROL_PREFIX, gate_name)
            control_gate = self._defgate(control_gate, control_name, controlled_gate)
            control_gate.inst((control_name, control_qubits[0], target_qubit))
            return control_gate

        else:
            control_sqrt_name = self.format_gate_name(CONTROL_PREFIX, sqrt_name)
            control_gate = self._defgate(control_gate, control_sqrt_name, controlled_root_gate)
            control_gate.inst((control_sqrt_name, control_qubits[-1], target_qubit))
            # Here we recurse to build a toffoli gate on n - 1 of the qubits.
            n_minus_one_toffoli = self._recursive_builder(NOT_GATE,
                                                          NOT_GATE_LABEL,
                                                          control_qubits[:-1],
                                                          control_qubits[-1])

            # We recurse to build a controlled sqrt of the target_gate, excluding the last control.
            n_minus_one_controlled_sqrt = self._recursive_builder(sqrtm(operation, disp=True),
                                                                  sqrt_name,
                                                                  control_qubits[:-1],
                                                                  target_qubit)
            controlled_subprogram += control_gate
            controlled_subprogram += n_minus_one_toffoli
            controlled_subprogram += control_gate.dagger()
            # We only add the instructions so that we don't redefine gates
            controlled_subprogram += n_minus_one_toffoli.instructions
            controlled_subprogram += n_minus_one_controlled_sqrt
            return controlled_subprogram