ABC matrices for multi-controlled unitary

qiskit/circuit/add_control.py

@@ -16,6 +16,17 @@
 
 from qiskit.circuit.exceptions import CircuitError
 from qiskit.extensions import UnitaryGate
+from qiskit.circuit.library.standard_gates import (
+    MCXRecursive,
+    RZGate,
+    RYGate,
+    CXGate,
+    CCXGate,
+    C3XGate,
+    C4XGate,
+    MCXGate,
+)
+from qiskit.quantum_info import OneQubitEulerDecomposer
 from . import ControlledGate, Gate, QuantumRegister, QuantumCircuit
 
 
@@ -188,19 +199,41 @@ def control(
                             lamb, q_control, q_target[bit_indices[qargs[0]]], use_basis_gates=True
                         )
                     else:
-                        controlled_circ.mcrz(
-                            lamb, q_control, q_target[bit_indices[qargs[0]]], use_basis_gates=True
-                        )
-                        controlled_circ.mcry(
-                            theta,
-                            q_control,
-                            q_target[bit_indices[qargs[0]]],
-                            q_ancillae,
-                            use_basis_gates=True,
-                        )
-                        controlled_circ.mcrz(
-                            phi, q_control, q_target[bit_indices[qargs[0]]], use_basis_gates=True
+                        mcx_control, abc_control = define_mcx_control_and_ancilla(q_control)
+
+                        mcxr_rule = define_mcx_rule(
+                            mcx_control, q_target[bit_indices[qargs[0]]], abc_control
                         )
+
+                        # Getting euler angles from zyz decomposition
+                        th, ph, lb, alpha = OneQubitEulerDecomposer._params_zyz(gate.to_matrix())
+
+                        local_phase = alpha
+                        if local_phase:
+                            theta = th
+                            phi = ph
+                            lamb = lb
+
+                        a, b, c = get_abc_matrices(phi, theta, lamb)
+
+                        controlled_circ.unitary(c, q_target[bit_indices[qargs[0]]])
+                        if abc_control is not None:
+                            controlled_circ.control(abc_control)
+
+                        controlled_circ.data += mcxr_rule
+
+                        controlled_circ.unitary(b, q_target[bit_indices[qargs[0]]])
+                        if abc_control is not None:
+                            controlled_circ.control(abc_control)
+
+                        controlled_circ.data += mcxr_rule
+
+                        controlled_circ.unitary(a, q_target[bit_indices[qargs[0]]])
+                        if abc_control is not None:
+                            controlled_circ.control(abc_control)
+
+                        controlled_circ.mcp(local_phase, mcx_control[:-1], mcx_control[-1])
+
             elif gate.name == "z":
                 controlled_circ.h(q_target[bit_indices[qargs[0]]])
                 controlled_circ.mcx(q_control, q_target[bit_indices[qargs[0]]], q_ancillae)
@@ -248,6 +281,65 @@ def control(
     return cgate
 
 
+def define_mcx_rule(q_controls, q_target, q_ancilla=None):
+    """
+    Applies a N-controlled cnot to the target qubit based on the number
+    of qubits in the q_controls
+    """
+    rule = []
+
+    controls_and_target = q_controls + [q_target]
+    if len(q_controls) == 1:
+        rule += [(CXGate(), controls_and_target, [])]
+    elif len(q_controls) == 2:
+        rule += [(CCXGate(), controls_and_target, [])]
+    elif len(q_controls) == 3:
+        rule += [(C3XGate(), controls_and_target, [])]
+    elif len(q_controls) == 4:
+        rule += [(C4XGate(), controls_and_target, [])]
+    elif len(q_controls) >= 5 and len(q_controls) < 7:
+        rule += [(MCXGate(len(q_controls)), controls_and_target, [])]
+    else:
+        rule = MCXRecursive(len(q_controls))._recurse(controls_and_target, q_ancilla)
+    return rule
+
+
+def define_mcx_control_and_ancilla(control_qubits):
+    """
+    Defines which qubits are to be used as control
+    and as an ancilla for the multicontrolled cnot
+    gate
+    """
+    if len(control_qubits) < 7:
+        mcx_controls = list(control_qubits)
+        abc_ancilla = None
+    else:
+        mcx_controls = control_qubits[1:]
+        abc_ancilla = control_qubits[0]
+    return mcx_controls, abc_ancilla
+
+
+def get_abc_matrices(beta, gamma, delta):
+    """
+    Creates A,B and C matrices such that
+    ABC = I
+    """
+    # A
+    a_rz = RZGate(beta).to_matrix()
+    a_ry = RYGate(gamma / 2).to_matrix()
+    a_matrix = a_rz.dot(a_ry)
+
+    # B
+    b_ry = RYGate(-gamma / 2).to_matrix()
+    b_rz = RZGate(-(delta + beta) / 2).to_matrix()
+    b_matrix = b_ry.dot(b_rz)
+
+    # C
+    c_matrix = RZGate((delta - beta) / 2).to_matrix()
+
+    return a_matrix, b_matrix, c_matrix
+
+
 def _gate_to_dag(operation):
     from qiskit.converters.circuit_to_dag import circuit_to_dag