import numpy as np
from itertools import product, combinations
# Pauli matrices
I = np.array([[1,0],[0,1]], dtype=complex)
X = np.array([[0,1],[1,0]], dtype=complex)
Y = np.array([[0,-1j],[1j,0]], dtype=complex)
Z = np.array([[1,0],[0,-1]], dtype=complex)
H = (1/np.sqrt(2))*np.array([[1,1],[1,-1]], dtype=complex)
paulis = {'I':I, 'X':X, 'Y':Y, 'Z':Z}
pauli_labels = ['I','X','Y','Z']
# Tensor product utility for 5 qubits
def kron_n(ops):
M = np.array([[1]], dtype=complex)
for op in ops:
M = np.kron(M, op)
return M
# 5-qubit identity
I5 = kron_n([I]*5)
# Build single-qubit gate on specific qubit index 0..4
def single_qubit_gate(U, q):
ops = [I]*5
ops[q] = U
return kron_n(ops)
# CNOT gate between control c and target t over 5 qubits
def CNOT(c, t):
# |0><0|_c ⊗ I + |1><1|_c ⊗ X_t, identity on others
P0 = np.array([[1,0],[0,0]], dtype=complex)
P1 = np.array([[0,0],[0,1]], dtype=complex)
ops0 = [I]*5
ops1 = [I]*5
ops0[c] = P0
ops1[c] = P1
# target X
Xt = [I]*5
Xt[t] = X
return kron_n(ops0) + kron_n(ops1) @ kron_n(Xt)
# two-qubit Pauli on qubits (a,b) with labels pa,b
def two_qubit_pauli(pa, pb, a, b):
ops = [I]*5
ops[a] = paulis[pa]
ops[b] = paulis[pb]
return kron_n(ops)
# Stabilizers and logicals for [[4,2,2]] on qubits 0..3
Sx = kron_n([X,X,X,X,I]) # XXXX on 0..3
Sz = kron_n([Z,Z,Z,Z,I]) # ZZZZ on 0..3
ZA = kron_n([Z,Z,I,I,I]) # Z_A = ZZII
ZB = kron_n([Z,I,Z,I,I]) # Z_B = ZIZI
XA = kron_n([X,I,X,I,I]) # X_A = XIXI
XB = kron_n([X,X,I,I,I]) # X_B = XXII
# Ideal |00000> state
zero = np.array([1,0], dtype=complex)
one = np.array([0,1], dtype=complex)
ket0_5 = kron_n([zero, zero, zero, zero, zero])
# Build the ideal unitary circuit U (without errors)
U = I5
# Right-to-left execution in circuit string: H1 then CNOT12, CNOT10, CNOT23, CNOT34, CNOT04
# We'll apply in that order
U = CNOT(0,4) @ CNOT(3,4) @ CNOT(2,3) @ CNOT(1,0) @ CNOT(1,2) @ single_qubit_gate(H,1)
# Prepare ideal state
psi_ideal = U @ ket0_5
# Ancilla Z-basis projectors
P_anc0 = single_qubit_gate(np.array([[1,0],[0,0]], dtype=complex), 4)
P_anc1 = single_qubit_gate(np.array([[0,0],[0,1]], dtype=complex), 4)
# Define projector to codespace (+1 eigenspace of XXXX and ZZZZ) on data qubits 0..3
P_code = (np.eye(32, dtype=complex) + Sx)/2
P_code = (P_code + Sz @ P_code)/2 # equivalent to (I+Sx)/2 (I+Sz)/2 but ensure Hermitian
# Extract data-only state vector after projecting ancilla to |0>
def postselect_anc0(state):
psi = P_anc0 @ state
p0 = float((psi.conj().T @ psi).real)
if p0 == 0:
return None, 0.0
psi = psi / np.sqrt(p0)
# trace out ancilla by reshaping
# Alternatively, project ancilla space to |0> and then pick data subspace components where ancilla=0 basis;
# We'll reshape as (2^4,2) and select column 0
psi_full = psi.reshape((16,2), order='C')
data = psi_full[:,0]
return data, p0
# Normalize data state
def normalize(v):
n = np.linalg.norm(v)
if n == 0:
return v, 0.0
return v/n, n**2
# Build ideal logical |00>_L on data qubits (should equal GHZ4)
# We'll explicitly construct GHZ4 as (|0000>+|1111>)/sqrt(2)
ket0000 = kron_n([zero, zero, zero, zero, zero])[:16] # first 4 qubits; ignore ancilla index
ket1111 = kron_n([one, one, one, one, zero])[:16]
phi_ghz = (ket0000 + ket1111)/np.sqrt(2)
# Verify it's stabilized
sx_val = (phi_ghz.conj().T @ Sx[:16,:16] @ phi_ghz).round(10)
sz_val = (phi_ghz.conj().T @ Sz[:16,:16] @ phi_ghz).round(10)
# Verify GHZ corresponds to +1 eigenstate of ZA and ZB
za_val = (phi_ghz.conj().T @ ZA[:16,:16] @ phi_ghz).round(10)
zb_val = (phi_ghz.conj().T @ ZB[:16,:16] @ phi_ghz).round(10)
print('Sx=', sx_val, 'Sz=', sz_val, 'ZA=', za_val, 'ZB=', zb_val)
# Prepare list of CNOTs in chronological application order, and track which qubits they act on
cnot_list = [(1,2), (1,0), (2,3), (3,4), (0,4),] # But we're missing H_1; we only inject errors after CNOTs; we also had 5 CNOTs? Wait there are 5 listed. Actually the circuit includes 5 CNOTs after H and then 2 ancilla CNOTs; No: the circuit prints CNOT_{23}, CNOT_{10}, CNOT_{12}, H_1 in right-to-left; plus CNOT_{34}, CNOT_{04}
# Let's ensure we have 5 CNOTs? Provided list includes CNOT_{23}, CNOT_{10}, CNOT_{12}, CNOT_{34}, CNOT_{04} -> that's 5.
print('Number of CNOTs we included:', len(cnot_list))
# But the problem text shows six operations: CNOT04, CNOT34, CNOT23, CNOT10, CNOT12: indeed 5 CNOTS, plus H1 and measurement. OK.
# We need the interleaved unitaries: After each CNOT, we can insert error.
# We'll precompute the prefix unitaries and suffix unitaries to efficiently apply errors.
# Actually U = CNOT(0,4) @ CNOT(3,4) @ CNOT(2,3) @ CNOT(1,0) @ CNOT(1,2) @ H1.
# We'll create the sequence
ops = [single_qubit_gate(H,1), CNOT(1,2), CNOT(1,0), CNOT(2,3), CNOT(3,4), CNOT(0,4)]
# Confirm that multiplying them equals U
UU = ops[5] @ ops[4] @ ops[3] @ ops[2] @ ops[1] @ ops[0]
print('U matches:', np.allclose(U, UU))
# Precompute prefixes after each step
prefix = [np.eye(32, dtype=complex)]
for k in range(len(ops)):
prefix.append(ops[k] @ prefix[-1])
# prefix[i] = ops[i-1] @ ... @ ops[0]
# For injecting error after k-th CNOT, we take state: prefix[k+1] (since ops[0] is H, ops[1] is first CNOT)
# We'll implement enumeration over pairs of distinct CNOT index positions among indices [1..5] corresponding to ops indices for CNOTs
cnot_indices = [1,2,3,4,5] # positions in ops of CNOTs
cnot_qubits = {1:(1,2),2:(1,0),3:(2,3),4:(3,4),5:(0,4)}
# Build suffix unitaries from the following operations after the injected error
suffix = [None]*(len(ops)+1)
# suffix[k] = ops[n-1] ... ops[k]
acc = np.eye(32,dtype=complex)
for k in range(len(ops), -1, -1):
suffix[k] = acc
if k>0:
acc = ops[k-1] @ acc
# Define 2-qubit Pauli label pairs excluding II
pauli_pairs = [(a,b) for a in pauli_labels for b in pauli_labels if not (a=='I' and b=='I')]
# Function to build error operator E_k acting after CNOT at index k
def error_after_k(k, pa, pb):
# Build Pauli on the two involved qubits of that CNOT
a,b = cnot_qubits[k]
return two_qubit_pauli(pa, pb, a, b)
# Given two errors at k1,k2 with (pa,pb) each, compute accepted probability contribution and whether final accepted state equals |00>_L
# We'll compute the net accepted state after projecting ancilla to 0 and projecting to codespace (+1,+1) and normalizing
phighz = phi_ghz # ideal |00>_L
# We'll accumulate counts for events (accepted? and fidelity correctness?)
accepted_pairs = 0
accepted_wrong_pairs = 0
accepted_correct_pairs = 0
# We also want to verify whether acceptance weights p0 and p_code are 0/1
partial_acceptance_cases = []
# enumerate unordered pairs (k1<k2)
for k1, k2 in combinations(cnot_indices, 2):
for pa1, pb1 in pauli_pairs:
E1 = error_after_k(k1, pa1, pb1)
# state after first error only: not required
for pa2, pb2 in pauli_pairs:
E2 = error_after_k(k2, pa2, pb2)
# Full state: apply the whole sequence with two errors inserted
# |ψ> = suffix[k2+1] @ E2 @ suffix[k1+1..k2] @ E1 @ prefix[k1+1] @ |0>
# But we can simply do: U with errors: ops[5] ... ops[k2+1] @ E2 @ ops[k2] ... ops[k1+1] @ E1 @ ops[k1] ... ops[0] |0>
# We have prefix and suffix: final = suffix[k2+1] @ E2 @ (ops[k2]...ops[k1+1]) @ E1 @ prefix[k1+1] |0>
# But we don't have mid-chunk handy; easier: directly multiply all ops sequentially with errors when needed
psi = ket0_5.copy()
for idx in range(len(ops)):
psi = ops[idx] @ psi
if idx+1 == k1:
psi = E1 @ psi
if idx+1 == k2:
psi = E2 @ psi
# Now we have final |ψ>
# Postselect ancilla=0
data, p0 = postselect_anc0(psi)
if p0 == 0:
continue # rejected by ancilla
# Normalize data
data = data / np.linalg.norm(data)
# Now project to codespace (+1,+1): apply P_code and renormalize
# Compute acceptance probability for code project
p_code = float((data.conj().T @ (P_code[:16,:16] @ data)).real)
if p_code == 0:
# Rejected by stabilizers
continue
# accepted
accepted_pairs += 1
# Form post-selected data state
psi_data = P_code[:16,:16] @ data
psi_data = psi_data / np.linalg.norm(psi_data)
# Fidelity with ideal codeword
F = float(np.abs(phighz.conj().T @ psi_data)**2)
if F > 1 - 1e-9:
accepted_correct_pairs += 1
elif F < 1e-9:
accepted_wrong_pairs += 1
else:
partial_acceptance_cases.append((k1,k2,pa1+pb1,pa2+pb2,F))
print('Accepted pair cases count:', accepted_pairs)
print('Accepted correct:', accepted_correct_pairs)
print('Accepted wrong:', accepted_wrong_pairs)
print('Partial F cases:', len(partial_acceptance_cases))
# Also compute total number of unordered pairs and error pair combos
total_unordered_pairs = len(list(combinations(cnot_indices,2)))
num_error_pairs_per_position = len(pauli_pairs)**2
print('Total unordered location pairs:', total_unordered_pairs)
print('Pauli-pair combos per location pair:', num_error_pairs_per_position)
print('Total combinations enumerated:', total_unordered_pairs * num_error_pairs_per_position)
# Compute coefficient C where P_bad_acc_2 = (accepted_wrong_pairs / total_combos) * p^2
# But be careful: acceptance is not necessarily all accepted; But we counted only accepted events; The probability weight is count of accepted wrong events divided by total combos (including unaccepted), because each event has probability (p/15)^2 * (1-p)^(others). To O(p^2), it's (p/15)^2 times the fraction of combos that are accepted and wrong.
all_combos = total_unordered_pairs * num_error_pairs_per_position
C = accepted_wrong_pairs / all_combos
print('Coefficient fraction C (pre-mult p^2):', C)
# For completeness, compute acceptance fraction among combos at order p^2
acceptance_fraction_order2 = accepted_pairs / all_combos
correct_fraction_order2 = accepted_correct_pairs / all_combos
print('Acceptance fraction among two-fault combos:', acceptance_fraction_order2)
print('Correct acceptance fraction among two-fault combos:', correct_fraction_order2)
# Let's also verify that single-fault events produce zero accepted wrong events
# Enumerate single-fault case
accepted_single = 0
wrong_single = 0
for k in cnot_indices:
for pa,pb in pauli_pairs:
psi = ket0_5.copy()
for idx in range(len(ops)):
psi = ops[idx] @ psi
if idx+1 == k:
psi = error_after_k(k, pa, pb) @ psi
data, p0 = postselect_anc0(psi)
if p0 == 0:
continue
data = data/np.linalg.norm(data)
p_code = float((data.conj().T @ (P_code[:16,:16] @ data)).real)
if p_code == 0:
continue
accepted_single += 1
psi_data = P_code[:16,:16] @ data
psi_data = psi_data / np.linalg.norm(psi_data)
F = float(np.abs(phighz.conj().T @ psi_data)**2)
if F < 1-1e-9:
wrong_single += 1
print('Single-fault accepted cases:', accepted_single, 'bad among them:', wrong_single)