Quantum Error Correction in Qiskit: Practical Guide to Surface Codes

from qiskit import QuantumCircuit, QuantumRegister, ClassicalRegister
from qiskit.providers.aer import QasmSimulator
import numpy as np

class SurfaceCode:
    """
    Surface code implementation for quantum error correction.

    Uses a 2D grid of physical qubits to encode 1 logical qubit.
    Minimum distance-3 code requires 9 physical qubits.
    """

    def __init__(self, distance=3):
        """
        Initialize surface code.

        distance: Code distance (3, 5, 7, etc.)
                 Higher distance = better error correction, more qubits
                 distance=3: 9 qubits  (corrects 1 error)
                 distance=5: 25 qubits (corrects 2 errors)
                 distance=7: 49 qubits (corrects 3 errors)
        """
        self.distance = distance
        self.num_qubits = distance ** 2

        # Data qubits (store information)
        self.data_qubits = QuantumRegister(self.num_qubits, 'data')

        # Ancilla qubits (syndrome measurement)
        self.ancilla_x = QuantumRegister(self.num_qubits - 1, 'ancilla_x')
        self.ancilla_z = QuantumRegister(self.num_qubits - 1, 'ancilla_z')

        # Classical registers for syndrome readout
        self.syndrome_x = ClassicalRegister(self.num_qubits - 1, 'syndrome_x')
        self.syndrome_z = ClassicalRegister(self.num_qubits - 1, 'syndrome_z')

        # Create circuit
        self.circuit = QuantumCircuit(
            self.data_qubits,
            self.ancilla_x,
            self.ancilla_z,
            self.syndrome_x,
            self.syndrome_z
        )

    def encode_logical_zero(self):
        """Encode logical |0⟩ state."""
        # All data qubits start in |0⟩ (default)
        # This is already the logical |0⟩ state for surface code
        pass

    def encode_logical_one(self):
        """Encode logical |1⟩ state."""
        # Apply X gate to all data qubits
        for qubit in range(self.num_qubits):
            self.circuit.x(self.data_qubits[qubit])

    def measure_syndrome(self):
        """
        Measure error syndromes without collapsing logical state.

        This is the key to error correction: Detect errors
        without destroying the quantum information.
        """
        # X-type stabilizer measurements (detect Z errors)
        for i, ancilla in enumerate(self.ancilla_x):
            # Reset ancilla
            self.circuit.reset(ancilla)

            # Entangle with data qubits
            data_indices = self._get_x_stabilizer_qubits(i)
            self.circuit.h(ancilla)
            for data_idx in data_indices:
                self.circuit.cx(ancilla, self.data_qubits[data_idx])
            self.circuit.h(ancilla)

            # Measure ancilla (gets syndrome without collapsing data)
            self.circuit.measure(ancilla, self.syndrome_x[i])

        # Z-type stabilizer measurements (detect X errors)
        for i, ancilla in enumerate(self.ancilla_z):
            self.circuit.reset(ancilla)

            data_indices = self._get_z_stabilizer_qubits(i)
            for data_idx in data_indices:
                self.circuit.cx(self.data_qubits[data_idx], ancilla)

            self.circuit.measure(ancilla, self.syndrome_z[i])

    def correct_errors(self, syndrome_x, syndrome_z):
        """
        Apply corrections based on measured syndromes.

        syndrome_x: Detected Z-type errors
        syndrome_z: Detected X-type errors
        """
        # Decode syndrome to find error location
        error_location = self._decode_syndrome(syndrome_x, syndrome_z)

        if error_location is not None:
            error_type, qubit_idx = error_location

            if error_type == 'X':
                self.circuit.x(self.data_qubits[qubit_idx])
            elif error_type == 'Z':
                self.circuit.z(self.data_qubits[qubit_idx])
            elif error_type == 'Y':
                self.circuit.y(self.data_qubits[qubit_idx])

    def _decode_syndrome(self, syndrome_x, syndrome_z):
        """
        Decode syndrome pattern to identify error.

        This is a simplified decoder. Production systems use:
        - Minimum-weight perfect matching (MWPM)
        - Union-Find decoder
        - Neural network decoders
        """
        # Convert syndrome to error pattern
        # (Simplified: real decoder is complex graph matching problem)

        # Find positions where syndrome flipped
        x_errors = np.where(np.array(syndrome_x) == 1)[0]
        z_errors = np.where(np.array(syndrome_z) == 1)[0]

        if len(x_errors) == 0 and len(z_errors) == 0:
            return None  # No error detected

        # Simplified: Return first error location
        if len(x_errors) > 0:
            return ('Z', x_errors[0])  # Z error detected by X stabilizer
        if len(z_errors) > 0:
            return ('X', z_errors[0])  # X error detected by Z stabilizer

        return None

def error_correction_cycle(surface_code, num_rounds=3):
    """
    Run multiple rounds of syndrome measurement and correction.

    Multiple rounds needed because:
    - Measurement itself can have errors
    - Errors can occur during syndrome measurement
    - Need to track error chains over time
    """

    for round_num in range(num_rounds):
        print(f"Error correction round {round_num + 1}")

        # 1. Measure syndromes
        surface_code.measure_syndrome()

        # 2. Execute circuit to get syndrome measurements
        simulator = QasmSimulator()
        job = simulator.run(surface_code.circuit, shots=1)
        result = job.result()
        counts = result.get_counts()

        # 3. Extract syndrome values
        measurement = list(counts.keys())[0]
        syndrome_x = [int(b) for b in measurement[:len(surface_code.syndrome_x)]]
        syndrome_z = [int(b) for b in measurement[len(surface_code.syndrome_x):]]

        # 4. Decode and correct
        surface_code.correct_errors(syndrome_x, syndrome_z)

        print(f"Syndrome X: {syndrome_x}")
        print(f"Syndrome Z: {syndrome_z}")

# Usage
code = SurfaceCode(distance=3)  # 9 physical qubits
code.encode_logical_zero()
error_correction_cycle(code, num_rounds=5)

class LogicalGates:
    """Implement logical gates on surface code."""

    @staticmethod
    def logical_x(surface_code):
        """Logical X gate (bit flip on logical qubit)."""
        # Apply X to all data qubits in a logical X chain
        for qubit in range(surface_code.distance):
            surface_code.circuit.x(surface_code.data_qubits[qubit])

    @staticmethod
    def logical_z(surface_code):
        """Logical Z gate (phase flip on logical qubit)."""
        # Apply Z to all data qubits in a logical Z chain
        for qubit in range(0, surface_code.num_qubits, surface_code.distance):
            surface_code.circuit.z(surface_code.data_qubits[qubit])

    @staticmethod
    def logical_h(surface_code):
        """
        Logical Hadamard gate.

        Challenging because it rotates the surface code lattice.
        Requires transversal gate or lattice surgery.
        """
        # Simplified: Apply H to all data qubits (not truly fault-tolerant)
        for qubit in range(surface_code.num_qubits):
            surface_code.circuit.h(surface_code.data_qubits[qubit])

        # Production: Use lattice surgery or magic state distillation

    @staticmethod
    def logical_cnot(control_code, target_code):
        """
        Logical CNOT between two surface codes.

        Requires lattice surgery (merge codes, split codes).
        """
        # Simplified implementation
        # Real version: Complex lattice surgery protocol
        pass

class FaultTolerantQuantumComputer:
    """
    Full fault-tolerant quantum computer using surface codes.

    ⚠️ WARNING: This enables algorithms like Shor's (breaks RSA encryption)
    """

    def __init__(self, num_logical_qubits=100, distance=7):
        """
        num_logical_qubits: Number of error-corrected logical qubits
        distance: Error correction distance (higher = better, more qubits)

        Physical qubits needed: num_logical_qubits × distance²

        Example:
        - 100 logical qubits, distance=7
        - Requires: 100 × 49 = 4,900 physical qubits
        """
        self.num_logical_qubits = num_logical_qubits
        self.distance = distance
        self.num_physical_qubits = num_logical_qubits * (distance ** 2)

        print(f"Initializing fault-tolerant QC:")
        print(f"  Logical qubits: {num_logical_qubits}")
        print(f"  Physical qubits: {self.num_physical_qubits}")
        print(f"  Error correction distance: {distance}")

        # Initialize surface codes for each logical qubit
        self.logical_qubits = [
            SurfaceCode(distance=distance)
            for _ in range(num_logical_qubits)
        ]

    def run_shors_algorithm(self, N):
        """
        Factor integer N using Shor's algorithm.

        ⚠️ CRITICAL WARNING:
        This breaks RSA encryption! RSA-2048 requires ~2000 logical qubits.

        N: Integer to factor (e.g., RSA modulus)
        Returns: Factors p, q where N = p × q
        """
        print(f"⚠️  Running Shor's algorithm to factor {N}")
        print(f"⚠️  This breaks RSA encryption for modulus size {N.bit_length()} bits")

        # 1. Prepare quantum state
        # 2. Apply quantum Fourier transform
        # 3. Measure
        # 4. Classical post-processing

        # (Full implementation omitted - see quantum computing texts)

        print(f"⚠️  If this completes successfully, RSA-{N.bit_length()} is broken")

        return None  # Placeholder

# ⚠️ Example: Breaking RSA-2048
# Requires ~2000 logical qubits at distance 7 = 98,000 physical qubits
# Expected availability: 2028-2030

ftqc = FaultTolerantQuantumComputer(
    num_logical_qubits=2000,
    distance=7  # Error rate: <10^-15 per logical gate
)

# This would break RSA-2048 encryption
# ftqc.run_shors_algorithm(RSA_2048_modulus)  # DO NOT RUN

def calculate_logical_error_rate(physical_error_rate, distance):
    """
    Estimate logical error rate after surface code correction.

    physical_error_rate: Error rate of physical qubits (e.g., 0.001)
    distance: Surface code distance

    Returns: Logical error rate (much lower)
    """
    # Simplified formula (real calculation more complex)
    # Logical error rate ≈ (physical_error_rate / threshold)^((distance+1)/2)

    threshold = 0.01  # Surface code threshold (~1%)

    if physical_error_rate > threshold:
        raise ValueError("Physical error rate too high for surface code")

    logical_error_rate = (physical_error_rate / threshold) ** ((distance + 1) / 2)

    return logical_error_rate

# Examples
print(f"Physical error rate: 0.1%")
print(f"Distance=3: Logical error rate: {calculate_logical_error_rate(0.001, 3):.2e}")
print(f"Distance=5: Logical error rate: {calculate_logical_error_rate(0.001, 5):.2e}")
print(f"Distance=7: Logical error rate: {calculate_logical_error_rate(0.001, 7):.2e}")

# Output:
# Physical error rate: 0.1%
# Distance=3: Logical error rate: 1.00e-06  (0.0001%)
# Distance=5: Logical error rate: 1.00e-09  (0.0000001%)
# Distance=7: Logical error rate: 1.00e-12  (0.0000000001%)

def estimate_shor_requirements(rsa_key_size):
    """
    Estimate quantum resources needed to break RSA encryption.

    rsa_key_size: RSA key size in bits (e.g., 2048, 4096)
    """
    # Shor's algorithm requires ~2n logical qubits for n-bit RSA
    logical_qubits_needed = 2 * rsa_key_size

    # With distance-7 surface codes (industry target)
    distance = 7
    physical_qubits_per_logical = distance ** 2

    total_physical_qubits = logical_qubits_needed * physical_qubits_per_logical

    # Estimated circuit depth (gates)
    circuit_depth = rsa_key_size ** 3  # Approximate

    return {
        'logical_qubits': logical_qubits_needed,
        'physical_qubits': total_physical_qubits,
        'circuit_depth': circuit_depth
    }

# RSA-2048 (current standard)
rsa_2048 = estimate_shor_requirements(2048)
print(f"To break RSA-2048:")
print(f"  Logical qubits: {rsa_2048['logical_qubits']:,}")
print(f"  Physical qubits: {rsa_2048['physical_qubits']:,}")
print(f"  Circuit depth: {rsa_2048['circuit_depth']:,}")

# RSA-4096 (post-quantum migration target)
rsa_4096 = estimate_shor_requirements(4096)
print(f"\nTo break RSA-4096:")
print(f"  Logical qubits: {rsa_4096['logical_qubits']:,}")
print(f"  Physical qubits: {rsa_4096['physical_qubits']:,}")
print(f"  Circuit_depth: {rsa_4096['circuit_depth']:,}")

# Output:
# To break RSA-2048:
#   Logical qubits: 4,096
#   Physical qubits: 200,704
#   Circuit depth: 8,589,934,592
#
# To break RSA-4096:
#   Logical qubits: 8,192
#   Physical qubits: 401,408
#   Circuit depth: 68,719,476,736

# Real quantum hardware (IBM, Google, etc.)
from qiskit import IBMQ

# Load account (requires token from IBM Quantum)
IBMQ.save_account('YOUR_API_TOKEN')
provider = IBMQ.load_account()

# Get fault-tolerant backend (when available ~2028)
backend = provider.get_backend('ibm_quantum_surface_code')  # Hypothetical

# Run error-corrected circuit
job = backend.run(surface_code.circuit, shots=1000)
result = job.result()