import numpy as np
from itertools import product, combinations
np.set_printoptions(linewidth=200, suppress=True)
# Define basic matrices
I = np.eye(2, 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}
# Tensor product on list of 5 qubits
def kron_n(ops):
out = np.array([[1]], dtype=complex)
for op in ops:
out = np.kron(out, op)
return out
# Build 5-qubit identity
I5 = kron_n([I]*5)
# Single-qubit gate on qubit q (0..4), 5-qubit system
def gate1(U, q, n=5):
ops = [I]*n
ops[q] = U
return kron_n(ops)
# CNOT with control c, target t on n-qubit system
def cnot(c, t, n=5):
# projectors on control
P0 = np.array([[1,0],[0,0]], dtype=complex)
P1 = np.array([[0,0],[0,1]], dtype=complex)
ops0 = [I]*n
ops1 = [I]*n
ops1[t] = X
U = kron_n([P0 if i==c else I for i in range(n)]) + kron_n([P1 if i==c else (X if i==t else I) for i in range(n)])
return U
# Two-qubit Pauli error on given pair (a,b) with labels p1,p2
def twoq_pauli_on_pair(p1, p2, a, b, n=5):
ops = [I]*n
ops[a] = paulis[p1]
ops[b] = paulis[p2]
return kron_n(ops)
# Projector for ancilla measurement |0> on qubit 4
P_a0 = gate1(np.array([[1,0],[0,0]], dtype=complex), 4)
# Stabilizer projectors on data qubits 0..3
# ZZZZ and XXXX on qubits 0..3 (qubits 0..3 correspond to indices 0..3)
Z4 = kron_n([Z,Z,Z,Z, I])
X4 = kron_n([X,X,X,X, I])
P_Zplus = (I5 + Z4)/2
P_Xplus = (I5 + X4)/2
# Target logical |00>_L = (|0000> + |1111>)/sqrt(2) on data and ancilla |0>
# We'll compute fidelity after projection to data-only subspace; but build GHZ+ on data and |0> on ancilla to compute overlap convenience
# Build GHZ+ on data (0..3)
zero = np.array([1,0], dtype=complex)
one = np.array([0,1], dtype=complex)
# |
psi0000 = kron_n([zero, zero, zero, zero, zero])
psi1111 = kron_n([one, one, one, one, zero])
# But ancilla zero for both; We'll define target state on all 5 qubits before ancilla measure: GHZ+ on data, ancilla |0>
GHZ_data = (kron_n([zero, zero, zero, zero, zero]) + kron_n([one, one, one, one, zero]))/np.sqrt(2)
# Build the circuit unitary without errors
U = gate1(H, 1)
U = cnot(1,2) @ U
U = cnot(1,0) @ U
U = cnot(2,3) @ U
U = cnot(3,4) @ U
U = cnot(0,4) @ U
# But according to the problem, the order is H1 then CNOT12, CNOT10, CNOT23, CNOT34, CNOT04 (right to left as in expression)
# We have built in natural left-to-right order for application to initial state |00000>.
# Verify ideal output state is GHZ on data and ancilla 0 deterministically
init = kron_n([zero, zero, zero, zero, zero])
psi_ideal = U @ init
# Verify ancilla amplitude in |0> equals 1 and data GHZ
# Let's project ancilla to 0
psi_a0 = P_a0 @ psi_ideal
p_a0 = np.vdot(psi_a0, psi_a0).real
psi_a0_norm = psi_a0 / np.sqrt(p_a0)
# Project onto code stabilizers +1
psi_code = P_Zplus @ psi_a0_norm
p_zplus = np.vdot(psi_code, psi_code).real
psi_code = psi_code/np.sqrt(p_zplus)
psi_code2 = P_Xplus @ psi_code
p_xplus = np.vdot(psi_code2, psi_code2).real
psi_code2 = psi_code2/np.sqrt(p_xplus)
# Extract the data-only vector by tracing out ancilla (ancilla is |0> after projection)
# We can build the 4-qubit GHZ vector from psi_code2: pick components where ancilla index is 0
# But fidelity is simply overlap with GHZ_data state on 5 qubits since ancilla is |0>
F0 = abs(np.vdot(GHZ_data, psi_code2))**2
print('Ideal acceptance probs:', p_a0, p_zplus, p_xplus, 'F0', F0)
# Now build the locations list for the five CNOTs for insertion of two-qubit errors after each
# Gate sequence list with pairs (control,target)
cnots = [(1,2), (1,0), (2,3), (3,4), (0,4)]
# Two-qubit error set: all tensor products of {I,X,Y,Z} x {I,X,Y,Z} minus II
single_qubit_labels = ['I','X','Y','Z']
err_pairs = [(a,b) for a in single_qubit_labels for b in single_qubit_labels if not (a=='I' and b=='I')]
# Precompute the unitaries for the circuit step by step to allow insertion of errors after each CNOT
Us = []
state = gate1(H,1) # after H
Us.append(state)
state = cnot(1,2) @ state
Us.append(state)
state = cnot(1,0) @ state
Us.append(state)
state = cnot(2,3) @ state
Us.append(state)
state = cnot(3,4) @ state
Us.append(state)
state = cnot(0,4) @ state
Us.append(state)
# Us[k] = unitary after applying up to gate k (with Us[0] after H, Us[5] after all gates)
# But inserting errors after each CNOT implies: after gate indices 1..5 (corresponding to those CNOTs), we insert a Pauli on that pair
# We'll create function to simulate for a set of errors across selected locations
def simulate_with_errors(error_list):
"""
error_list: dict mapping loc index in [1..5] to a tuple (p1,p2) of Pauli labels on that CNOT's (control,target)
Returns: acceptance probability p_acc, fidelity to target after post-selection (if accepted; if p_acc=0, fidelity irrelevant), and flags about acceptance
"""
# Build the total operator as product of gates with error channels inserted after each relevant CNOT
# Alternatively, simulate state progression stepwise for efficiency
psi = init.copy()
# Apply H
psi = gate1(H,1) @ psi
# Apply each CNOT, then error if present
for idx,(c,t) in enumerate(cnots, start=1):
psi = cnot(c,t) @ psi
if idx in error_list:
p1,p2 = error_list[idx]
E = twoq_pauli_on_pair(p1,p2,c,t)
psi = E @ psi
# Measure ancilla |0>
psi_a0 = P_a0 @ psi
p_a0 = float(np.vdot(psi_a0, psi_a0).real)
if p_a0 <= 1e-15:
return 0.0, None, {'a0':0.0,'z+':None,'x+':None}
psi = psi_a0/np.sqrt(p_a0)
# Measure ZZZZ +1
psi_z = P_Zplus @ psi
p_z = float(np.vdot(psi_z, psi_z).real)
if p_z <= 1e-15:
return 0.0, None, {'a0':p_a0,'z+':0.0,'x+':None}
psi = psi_z/np.sqrt(p_z)
# Measure XXXX +1
psi_x = P_Xplus @ psi
p_x = float(np.vdot(psi_x, psi_x).real)
if p_x <= 1e-15:
return 0.0, None, {'a0':p_a0,'z+':p_z,'x+':0.0}
psi = psi_x/np.sqrt(p_x)
# Fidelity with GHZ_data (on data, ancilla in |0>)
F = float(abs(np.vdot(GHZ_data, psi))**2)
p_acc = p_a0*p_z*p_x
return p_acc, F, {'a0':p_a0,'z+':p_z,'x+':p_x}
# Test no errors case
pacc0,F0,probs= simulate_with_errors({})
print('No error: pacc', pacc0, 'F',F0, 'probs', probs)
# Enumerate single error patterns
single_results = []
for loc in range(1,6):
c,t = cnots[loc-1]
for e in err_pairs:
pacc, F, _ = simulate_with_errors({loc:e})
single_results.append((loc, e, pacc, F))
num_acc_single = sum(1 for r in single_results if r[2] > 1e-12)
num_total_single = len(single_results)
print('Single-fault accepted count:', num_acc_single, 'of', num_total_single)
# Summarize single accepted faults: check if F==1 always
fa = [r for r in single_results if r[2]>1e-12]
print('Unique fidelities among accepted single faults:', sorted(set([round(r[3],12) for r in fa])))
# Enumerate double-fault patterns across distinct locations
from collections import defaultdict
pair_results = []
for (l1,l2) in combinations(range(1,6), 2):
c1,t1 = cnots[l1-1]
c2,t2 = cnots[l2-1]
for e1 in err_pairs:
for e2 in err_pairs:
pacc, F, _ = simulate_with_errors({l1:e1, l2:e2})
pair_results.append((l1,l2,e1,e2,pacc,F))
num_acc_pairs = sum(1 for r in pair_results if r[4] > 1e-12)
print('Double-fault accepted count:', num_acc_pairs, 'of', len(pair_results))
# Compute totals to O(p^2)
p = 1e-6 # symbolic not used; we'll compute coefficients explicitly
# Zero-fault contribution
P0 = 1 # to O(p^2)
N0 = P0 * pacc0 * F0
D0 = P0 * pacc0
# Single-fault contribution at O(p): probability per pattern is p/15; We need sum of pacc*F over patterns divided by 15? precisely: sum over loc and error of (p/15) * pacc(F)
S1_sum_N = sum(r[2]*r[3] for r in single_results) # sum over pacc*F
S1_sum_D = sum(r[2] for r in single_results) # sum over pacc
# Multiply by p/15 per pattern
# We'll keep symbolic: coefficients A1 = (1/15)*S1_sum_N, B1=(1/15)*S1_sum_D
A1 = S1_sum_N/15.0
B1 = S1_sum_D/15.0
# Double-fault contribution at O(p^2): each pattern weight (p/15)^2
S2_sum_N = sum(r[4]*r[5] for r in pair_results)
S2_sum_D = sum(r[4] for r in pair_results)
A2 = S2_sum_N/(15.0**2)
B2 = S2_sum_D/(15.0**2)
print('Coefficients: N0,D0=', N0,D0)
print('Sums: S1_N=', S1_sum_N, 'S1_D=', S1_sum_D, 'S2_N=', S2_sum_N, 'S2_D=', S2_sum_D)
print('A1,B1=', A1,B1, 'A2,B2=', A2,B2)
# Expand F_cond = (N0 + A1 p + A2 p^2) / (D0 + B1 p + B2 p^2) + O(p^3)
# Using series expansion: (a0 + a1 p + a2 p^2)/(b0 + b1 p + b2 p^2) = (a0/b0) + ((a1 b0 - a0 b1)/b0^2) p + ... second order term is (a2 b0^2 - a1 b0 b1 - a0 b0 b2 + a0 b1^2)/b0^3
a0 = N0; b0 = D0; a1 = A1; b1 = B1; a2 = A2; b2 = B2
c1 = (a1*b0 - a0*b1)/(b0**2)
c2 = (a2*(b0**2) - a1*b0*b1 - a0*b0*b2 + a0*(b1**2))/(b0**3)
print('Series coeffs: c1=', c1, 'c2=', c2)
# Let's also compute acceptance probability series for reference: D(p) = D0 + B1 p + B2 p^2
print('Acceptance series D(p) ~', D0, '+', B1, 'p +', B2, 'p^2')