import itertools
import numpy as np
from fractions import Fraction
# We'll compute exact rational arithmetic using Fractions where possible, but matrices will be complex floats.
# For fidelity polynomials, we will accumulate exact Fractions for A_k and G_k (sums of acceptance probabilities and weighted fidelities)
# Define single-qubit Paulis
I = np.array([[1,0],[0,1]], 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)
paulis = {'I': I, 'X': X, 'Y': Y, 'Z': Z}
pauli_labels = ['I','X','Y','Z']
# Construct n-qubit kron
def kron_n(ops):
out = np.array([[1]], dtype=complex)
for op in ops:
out = np.kron(out, op)
return out
# Build an n-qubit operator that applies 'op' on qubit q (0-based) and identity elsewhere; qubit order 0..n-1
def op_on(op, q, n):
ops = [I]*n
ops[q] = op
return kron_n(ops)
# Build two-qubit CNOT i->j on n qubits
def CNOT_ij(i, j, n):
# construct using projector |0><0| on control plus |1><1| with X on target
P0 = np.array([[1,0],[0,0]], dtype=complex)
P1 = np.array([[0,0],[0,1]], dtype=complex)
return (op_on(P0, i, n) @ np.eye(2**n) + op_on(P1, i, n) @ op_on(X, j, n))
# Two-qubit Pauli on qubits (i,j) with labels a,b in {I,X,Y,Z}
def two_pauli(a, i, b, j, n):
ops = [I]*n
ops[i] = paulis[a]
ops[j] = paulis[b]
return kron_n(ops)
# Projector onto ancilla qubit (index 4) state |0>
def P_anc0(n):
P0 = np.array([[1,0],[0,0]], dtype=complex)
return op_on(P0, 4, n)
# Stabilizer projectors on qubits 0..3
def projector_codespace(n):
# P = (I + XXXX)/2 * (I + ZZZZ)/2 on qubits 0..3
X4 = op_on(X,0,n) @ op_on(X,1,n) @ op_on(X,2,n) @ op_on(X,3,n)
Z4 = op_on(Z,0,n) @ op_on(Z,1,n) @ op_on(Z,2,n) @ op_on(Z,3,n)
P = 0.5*(np.eye(2**n) + X4) @ 0.5*(np.eye(2**n) + Z4)
return P
# Build the circuit U (without noise), on n=5 qubits, order: H1, CNOT12, CNOT10, CNOT23, CNOT34, CNOT04
n=5
U = np.eye(2**n, dtype=complex)
# Apply gates in chronological order (right-to-left order in math expression)
# 1) H_1
U = op_on(H, 1, n) @ U
# 2) CNOT12
U_cnot12 = CNOT_ij(1,2,n)
U = U_cnot12 @ U
# 3) CNOT10
U_cnot10 = CNOT_ij(1,0,n)
U = U_cnot10 @ U
# 4) CNOT23
U_cnot23 = CNOT_ij(2,3,n)
U = U_cnot23 @ U
# 5) CNOT34
U_cnot34 = CNOT_ij(3,4,n)
U = U_cnot34 @ U
# 6) CNOT04
U_cnot04 = CNOT_ij(0,4,n)
U = U_cnot04 @ U
# Define error locations: after each CNOT, the noise acts on the pair of those qubits
error_qubit_pairs = [(1,2),(1,0),(2,3),(3,4),(0,4)]
# Generate the set of 15 non-identity two-qubit Paulis as label pairs
non_id_two_paulis = []
for a in pauli_labels:
for b in pauli_labels:
if (a=='I' and b=='I'): continue
non_id_two_paulis.append((a,b))
# Verify length 15
print('num non-id 2q Paulis:', len(non_id_two_paulis))
# Initial state |00000>
zero = np.array([[1],[0]], dtype=complex)
psi0 = kron_n([zero]*n)
# Ideal reference state after circuit with projection onto ancilla |0> and codespace
rho_ideal = U @ psi0 @ (U @ psi0).conj().T
P0_anc = P_anc0(n)
P_code = projector_codespace(n)
P_total = P_code @ P0_anc
rho_ideal_acc = P_total @ rho_ideal @ P_total.conj().T
p_acc_ideal = np.trace(rho_ideal_acc).real
rho_ideal_acc = rho_ideal_acc / p_acc_ideal
# Reduce to code qubits (trace out ancilla qubit 4)
# We'll implement partial trace over qubit 4 for 5-qubit density matrix
def partial_trace_last_qubit(rho, n):
# trace out last qubit (index n-1)
dim = 2**(n-1)
rho = rho.reshape((2, dim, 2, dim))
# indices: (a, i, b, j), trace over a=b
red = rho[0,:,0,:] + rho[1,:,1,:]
return red
rho_ideal_code = partial_trace_last_qubit(rho_ideal_acc, n)
# Normalize: it should already be normalized (trace 1)
print('trace ideal code:', np.trace(rho_ideal_code).real)
# Check that rho_ideal_code is pure and stabilized by XXXX and ZZZZ
# Build 4-qubit operators embedded into 4-qubit space (we already traced out ancilla)
I4 = np.eye(2**4)
X4 = kron_n([X,X,X,X])
Z4 = kron_n([Z,Z,Z,Z])
print('Stab XXXX expectation:', np.trace(rho_ideal_code @ X4).real)
print('Stab ZZZZ expectation:', np.trace(rho_ideal_code @ Z4).real)
# We also want the ideal pure state vector |psi_ref> on code qubits by eigen-decomposition
val, vec = np.linalg.eigh(rho_ideal_code)
idx = np.argmax(val.real)
psi_ref = vec[:, idx]
# Normalize
psi_ref = psi_ref / np.linalg.norm(psi_ref)
# Function to compute accept probability and fidelity contribution for a given error pattern
# Precompute unitary for the fixed gates separated by error insertions? We can build step-by-step to insert errors interleaved.
# We'll implement the full process for each pattern: start with psi0, apply H1, then after each CNOT insert error operator (two-qubit Pauli on that pair), then continue
# Pre-build list of gate unitaries in sequence for convenience
gate_sequence = [
('single', 1, H),
('cnot', (1,2)),
('cnot', (1,0)),
('cnot', (2,3)),
('cnot', (3,4)),
('cnot', (0,4)),
]
cnot_map = { (1,2): U_cnot12, (1,0): U_cnot10, (2,3): U_cnot23, (3,4): U_cnot34, (0,4): U_cnot04 }
# Transform error pattern representation: a list of 5 entries; each either None for identity or a tuple (a,b) Pauli labels
from collections import defaultdict
A_k = defaultdict(Fraction) # sum of acceptance probabilities over patterns of weight k
G_k = defaultdict(Fraction) # sum of acceptance_prob * fidelity over patterns of weight k
# We'll iterate over all patterns. That's 16^5 ~ 1,048,576; heavy but doable.
# To speed up, precompute two-qubit Pauli operators on specific pairs for all 15 types at each error location
err_ops_by_loc = []
for (i,j) in error_qubit_pairs:
ops = []
for (a,b) in non_id_two_paulis:
ops.append(two_pauli(a,i,b,j,n))
err_ops_by_loc.append(ops)
# Identity operator for no error
Id = np.eye(2**n, dtype=complex)
# We'll represent patterns as indices 0..15 at each loc, where 0 means identity, 1..15 correspond to non-id list order
# Let's enumerate efficiently using integer base-16 representation
num_patterns = 16**5
print('Enumerating patterns:', num_patterns)
# We'll store G_k and A_k as Fractions by rounding small numerical errors to Fractions with limited denominator
def to_fraction(x, tol=1e-12):
# x is real close to rational with denominator 2^m; we convert using Fraction
return Fraction.from_float(float(x)).limit_denominator(2**20)
for pat_idx in range(num_patterns):
# Decode base-16 digits for 5 positions
digits = []
x = pat_idx
weight = 0
for _ in range(5):
d = x % 16
x //= 16
digits.append(d)
if d != 0:
weight += 1
# digits list corresponds to positions [0..4]
# Apply circuit with errors
psi = psi0
# Step 1: H1
psi = op_on(H, 1, n) @ psi
# After CNOT12
psi = U_cnot12 @ psi
if digits[0] != 0:
E = err_ops_by_loc[0][digits[0]-1]
psi = E @ psi
# After CNOT10
psi = U_cnot10 @ psi
if digits[1] != 0:
E = err_ops_by_loc[1][digits[1]-1]
psi = E @ psi
# After CNOT23
psi = U_cnot23 @ psi
if digits[2] != 0:
E = err_ops_by_loc[2][digits[2]-1]
psi = E @ psi
# After CNOT34
psi = U_cnot34 @ psi
if digits[3] != 0:
E = err_ops_by_loc[3][digits[3]-1]
psi = E @ psi
# After CNOT04
psi = U_cnot04 @ psi
if digits[4] != 0:
E = err_ops_by_loc[4][digits[4]-1]
psi = E @ psi
rho = psi @ psi.conj().T
rho_acc = P_total @ rho @ P_total.conj().T
p_acc = np.trace(rho_acc).real
if p_acc < 1e-15:
# no acceptance
continue
# Normalized accepted state
rho_acc_norm = rho_acc / p_acc
# Reduce to code qubits by tracing out ancilla
rho_code = partial_trace_last_qubit(rho_acc_norm, n)
# Fidelity with ideal state |psi_ref>
F = np.vdot(psi_ref, rho_code @ psi_ref).real
# Accumulate
A_k[weight] += to_fraction(p_acc)
G_k[weight] += to_fraction(p_acc * F)
# Convert dictionaries to lists for k=0..5
A_list = [A_k[k] for k in range(6)]
G_list = [G_k[k] for k in range(6)]
print('A_k (acceptance mass sums):', A_list)
print('G_k (good mass sums):', G_list)
# Also we can print total acceptance probability polynomial coefficients
# Den(p) = sum_k A_k * (1-p)^(5-k) * (p/15)^k
# Num(p) = sum_k G_k * (1-p)^(5-k) * (p/15)^k
# We'll expand as polynomials in p
from math import comb
# Compute coefficients of (1-p)^(5-k) * (p/15)^k = sum_{m=0}^{5-k} (-1)^m comb(5-k, m) p^m * (p/15)^k
# => sum_{m=0}^{5-k} (-1)^m comb(5-k,m) p^{m+k} / 15^k
max_deg = 5
coeff_den = [Fraction(0,1) for _ in range(max_deg+1)]
coeff_num = [Fraction(0,1) for _ in range(max_deg+1)]
for k in range(6):
A = A_list[k]
G = G_list[k]
for m in range(0, 5-k+1):
coeff = Fraction((-1)**m * comb(5-k, m), 15**k)
deg = m + k
coeff_den[deg] += A * coeff
coeff_num[deg] += G * coeff
print('Denominator polynomial coefficients (deg 0..5):', coeff_den)
print('Numerator polynomial coefficients (deg 0..5):', coeff_num)
# Let's also form floating representations and series expansion up to p^2
coeff_den_f = [float(c) for c in coeff_den]
coeff_num_f = [float(c) for c in coeff_num]
print('Den floats:', coeff_den_f)
print('Num floats:', coeff_num_f)
# Compute series of F(p) = Num/Den around p=0 up to p^2
import sympy as sp
p = sp.symbols('p')
Den_poly = sum(sp.Rational(c.numerator, c.denominator) * p**i for i,c in enumerate(coeff_den))
Num_poly = sum(sp.Rational(c.numerator, c.denominator) * p**i for i,c in enumerate(coeff_num))
F_series = sp.series(Num_poly/Den_poly, p, 0, 3) # up to p^2
print('F series up to p^2:', F_series)
# Also perhaps simplify the exact rational function
F_simplified = sp.simplify(Num_poly/Den_poly)
print('F rational simplified:', F_simplified)