QCQI — Chapter 12: Quantum Cryptography & Communication

# Backend selector
backend = "cirq"  # change: "pennylane", "braket", "pyquil", "qiskit"
print("Selected backend:", backend)
# Optional installs (if needed)
# !pip install numpy scipy cirq pennylane amazon-braket-sdk pyquil qiskit qiskit-aer

import numpy as np, math, random
rng = np.random.default_rng(1234)
def binent(p):
    return 0 if p==0 or p==1 else -(p*math.log2(p)+(1-p)*math.log2(1-p))

def bb84_run(N=1000, eavesdrop_frac=0.0, physical_error=0.0):
    # Alice bits & bases
    a = rng.integers(0,2,size=N)
    b = rng.integers(0,2,size=N)  # 0->Z, 1->X
    # Eve: with prob eavesdrop_frac, measure in random basis then resend
    e_do = rng.random(N) < eavesdrop_frac
    e_basis = rng.integers(0,2,size=N)
    # Bob bases
    bprime = rng.integers(0,2,size=N)
    # Outcomes
    x = np.zeros(N, dtype=int)
    for i in range(N):
        # ideal outcome if Bob's basis = Alice's basis (without Eve/physical errors)
        if e_do[i]:
            # Eve measurement
            # if e_basis matches Alice, she learns bit; else random
            if e_basis[i]==b[i]:
                ebit = a[i]
            else:
                ebit = rng.integers(0,2)
            # resend eigenstate of e_basis with ebit
            # Bob measures:
            if bprime[i]==e_basis[i]:
                x[i] = ebit
            else:
                x[i] = rng.integers(0,2)
        else:
            if bprime[i]==b[i]:
                x[i] = a[i]
            else:
                x[i] = rng.integers(0,2)
        # physical error flips
        if rng.random() < physical_error:
            x[i] ^= 1
    sift = (b==bprime)
    a_s = a[sift]; x_s = x[sift]
    if len(a_s)==0:
        return dict(sifted_bits=0, qber=0, R=0, a_s=a_s, x_s=x_s)
    qber = np.mean(a_s != x_s)
    # crude IR leakage: parity of blocks of size 32
    blocks = max(1, len(a_s)//32)
    ir_leak = blocks  # publish one parity per block (toy)
    key_len = len(a_s) - ir_leak
    # PA compression length estimate: key_len*(1 - 2*h(qber)) (toy asymptotic)
    R = max(0, int(key_len * max(0, 1 - 2*binent(qber))))
    return dict(sifted_bits=len(a_s), qber=qber, R=R, a_s=a_s, x_s=x_s)
for f in [0.0, 0.1, 0.25]:
    res = bb84_run(4000, eavesdrop_frac=f, physical_error=0.0)
    print(f"Eavesdrop {int(100*f)}% -> sifted {res['sifted_bits']}, QBER={res['qber']:.3f}, toy R≈{res['R']}")

def toy_ir(sender, receiver, block=32):
    # return corrected receiver string and number of parity bits leaked
    n = len(sender); r = receiver.copy()
    leak = 0
    for start in range(0, n, block):
        end = min(n, start+block)
        s_par = int(np.sum(sender[start:end]) % 2)
        r_par = int(np.sum(r[start:end]) % 2)
        leak += 1
        if s_par != r_par:
            # binary search error location
            L, R = start, end
            while R-L > 1:
                mid = (L+R)//2
                s_par = int(np.sum(sender[L:mid]) % 2)
                r_par = int(np.sum(r[L:mid]) % 2)
                leak += 1
                if s_par != r_par:
                    R = mid
                else:
                    L = mid
            r[L] ^= 1
    return r, leak
# demo
res = bb84_run(1024, eavesdrop_frac=0.05, physical_error=0.02)
corrected, leak = toy_ir(res['a_s'], res['x_s'])
print("IR: parity bits leaked =", leak, " residual errors =", int(np.sum(corrected != res['a_s'])))

def toeplitz_hash(bits, m):
    n = len(bits)
    # Build Toeplitz matrix via random first row+col
    tcol = rng.integers(0,2,size=n)
    trow = rng.integers(0,2,size=m)
    # Convolution-style multiply
    out = np.zeros(m, dtype=int)
    for i in range(m):
        # row i uses [tcol[i], tcol[i-1], ..., tcol[0], trow[1], ..., trow[i]] appropriately
        acc = 0
        for j in range(n):
            idx = i - j
            tij = tcol[j] if idx==0 else (trow[idx] if idx>0 else tcol[-idx] if -idx < n else 0)
            acc ^= (tij & bits[j])
        out[i] = acc
    return out
# demo compression
key = corrected
target_len = max(0, len(key) - leak - 64)
pa_key = toeplitz_hash(key, target_len) if target_len>0 else np.array([], dtype=int)
print("PA: compressed length =", len(pa_key))

if backend == "cirq":
    import cirq
    q = cirq.LineQubit.range(1)
    c = cirq.Circuit()
    # Prepare |+> or |0> at random and measure in random basis
    c.append(cirq.H(q[0]))
    c.append(cirq.measure(q[0], key='m'))
    print("Cirq: built minimal demo circuit.")
elif backend == "pennylane":
    import pennylane as qml
    dev = qml.device("default.qubit", wires=1, shots=10)
    @qml.qnode(dev)
    def measX():
        qml.Hadamard(0); return qml.sample(qml.PauliZ(0))
    print("PennyLane: sample X-basis by H+Z-measure:", measX())
elif backend == "braket":
    from braket.circuits import Circuit
    from braket.devices import LocalSimulator
    circ = Circuit().h(0).measure(0,0)
    _ = LocalSimulator().run(circ, shots=5).result()
    print("Braket: ran minimal H+measure.")
elif backend == "pyquil":
    from pyquil import Program
    from pyquil.gates import H, MEASURE
    from pyquil.api import get_qc
    p = Program(); ro = p.declare('ro','BIT',1); p += H(0); p += MEASURE(0,ro[0])
    _ = get_qc('1q-qvm').run(p.wrap_in_numshots_loop(5))
    print("PyQuil: ran minimal H+measure.")
elif backend == "qiskit":
    from qiskit import QuantumCircuit
    from qiskit_aer import AerSimulator
    qc = QuantumCircuit(1,1); qc.h(0); qc.measure(0,0)
    _ = AerSimulator().run(qc, shots=8).result()
    print("Qiskit: ran minimal H+measure.")

def b92_sweep(theta_deg=45, N=10000):
    # two pure states separated by angle theta on the Bloch sphere equator
    theta = math.radians(theta_deg)
    # Conclusive probability roughly ~ 1 - |<ψ0|ψ1>| = 1 - cos(theta/2)^2*2? (toy Monte Carlo below)
    # We'll simulate a simple USD-like POVM with an inconclusive outcome probability dependent on overlap.
    # Here we just use the Helstrom bound-inspired random decision (toy).
    overlap = abs(math.cos(theta/2))**2 + abs(math.sin(theta/2))**2 - 1  # bogus placeholder; use MC
    # Instead do MC by projecting onto orthogonal complements:
    # |ψ0>=|0>, |ψ1>=cosθ|0>+sinθ|1>; conclusive when outcome orthogonal to the other state's support.
    a = rng.integers(0,2,size=N)
    conclusive = 0
    for i in range(N):
        if a[i]==0:
            # measure {Π_perp(ψ1), I-Π_perp(ψ1)}; conclusive if Π_perp(ψ1)
            v1 = np.array([math.cos(theta), math.sin(theta)])
            perp = np.array([-v1[1], v1[0]]); perp = perp/np.linalg.norm(perp)
            prob = abs(perp[0])**2  # since ψ0=|0>
        else:
            prob = math.sin(theta)**2  # symmetric
        if rng.random() < prob: conclusive += 1
    return conclusive/N
for ang in [20,40,60]:
    print(f"B92: angle {ang}° -> conclusive rate ~ {b92_sweep(ang):.3f}")

# Sample CHSH with ideal correlations E(a,b) = -cos(theta_ab) using local-hidden-variable-violating settings
def chsh_value():
    import numpy as np, math
    # angles (radians) for A0,A1,B0,B1 on equator
    A0, A1 = 0, math.pi/2
    B0, B1 = math.pi/4, -math.pi/4
    def E(a,b): return -math.cos(a-b)
    S = E(A0,B0) + E(A0,B1) + E(A1,B0) - E(A1,B1)
    return S
print("Ideal singlet CHSH S =", round(chsh_value(),3), "(>2 violates Bell)")