Quantum Supremacy Explained — Milestones, Experiments and Debate

# supremacy_criteria.py
import numpy as np
import math

class SupremacyAnalysis:
    """Analyze the criteria for quantum supremacy experiments."""

    @staticmethod
    def classical_simulation_cost(n_qubits, depth):
        """
        Estimate the classical cost of simulating a quantum circuit.
        Uses tensor network contraction cost estimation.
        """
        # Statevector simulation: O(2^n * depth)
        sv_cost = (2 ** n_qubits) * depth

        # Tensor network simulation (approximate)
        # For random circuits, the cost scales as ~2^(c*n) for some c < 1
        tn_cost = 2 ** (0.5 * n_qubits)  # Simplified

        return sv_cost, tn_cost

    @staticmethod
    def supremacy_frontier(available_flops=1e17, time_limit=3600):
        """
        Find the frontier where quantum circuits become classically hard.
        available_flops: Summit supercomputer ~200 PFLOPS = 2e17 FLOPS
        time_limit: maximum acceptable runtime in seconds (1 hour)
        """
        max_classical_ops = available_flops * time_limit

        print("=== Supremacy Frontier ===")
        print(f"Classical budget: {max_classical_ops:.2e} operations")

        for n in [20, 30, 40, 50, 53, 60, 70]:
            for d in [10, 20, 40]:
                sv_cost, tn_cost = SupremacyAnalysis.classical_simulation_cost(n, d)

                feasible_sv = sv_cost < max_classical_ops
                feasible_tn = tn_cost < max_classical_ops

                if not feasible_sv:
                    print(f"n={n:2d}, d={d:2d}: "
                          f"SV={sv_cost:.2e} {'[INFEASIBLE]':>12} "
                          f"TN={tn_cost:.2e} {'[OK]' if feasible_tn else '[INFEASIBLE]':>12}")

# Analysis
SupremacyAnalysis.supremacy_frontier()

# Key milestone timeline
print(f"\n=== Supremacy Milestones ===")
milestones = [
    (2010, "1-2 qubits", "Basic gate demonstrations"),
    (2015, "5-10 qubits", "Small algorithms (Deutsch-Jozsa, Grover)"),
    (2017, "20 qubits", "Quantum chemistry simulations"),
    (2019, "53 qubits", "Google Sycamore supremacy experiment"),
    (2020, "60 qubits", "Chinese Zuchongzhi 2.1 (2-3 orders faster)"),
    (2023, "100+ qubits", "Quantum utility demonstrations"),
    (2026, "1000+ qubits", "Error correction below threshold at scale"),
]

for year, qubits, event in milestones:
    print(f"  {year}: {qubits:>10} — {event}")

# sycamore_simulation.py
import numpy as np

class SycamoreExperiment:
    """Simulate and analyze the Google Sycamore supremacy experiment."""

    def __init__(self, n_qubits=53, n_cycles=20):
        self.n = n_qubits
        self.cycles = n_cycles

    def random_circuit(self):
        """
        Generate a random circuit similar to Sycamore's supremacy circuits.
        Each cycle consists of single-qubit gates + two-qubit gates on a grid.
        """
        # Sycamore used a 2D grid with specific connectivity
        # Each cycle: random single-qubit gates on all qubits,
        # followed by CZ gates on a subset of edges

        circuit = []
        gates_pool = [
            lambda: ('X^1/2', 0),
            lambda: ('Y^1/2', 0),
            lambda: ('W', 0),
        ]

        for cycle in range(self.cycles):
            # Single qubit gates (random choice)
            layer1 = []
            for q in range(min(10, self.n)):  # Show first 10
                gate = np.random.choice(gates_pool)()
                layer1.append((q, gate))
            circuit.append(layer1)

            # Two qubit gates (CZ on grid edges)
            layer2 = []
            for q in range(0, min(10, self.n) - 1, 2):
                layer2.append((q, q + 1, 'CZ'))
            circuit.append(layer2)

        return circuit

    def cross_entropy_fidelity(self, ideal_probs, measured_probs):
        """
        Compute the cross-entropy fidelity benchmark.
        F = sum(p_ideal * log(p_measured)) for heavy outputs
        """
        # Sort by ideal probability (descending)
        sorted_indices = np.argsort(ideal_probs)[::-1]

        # Heavy outputs: bitstrings with probability above median
        median = np.median(ideal_probs)
        heavy = [i for i in range(len(ideal_probs)) if ideal_probs[i] > median]

        # XEB fidelity
        F = 0
        for h in heavy:
            if ideal_probs[h] > 0:
                F += measured_probs[h] / ideal_probs[h]
        F = F / len(heavy) - 1

        return F

    def analyze(self):
        """Analyze the supremacy experiment parameters."""
        # State space size
        state_space = 2 ** self.n
        print(f"=== Sycamore Experiment Analysis ===")
        print(f"Qubits: {self.n}")
        print(f"Cycles: {self.cycles}")
        print(f"State space: 2^{self.n} = {state_space:.2e}")

        # Circuit statistics
        n_gates = self.n * self.cycles * 2  # Approximate
        circuit = self.random_circuit()
        print(f"Approximate gates: {n_gates:,}")
        print(f"Circuit depth: {self.cycles * 2}")

        # Classical simulation difficulty
        sv_memory = (2 ** self.n) * 16  # bytes (complex128)
        print(f"\nStatevector memory: {sv_memory:.2e} bytes "
              f"({sv_memory/1e12:.2f} TB)")

        # Sycamore's claim: 200 seconds vs 10,000 years
        speedup = 10000 * 365.25 * 24 * 3600 / 200
        print(f"\nClaimed speedup: ~{speedup:.2e}x")

        # Cross-entropy benchmarking
        print(f"\nCross-Entropy Benchmarking:")
        print(f"  Sycamore achieved F_XEB ≈ 0.002 (0.2%)")
        print(f"  Verified using 1M samples (bitstrings)")
        print(f"  Classical verification cost: ~days on supercomputer")

# Analysis
sycamore = SycamoreExperiment(53, 20)
sycamore.analyze()

# xeb_verification.py
import numpy as np

class CrossEntropyBenchmarking:
    """Implement cross-entropy benchmarking for quantum circuit verification."""

    def __init__(self, n_qubits=12):  # Small for simulation
        self.n = n_qubits
        self.dim = 2 ** n_qubits

    def ideal_simulation(self, depth=10):
        """Simulate the ideal circuit output probabilities."""
        # Random unitary (simplified: random orthogonal matrix)
        np.random.seed(42)
        random_matrix = np.random.randn(self.dim, self.dim) + 1j * np.random.randn(self.dim, self.dim)
        Q, R = np.linalg.qr(random_matrix)
        U = Q @ np.diag(np.diag(R) / np.abs(np.diag(R)))

        # Start in |0...0⟩
        state = np.zeros(self.dim, dtype=complex)
        state[0] = 1.0

        # Apply circuit (simplified)
        for d in range(depth):
            state = U @ state

        # Compute probabilities
        probs = np.abs(state) ** 2
        return probs

    def noisy_simulation(self, ideal_probs, noise_level=0.01):
        """Simulate noisy measurement."""
        # Add depolarizing noise to output probabilities
        noisy = (1 - noise_level) * ideal_probs + noise_level / self.dim
        return noisy

    def compute_xeb(self, ideal_probs, measured_probs, n_samples=10000):
        """Compute cross-entropy fidelity."""
        # Draw samples from measured distribution
        rng = np.random.default_rng(42)
        samples = rng.choice(self.dim, size=n_samples, p=measured_probs)

        # Compute linear XEB
        sum_log_p = 0
        for s in samples:
            if ideal_probs[s] > 0:
                sum_log_p += np.log(ideal_probs[s])

        # Average log probability
        avg_log_p = sum_log_p / n_samples

        # Expected value for uniform distribution
        uniform_log_p = np.log(1 / self.dim)

        # XEB fidelity
        # F_XEB = (avg_log_p - uniform_log_p) / (max_log_p - uniform_log_p)
        # where max_log_p = 0 (probabilities ≤ 1)

        # Simpler: compute linear XEB
        # F = D * sum(p_ideal * p_measured) - 1
        F = self.dim * np.sum(ideal_probs * measured_probs) - 1

        return F

    def run_benchmark(self):
        print("=== Cross-Entropy Benchmarking ===")
        print(f"Qubits: {self.n}, State space: {self.dim}")

        # Ideal run
        ideal = self.ideal_simulation(10)

        # Noisy runs with different noise levels
        for noise in [0.0, 0.01, 0.05, 0.1, 0.2]:
            noisy = self.noisy_simulation(ideal, noise)
            F = self.compute_xeb(ideal, noisy, 50000)

            # Determine if supremacy is demonstrated
            if F > 0:
                bits_per_sample = -np.log2(noise) if noise > 0 else self.n
                print(f"  Noise={noise:.0%}: F_XEB={F:.6f} "
                      f"{'[Supremacy threshold]' if F > 0.01 else ''}")

        # Statistical significance
        print(f"\n=== Verification Requirements ===")
        n_qubits_range = [12, 20, 30, 40, 50, 53]
        target_F = 0.002  # Sycamore's achieved XEB

        for n in n_qubits_range:
            dim = 2 ** n
            # To verify F_XEB with 3 sigma confidence
            n_samples = int(3**2 / (target_F**2 * dim))
            n_samples = min(max(n_samples, 1000), 10_000_000)
            print(f"  n={n:2d}: {n_samples:>8,} samples needed for F={target_F}")

xeb = CrossEntropyBenchmarking()
xeb.run_benchmark()

# classical_simulation.py
import numpy as np

class ClassicalSimulationChallenge:
    """Analyze why classical simulation of quantum circuits is hard."""

    def __init__(self):
        pass

    def tensor_network_cost(self, n_qubits, depth, bond_dim=2):
        """Estimate tensor network contraction cost."""
        # For a 1D chain, contraction cost is O(n * d * bond^3)
        # For a 2D grid, cost is O(2^(sqrt(n)))
        cost_1d = n_qubits * depth * bond_dim ** 3
        cost_2d = 2 ** (0.5 * np.sqrt(n_qubits))
        cost_statevector = 2 ** n_qubits

        return cost_1d, cost_2d, cost_statevector

    def compare_methods(self):
        print("=== Classical Simulation Methods Comparison ===")
        methods = [
            ("Statevector", lambda n, d: 2**n, "Exponential memory"),
            ("Feynman paths", lambda n, d: 4**(n*d), "Doubly exponential"),
            ("Tensor network (1D)", lambda n, d: n*d*8, "Efficient for 1D"),
            ("Tensor network (2D)", lambda n, d: 2**(0.5*np.sqrt(n)), "Hard for 2D"),
            ("Clifford (stabilizer)", lambda n, d: n**2, "Only Clifford gates"),
            ("MPS (low entanglement)", lambda n, d: n*d*100, "Fails for deep circuits"),
        ]

        for name, cost_fn, desc in methods:
            print(f"\n{name}:")
            for n, d in [(10, 10), (20, 10), (30, 20), (53, 20), (70, 20)]:
                cost = cost_fn(n, d)
                if cost < 1e3:
                    print(f"  n={n:2d}, d={d:2d}: {cost:.0f} ops [EASY]")
                elif cost < 1e12:
                    print(f"  n={n:2d}, d={d:2d}: {cost:.2e} ops [POSSIBLE]")
                elif cost < 1e18:
                    print(f"  n={n:2d}, d={d:2d}: {cost:.2e} ops [HARD]")
                else:
                    print(f"  n={n:2d}, d={d:2d}: {cost:.2e} ops [INFEASIBLE]")

    def chinese_supercomputer_challenge(self):
        """
        In 2020, Chinese researchers claimed classical simulation
        of Sycamore's circuit using a supercomputer.
        """
        print(f"\n=== The Simulation Debate ===")
        print(f"2020: Chinese team simulates Sycamore circuit")
        print(f"  Method: Tensor network on Sunway TaihuLight")
        print(f"  Time: ~20 days (claimed)")
        print(f"  2021: Improved to ~304 seconds on 512 GPUs")
        print(f"  Counter-claim: Still uses approximations")
        print(f"  Google response: Full verification still impossible")
        print(f"")
        print(f"Key takeaway: Supremacy is not a fixed boundary.")
        print(f"As classical algorithms improve, the frontier shifts.")
        print(f"Quantum hardware must stay ahead of classical simulation.")

sim = ClassicalSimulationChallenge()
sim.compare_methods()
sim.chinese_supercomputer_challenge()

# zuchongzhi_simulation.py
import numpy as np

print("=== Zuchongzhi 2.1 (2021) ===")

# Specifications
specs = {
    "Processor": "Zuchongzhi 2.1",
    "Qubits": 66,
    "Used for experiment": 56,
    "Gate type": "Superconducting (transmon)",
    "Connectivity": "2D grid",
    "Chip area": "~10 mm x 10 mm",
}

for k, v in specs.items():
    print(f"  {k}: {v}")

# Performance comparison
print(f"\n=== Performance Comparison ===")
experiments = [
    ("Google Sycamore (2019)", 53, 200, 10000*365.25*24*3600),
    ("Zuchongzhi 2.0 (2021)", 56, 72, 8*24*3600),
    ("Zuchongzhi 2.1 (2021)", 60, 84, 37.5*3600),
    ("IBM Quantum Utility (2023)", 127, 5*3600, 15*3600),
]

print(f"{'Experiment':<30} {'Qubits':>8} {'Time(s)':>10} {'Classical(s)':>15} {'Speedup':>10}")
print("-" * 75)

for name, qubits, q_time, c_time in experiments:
    speedup = c_time / q_time if q_time > 0 else float('inf')
    c_str = f"{c_time:.1e}" if c_time > 1e6 else f"{c_time:.0f}"
    print(f"{name:<30} {qubits:>8} {q_time:>10.1f} {c_str:>15} {speedup:>10.2e}")

# Scaling analysis
print(f"\n=== Scaling Analysis ===")
n_qubits = np.array([53, 56, 60, 67, 70, 80])
classical_time = lambda n: 2 ** (n - 53) * 10000 * 365.25 * 24 * 3600

for n in [53, 56, 60, 70, 80]:
    ct = classical_time(n)
    days = ct / (24 * 3600)
    years = ct / (365.25 * 24 * 3600)
    if days < 1:
        print(f"  n={n}: {ct:.1e} seconds")
    elif years < 1:
        print(f"  n={n}: {days:.1f} days")
    else:
        print(f"  n={n}: {years:.1e} years")

def is_supremacy_circuit(n_qubits, depth, gate_set, topology):
    """
    Determine if a circuit is likely classically hard.
    
    Returns: (bool, dict) where bool indicates supremacy potential
             and dict contains analysis metrics.
    """
    pass