from itertools import product
from fractions import Fraction
# Qubit indexing: 0,1,2,3 are data; 4 is ancilla
n = 5
# Helper: Pauli as (x,z) bitvectors lengths n
# Multiplication: XOR x and z (phase ignored)
def pauli_mul(p, q):
x1, z1 = p
x2, z2 = q
return [x1 ^ x2, z1 ^ z2]
# Conjugation by CNOT(c,t)
def conj_cnot(p, c, t):
x, z = p
# x_t <- x_t XOR x_c
xt = (x >> t) & 1
xc = (x >> c) & 1
x ^= (xc << t)
# z_c <- z_c XOR z_t
zt = (z >> t) & 1
zc = (z >> c) & 1
z ^= (zt << c)
return [x, z]
# Build single-qubit Paulis
label_to_bits = {
'I': (0,0),
'X': (1,0),
'Z': (0,1),
'Y': (1,1),
}
# Two-qubit Paulis on a given pair (a,b)
def two_qubit_paulis_on_pair(a, b):
ops = []
labels = ['I','X','Y','Z']
for la in labels:
for lb in labels:
if la=='I' and lb=='I':
lab = 'II'
else:
lab = la+lb
# Build 5-qubit Pauli
x = 0
z = 0
xa, za = label_to_bits[la]
xb, zb = label_to_bits[lb]
if xa: x |= (1<<a)
if xb: x |= (1<<b)
if za: z |= (1<<a)
if zb: z |= (1<<b)
ops.append((lab, [x, z]))
return ops
# Gates in temporal order with error inserted after each CNOT
# 0: H1 (no error here)
# 1: CNOT12 -> error E1 on (1,2)
# 2: CNOT10 -> error E2 on (1,0)
# 3: CNOT23 -> error E3 on (2,3)
# 4: CNOT34 -> error E4 on (3,4)
# 5: CNOT04 -> error E5 on (0,4)
gates = [
('H', 1, None),
('CNOT', 1, 2),
('CNOT', 1, 0),
('CNOT', 2, 3),
('CNOT', 3, 4),
('CNOT', 0, 4),
]
# Precompute conjugated error for each insertion point and each 2-qubit Pauli
insertion_points = [1,2,3,4,5] # after these gates
pairs = {1:(1,2), 2:(1,0), 3:(2,3), 4:(3,4), 5:(0,4)}
ops_by_point = {}
for s in insertion_points:
a,b = pairs[s]
ops = two_qubit_paulis_on_pair(a,b)
# For each Pauli, propagate through subsequent gates s+1..5
propagated = []
for lab, p in ops:
x,z = p
p_cur = [x,z]
for idx in range(s+1, 6):
g, c, t = gates[idx]
if g=='CNOT':
p_cur = conj_cnot(p_cur, c, t)
else:
# No more H gates after s>=1
pass
propagated.append((lab, p_cur))
ops_by_point[s] = propagated
# Stabilizers Sx=XXXX, Sz=ZZZZ on data qubits
Sx = [0,0]
for q in [0,1,2,3]:
Sx[0] |= (1<<q) # x bits
Sz = [0,0]
for q in [0,1,2,3]:
Sz[1] |= (1<<q) # z bits
# Logical operators
# X_A = X I X I (0,2)
XA = [0,0]
XA[0] |= (1<<0) | (1<<2)
# Z_A = Z Z I I (0,1)
ZA = [0,0]
ZA[1] |= (1<<0) | (1<<1)
# X_B = X X I I (0,1)
XB = [0,0]
XB[0] |= (1<<0) | (1<<1)
# Z_B = Z I Z I (0,2)
ZB = [0,0]
ZB[1] |= (1<<0) | (1<<2)
# Commutation check: returns 0 if commute, 1 if anticommute
def anticomm(p, q):
x1,z1 = p
x2,z2 = q
# (-1)^{x1·z2 + z1·x2}
s = ((x1 & z2).bit_count() + (z1 & x2).bit_count()) & 1
return s
# Restrict pauli to data (zero ancilla components)
def restrict_data(p):
x,z = p
# zero out ancilla bit 4
mask_data = ((1<<4)-1) # bits 0..3 set
return [x & mask_data, z & mask_data]
# acceptance checks
def accept_code(p):
pd = restrict_data(p)
return (anticomm(pd, Sx)==0) and (anticomm(pd, Sz)==0)
def accept_anc(p):
x,z = p
# accept iff no X/Y on ancilla -> x_4 == 0
return ((x>>4)&1)==0
# classify logical good: accepted and no logical X on A or B
def logical_good(p):
if not (accept_code(p) and accept_anc(p)):
return False
pd = restrict_data(p)
xA = anticomm(pd, ZA) # anticomm with Z_A indicates X_A component
xB = anticomm(pd, ZB)
return (xA==0) and (xB==0)
# Iterate all combinations across 5 insertion points, each with 16 options. Compute counts by k (# non-identity)
from collections import defaultdict
acc_counts = [0]*6
good_counts = [0]*6
all_counts = [0]*6
num_total = 0
# Precompute length 16 for each point
ops_per_point = []
labels_per_point = []
for s in insertion_points:
ops = ops_by_point[s]
# They are currently 16 ops including 'II' label
labels_per_point.append([lab for lab,_ in ops])
ops_per_point.append([p for _,p in ops])
# iterate
indices = [range(16) for _ in range(5)]
for i1 in indices[0]:
p1 = ops_per_point[0][i1]
k1 = 0 if labels_per_point[0][i1]=='II' else 1
for i2 in indices[1]:
p2 = ops_per_point[1][i2]
k2 = k1 + (0 if labels_per_point[1][i2]=='II' else 1)
p12 = pauli_mul(p1, p2)
for i3 in indices[2]:
p3 = ops_per_point[2][i3]
k3 = k2 + (0 if labels_per_point[2][i3]=='II' else 1)
p123 = pauli_mul(p12, p3)
for i4 in indices[3]:
p4 = ops_per_point[3][i4]
k4 = k3 + (0 if labels_per_point[3][i4]=='II' else 1)
p1234 = pauli_mul(p123, p4)
for i5 in indices[4]:
p5 = ops_per_point[4][i5]
k = k4 + (0 if labels_per_point[4][i5]=='II' else 1)
p_final = pauli_mul(p1234, p5)
all_counts[k] += 1
if accept_code(p_final) and accept_anc(p_final):
acc_counts[k] += 1
if logical_good(p_final):
good_counts[k] += 1
print('all_counts', all_counts)
print('acc_counts', acc_counts)
print('good_counts', good_counts)
# Sanity: total combos
print('Total combos', sum(all_counts))
# Build polynomials Numerator (good) and Denominator (acc) as exact Fractions
from math import comb
maxdeg = 5
num_coeff = [Fraction(0,1) for _ in range(maxdeg+1)]
den_coeff = [Fraction(0,1) for _ in range(maxdeg+1)]
for k in range(0,6):
# weight factor (p/15)^k * (1-p)^{5-k}
coef_base = Fraction(1, 15**k)
# expand (1-p)^{5-k}
for m in range(0, 6-k):
coeff = coef_base * Fraction(comb(5-k, m), 1) * ((-1)**m)
power = k + m
den_coeff[power] += coeff * acc_counts[k]
num_coeff[power] += coeff * good_counts[k]
print('Denominator coefficients (lowest to highest degree):')
print([f'{c.numerator}/{c.denominator}' for c in den_coeff])
print('Numerator coefficients:')
print([f'{c.numerator}/{c.denominator}' for c in num_coeff])
# Also print small-p expansion ratio up to order 2
# Compute ratio series: num/den = a0 + a1 p + a2 p^2 + ... using series division
# We'll compute to order 2
from copy import deepcopy
an = [Fraction(c) for c in num_coeff]
ad = [Fraction(c) for c in den_coeff]
# ensure ad[0] != 0
# series division
F0 = an[0]/ad[0]
# Next coefficient
F1 = (an[1] - F0*ad[1]) / ad[0]
F2 = (an[2] - F0*ad[2] - F1*ad[1]) / ad[0]
print('Series coefficients F(p) = F0 + F1 p + F2 p^2 + ...:')
print(F0, F1, F2)