"""Optimal CX-cost decomposition of a :class:`.Clifford` operator on 2 qubits

or 3 qubits into a :class:`.QuantumCircuit` based on the Bravyi-Maslov method [1].

Args:

clifford: A Clifford operator.

Returns:

A circuit implementation of the Clifford.

Raises:

QiskitError: if Clifford is on more than 3 qubits.

References:

1. S. Bravyi, D. Maslov, *Hadamard-free circuits expose the

structure of the Clifford group*,

`arXiv:2003.09412 [quant-ph] <https://arxiv.org/abs/2003.09412>`_

"""

num_qubits = clifford.num_qubits

if num_qubits > 3:

raise QiskitError("Can only decompose up to 3-qubit Clifford circuits.")

if num_qubits == 1:

return _decompose_clifford_1q(clifford.tableau)

clifford_name = str(clifford)

# Inverse of final decomposed circuit

inv_circuit = QuantumCircuit(num_qubits, name="inv_circ")

# CNOT cost of clifford

cost = _cx_cost(clifford)

# Find composition of circuits with CX and (H.S)^a gates to reduce CNOT count

while cost > 0:

clifford, inv_circuit, cost = _reduce_cost(clifford, inv_circuit, cost)

# Decompose the remaining product of 1-qubit cliffords

ret_circ = QuantumCircuit(num_qubits, name=clifford_name)

for qubit in range(num_qubits):

pos = [qubit, qubit + num_qubits]

circ = _decompose_clifford_1q(clifford.tableau[pos][:, pos + [-1]])

if len(circ) > 0:

ret_circ.append(circ, [qubit])

# Add the inverse of the 2-qubit reductions circuit

if len(inv_circuit) > 0:

ret_circ.append(inv_circuit.inverse(), range(num_qubits))

"""Decompose a single-qubit clifford"""

circuit = QuantumCircuit(1, name="temp")

# Add phase correction

destab_phase, stab_phase = tableau[:, 2]

if destab_phase and not stab_phase:

circuit.z(0)

elif not destab_phase and stab_phase:

circuit.x(0)

elif destab_phase and stab_phase:

circuit.y(0)

destab_phase_label = "-" if destab_phase else "+"

stab_phase_label = "-" if stab_phase else "+"

destab_x, destab_z = tableau[0, 0], tableau[0, 1]

stab_x, stab_z = tableau[1, 0], tableau[1, 1]

# Z-stabilizer

if stab_z and not stab_x:

stab_label = "Z"

if destab_z:

destab_label = "Y"

circuit.s(0)

else:

destab_label = "X"

# X-stabilizer

elif not stab_z and stab_x:

stab_label = "X"

if destab_x:

destab_label = "Y"

circuit.sdg(0)

else:

destab_label = "Z"

circuit.h(0)

# Y-stabilizer

else:

stab_label = "Y"

if destab_z:

destab_label = "Z"

else:

destab_label = "X"

circuit.s(0)

circuit.h(0)

circuit.s(0)

# Add circuit name

name_destab = f"Destabilizer = ['{destab_phase_label}{destab_label}']"

name_stab = f"Stabilizer = ['{stab_phase_label}{stab_label}']"

circuit.name = f"Clifford: {name_stab}, {name_destab}"

"""Two-qubit cost reduction step"""

num_qubits = clifford.num_qubits

for qubit0 in range(num_qubits):

for qubit1 in range(qubit0 + 1, num_qubits):

for n0, n1 in product(range(3), repeat=2):

# Apply a 2-qubit block

reduced = clifford.copy()

for qubit, n in [(qubit0, n0), (qubit1, n1)]:

if n == 1:

_append_v(reduced, qubit)

elif n == 2:

_append_w(reduced, qubit)

_append_cx(reduced, qubit0, qubit1)

# Compute new cost

new_cost = _cx_cost(reduced)

if new_cost == cost - 1:

# Add decomposition to inverse circuit

for qubit, n in [(qubit0, n0), (qubit1, n1)]:

if n == 1:

inv_circuit.sdg(qubit)

inv_circuit.h(qubit)

elif n == 2:

inv_circuit.h(qubit)

inv_circuit.s(qubit)

inv_circuit.cx(qubit0, qubit1)

return reduced, inv_circuit, new_cost

# If we didn't reduce cost

"""Return the number of CX gates required for Clifford decomposition."""

if clifford.num_qubits == 2:

return _cx_cost2(clifford)

if clifford.num_qubits == 3:

return _cx_cost3(clifford)

"""Return rank of 2x2 boolean matrix."""

if (a & d) ^ (b & c):

return 2

if a or b or c or d:

return 1

"""Return CX cost of a 2-qubit clifford."""

U = clifford.tableau[:, :-1]

r00 = _rank2(U[0, 0], U[0, 2], U[2, 0], U[2, 2])

r01 = _rank2(U[0, 1], U[0, 3], U[2, 1], U[2, 3])

if r00 == 2:

return r01

"""Return CX cost of a 3-qubit clifford."""

# pylint: disable=too-many-return-statements,too-many-boolean-expressions

U = clifford.tableau[:, :-1]

n = 3

# create information transfer matrices R1, R2

R1 = np.zeros((n, n), dtype=int)

R2 = np.zeros((n, n), dtype=int)

for q1 in range(n):

for q2 in range(n):

R2[q1, q2] = _rank2(U[q1, q2], U[q1, q2 + n], U[q1 + n, q2], U[q1 + n, q2 + n])

mask = np.zeros(2 * n, dtype=int)

mask[[q2, q2 + n]] = 1

isLocX = np.array_equal(U[q1, :] & mask, U[q1, :])

isLocZ = np.array_equal(U[q1 + n, :] & mask, U[q1 + n, :])

isLocY = np.array_equal((U[q1, :] ^ U[q1 + n, :]) & mask, (U[q1, :] ^ U[q1 + n, :]))

R1[q1, q2] = 1 * (isLocX or isLocZ or isLocY) + 1 * (isLocX and isLocZ and isLocY)

diag1 = np.sort(np.diag(R1)).tolist()

diag2 = np.sort(np.diag(R2)).tolist()

nz1 = np.count_nonzero(R1)

nz2 = np.count_nonzero(R2)

if diag1 == [2, 2, 2]:

return 0

if diag1 == [1, 1, 2]:

return 1

if (

diag1 == [0, 1, 1]

or (diag1 == [1, 1, 1] and nz2 < 9)

or (diag1 == [0, 0, 2] and diag2 == [1, 1, 2])

):

return 2

if (

(diag1 == [1, 1, 1] and nz2 == 9)

or (

diag1 == [0, 0, 1]

and (nz1 == 1 or diag2 == [2, 2, 2] or (diag2 == [1, 1, 2] and nz2 < 9))

)

or (diag1 == [0, 0, 2] and diag2 == [0, 0, 2])

or (diag2 == [1, 2, 2] and nz1 == 0)

):

return 3

if diag2 == [0, 0, 1] or (

diag1 == [0, 0, 0]

and (

(diag2 == [1, 1, 1] and nz2 == 9 and nz1 == 3)

or (diag2 == [0, 1, 1] and nz2 == 8 and nz1 == 2)

)

):

return 5

if nz1 == 3 and nz2 == 3:

return 6