# We'll enumerate error patterns for the given circuit and compute conditional logical fidelity.
# 5 data qubits (0..3) + ancilla (4)
# Circuit order (executed right-to-left): H1, CNOT12, CNOT10, CNOT23, CNOT34, CNOT04, M4
# Errors after each CNOT: E1 after CNOT12, E2 after CNOT10, E3 after CNOT23, E4 after CNOT34, E5 after CNOT04.
from itertools import product
n=5
# Bit masks
DATA_MASK = (1<<4)-1 # qubits 0..3
ANC_BIT = 1<<4
# Symplectic representation helpers
# Pauli mapping: 0:I, 1:X, 2:Y, 3:Z
px = [0,1,1,0]
pz = [0,0,1,1]
# Conjugation by CNOT c->t on symplectic bits (x,z) integers
def conj_cnot(x,z,c,t):
# x_t' = x_t XOR x_c
# z_c' = z_c XOR z_t
xc = (x>>c)&1
xt = (x>>t)&1
zc = (z>>c)&1
zt = (z>>t)&1
# update bits
xt_new = xt ^ xc
zc_new = zc ^ zt
# apply updates
if xt_new != xt:
x ^= (1<<t)
if zc_new != zc:
z ^= (1<<c)
return x,z
# Build suffix gate lists for each error location (after that gate)
# order: after CNOT12 -> gates: CNOT10, CNOT23, CNOT34, CNOT04
suffixes = [
[(1,0),(2,3),(3,4),(0,4)],
[(2,3),(3,4),(0,4)],
[(3,4),(0,4)],
[(0,4)],
[]
]
# For each location, map 16 2-qubit Paulis on its gate qubits to final (x,z)
# gate qubit pairs in order
pairs = [(1,2),(1,0),(2,3),(3,4),(0,4)]
# Precompute mapping arrays E_maps[loc][p_idx] -> (x,z)
E_maps = []
for loc,(c,t) in enumerate(pairs):
arr = []
for a in range(4):
for b in range(4):
# Build initial x,z
x = 0
z = 0
if a:
x |= px[a]<<c
z |= pz[a]<<c
if b:
x |= px[b]<<t
z |= pz[b]<<t
# propagate through suffix
for (cc,tt) in suffixes[loc]:
x,z = conj_cnot(x,z,cc,tt)
arr.append((x,z))
E_maps.append(arr)
# Build stabilizer state generators for |00>_L on data qubits: S_X=XXXX, S_Z=ZZZZ, Z_A=ZZII, Z_B=ZIZI
# Represented as (x,z) on 4 data qubits; ancilla qubit ignored here
SX_x = 0b1111; SX_z = 0
SZ_x = 0; SZ_z = 0b1111
ZA_x = 0; ZA_z = (1<<0)|(1<<1) # Z0 Z1
ZB_x = 0; ZB_z = (1<<0)|(1<<2) # Z0 Z2
# Generate the 16-element stabilizer group of the GHZ |00>_L state on 4 qubits
stab_set = set()
# We can generate by combining the four independent generators
gens = [(SX_x,SX_z),(SZ_x,SZ_z),(ZA_x,ZA_z),(ZB_x,ZB_z)]
# enumerate all 16 products ignoring phase; XOR for bits
for s0 in [0,1]:
for s1 in [0,1]:
for s2 in [0,1]:
for s3 in [0,1]:
x=z=0
if s0:
x ^= gens[0][0]; z ^= gens[0][1]
if s1:
x ^= gens[1][0]; z ^= gens[1][1]
if s2:
x ^= gens[2][0]; z ^= gens[2][1]
if s3:
x ^= gens[3][0]; z ^= gens[3][1]
stab_set.add((x,z))
# Now enumerate all 16^5 patterns
from collections import defaultdict
N_acc = [0]*6 # counts of accepted patterns by k errors
N_good = [0]*6
N_total_by_k = [0]*6
N_acc_total = 0
N_good_total = 0
# We'll also verify single-fault detection property by counting accepted k=1 events
for a0 in range(16):
x0,z0 = E_maps[0][a0]
for a1 in range(16):
x1,z1 = E_maps[1][a1]
x01 = x0 ^ x1
z01 = z0 ^ z1
for a2 in range(16):
x2,z2 = E_maps[2][a2]
x012 = x01 ^ x2
z012 = z01 ^ z2
for a3 in range(16):
x3,z3 = E_maps[3][a3]
x0123 = x012 ^ x3
z0123 = z012 ^ z3
for a4 in range(16):
x4,z4 = E_maps[4][a4]
xfin = x0123 ^ x4
zfin = z0123 ^ z4
# number of non-identity errors k
k = (a0!=0) + (a1!=0) + (a2!=0) + (a3!=0) + (a4!=0)
N_total_by_k[k] += 1
# acceptance: ancilla X must be 0, and commute with SX,SZ on data
anc_x = (xfin>>4)&1
if anc_x:
continue
# data bits
xdat = xfin & DATA_MASK
zdat = zfin & DATA_MASK
# commute with SX: parity of zdat must be 0
# commute with SZ: parity of xdat must be 0
if (bin(zdat).count('1') & 1) != 0:
continue
if (bin(xdat).count('1') & 1) != 0:
continue
# accepted
N_acc[k] += 1
# good if (xdat,zdat) in stab_set
if (xdat,zdat) in stab_set:
N_good[k] += 1
# Summaries
print("N_total_by_k:", N_total_by_k)
print("N_acc by k:", N_acc)
print("N_good by k:", N_good)
# Now produce symbolic numerator and denominator coefficients for F(p)
# Prob per pattern with k errors: (1-p)^(5-k)*(p/15)^k
# So N_good(p) = sum_k N_good[k] * (1-p)^(5-k) * (p/15)^k
# P_acc(p) = sum_k N_acc[k] * (1-p)^(5-k) * (p/15)^k
# We'll express as polynomials in p up to degree 5. Let q=1-p.
from fractions import Fraction
# We'll compute exact rational coefficients for numerator and denominator in terms of p^k.
# Expand (1-p)^(5-k) p^k * (1/15)^k = sum_{j=0}^{5-k} C(5-k,j) (-p)^j p^k * 15^{-k}
# = sum_{t=k}^{5} C(5-k, t-k) (-1)^{t-k} p^t * 15^{-k}
# We'll compute arrays coeff_num[t], coeff_den[t] rational numbers.
import math
coeff_num = [Fraction(0,1) for _ in range(6)]
coeff_den = [Fraction(0,1) for _ in range(6)]
for k in range(6):
for t in range(k,6):
comb = math.comb(5-k, t-k)
factor = Fraction(((-1)**(t-k)) * comb, 15**k)
coeff_num[t] += Fraction(N_good[k],1) * factor
coeff_den[t] += Fraction(N_acc[k],1) * factor
print("Numerator coeffs (p^t):", coeff_num)
print("Denominator coeffs (p^t):", coeff_den)
# Optionally, display F(p) rational polynomial: [We'll present as quotient of these polynomials]
# Let's also produce simplified rational function as a ratio of polynomials with integer coefficients by multiplying out denominators gcd.
def lcm(a,b):
from math import gcd
return a*b//gcd(a,b)
# Find common denominator for numerator and denominator polynomials
from math import gcd
D_num = 1
D_den = 1
for c in coeff_num:
D_num = lcm(D_num, c.denominator)
for c in coeff_den:
D_den = lcm(D_den, c.denominator)
int_num = [int(c*D_num) for c in coeff_num]
int_den = [int(c*D_den) for c in coeff_den]
print("Integer-coeff numerator (over {}):".format(D_num), int_num)
print("Integer-coeff denominator (over {}):".format(D_den), int_den)
# For sanity, compute small-p expansion of F(p) to p^2 or higher by series division
import decimal
def series_div(a,b):
# compute series c such that sum a_i p^i / sum b_i p^i = sum c_i p^i up to len(a)-1
n = len(a)
c = [0.0]*n
c[0] = a[0]/b[0]
for i in range(1,n):
s = a[i]
for j in range(i):
s -= c[j]*b[i-j]
c[i] = s/b[0]
return c
# Convert to floats for series
an = [float(c) for c in coeff_num]
bd = [float(c) for c in coeff_den]
ser = series_div(an,bd)
print("Series coefficients (to p^5):", ser)