import numpy as np
from itertools import product, combinations
# Define single-qubit 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}
labels = ['I','X','Y','Z']
# Tensor product utility for applying single-qubit unitary on k-th qubit (0..4, 0 is left-most)
def apply_1q(U, k, num=5):
ops = [I]*num
ops[k] = U
M = ops[0]
for i in range(1,num):
M = np.kron(M, ops[i])
return M
# Build CNOT matrix on n-qubit system, control c, target t (0..n-1)
def CNOT(n, c, t):
dim = 2**n
U = np.zeros((dim,dim), dtype=complex)
# basis |q0 q1 ... q_{n-1}>
for idx in range(dim):
# compute bit values
bits = [(idx >> (n-1-k)) & 1 for k in range(n)] # q0 is most significant
new_bits = bits.copy()
if bits[c] == 1:
new_bits[t] ^= 1
# map to new index
new_idx = 0
for k in range(n):
new_idx = (new_idx << 1) | new_bits[k]
U[new_idx, idx] = 1.0
return U
# Apply 2-qubit Pauli P_a \otimes P_b on qubits (qa, qb) in n qubits
# qa and qb are distinct
def two_qubit_pauli(n, qa, qb, label_a, label_b):
ops = [I]*n
ops[qa] = paulis[label_a]
ops[qb] = paulis[label_b]
M = ops[0]
for i in range(1,n):
M = np.kron(M, ops[i])
return M
# Projector onto ancilla |0> on qubit index 4
n = 5
anc = 4
# Build projectors for code stabilizers on 4 data qubits (0..3)
# S_X = X X X X; S_Z = Z Z Z Z
Sx = np.kron(np.kron(np.kron(X,X),X),X)
Sz = np.kron(np.kron(np.kron(Z,Z),Z),Z)
P_code = 0.25*(np.eye(16) + Sx) @ (np.eye(16) + Sz)
# Target codeword |00>_AB physical representative: take GHZ+ = (|0000>+|1111>)/sqrt2
v0000 = np.zeros(16, dtype=complex); v0000[0] = 1.0
v1111 = np.zeros(16, dtype=complex); v1111[-1] = 1.0
code00 = (v0000 + v1111)/np.sqrt(2)
# Normalize (should already be normalized)
code00 = code00 / np.linalg.norm(code00)
# Build unitaries for the ideal gates: H1, CNOTs in order
U_H1 = apply_1q(H, 1, num=5)
U_c12 = CNOT(5, 1, 2)
U_c10 = CNOT(5, 1, 0)
U_c23 = CNOT(5, 2, 3)
U_c34 = CNOT(5, 3, 4)
U_c04 = CNOT(5, 0, 4)
# The circuit unitary with error placeholders: We'll apply sequentially
# Initial state |00000>
psi0 = np.zeros(32, dtype=complex)
psi0[0] = 1.0
# Enumerate single- and double-error cases
# List gates and their qubit pairs for errors
# Order: after each CNOT: (12), (10), (23), (34), (04)
gate_pairs = [(1,2),(1,0),(2,3),(3,4),(0,4)]
U_gates = [U_H1, U_c12, U_c10, U_c23, U_c34, U_c04] # include H as first then CNOTs
# Precompute sequence unitaries for speed? We'll just apply directly each time.
# Generate non-identity two-qubit Pauli labels
pauli2 = [(a,b) for a in labels for b in labels]
pauli2_nonid = [(a,b) for (a,b) in pauli2 if not (a=='I' and b=='I')]
# Function to run the circuit with a list of 5 errors (each either None or (label_a,label_b) for that CNOT)
def run_with_errors(err_list):
psi = psi0.copy()
# H1
psi = U_H1 @ psi
# CNOT12 then error 0
psi = U_c12 @ psi
if err_list[0] is not None:
a,b = err_list[0]
E = two_qubit_pauli(5, 1, 2, a, b)
psi = E @ psi
# CNOT10 then error 1
psi = U_c10 @ psi
if err_list[1] is not None:
a,b = err_list[1]
E = two_qubit_pauli(5, 1, 0, a, b)
psi = E @ psi
# CNOT23 then error 2
psi = U_c23 @ psi
if err_list[2] is not None:
a,b = err_list[2]
E = two_qubit_pauli(5, 2, 3, a, b)
psi = E @ psi
# CNOT34 then error 3
psi = U_c34 @ psi
if err_list[3] is not None:
a,b = err_list[3]
E = two_qubit_pauli(5, 3, 4, a, b)
psi = E @ psi
# CNOT04 then error 4
psi = U_c04 @ psi
if err_list[4] is not None:
a,b = err_list[4]
E = two_qubit_pauli(5, 0, 4, a, b)
psi = E @ psi
return psi
# Function to compute acceptance probabilities and fidelity to code00
def acceptance_and_fidelity(psi):
# Project ancilla to |0>
psi_tensor = psi.reshape(2,2,2,2,2)
# Post-select ancilla = 0 component
phi_data = psi_tensor[:,:,:,:,0].reshape(16)
p_anc0 = np.vdot(phi_data, phi_data).real
if p_anc0 == 0:
return 0.0, 0.0
# Normalize after ancilla measurement
phi_data = phi_data / np.sqrt(p_anc0)
# Project onto codespace
chi = P_code @ phi_data
p_code = np.vdot(chi, chi).real
if p_code == 0:
return 0.0, 0.0
chi = chi / np.sqrt(p_code)
# Fidelity to code00
F = abs(np.vdot(code00, chi))**2
return p_anc0 * p_code, p_anc0 * p_code * F
# Baseline (no error)
psi_clean = run_with_errors([None]*5)
p_acc0, p_good0 = acceptance_and_fidelity(psi_clean)
# Sanity checks
print('Clean acceptance prob:', p_acc0, 'Clean fidelity contribution:', p_good0)
# Single-error enumeration
acc1_sum = 0.0
num1_sum = 0.0
for loc in range(5):
for (a,b) in pauli2_nonid:
errs = [None]*5
errs[loc] = (a,b)
psi = run_with_errors(errs)
p_acc, p_good = acceptance_and_fidelity(psi)
acc1_sum += p_acc
num1_sum += p_good
print('Single-error accepted total (sum over all 5*15 cases):', acc1_sum)
print('Single-error good total:', num1_sum)
# Double-error enumeration
acc2_sum = 0.0
num2_sum = 0.0
for loc1, loc2 in combinations(range(5), 2):
for (a1,b1) in pauli2_nonid:
for (a2,b2) in pauli2_nonid:
errs = [None]*5
errs[loc1] = (a1,b1)
errs[loc2] = (a2,b2)
psi = run_with_errors(errs)
p_acc, p_good = acceptance_and_fidelity(psi)
acc2_sum += p_acc
num2_sum += p_good
print('Double-error accepted sum (over 10*225 cases):', acc2_sum)
print('Double-error good sum:', num2_sum)
# Now compute series coefficients up to O(p^2)
# D = (1-p)^5 + (p/15)(1-p)^4 * acc1_sum + (p/15)^2 (1-p)^3 * acc2_sum + O(p^3)
# N = (1-p)^5 + (p/15)(1-p)^4 * num1_sum + (p/15)^2 (1-p)^3 * num2_sum + O(p^3)
# Expand to O(p^2):
# (1-p)^5 = 1 - 5p + 10 p^2
# (p/15)(1-p)^4 = (p/15)(1 - 4p)
# (p/15)^2 (1-p)^3 = (p^2/225)
# Coefficients d1, d2, n1, n2 such that D = 1 - d1 p + d2 p^2, N = 1 - n1 p + n2 p^2
acc1 = acc1_sum
acc2 = acc2_sum
num1 = num1_sum
num2 = num2_sum
# Compute d1, d2, n1, n2
# D = 1 + [ -5 + acc1/15 ] p + [ 10 - 4*acc1/15 + acc2/225 ] p^2
# N = 1 + [ -5 + num1/15 ] p + [ 10 - 4*num1/15 + num2/225 ] p^2
d1 = 5 - acc1/15.0
n1 = 5 - num1/15.0
d2 = 10 - (4*acc1)/15.0 + acc2/225.0
n2 = 10 - (4*num1)/15.0 + num2/225.0
print('d1, d2 =', d1, d2)
print('n1, n2 =', n1, n2)
# Provide rational approximations maybe