from itertools import product
from fractions import Fraction
# Number of qubits: 5 (0..3 data, 4 ancilla)
N = 5
def pauli_vec_from_label_on_pair(q1, q2, label):
# label is a pair like 'IX','ZY', etc on the ordered tuple (q1,q2)
# returns (x,z) as lists of 0/1 length N
x = [0]*N
z = [0]*N
for idx, q in enumerate((q1, q2)):
ch = label[idx]
if ch == 'I':
pass
elif ch == 'X':
x[q] ^= 1
elif ch == 'Z':
z[q] ^= 1
elif ch == 'Y':
x[q] ^= 1
z[q] ^= 1
else:
raise ValueError('bad')
return tuple(x), tuple(z)
# Conjugation by CNOT(c->t):
# x_c' = x_c; x_t' = x_t ^ x_c
# z_c' = z_c ^ z_t; z_t' = z_t
def conj_cnot(x, z, c, t):
x = list(x)
z = list(z)
# x_t ^= x_c
x[t] ^= x[c]
# z_c ^= z_t
z[c] ^= z[t]
return tuple(x), tuple(z)
# Multiplying two Paulis ignoring phase: XOR the x,z bits
def mul_paulixz(a, b):
x1, z1 = a
x2, z2 = b
return tuple((x1[i]^x2[i]) for i in range(N)), tuple((z1[i]^z2[i]) for i in range(N))
# Build all 2-qubit Pauli labels on a pair
singles = ['I','X','Y','Z']
labels_pair = [a+b for a in singles for b in singles]
# Remove 'II' later when needed
# Gates in chronological order (right-to-left execution):
# After each of these CNOTs there is a noise channel
cnot_list = [(1,2), (1,0), (2,3), (3,4), (0,4)]
num_gates = len(cnot_list)
# Precompute pushforward of each pair-Pauli at each location to the end
pushforward = [] # list of dict: label->(x,z)
for i, (c,t) in enumerate(cnot_list):
mp = {}
for lab in labels_pair:
# Build Pauli supported on the pair
x,z = pauli_vec_from_label_on_pair(c,t,lab)
# Conjugate by subsequent CNOTs (i+1..end)
for (c2,t2) in cnot_list[i+1:]:
x,z = conj_cnot(x,z,c2,t2)
mp[lab] = (x,z)
pushforward.append(mp)
# Build the 16-element stabilizer group of the target logical |00>_AB on data qubits
# Generators: XXXX, ZZZZ, ZZII, ZIZI on qubits 0..3
def pauli_from_str_data(s):
# s is 4-char over IXYZ for data 0..3
x = [0]*N
z = [0]*N
for q,ch in enumerate(s):
if ch=='I':
pass
elif ch=='X':
x[q]=1
elif ch=='Z':
z[q]=1
elif ch=='Y':
x[q]=1; z[q]=1
return (tuple(x), tuple(z))
# Data-only Paulis (ancilla bits zero)
S1 = pauli_from_str_data('XXXX')
S2 = pauli_from_str_data('ZZZZ')
S3 = pauli_from_str_data('ZZII')
S4 = pauli_from_str_data('ZIZI')
# Generate group
T = set()
for a in [0,1]:
for b in [0,1]:
for c in [0,1]:
for d in [0,1]:
x = [0]*N; z = [0]*N
for (flag, P) in ((a,S1),(b,S2),(c,S3),(d,S4)):
if flag:
x = [x[i]^P[0][i] for i in range(N)]
z = [z[i]^P[1][i] for i in range(N)]
T.add((tuple(x),tuple(z)))
# Sanity: should be 16 elements
# print('Size T:', len(T))
# Functions to test acceptance
def ancilla_accept(x,z):
# accept if no X (or Y) on qubit 4, i.e., x[4]==0
return (x[4] == 0)
def code_accept(x,z):
# accept if data X parity even and data Z parity even
sx = sum(x[:4]) % 2
sz = sum(z[:4]) % 2
return (sx==0 and sz==0)
def good_data(x,z):
# Check if data part (qubits 0..3) lies in T
# Zero out ancilla bits
x0 = list(x); z0 = list(z)
x0[4]=0; z0[4]=0
return (tuple(x0),tuple(z0)) in T
# Now enumerate all error combinations. We will accumulate coefficients of polynomial in p
# Probability weight for each per-gate error choice is (1-p) if 'II', else (p/15)
# The total accepted probability is sum_k coef[k] * p^k, where coef[k] rational numbers
from collections import defaultdict
coef_pass = defaultdict(Fraction)
coef_good = defaultdict(Fraction)
# Precompute list of labels and their (k, weight_factor)
# Represent (k, c) where k = number of non-identity errors contributed (0 or 1 per gate), and c = fractional coefficient (1 for identity, 1/15 for non-identity)
choices = []
for lab in labels_pair:
if lab == 'II':
choices.append((lab, 0, Fraction(1,1)))
else:
choices.append((lab, 1, Fraction(1,15)))
# For speed, we can precompute for each gate i and each label the pushed (x,z)
# Already in pushforward
# Enumerate
cnt=0
for labs in product(range(len(choices)), repeat=num_gates):
# Build final E_total
x = tuple([0]*N)
z = tuple([0]*N)
k_total = 0
weight = Fraction(1,1)
for i, choice_idx in enumerate(labs):
lab, k, fac = choices[choice_idx]
k_total += k
weight *= fac
# pushforward to end
x_i, z_i = pushforward[i][lab]
# multiply
x = tuple((x[j]^x_i[j]) for j in range(N))
z = tuple((z[j]^z_i[j]) for j in range(N))
# Now acceptance tests
if ancilla_accept(x,z) and code_accept(x,z):
coef_pass[k_total] += weight
if good_data(x,z):
coef_good[k_total] += weight
cnt += 1
# Now build polynomials P_pass(p) = sum_k coef_pass[k]*(1-p)^{5-k}*(p)^k factoring 1/15^k; But we already included 1/15^k in weight.
# Actually our weight is product over gates of either 1 for identity or 1/15 for non-identity, but missing the (1-p) and p factors.
# The total probability contribution of all cases with exactly k non-identity errors is weight * (1-p)^{5-k} * p^k.
# Let's express polynomials as rational coefficients for powers of p
# Expand sum_k A_k * (1-p)^{5-k} * p^k where A_k = coef_pass[k]
from math import comb
def expand_poly(A):
# A: dict k->Fraction coefficient
# Returns dict deg->Fraction for polynomial in p up to degree 5
P = defaultdict(Fraction)
for k, Ak in A.items():
# (1-p)^{5-k} * p^k = sum_{m=0}^{5-k} C(5-k,m)*(-p)^m * p^k = sum_{m} C(5-k,m)*(-1)^m * p^{k+m}
for m in range(0, 5-k+1):
deg = k + m
coeff = Ak * Fraction(comb(5-k, m), 1) * (Fraction(-1,1) if (m%2==1) else Fraction(1,1))
P[deg] += coeff
return P
P_pass = expand_poly(coef_pass)
P_good = expand_poly(coef_good)
# Convert to simplified Fractions
# Show polynomials
maxdeg = 5
#print('P_pass coeffs by degree:')
#for d in range(maxdeg+1):
# print(d, P_pass.get(d, Fraction(0,1)))
#print('P_good coeffs by degree:')
#for d in range(maxdeg+1):
# print(d, P_good.get(d, Fraction(0,1)))
# Let's format as rational numbers
P_pass_list = [P_pass.get(d, Fraction(0,1)) for d in range(maxdeg+1)]
P_good_list = [P_good.get(d, Fraction(0,1)) for d in range(maxdeg+1)]
P_pass_list, P_good_list, sum(P_pass_list)