import numpy as np
from itertools import product, combinations
# Define basic 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)
# Pauli dict
paulis = {'I':I,'X':X,'Y':Y,'Z':Z}
pauli_labels = ['I','X','Y','Z']
# Utilities to build n-qubit ops
def kron_n(ops):
out = np.array([[1]], dtype=complex)
for op in ops:
out = np.kron(out, op)
return out
def embed_single(op, n, q):
# op acts on qubit q (0-indexed, with 0 as leftmost? We'll choose qubit ordering 0..n-1 left-to-right.)
# We'll map qubit index 0 as most significant bit (leftmost in tensor ordering)
ops = []
for i in range(n):
ops.append(op if i==q else I)
return kron_n(ops)
def embed_two(opA, qA, opB, qB, n):
# return operator acting with opA on qA and opB on qB, identity elsewhere
# qubits are 0..n-1 left-to-right
ops = []
for i in range(n):
if i==qA:
ops.append(opA)
elif i==qB:
ops.append(opB)
else:
ops.append(I)
return kron_n(ops)
def cnot(n, control, target):
# Build full matrix of CNOT on n qubits
# Using projector method: |0><0|_c ⊗ I + |1><1|_c ⊗ X_t
P0 = np.array([[1,0],[0,0]], dtype=complex)
P1 = np.array([[0,0],[0,1]], dtype=complex)
term0 = embed_single(P0, n, control)
term1 = embed_single(P1, n, control)
Xt = embed_single(X, n, target)
return term0 + term1 @ Xt
# Stabilizers for data qubits 0..3 (ancilla is qubit 4)
def stabilizer_projector(n):
# n includes ancilla qubit as highest index (4)
# Build XXXX and ZZZZ on data qubits only (0..3)
X4 = embed_two(X,0,X,1,n) @ embed_two(X,2,X,3,n) # X on 0,1,2,3
# But embed_two only handles two qubits; we build product of XX on (0,1) and (2,3)
# Actually above line is correct: apply XX on 0-1 then on 2-3; they commute
Z4 = embed_two(Z,0,Z,1,n) @ embed_two(Z,2,Z,3,n)
Iall = np.eye(2**n, dtype=complex)
P = (Iall + X4)/2 @ (Iall + Z4)/2
return P
# Target logical |00> state on data qubits is (|0000>+|1111>)/sqrt(2)
def logical00_state(n):
# n=5 qubits, with ancilla 4; we return vector in full Hilbert space with ancilla |0>
ket0000 = np.zeros(2**n, dtype=complex)
ket1111 = np.zeros(2**n, dtype=complex)
# index mapping: basis index with qubits 0..4: basis bits b0 b1 b2 b3 b4; state index int(b0...b4,2)
def idx(bits):
val=0
for b in bits:
val = (val<<1) | b
return val
# ancilla fixed 0
ket0000[idx([0,0,0,0,0])] = 1.0
ket1111[idx([1,1,1,1,0])] = 1.0
state = (ket0000 + ket1111)/np.sqrt(2)
return state
# Build the noiseless state prep unitary
n=5
U = np.eye(2**n, dtype=complex)
# Apply operations right-to-left: H1, CNOT12, CNOT10, CNOT23, CNOT34, CNOT04
U = cnot(n,0,4) @ U
U = cnot(n,3,4) @ U
U = cnot(n,2,3) @ U
U = cnot(n,1,0) @ U
U = cnot(n,1,2) @ U
U = embed_single(H, n, 1) @ U
# Verify that U |00000> yields |overline{00} on data, ancilla 0
psi0 = np.zeros(2**n, dtype=complex)
psi0[0] = 1.0
psi = U @ psi0
# compute overlap with logical00_state
phi = logical00_state(n)
ov = np.vdot(phi, psi)
print('Noiseless overlap amplitude with |L00, a0>:', ov)
print('Squared magnitude:', abs(ov)**2)
# Build projector onto ancilla |0>
P_a0 = embed_single(np.array([[1,0],[0,0]], dtype=complex), n, 4)
P_code = stabilizer_projector(n)
# Confirm acceptance prob noiseless
acc = np.vdot(psi, P_a0 @ P_code @ psi).real
print('Noiseless acceptance probability:', acc)
# Generate list of 15 non-identity two-qubit Paulis on a given pair as label tuples
nonId_two_paulis = []
for a in pauli_labels:
for b in pauli_labels:
if a=='I' and b=='I':
continue
nonId_two_paulis.append((a,b))
# Define function to embed two-qubit Pauli on specific qubit pair
def two_qubit_pauli_op(n, q1, q2, label_pair):
A = paulis[label_pair[0]]
B = paulis[label_pair[1]]
return embed_two(A, q1, B, q2, n)
# Function to simulate with specific error tuple for up to two faults
def simulate_with_errors(errors):
# errors is dict mapping gate index to label pair, gate indices 1..5 for the five CNOTs in order [12,10,23,34,04]
# Build the unitary with errors inserted after each CNOT
ops = []
# We'll apply sequentially starting from |000..0>
state = psi0.copy()
# define sequence: H1, CNOT12, E1, CNOT10, E2, CNOT23, E3, CNOT34, E4, CNOT04, E5
# H1
state = embed_single(H, n, 1) @ state
# CNOT12
state = cnot(n,1,2) @ state
if 1 in errors:
state = two_qubit_pauli_op(n,1,2,errors[1]) @ state
# CNOT10
state = cnot(n,1,0) @ state
if 2 in errors:
state = two_qubit_pauli_op(n,1,0,errors[2]) @ state
# CNOT23
state = cnot(n,2,3) @ state
if 3 in errors:
state = two_qubit_pauli_op(n,2,3,errors[3]) @ state
# CNOT34
state = cnot(n,3,4) @ state
if 4 in errors:
state = two_qubit_pauli_op(n,3,4,errors[4]) @ state
# CNOT04
state = cnot(n,0,4) @ state
if 5 in errors:
state = two_qubit_pauli_op(n,0,4,errors[5]) @ state
# Compute acceptance probability and correct-target probability
# acceptance projector P = P_code ⊗ |0><0|_a
P_acc = P_a0 @ P_code
psi_acc = P_acc @ state
p_acc = np.vdot(state, psi_acc).real # ||P psi||^2, since P is projector but order matters? Here it's <psi|P|psi>
# Correct-target projector
phi = logical00_state(n) # includes ancilla 0
P_target = np.outer(phi, np.conj(phi))
p_target = np.vdot(state, P_target @ state).real
# But p_target includes ancilla=0 implicitly and no stabilizer measurement; Actually phi lies in code subspace, so P_target <= P_acc
return p_acc, p_target
# Enumerate 0,1,2-fault cases and compute coefficients
from collections import defaultdict
# Gate mapping indices and qubit pairs
pairs = {1:(1,2), 2:(1,0), 3:(2,3), 4:(3,4), 5:(0,4)}
# 0 faults
pacc0, ptarget0 = simulate_with_errors({})
print('0 faults: P_acc=', pacc0, 'P_target=', ptarget0)
# 1 fault cases
acc1 = []
cor1 = []
for gi in range(1,6):
for lab in nonId_two_paulis:
pacc, ptgt = simulate_with_errors({gi: lab})
acc1.append(pacc)
cor1.append(ptgt)
acc1_sum = float(np.sum(acc1))
cor1_sum = float(np.sum(cor1))
print('Sum over single fault P_acc contributions (unweighted):', acc1_sum)
print('Sum over single fault P_target contributions (unweighted):', cor1_sum)
# 2 fault cases
acc2 = []
cor2 = []
for gi, gj in combinations(range(1,6), 2):
for lab1 in nonId_two_paulis:
for lab2 in nonId_two_paulis:
pacc, ptgt = simulate_with_errors({gi: lab1, gj: lab2})
acc2.append(pacc)
cor2.append(ptgt)
acc2_sum = float(np.sum(acc2))
cor2_sum = float(np.sum(cor2))
print('Sum over two-fault P_acc contributions (unweighted):', acc2_sum)
print('Sum over two-fault P_target contributions (unweighted):', cor2_sum)
# Now assemble series up to O(p^2)
# Using weights: k=0: (1 - 5p + 10 p^2), k=1: (p/15 - 4 p^2/15), k=2: (p^2/225)
import sympy as sp
p = sp.symbols('p', real=True)
Pacc_series = (pacc0*(1 - 5*p + 10*p**2)
+ (acc1_sum)*(p/15 - 4*p**2/15)
+ (acc2_sum)*(p**2/225))
Pcor_series = (ptarget0*(1 - 5*p + 10*p**2)
+ (cor1_sum)*(p/15 - 4*p**2/15)
+ (cor2_sum)*(p**2/225))
Pacc_series = sp.simplify(sp.expand(Pacc_series))
Pcor_series = sp.simplify(sp.expand(Pcor_series))
print('Pacc_series:', Pacc_series)
print('Pcor_series:', Pcor_series)
# Fidelity F = Pcor_series / Pacc_series up to O(p^2)
F_series = sp.series(Pcor_series / Pacc_series, p, 0, 3) # series up to p^2
print('Fidelity series up to O(p^2):', F_series)
# Also give exact rational coefficients numerically
F2_coeff = sp.expand(F_series.removeO()).coeff(p,2)
print('Quadratic coefficient (p^2 term) in infidelity 1-F:', sp.nsimplify(1 - (1 + 0*p + F2_coeff*p**2).subs(p,0)))
print('F_series no O():', sp.expand(F_series.removeO()))