import numpy as np
from itertools import product
# Define single-qubit Pauli matrices
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)
pauli_dict = {'I':I, 'X':X, 'Y':Y, 'Z':Z}
paulis = ['I','X','Y','Z']
# Two-qubit CNOT gate (control c, target t) acting on two-qubit space (control first)
CNOT = np.array([[1,0,0,0],
[0,1,0,0],
[0,0,0,1],
[0,0,1,0]], dtype=complex)
H = (1/np.sqrt(2)) * np.array([[1,1],[1,-1]], dtype=complex)
# Build mapping tables by explicit conjugation for 1-qubit H and 2-qubit CNOT
def single_qubit_conj_map(U):
# returns dict mapping 'I','X','Y','Z' -> (phase, pauli)
M = {}
for p in paulis:
Pin = pauli_dict[p]
Pout = U @ Pin @ U.conj().T
# Find which Pauli (up to phase) equals Pout
found = False
for q in paulis:
Q = pauli_dict[q]
# Try factors 1, -1, 1j, -1j
for c in [1, -1, 1j, -1j]:
if np.allclose(Pout, c*Q):
M[p] = (c, q)
found = True
break
if found:
break
if not found:
raise RuntimeError('No match for H conjugation of', p)
return M
H_map = single_qubit_conj_map(H)
# Two-qubit conjugation map for CNOT (control first, target second)
def two_qubit_conj_map(U):
M = {}
for pc in paulis:
for pt in paulis:
Pin = np.kron(pauli_dict[pc], pauli_dict[pt])
Pout = U @ Pin @ U.conj().T
found = False
for qc in paulis:
for qt in paulis:
Q = np.kron(pauli_dict[qc], pauli_dict[qt])
for c in [1, -1, 1j, -1j]:
if np.allclose(Pout, c*Q):
M[(pc,pt)] = (c, qc, qt)
found = True
break
if found:
break
if found:
break
if not found:
raise RuntimeError('No match for CNOT conjugation of', pc, pt)
return M
CNOT_map = two_qubit_conj_map(CNOT)
# Utility: multiply two single-qubit Paulis (with phases)
# returns (phase, resulting_pauli)
mult_table = {}
for a in paulis:
for b in paulis:
A = pauli_dict[a]
B = pauli_dict[b]
P = A @ B
found=False
for c in [1, -1, 1j, -1j]:
for r in paulis:
if np.allclose(P, c*pauli_dict[r]):
mult_table[(a,b)] = (c, r)
found=True
break
if found: break
if not found:
raise RuntimeError('No match for mult', a, b)
# Multiply multi-qubit Pauli strings with phase tracking
# strings are like ['X','I','Z','Y','I'] for 5 qubits
def multiply_strings(s1, s2):
c = 1+0j
out = []
for a,b in zip(s1,s2):
cb, r = mult_table[(a,b)]
c *= cb
out.append(r)
return c, out
# Build projectors as Pauli sums
# P_code on data qubits 0..3: (I + XXXX)(I + ZZZZ)/4 = (I + XXXX + ZZZZ - YYYY)/4
def pauli_string(qubits5, letters, coeff=1+0j):
# letters length equals number of active qubits; we will place on 5 qubits with last being ancilla optionally
return (coeff, letters)
# Generate P_accept terms on 5 qubits: (1/8)*(I+Z4) tensor (I + XXXX + ZZZZ - YYYY)
Pacc_terms = []
for anc in ['I','Z']:
for data in [ ['I','I','I','I'], ['X','X','X','X'], ['Z','Z','Z','Z'], ['Y','Y','Y','Y'] ]:
coeff = 1/8
if data == ['Y','Y','Y','Y']:
coeff *= -1 # because of minus sign
letters5 = data + [anc]
Pacc_terms.append((coeff+0j, letters5))
# Build target projector on data qubits: stabilized by XXXX, ZZZZ, Z_A=ZZII, Z_B=ZIZI with +1 eigenvalues
# P_target = (1/16) sum over group generated by these four independent commuting Paulis
Ggens = [
['X','X','X','X'],
['Z','Z','Z','Z'],
['Z','Z','I','I'], # Z_A
['Z','I','Z','I'], # Z_B
]
# Enumerate group elements
group_elems = {}
for a,b,c,d in product([0,1],[0,1],[0,1],[0,1]):
coeff = 1+0j
s = ['I','I','I','I']
for e, g in zip([a,b,c,d], Ggens):
if e==1:
c2, s = multiply_strings(s, g)
coeff *= c2
key = tuple(s)
# Sum coefficients if same Pauli appears (shouldn't happen often)
group_elems[key] = group_elems.get(key, 0+0j) + coeff
# Convert to list of (coeff, letters)
Ptarget_data_terms = []
for letters, coeff in group_elems.items():
# Each term has coefficient 1/16 times its coeff from multiplication
Ptarget_data_terms.append(((coeff/16)+0j, list(letters)))
# Build P_good terms on 5 qubits: (1/2)*(I+Z4) tensor P_target_data
Pgood_terms = []
for anc in ['I','Z']:
for coeff, letters4 in Ptarget_data_terms:
Pgood_terms.append(((coeff/2)+0j, letters4 + [anc]))
# Define Heisenberg propagation through the full circuit
# Steps: [H1], [CNOT12 + E1], [CNOT10 + E2], [CNOT23 + E3], [CNOT34 + E4], [CNOT04 + E5]
# We implement backward propagation: for k=last..first: apply E_k^† then conjugate by gate^
# Conjugation helpers
def conj_H_on_qubit(term_letters, q):
letters = term_letters.copy()
c = 1+0j
p = letters[q]
c1, q1 = H_map[p]
letters[q] = q1
c *= c1
return c, letters
def conj_CNOT_on_pair(term_letters, c_idx, t_idx):
letters = term_letters.copy()
pc = letters[c_idx]
pt = letters[t_idx]
c1, qc, qt = CNOT_map[(pc,pt)]
letters[c_idx] = qc
letters[t_idx] = qt
return c1, letters
# For depolarizing adjoint factor: if either of the two qubits at indices (a,b) is non-identity
def E_dag_factor(term_letters, a, b, lam):
return lam if (term_letters[a] != 'I' or term_letters[b] != 'I') else 1.0
# Propagate a term and return its contribution factor to expectation on |00000> (real scalar) as a function of lam
def propagate_term(term_coeff, term_letters, lam):
c = term_coeff
letters = term_letters.copy()
# k=6: E5 on (0,4) then CNOT04
c *= E_dag_factor(letters, 0, 4, lam)
c1, letters = conj_CNOT_on_pair(letters, 0, 4)
c *= c1
# k=5: E4 on (3,4) then CNOT34
c *= E_dag_factor(letters, 3, 4, lam)
c1, letters = conj_CNOT_on_pair(letters, 3, 4)
c *= c1
# k=4: E3 on (2,3) then CNOT23
c *= E_dag_factor(letters, 2, 3, lam)
c1, letters = conj_CNOT_on_pair(letters, 2, 3)
c *= c1
# k=3: E2 on (1,0) then CNOT10
c *= E_dag_factor(letters, 1, 0, lam)
c1, letters = conj_CNOT_on_pair(letters, 1, 0)
c *= c1
# k=2: E1 on (1,2) then CNOT12
c *= E_dag_factor(letters, 1, 2, lam)
c1, letters = conj_CNOT_on_pair(letters, 1, 2)
c *= c1
# k=1: H1
c1, letters = conj_H_on_qubit(letters, 1)
c *= c1
# Now compute expectation on |00000>: nonzero only if letters are only I or Z
for L in letters:
if L in ('X','Y'):
return 0.0
# For I/Z-only string, expectation equals real part of global coefficient c (since <0|Z|0>=1)
return float(np.real_if_close(c).real)
# Compute A(lam) and G(lam)
def compute_A_G(lam):
A = 0.0
for coeff, letters in Pacc_terms:
A += propagate_term(coeff, letters, lam)
G = 0.0
for coeff, letters in Pgood_terms:
G += propagate_term(coeff, letters, lam)
return A, G
# Sanity checks: p=0 -> lam=1; both A and G should be 1
for lam in [1.0, 0.0, 0.5]:
A, G = compute_A_G(lam)
print('lam=', lam, 'A=', A, 'G=', G, 'F=G/A' if A>0 else 'A=0')
# Now, let's try to express A and G as polynomials in lam
# We can fit polynomials up to degree 5 since there are 5 noise channels
L_vals = np.linspace(0,1,6)
A_coeffs = np.polyfit(L_vals, [compute_A_G(l)[0] for l in L_vals], 5)
G_coeffs = np.polyfit(L_vals, [compute_A_G(l)[1] for l in L_vals], 5)
print('A(lam) poly coeffs (deg 5 -> deg 0):', A_coeffs)
print('G(lam) poly coeffs (deg 5 -> deg 0):', G_coeffs)
# For readability, round coefficients to rationals with small denominators
from fractions import Fraction
def rationalize_coeffs(coeffs):
rats = []
for c in coeffs:
f = Fraction(float(np.real_if_close(c).real)).limit_denominator(10**6)
rats.append(f)
return rats
print('A rational coeffs:', rationalize_coeffs(A_coeffs))
print('G rational coeffs:', rationalize_coeffs(G_coeffs))
# Convert lam to p: lam = 1 - 16p/15
# We can express A(p) and G(p) by substituting lam and expanding symbolic polynomial using numpy.poly1d
A_poly_lam = np.poly1d(A_coeffs)
G_poly_lam = np.poly1d(G_coeffs)
import sympy as sp
p = sp.symbols('p', real=True)
lam = 1 - sp.Rational(16,15)*p
# Convert numpy poly1d to sympy poly
A_sym_lam = sum(sp.Rational(Fraction(float(c)).limit_denominator(10**6)) * lam**i for i,c in enumerate(A_poly_lam.c[::-1]))
G_sym_lam = sum(sp.Rational(Fraction(float(c)).limit_denominator(10**6)) * lam**i for i,c in enumerate(G_poly_lam.c[::-1]))
A_sym = sp.simplify(sp.expand(A_sym_lam))
G_sym = sp.simplify(sp.expand(G_sym_lam))
F_sym = sp.simplify(sp.factor(G_sym / A_sym))
print('A_sym(p)=', A_sym)
print('G_sym(p)=', G_sym)
print('F_sym(p)=', F_sym)
# Also expand series up to p^2
print('Series F up to p^3:', sp.series(F_sym, p, 0, 4))