import numpy as np
from itertools import product, combinations
# Define 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)
paulis = {'I': I, 'X': X, 'Y': Y, 'Z': Z}
# Single-qubit gates
H = (1/np.sqrt(2))*np.array([[1,1],[1,-1]], dtype=complex)
# Tensor utility: build n-qubit operator from list of single-qubit ops (leftmost is qubit 0)
def kronN(ops):
out = np.array([[1]], dtype=complex)
for op in ops:
out = np.kron(out, op)
return out
# n-qubit identity
def In(n):
out = np.array([[1]], dtype=complex)
for _ in range(n):
out = np.kron(out, I)
return out
nq = 5 # qubits 0..4 (4 data + ancilla)
# Build CNOT for arbitrary control c, target t on nq qubits
def CNOT(n, c, t):
# 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)
ops0 = [I]*n
ops1 = [I]*n
ops0[c] = P0
ops1[c] = P1
term0 = kronN(ops0)
ops1[t] = X
term1 = kronN(ops1)
return term0 + term1
# Projectors for ancilla measurement on qubit 4
Pi_a0 = kronN([I,I,I,I, np.array([[1,0],[0,0]], dtype=complex)])
Pi_a1 = kronN([I,I,I,I, np.array([[0,0],[0,1]], dtype=complex)])
# Code states on data qubits 0..3
# |00>_L = (|0000> + |1111>)/sqrt2
basis_4 = np.eye(16, dtype=complex)
ket0000 = np.zeros((16,1), dtype=complex); ket0000[0,0] = 1
ket1111 = np.zeros((16,1), dtype=complex); ket1111[-1,0] = 1
ketGHZp = (ket0000 + ket1111)/np.sqrt(2)
# Logical X operators per problem
# X_A = X I X I on qubits 0..3
# X_B = X X I I on qubits 0..3
XA = kronN([X,I,X,I])
XB = kronN([X,X,I,I])
# Generate code basis |00>, |10>=XA|00>, |01>=XB|00>, |11>=XA XB |00>
ket00 = ketGHZp
ket10 = XA @ ket00
ket01 = XB @ ket00
ket11 = XA @ ket01
# Normalize (should be normalized already)
for k in [ket00, ket10, ket01, ket11]:
k /= np.linalg.norm(k)
# Projector onto code space P_code = sum |ij><ij|
P_code = np.zeros((16,16), dtype=complex)
for k in [ket00, ket10, ket01, ket11]:
P_code += k @ k.conj().T
# Projector onto |00> logical
P_00 = ket00 @ ket00.conj().T
# Expand to 5-qubit operators by tensoring with ancilla identity where necessary
# We'll need to apply gates and errors on full 5-qubit space; For code projectors we act only on 0..3
# So to compute probabilities after ancilla projection, we will extract the data state by partial projection.
# Build the noiseless unitary for the circuit (up to just before measuring ancilla)
U = np.eye(2**nq, dtype=complex)
# Apply H on qubit 1
ops = [I]*nq; ops[1] = H
U = kronN(ops) @ U
# Apply CNOTs in order: 1->2, 1->0, 2->3, 3->4, 0->4
cnots = [(1,2), (1,0), (2,3), (3,4), (0,4)]
CNOTs = [CNOT(nq, c, t) for (c,t) in cnots]
# We'll build functions to apply errors after each CNOT
# List all two-qubit Pauli pairs (including II); return as 5-qubit operator placed on the two qubits of given pair
labels = ['I','X','Y','Z']
pauli_pairs = [] # tuples (label_c, label_t, matrix 4x4)
for a in labels:
for b in labels:
pauli_pairs.append((a,b, np.kron(paulis[a], paulis[b])))
# Map location index -> pair of qubit indices
loc_qubits = {i: cnots[i] for i in range(5)}
# For each pair (a,b), build the 5-qubit operator at location i
from functools import lru_cache
@lru_cache(None)
def pair_op_on_loc(i, a, b):
q1, q2 = loc_qubits[i]
ops = [I]*nq
ops[q1] = paulis[a]
ops[q2] = paulis[b]
return kronN(ops)
# Initial state |00000>
ket0_5 = np.zeros((2**nq,1), dtype=complex); ket0_5[0,0] = 1
# Build a function to apply the circuit with a specified error pattern up to 2 faults
# error_pattern: dict mapping location index i (0..4) -> (a, b) Pauli labels (non-identity permitted; 'I','I' means no error at that site)
def run_with_errors(error_pattern):
state = ket0_5.copy()
# H
ops = [I]*nq; ops[1] = H
state = kronN(ops) @ state
# For each CNOT + error after it
for i, (c,t) in enumerate(cnots):
state = CNOTs[i] @ state
a,b = error_pattern.get(i, ('I','I'))
if not (a=='I' and b=='I'):
state = pair_op_on_loc(i,a,b) @ state
return state
# Build projection operators for data (0..3) only
# To compute <ψ| Π_a0 ⊗ P_code |ψ> etc, we'll expand them into full 5-qubit space by tensoring with ancilla projectors.
P_code_5 = kronN([None,None,None,None,None]) # placeholder
# Build tensored operators conveniently
def lift_data_op_to_5(op4):
# op4 acts on qubits 0..3, ancilla identity on qubit 4
return np.kron(op4, I)
P_code_5 = lift_data_op_to_5(P_code)
P_00_5 = lift_data_op_to_5(P_00)
# Now enumerate 0,1,2 faults, compute contributions to numerator and denominator as functions of p
# Two-qubit non-identity pairs
nonid_pairs = [(a,b) for (a,b,_) in pauli_pairs if not (a=='I' and b=='I')]
# For clarity, precompute the weights coefficients for each k when expanding to O(p^2)
# For a given pattern with k faults, its weight = (1-p)^(5-k) * (p/15)^k.
# We'll accumulate contributions for polynomial up to p^2: we represent as coeffs [c0, c1, c2] such that sum = c0 + c1 p + c2 p^2
from collections import defaultdict
num_coeff = np.array([0.0, 0.0, 0.0]) # numerator coefficients up to p^2
den_coeff = np.array([0.0, 0.0, 0.0]) # denominator coefficients up to p^2