from fractions import Fraction
import numpy as np
# 1) Build circuit state and code projector
I = np.eye(2, 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)
H = (1/np.sqrt(2))*np.array([[1,1],[1,-1]], dtype=complex)
n=5
def kronN(ops):
out = ops[0]
for op in ops[1:]:
out = np.kron(out, op)
return out
def gate_single(U, q):
ops = [I]*n
ops[q] = U
return kronN(ops)
def gate_cnot(c,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
ops1[t] = X
return kronN(ops0) + kronN(ops1)
U = np.eye(2**n, dtype=complex)
U = gate_single(H,1) @ U
U = gate_cnot(1,2) @ U
U = gate_cnot(1,0) @ U
U = gate_cnot(2,3) @ U
U = gate_cnot(3,4) @ U
U = gate_cnot(0,4) @ U
psi0 = np.zeros((2**n,), dtype=complex)
psi0[0] = 1
psi = U @ psi0
phi = psi.reshape((16,2))[:,0]
I4 = np.eye(16, dtype=complex)
X1 = np.array([[0,1],[1,0]], dtype=complex)
Z1 = np.array([[1,0],[0,-1]], dtype=complex)
X4 = kronN([X1,X1,X1,X1])
Z4 = kronN([Z1,Z1,Z1,Z1])
P_code = 0.5*(I4 + X4) @ (0.5*(I4 + Z4))
phi_code = P_code @ phi
phi_code /= np.sqrt((phi_code.conj()@phi_code).real)
# 2) Precompute s_E and t_E for all 1024 Paulis
single_paulis = [np.eye(2, dtype=complex), X, Y, Z]
pauli_mats = []
for a in range(4):
for b in range(4):
for c in range(4):
for d in range(4):
for e in range(4):
M = kronN([single_paulis[a], single_paulis[b], single_paulis[c], single_paulis[d], single_paulis[e]])
pauli_mats.append(M)
s = np.zeros(4**5)
t = np.zeros(4**5)
for idx, E in enumerate(pauli_mats):
state = E @ psi
state_a0 = state.reshape((16,2))[:,0]
state_code = P_code @ state_a0
t[idx] = (state_code.conj()@state_code).real
if t[idx] > 1e-14:
vec = state_code/np.sqrt(t[idx])
s[idx] = np.abs((phi_code.conj()@vec))**2
else:
s[idx] = 0.0
# 3) Build Pauli bit representation (x,z) and index mapping
# mapping digits: 0->I, 1->X, 2->Y, 3->Z
pow4 = [4**k for k in range(5)][::-1]
# order of qubits: q0,q1,q2,q3,q4
def idx_to_digits(idx):
d=[]
for p in pow4:
d.append(idx//p)
idx%=p
return d
def digits_to_idx(d):
idx=0
for i in range(5):
idx = idx*4 + d[i]
return idx
# maps per digit to (x,z) bits
x_per = [0,1,1,0]
z_per = [0,0,1,1]
def idx_to_bits(idx):
d = idx_to_digits(idx)
xb=0; zb=0
for i in range(5):
if x_per[d[i]]: xb |= (1<<(4-i))
if z_per[d[i]]: zb |= (1<<(4-i))
return xb, zb
# idx from bits
def bits_to_idx(xb,zb):
d=[0]*5
for i in range(5):
bitpos = 4-i
x = (xb>>bitpos)&1
z = (zb>>bitpos)&1
if x==0 and z==0: d[i]=0
elif x==1 and z==0: d[i]=1
elif x==1 and z==1: d[i]=2
elif x==0 and z==1: d[i]=3
return digits_to_idx(d)
# 4) Define CNOT conjugation on bits (control c, target t)
def conj_cnot_bits(xb,zb,c,t):
# Toggle x target by x control
xc = (xb>> (4-c)) & 1
if xc:
xb ^= (1<< (4-t))
# Toggle z control by z target
zt = (zb>> (4-t)) & 1
if zt:
zb ^= (1<< (4-c))
return xb,zb
# 5) Build per-gate error sets (E_final idx, poly)
# Gates: after CNOT12, after CNOT10, after CNOT23, after CNOT34, after CNOT04
cnots = [(1,2),(1,0),(2,3),(3,4),(0,4)]
# For 2-qubit Paulis on (i,j), generate 16 combinations
single_digits = [0,1,2,3] # I,X,Y,Z
def pauli_bits_for_two_qubit(i,j, di, dj):
# returns initial (xb,zb) having di at qubit i and dj at qubit j
xb=0; zb=0
x_i = x_per[di]; z_i = z_per[di]
x_j = x_per[dj]; z_j = z_per[dj]
if x_i: xb |= (1<<(4-i))
if z_i: zb |= (1<<(4-i))
if x_j: xb |= (1<<(4-j))
if z_j: zb |= (1<<(4-j))
return xb,zb
# Build subsequent gates lists for each insertion point
subseq = {
0: [1,2,3,4],
1: [2,3,4],
2: [3,4],
3: [4],
4: [],
}
# polynomial operations with Fraction
def poly_zero():
return [Fraction(0,1)]*6
def poly_add(a,b):
return [ai+bi for ai,bi in zip(a,b)]
def poly_mul(a,b):
c = [Fraction(0,1)]*6
for i in range(6):
if a[i]==0: continue
for j in range(6-i):
if b[j]==0: continue
c[i+j] += a[i]*b[j]
return c
# Polynomials for identity and error per gate
poly_id = [Fraction(1,1), Fraction(-1,1)] + [Fraction(0,1)]*4 # 1 - p
poly_err = [Fraction(0,1), Fraction(1,15)] + [Fraction(0,1)]*4 # p/15
# For each gate g (0..4), build list of (idx_final, poly)
per_gate_sets = []
for g, (c,tg) in enumerate(cnots):
items = []
# identity
items.append((0, poly_id))
for di in single_digits:
for dj in single_digits:
if di==0 and dj==0: continue
xb,zb = pauli_bits_for_two_qubit(c,tg, di, dj)
# conjugate through subsequent CNOTs
for h in subseq[g]:
cc,tt = cnots[h]
xb,zb = conj_cnot_bits(xb,zb, cc,tt)
idx_final = bits_to_idx(xb,zb)
items.append((idx_final, poly_err))
per_gate_sets.append(items)
# Convolve distributions across gates
from collections import defaultdict
dist = {0: [Fraction(1,1)] + [Fraction(0,1)]*5} # start with identity with weight 1
for g in range(5):
newdist = defaultdict(lambda: [Fraction(0,1)]*6)
for idx_prev, poly_prev in dist.items():
for idx_g, poly_g in per_gate_sets[g]:
# combine Paul's: XOR in bitspace -> compute idx_total by mapping bits; we can derive via digits: using mapping idx XOR function
# We'll use bits for idx_prev and idx_g
# Precompute bits for these? compute on the fly
# Convert idx to bits
# We'll compute using helper functions
# But we need a fast XOR composition; we implement idx_to_bits for both and compute bits_to_idx of XOR
xb_prev, zb_prev = 0,0
# We'll use cached bit maps to avoid recomputing each time; Let's build caches outside loops
# We'll rebuild with caches outside main loop
break
# Build caches for idx->bits and bits->idx mapping arrays
idx_to_bits_cache = {}
for idx in range(4**5):
idx_to_bits_cache[idx] = idx_to_bits(idx)
# Now convolve properly
dist = {0: [Fraction(1,1)] + [Fraction(0,1)]*5}
for g in range(5):
newdist = defaultdict(lambda: [Fraction(0,1)]*6)
for idx_prev, poly_prev in dist.items():
xb_prev, zb_prev = idx_to_bits_cache[idx_prev]
for idx_g, poly_g in per_gate_sets[g]:
xb_g, zb_g = idx_to_bits_cache[idx_g]
xb_tot = xb_prev ^ xb_g
zb_tot = zb_prev ^ zb_g
idx_tot = bits_to_idx(xb_tot, zb_tot)
poly_tot = poly_mul(poly_prev, poly_g)
cur = newdist[idx_tot]
newdist[idx_tot] = poly_add(cur, poly_tot)
dist = newdist
# Verify normalization: sum of polynomials over all idx equals 1
sum_poly = [Fraction(0,1)]*6
for poly in dist.values():
sum_poly = poly_add(sum_poly, poly)
# Evaluate sum_poly at some p to check equals 1 for all p: coefficients should be [1,0,0,0,0,0]
print('sum_poly coefficients:', sum_poly)
# 6) Compute numerator and denominator polynomials
import math
# Build s,t arrays again to avoid missing
# (But they may not require exact coefficients; they are 0/1 but we computed numerically earlier; We will now compute them logically using phi vector approach again)
# We'll recompute numeric s and t within this run to use them
single_paulis2 = [np.eye(2, dtype=complex), np.array([[0,1],[1,0]], dtype=complex), np.array([[0,-1j],[1j,0]], dtype=complex), np.array([[1,0],[0,-1]], dtype=complex)]
# Build psi again to have it in this execution (already above). We'll reuse 'psi', 'phi_code', 'P_code'
s_arr = np.zeros(4**5)
t_arr = np.zeros(4**5)
for idx in range(4**5):
# reconstruct matrix for idx
d = idx_to_digits(idx)
E = kronN([single_paulis2[d[0]], single_paulis2[d[1]], single_paulis2[d[2]], single_paulis2[d[3]], single_paulis2[d[4]]])
state = E @ psi
state_a0 = state.reshape((16,2))[:,0]
state_code = P_code @ state_a0
t_arr[idx] = (state_code.conj()@state_code).real
if t_arr[idx] > 1e-12:
vec = state_code/np.sqrt(t_arr[idx])
s_arr[idx] = np.abs((phi_code.conj()@vec))**2
else:
s_arr[idx] = 0.0
# Identify accepted and correct sets
accepted_idxs = [idx for idx in range(4**5) if t_arr[idx] > 1e-12]
correct_idxs = [idx for idx in range(4**5) if t_arr[idx] > 1e-12 and abs(s_arr[idx]-1.0)<1e-12]
print('counts: accepted', len(accepted_idxs), 'correct', len(correct_idxs))
# Compute numerator/denominator polynomials
num_poly = [Fraction(0,1)]*6
den_poly = [Fraction(0,1)]*6
for idx, poly in dist.items():
if t_arr[idx] > 1e-12:
den_poly = poly_add(den_poly, poly)
if abs(s_arr[idx]-1.0) < 1e-12:
num_poly = poly_add(num_poly, poly)
print('den_poly:', den_poly)
print('num_poly:', num_poly)
# Simplify: print as fractions and also as decimals
def poly_to_str(poly):
terms=[]
for k,coef in enumerate(poly):
if coef == 0:
continue
# format fraction
if coef.denominator==1:
coeff_str = f"{coef.numerator}"
else:
coeff_str = f"{coef.numerator}/{coef.denominator}"
if k==0:
terms.append(coeff_str)
elif k==1:
terms.append(f"({coeff_str}) p")
else:
terms.append(f"({coeff_str}) p^{k}")
return ' + '.join(terms) if terms else '0'
print('Denominator polynomial D(p) =', poly_to_str(den_poly))
print('Numerator polynomial N(p) =', poly_to_str(num_poly))
# Optionally, compute series expansion for F(p) = N(p)/D(p) to p^3
from fractions import Fraction
# We'll compute rational series expansion up to order 3 via formal power series
maxdeg=5
# F = a/b; build series by polynomial long division
def poly_series_div(num, den, order):
# ensure den[0] != 0; here den[0]=1
res=[Fraction(0,1)]*(order+1)
rem=num.copy()
for k in range(order+1):
res[k] = rem[k]/den[0]
# subtract res[k]*den shifted by k
for j in range(1, order-k+1):
rem[k+j] -= res[k]*den[j]
return res
series = poly_series_div(num_poly, den_poly, 5)
print('Series coefficients for F(p):', series)