import numpy as np
from functools import reduce
# Single-qubit Paulis
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)
Psingle = [I,X,Y,Z]
# Kronecker helper
kron = lambda *ops: reduce(np.kron, ops)
# Build circuit
def CNOT_ij(i,j):
P0 = np.array([[1,0],[0,0]], complex)
P1 = np.array([[0,0],[0,1]], complex)
ops0 = [I,I,I,I]
ops1 = [I,I,I,I]
ops0[i] = P0
ops1[i] = P1
ops1[j] = X
return kron(*ops0) + kron(*ops1)
def H_q(q):
ops = [I,I,I,I]
ops[q] = H
return kron(*ops)
# Circuit: U = CNOT_03 * H0 * CNOT_21 * H2
U = CNOT_ij(0,3) @ H_q(0) @ CNOT_ij(2,1) @ H_q(2)
zero = np.array([[1],[0]], complex)
psi0 = kron(zero,zero,zero,zero)
psi_ideal = U @ psi0
# Stabilizers S1=XXXX, S2=ZZZZ
S1 = kron(X,X,X,X)
S2 = kron(Z,Z,Z,Z)
# Logical operators (one representative choice matching problem statement)
XA = kron(X, I, X, I)
XB = kron(X, X, I, I)
ZA = kron(Z, Z, I, I)
ZB = kron(Z, I, Z, I)
# Logical stabilizers of |Phi+>_AB in code space
Lstab = [np.eye(16), XA@XB, ZA@ZB, (XA@XB)@(ZA@ZB)]
# Enumerate all two-qubit Pauli on pairs (0,3) and (2,1)
# Build E_total(P03, P21)
Ops = [I,X,Y,Z]
# Function to classify a 4-qubit Pauli E as accepted and as good (leaves psi_ideal up to phase)
def is_in_normalizer(E):
c1 = np.allclose(E@S1 - S1@E, 0)
c2 = np.allclose(E@S2 - S2@E, 0)
return c1 and c2
def leaves_state(E):
psi = E @ psi_ideal
ov = (psi_ideal.conj().T @ psi)[0,0]
return np.allclose(abs(ov), 1.0, atol=1e-10)
# Precompute acceptance and good masks for each pair of local paulis
accepted = {}
leaves = {}
for a03 in range(4):
for b03 in range(4):
E03 = kron(Ops[a03], I, I, Ops[b03])
for a21 in range(4):
for b21 in range(4):
E21 = kron(I, Ops[b21], Ops[a21], I)
E = E03 @ E21
key = (a03,b03,a21,b21)
accepted[key] = is_in_normalizer(E)
leaves[key] = leaves_state(E)
# Count categories for sanity
n_acc = sum(accepted.values())
n_leave = sum(v and accepted[k] for k,v in leaves.items())
print('Total combos:', 16*16, 'Accepted (normalizer):', n_acc, 'Good within accepted:', n_leave)
# As expected: normalizer has size 64; E_total restricted to product of two local 2-qubit pauli sets gives 256 combos; acceptance 64
# Build weights for each key given depolarizing probability p
def weights_of_p(p):
# weights for each local 2-qubit pauli on a pair
w = np.zeros((4,4))
for a in range(4):
for b in range(4):
if a==0 and b==0:
w[a,b] = 1-p
else:
w[a,b] = p/15
return w
# Compute acceptance probability and conditional fidelity
def acceptance_and_fidelity(p):
w03 = weights_of_p(p)
w21 = weights_of_p(p)
Pacc = 0.0
Pgood = 0.0
for a03 in range(4):
for b03 in range(4):
for a21 in range(4):
for b21 in range(4):
key = (a03,b03,a21,b21)
w = w03[a03,b03]*w21[a21,b21]
if accepted[key]:
Pacc += w
if leaves[key]:
Pgood += w
Fcond = Pgood/Pacc if Pacc>0 else 1.0
return Pacc, Fcond
for p in [0.0, 0.01, 0.1, 0.2, 0.5, 1.0]:
Pacc, Fcond = acceptance_and_fidelity(p)
print(p, 'Pacc=', Pacc, 'F_logical=', Fcond)