import numpy as np
from itertools import product, combinations
# Qubit indexing: 0..4, with 4 the ancilla
nq = 5
N = 2**nq
def kron_n(ops):
out = ops[0]
for m in ops[1:]:
out = np.kron(out, m)
return out
# Single-qubit Paulis
I = np.array([[1,0],[0,1]], complex)
X = np.array([[0,1],[1,0]], complex)
Y = np.array([[0,-1j],[1j,0]], complex)
Z = np.array([[1,0],[0,-1]], complex)
H = (1/np.sqrt(2))*np.array([[1,1],[1,-1]], complex)
paulis = {'I':I,'X':X,'Y':Y,'Z':Z}
# Build n-qubit op acting on specified qubits for 1q or 2q gate
def one_qubit_op(U, q, nq=5):
ops = []
for i in range(nq):
ops.append(U if i==q else I)
return kron_n(ops)
def two_qubit_op(U2, q1, q2, nq=5):
# U2 acts on (q1,q2) in that order in the tensor product
ops = []
for i in range(nq):
if i==q1:
ops.append(None) # placeholder
elif i==q2:
ops.append(None)
else:
ops.append(I)
# Build kron in order 0..nq-1
out = None
for i in range(nq):
if ops[i] is None:
# when hitting first of q1 or q2, insert U2 appropriately
if i==q1:
block = U2
elif i==q2:
# if q2 comes before q1 in ordering; then need to swap order of kron product inputs
# Our U2 is defined in order (q1,q2), but if q2<q1, we need to account by conjugating with SWAP
# Simpler: build embedding that respects placement by using identities and then permute?
# To avoid complexity, reorder so that q1<q2 always when building U2 embedding
raise RuntimeError("two_qubit_op called with q2 before q1; provide q1<q2")
# now append block and skip q2 when we reach it
if out is None:
out = block
else:
out = np.kron(out, block)
# skip the partner qubit in the iteration by marking it as used
# But we can't skip in for loop; We'll mark ops[q2] to identity after adding block
ops[q2] = 'DONE'
elif ops[i] == 'DONE':
continue
else:
out = ops[i] if out is None else np.kron(out, ops[i])
return out
# Let's avoid the above complexity: instead construct general n-qubit op via tensor product building, not using placeholders.
# Simpler approach: Use index mapping to place U2 at positions q1 and q2 by building kron product of 5 factors, where the factor is a 2x2 for other qubits,
# and for q1 use basis e00 T e00 etc is complicated. Alternative: Build U as 32x32 using computational basis mapping by applying U2 to bits of indices.
# Let's implement a general function apply_two_qubit_gate(state, U2, q1, q2) acting on state vector without forming full matrix.
def apply_two_qubit_gate_statevec(state, U2, q1, q2, nq=5):
# state is vector of length 2**nq
# U2 acts on qubits q1 (most significant within the pair) and q2 (least significant within the pair) as given
# We'll loop over basis states and update amplitudes in blocks of size 4 corresponding to q1,q2 bits.
# For efficiency, precompute bit masks
if q1==q2:
raise ValueError
# Ensure q1<q2 for block enumeration ease
if q1>q2:
q1, q2 = q2, q1
# but then U2 should be swapped: i.e., we need to apply U2' that corresponds to SWAP-conjugated U2
# For Clifford gates and Paulis it'll be okay if we pass the swapped version.
# Let's build swap matrix S and compute U2_swapped = S U2 S
S = np.array([[1,0,0,0],
[0,0,1,0],
[0,1,0,0],
[0,0,0,1]], complex)
U2 = S @ U2 @ S
out = np.zeros_like(state)
bit1 = 1 << (nq-1-q1)
bit2 = 1 << (nq-1-q2)
# Iterate over all indices, group them by fixing other bits and the two bits of interest vary over 00,01,10,11
N = state.shape[0]
for base in range(0, N, 1):
# We'll only process when the bits at q1 and q2 are 00; then build the 4-vector with other bits same
if (base & bit1) or (base & bit2):
continue
# collect indices for 00,01,10,11
idx00 = base
idx01 = base | bit2
idx10 = base | bit1
idx11 = base | bit1 | bit2
v = np.array([state[idx00], state[idx01], state[idx10], state[idx11]], complex)
v2 = U2 @ v
out[idx00] += v2[0]
out[idx01] += v2[1]
out[idx10] += v2[2]
out[idx11] += v2[3]
return out
# Similarly for one-qubit gate
def apply_one_qubit_gate_statevec(state, U, q, nq=5):
out = np.zeros_like(state)
bit = 1 << (nq-1-q)
N = state.shape[0]
for base in range(0, N, 1):
if base & bit:
continue
idx0 = base
idx1 = base | bit
v = np.array([state[idx0], state[idx1]], complex)
v2 = U @ v
out[idx0] += v2[0]
out[idx1] += v2[1]
return out
# Build CNOT 2-qubit unitary (4x4) for control->target on those two local qubits
CNOT_2 = np.array([[1,0,0,0],
[0,1,0,0],
[0,0,0,1],
[0,0,1,0]], complex) # in basis |00>,|01>,|10>,|11> with first = control, second = target
# 2-qubit Pauli errors (on the pair) as 4x4 matrices in order (q1,q2)
pauli_labels = ['I','X','Y','Z']
P2 = {}
for a in pauli_labels:
for b in pauli_labels:
lab = a+b
P2[lab] = np.kron(paulis[a], paulis[b])
# GHZ+ on 4 qubits (data), with ordering qubits 0,1,2,3 as in our labeling
psi_ghz4 = np.zeros(16, complex)
psi_ghz4[0] = 1/np.sqrt(2) # |0000>
psi_ghz4[-1] = 1/np.sqrt(2) # |1111>
# Stabilizers on 4 qubits
Z4 = np.kron(np.kron(np.kron(Z,Z),Z),Z)
X4 = np.kron(np.kron(np.kron(X,X),X),X)
# Helper: extract data substate conditioned on ancilla = b
def measure_ancilla_and_postselect(psi5, b, nq=5):
# Return conditional state on data (4-qubit) if measurement outcome b occurs with prob 1, else return None
# Compute probability of outcome b
# ancilla is qubit 4 (index 4), as LSB? In our indexing, q index increases left in bit significance; We used bit mask bit = 1 << (nq-1-q)
# So ancilla bit is bit = 1 << (nq-1-4) = 1 << 0 = 1
bit = 1
N = psi5.shape[0]
tot0 = 0.0
tot1 = 0.0
# Gather data subvector for ancilla=b
data = np.zeros(16, complex)
for idx in range(N):
amp = psi5[idx]
if idx & bit:
tot1 += (abs(amp)**2)
else:
tot0 += (abs(amp)**2)
# Determine if deterministic
if b==0:
if abs(tot0-1.0) < 1e-12:
# build data state for ancilla=0 branch
# mapping: index 4 is LSB; So data index corresponds to (idx >> 1)
for idx in range(N):
if (idx & bit)==0:
data_idx = idx >> 1
data[data_idx] = psi5[idx]
# normalize
data = data / np.sqrt(tot0)
return data
else:
return None
else:
if abs(tot1-1.0) < 1e-12:
for idx in range(N):
if (idx & bit):
data_idx = (idx ^ bit) >> 1 # remove LSB
data[data_idx] = psi5[idx]
data = data / np.sqrt(tot1)
return data
else:
return None
# Apply stabilizer acceptance check: require +1 for both ZZZZ and XXXX
def stabilizer_accepts(psi4):
vZ = (psi4.conj().T @ (Z4 @ psi4)).real
vX = (psi4.conj().T @ (X4 @ psi4)).real
# round near integers
return (abs(vZ-1.0) < 1e-9) and (abs(vX-1.0) < 1e-9)
# Fidelity with GHZ+
def fid_ghz(psi4):
amp = np.vdot(psi_ghz4, psi4)
return (abs(amp)**2).real
# Build function to run the circuit with two specified errors at gates i<j and labels lab_i, lab_j
# Gates: 1:H1; 2:C12; 3:C10; 4:C23; 5:C34; 6:C04
# We'll label CNOT gates indices 1..5 for errors
# Precompute nothing; We'll apply stepwise
def run_with_two_errors(lab_i, gate_i, lab_j, gate_j):
# gate indices start 1: C12, 2: C10, 3: C23, 4: C34, 5: C04
# lab is two-letter string for two-qubit Pauli on that gate's qubits
# Build initial state
psi = np.zeros(N, complex)
psi[0] = 1.0
# H on qubit 1
psi = apply_one_qubit_gate_statevec(psi, H, 1, nq)
# gate 1: CNOT 1->2; then error if gate_i==1 or gate_j==1
psi = apply_two_qubit_gate_statevec(psi, CNOT_2, 1, 2, nq)
if gate_i==1:
E = P2[lab_i]
psi = apply_two_qubit_gate_statevec(psi, E, 1, 2, nq)
if gate_j==1:
E = P2[lab_j]
psi = apply_two_qubit_gate_statevec(psi, E, 1, 2, nq)
# gate 2: CNOT 1->0
psi = apply_two_qubit_gate_statevec(psi, CNOT_2, 1, 0, nq)
if gate_i==2:
E = P2[lab_i]
psi = apply_two_qubit_gate_statevec(psi, E, 1, 0, nq)
if gate_j==2:
E = P2[lab_j]
psi = apply_two_qubit_gate_statevec(psi, E, 1, 0, nq)
# gate 3: CNOT 2->3
psi = apply_two_qubit_gate_statevec(psi, CNOT_2, 2, 3, nq)
if gate_i==3:
E = P2[lab_i]
psi = apply_two_qubit_gate_statevec(psi, E, 2, 3, nq)
if gate_j==3:
E = P2[lab_j]
psi = apply_two_qubit_gate_statevec(psi, E, 2, 3, nq)
# gate 4: CNOT 3->4
psi = apply_two_qubit_gate_statevec(psi, CNOT_2, 3, 4, nq)
if gate_i==4:
E = P2[lab_i]
psi = apply_two_qubit_gate_statevec(psi, E, 3, 4, nq)
if gate_j==4:
E = P2[lab_j]
psi = apply_two_qubit_gate_statevec(psi, E, 3, 4, nq)
# gate 5: CNOT 0->4
psi = apply_two_qubit_gate_statevec(psi, CNOT_2, 0, 4, nq)
if gate_i==5:
E = P2[lab_i]
psi = apply_two_qubit_gate_statevec(psi, E, 0, 4, nq)
if gate_j==5:
E = P2[lab_j]
psi = apply_two_qubit_gate_statevec(psi, E, 0, 4, nq)
# ancilla measurement postselect on 0 deterministically
psi4 = measure_ancilla_and_postselect(psi, 0, nq)
if psi4 is None:
return False, False # not accepted by ancilla (or not deterministic 0)
# stabilizer postselection
if not stabilizer_accepts(psi4):
return False, False
# accepted; now fidelity
F = fid_ghz(psi4)
good = abs(F-1.0) < 1e-9
return True, good
# Now enumerate all pairs of gates and 2q Pauli labels (excluding II)
labels = [a+b for a in pauli_labels for b in pauli_labels if not (a=='I' and b=='I')]
accepted_count = 0
accepted_good = 0
accepted_bad = 0
# We'll accumulate by gate pair
malignant_by_pair = {}
acc_by_pair = {}
for (i,j) in combinations(range(1,6), 2):
mal = 0
acc = 0
for lab_i in labels:
for lab_j in labels:
accepted, good = run_with_two_errors(lab_i, i, lab_j, j)
if accepted:
acc += 1
if not good:
mal += 1
malignant_by_pair[(i,j)] = mal
acc_by_pair[(i,j)] = acc
accepted_count += acc
accepted_bad += mal
accepted_good = accepted_count - accepted_bad
accepted_count, accepted_bad, malignant_by_pair