import itertools
from math import comb
# Pauli binary symplectic representation utilities
# Represent Pauli on n qubits as (x,z) binary lists length n, ignoring global phase.
I = (0,0)
X = (1,0)
Z = (0,1)
Y = (1,1)
pauli_letter_to_xz = {'I':I,'X':X,'Y':Y,'Z':Z}
letters = ['I','X','Y','Z']
n=5 # qubits 0..3 data, 4 ancilla
# Conjugation by CNOT(c,t):
def conj_cnot(x,z,c,t):
# x,z are lists of length n
x = x.copy(); z = z.copy()
# Update: x_t' = x_t XOR x_c; z_c' = z_c XOR z_t
x[t] ^= x[c]
z[c] ^= z[t]
return x,z
# Multiply two Paulis: XOR x and z
def multiply_pauli(x1,z1,x2,z2):
x = [(a^b) for a,b in zip(x1,x2)]
z = [(a^b) for a,b in zip(z1,z2)]
return x,z
# Build two-qubit Pauli on pair (i,j) given letters li, lj, embedded in n qubits
def twoq_pauli_xz(n, i, j, li, lj):
x = [0]*n
z = [0]*n
xi,zi = pauli_letter_to_xz[li]
xj,zj = pauli_letter_to_xz[lj]
x[i]=xi; z[i]=zi
x[j]=xj; z[j]=zj
return x,z
# Define gates sequence (CNOTs) indices in execution order (right to left in provided formula)
# The list of (control, target)
sequence = [(1,2),(1,0),(2,3),(3,4),(0,4)]
# Precompute propagated error for each gate and each 2q Pauli (including identity)
propagated = [] # list for each gate: dict mapping (li, lj) -> (x,z)
nonid_flag = [] # list for each gate: dict mapping (li, lj) -> 0/1
for g_idx,(c,t) in enumerate(sequence):
# subsequent gates to conjugate through are the ones after g_idx
subs = sequence[g_idx+1:]
table = {}
flags = {}
for li in letters:
for lj in letters:
x,z = twoq_pauli_xz(n, c if c< t else t, t if c< t else c, 'I','I')
# Actually need to place on qubits c and t specifically
x,z = twoq_pauli_xz(n, c, t, li, lj)
# Conjugate forward through subsequent gates
for (cc,tt) in subs:
x,z = conj_cnot(x,z,cc,tt)
table[(li,lj)] = (tuple(x), tuple(z))
flags[(li,lj)] = 0 if (li=='I' and lj=='I') else 1
propagated.append(table)
nonid_flag.append(flags)
# Build GHZ_4 stabilizer group on data qubits 0..3
# Generators g1=XXXX, g2=Z0Z1, g3=Z1Z2, g4=Z2Z3
def pauli_from_gen_list(gens):
# gens is list of (x,z)
# return set of (x_data, z_data) tuples for 4-qubit data
m = len(gens)
elems = set()
for mask in range(1<<m):
x = [0]*4
z = [0]*4
for i in range(m):
if (mask>>i)&1:
gx, gz = gens[i]
x = [a^b for a,b in zip(x,gx)]
z = [a^b for a,b in zip(z,gz)]
elems.add((tuple(x), tuple(z)))
return elems
# Generators in 4-qubit space
x_g1 = ([1,1,1,1],[0,0,0,0])
z_g2 = ([0,0,0,0],[1,1,0,0])
z_g3 = ([0,0,0,0],[0,1,1,0])
z_g4 = ([0,0,0,0],[0,0,1,1])
GHZ_stab = pauli_from_gen_list([x_g1, z_g2, z_g3, z_g4])
# For sanity, size should be 16
assert len(GHZ_stab)==16
# Acceptance checks
# Code acceptance: parity of x and z on data are even.
# Ancilla acceptance: x[4]==0
def accept(x,z):
# data parity even conditions
if (sum(z[q] for q in range(4))%2)!=0: return False
if (sum(x[q] for q in range(4))%2)!=0: return False
# ancilla
if x[4]!=0: return False
return True
# Goodness: data part in GHZ stabilizer
def good(x,z):
x_d = tuple(x[:4])
z_d = tuple(z[:4])
return (x_d, z_d) in GHZ_stab
# Enumerate over all 16^5 combinations and count
from collections import defaultdict
A_counts = [0]*6 # good accepted counts by k non-identity errors
B_counts = [0]*6 # accepted counts by k
# We'll also compute exact polynomials if needed later
choices_per_gate = [list(propagated[g].keys()) for g in range(5)]
# We'll pre-extract (x,z) and flag for speed
prop_lists = []
flag_lists = []
for g in range(5):
prop_lists.append([propagated[g][key] for key in choices_per_gate[g]])
flag_lists.append([nonid_flag[g][key] for key in choices_per_gate[g]])
# For speed: number of combinations is 1,048,576 loops; should be OK
count_total = 0
for idx0 in range(16):
x0,z0 = prop_lists[0][idx0]
f0 = flag_lists[0][idx0]
for idx1 in range(16):
x1,z1 = prop_lists[1][idx1]
f1 = flag_lists[1][idx1]
# mul 0*1
x01 = [a^b for a,b in zip(x0,x1)]
z01 = [a^b for a,b in zip(z0,z1)]
k01 = f0+f1
for idx2 in range(16):
x2,z2 = prop_lists[2][idx2]
f2 = flag_lists[2][idx2]
x012 = [a^b for a,b in zip(x01,x2)]
z012 = [a^b for a,b in zip(z01,z2)]
k012 = k01+f2
for idx3 in range(16):
x3,z3 = prop_lists[3][idx3]
f3 = flag_lists[3][idx3]
x0123 = [a^b for a,b in zip(x012,x3)]
z0123 = [a^b for a,b in zip(z012,z3)]
k0123 = k012+f3
for idx4 in range(16):
x4,z4 = prop_lists[4][idx4]
f4 = flag_lists[4][idx4]
x_final = [a^b for a,b in zip(x0123,x4)]
z_final = [a^b for a,b in zip(z0123,z4)]
k = k0123+f4
if accept(x_final, z_final):
B_counts[k]+=1
if good(x_final, z_final):
A_counts[k]+=1
count_total+=1
print("Total combos:", count_total)
print("Accepted counts B_k:", B_counts)
print("Good accepted counts A_k:", A_counts)
# Now construct P_acc(p) and P_good(p) polynomials: sum_k B_k (1-p)^{5-k} (p/15)^k
# We'll expand into polynomial in p up to degree 5: (1-p)^{5-k} * p^k = sum_{t=0}^{5-k} (-1)^t C(5-k, t) p^{t+k}
# Thus P_acc(p) = sum_{m=0}^{5} c_m p^m with coefficients rational numbers with denominator 15^{m'} ??? We'll keep as fractions of 15^k.
from fractions import Fraction
# Compute coefficients as Fractions
# Denominator due to 15^k; We'll keep them separate or convert to rational function form later
# We'll compute P(p) = sum_m coeff[m] * p^m, where coeff[m] is Fraction
def poly_coeffs_from_counts(counts):
coeff = [Fraction(0,1) for _ in range(6)]
for k in range(6):
Bk = counts[k]
# weight factor: (1)^(5-k) for constant; then expansion:
for t in range(0, 6-k):
m = t + k
# contribution: Bk * C(5-k, t) * (-1)^t * 1 * (1) / 15^k
val = Fraction(Bk * comb(5-k, t) * ((-1)**t), 15**k)
coeff[m] += val
return coeff
coeff_acc = poly_coeffs_from_counts(B_counts)
coeff_good = poly_coeffs_from_counts(A_counts)
print("P_acc polynomial coeffs (p^0..p^5):", coeff_acc)
print("P_good polynomial coeffs (p^0..p^5):", coeff_good)
# Simplify to decimal or rational forms; Also compute F_logical(p) = P_good / P_acc as rational function.
# For presentation, we can give both series expansion up to p^2 and exact rational function via quotient of polynomials.
# Let's also compute exact rational function numerator and denominator polynomials with integer coefficients by clearing denominators.
def scale_to_common_denom(coeffs):
# coeffs is list of Fractions
from math import lcm
denoms = [c.denominator for c in coeffs]
D = 1
for d in denoms:
from math import gcd
D = (D*d)//__import__('math').gcd(D,d)
ints = [int(c*D) for c in coeffs]
return D, ints
D_acc, poly_acc = scale_to_common_denom(coeff_acc)
D_good, poly_good = scale_to_common_denom(coeff_good)
print("Acc denom:", D_acc)
print("Acc integer poly:", poly_acc)
print("Good denom:", D_good)
print("Good integer poly:", poly_good)
# Now F(p) = (poly_good/D_good)/(poly_acc/D_acc) = (poly_good * D_acc) / (poly_acc * D_good)
# We'll try to reduce common factor if any
from math import gcd
def poly_gcd_int(a,b):
# naive gcd of integer polynomials by content GCD (not polynomial gcd), we'll just reduce by gcd of coefficients
g = 0
for v in a:
g = gcd(g, abs(v))
h = 0
for v in b:
h = gcd(h, abs(v))
return g,h
# We'll also compute small-p expansion of F(p) up to p^3 by series division
import sympy as sp
p=sp.symbols('p')
Poly_acc = sum(sp.Rational(c.numerator,c.denominator)*p**i for i,c in enumerate(coeff_acc))
Poly_good = sum(sp.Rational(c.numerator,c.denominator)*p**i for i,c in enumerate(coeff_good))
F_series = sp.series(Poly_good/Poly_acc, p, 0, 4) # up to p^3
print("F_series:", F_series)
# Try to express F as rational function with (small) integer coefficients by multiplying numerator and denominator by LCM denominators
F_simplified = sp.simplify(Poly_good/Poly_acc)
print("F_simplified rational:", sp.together(F_simplified))