Breaking a 4-Bit Elliptic Curve Key with a 133-Qubit Quantum Computer

Code:

# Main circuit
# Imports
import logging, json
from math import gcd
import numpy as np
from qiskit import QuantumCircuit, QuantumRegister, ClassicalRegister, transpile
from qiskit.circuit.library import UnitaryGate, QFT
from qiskit_ibm_runtime import QiskitRuntimeService, SamplerV2
from qiskit.visualization import plot_histogram
import matplotlib.pyplot as plt
import pandas as pd

# IBMQ
TOKEN = "YOUR_IBMQ_API_KEY"
INSTANCE = "YOUR_IBMQ_CRN"
BACKEND = "ibm_torino"
CAL_CSV = "/Users/steventippeconnic/Downloads/ibm_torino_calibrations_2025-06-24T17_23_41Z.csv"
SHOTS = 16384

# Toy‑curve parameters  (order‑16 subgroup of E(F_p))
ORDER = 16 # |E(F_p)| = 16
P_IDX = 1  # generator P  ->  index 1
Q_IDX = 7  # public point Q = kP, here “7 mod 16”

# Logging helper
logging.basicConfig(level=logging .INFO,
                    format="%(asctime)s | %(levelname)s | %(message)s")
log = logging.getLogger(__name__)

# Calibration‑based qubit pick
def best_qubits(csv_path: str, n: int) -> list[int]:
    df = pd .read_csv(csv_path)
    df.columns = df.columns.str.strip()
    winners = (
        df.sort_values(["√x (sx) error", "T1 (us)", "T2 (us)"],
                       ascending=[True, False, False])
        ["Qubit"].head(n).tolist()
    )
    log .info("Best physical qubits: %s", winners)
    return winners

N_Q_TOTAL = 4 + 4 + 4          # a, b, point registers
PHYSICAL  = best_qubits(CAL_CSV, N_Q_TOTAL)

# Constant‑adder modulo 16 as a reusable gate
def add_const_mod16_gate(c: int) -> UnitaryGate:
    """Return a 4‑qubit gate that maps |x⟩ ↦ |x+c (mod 16)⟩."""
    mat = np.zeros((16, 16))
    for x in range(16):
        mat[(x + c) % 16, x] = 1
    return UnitaryGate(mat, label=f"+{c}")

ADDERS = {c: add_const_mod16_gate(c) for c in range(1, 16)}

def controlled_add(qc: QuantumCircuit, ctrl_qubit, point_reg, constant):
    """Apply |x⟩ → |x+constant (mod 16)⟩ controlled by one qubit."""
    qc.append(ADDERS[constant].control(), [ctrl_qubit, *point_reg])

# Oracle  U_f : |a⟩|b⟩|0⟩ ⟶ |a⟩|b⟩|aP + bQ⟩   (index arithmetic mod 16)
def ecdlp_oracle(qc, a_reg, b_reg, point_reg):
    for i in range(4):                          # bit 2^i of a
        constant = (P_IDX * (1 << i)) % ORDER
        if constant:
            controlled_add(qc, a_reg[i], point_reg, constant)

    for i in range(4):
        constant = (Q_IDX * (1 << i)) % ORDER
        if constant:
            controlled_add(qc, b_reg[i], point_reg, constant)

# Build the full Shor circuit
def shor_ecdlp_circuit() -> QuantumCircuit:
    a = QuantumRegister(4, "a")
    b = QuantumRegister(4, "b")
    p = QuantumRegister(4, "p")
    c = ClassicalRegister(8, "c")
    qc = QuantumCircuit(a, b, p, c, name="ECDLP_16pts")

    qc.h(a); qc.h(b)
    ecdlp_oracle(qc, a, b, p)
    qc.barrier()

    qc.append(QFT(4, do_swaps=False), a)
    qc.append(QFT(4, do_swaps=False), b)

    qc.measure(a, c[:4])
    qc.measure(b, c[4:])

    return qc

# IBM Runtime execution
service = QiskitRuntimeService(channel="ibm_cloud",
                               token=TOKEN,
                               instance=INSTANCE)
backend = service.backend(BACKEND)
log .info("Backend → %s", backend .name)

qc_raw = shor_ecdlp_circuit()
trans   = transpile(qc_raw,
                    backend         = backend,
                    initial_layout  = PHYSICAL,
                    optimization_level = 3)
log .info("Circuit depth %d, gate counts %s", trans.depth(), trans.count_ops())

sampler = SamplerV2(mode=backend)
job     = sampler .run([trans], shots=SHOTS)
result  = job.result()

# Classical post‑processing
creg_name = trans.cregs[0].name
counts_raw = result[0].data.__getattribute__(creg_name).get_counts()

def bits_to_int(bs): return int(bs[::-1], 2)

counts = {(bits_to_int(k[4:]), bits_to_int(k[:4])): v
          for k, v in counts_raw.items()}
top = sorted(counts.items(), key=lambda kv: kv[1], reverse=True)

# Success criteria. Check top 10 invertible rows for k = 7
top_invertibles = []
for (a_val, b_val), freq in top:
    if gcd(b_val, ORDER) != 1:
        continue
    inv_b       = pow(b_val, -1, ORDER)
    k_candidate = (-a_val * inv_b) % ORDER
    top_invertibles.append(((a_val, b_val), k_candidate, freq))
    if len(top_invertibles) == 10:
        break

# Check for success and print results
found_k7 = any(k == 7 for (_, k, _) in top_invertibles)

if found_k7:
    print("\nSUCCESS — k = 7 found in top 10 results\n")
else:
    print("\nWARNING — k = 7 NOT found in top 10 results\n")

print("Top 10 invertible (a, b) pairs and recovered k:")
for (a, b), k, count in top_invertibles:
    tag = " <<<" if k == 7 else ""
    print(f"  (a={a:2}, b={b:2})  →  k = {k:2}   (count = {count}){tag}")

# Save raw data
out = {
    "experiment": "ECDLP_16pts_Shors",
    "backend": backend .name,
    "physical_qubits": PHYSICAL,
    "shots": SHOTS,
    "counts": counts_raw
}
JSON_PATH = "/Users/steventippeconnic/Documents/QC/Shors_ECC_4_Bit_Key_0.json"
with open(JSON_PATH, "w") as fp:
    json.dump(out, fp, indent=4)
log .info("Results saved → %s", JSON_PATH)

# End


/////////////////////////////////////////////////////////////////

          
# Code for all visuals from experiment JSON
# imports
import json
import matplotlib.pyplot as plt
import numpy as np
from math import gcd
from collections import Counter
import math

# Load result data
path = '/Users/steventippeconnic/Documents/QC/Shors_ECC_4_Bit_Key_0.json'
with open(path, 'r') as f:
    data = json.load(f)

raw_counts = data['counts']
ORDER = 16

# Helper: convert bitstring to integers (a, b)
def decode(bitstr):
    a = int(bitstr[0:4][::-1], 2)
    b = int(bitstr[4:8][::-1], 2)
    return (a, b)

counts = {}
for bitstr, count in raw_counts.items():
    a, b = decode(bitstr)
    counts[(a, b)] = counts.get((a, b), 0) + count

# Sort for top list
sorted_counts = sorted(counts.items(), key=lambda kv: kv[1], reverse=True)
top10 = sorted_counts[:10]

# Helper: compute recovered k for valid (a, b)
def recover_k(a, b):
    if gcd(b, ORDER) != 1:
        return None
    return (-a * pow(b, -1, ORDER)) % ORDER

# Interference Ridge - Full 2D Grid 

grid = np.zeros((ORDER, ORDER))
for (a, b), c in counts.items():
    grid[a, b] = c

plt.figure(figsize=(8, 6))
plt.imshow(grid, origin='lower', cmap='viridis')
plt.colorbar(label='Shot Count')
plt.title("2D Interference Pattern over ℤ₁₆ × ℤ₁₆")
plt.xlabel('b register')
plt.ylabel('a register')
plt.grid(False)
plt.show()

# Modular Line Geometry - Plot of Recovered k 
k_plane = np.full((ORDER, ORDER), np.nan)
for (a, b), c in counts.items():
    if gcd(b, ORDER) == 1:
        k = recover_k(a, b)
        k_plane[a, b] = k

plt.figure(figsize=(8, 6))
plt.imshow(k_plane, origin='lower', cmap='twilight', vmin=0, vmax=ORDER-1)
plt.colorbar(label='Recovered k value')
plt.title("Modular Key Recovery Landscape")
plt.xlabel('b register')
plt.ylabel('a register')
plt.grid(False)
plt.show()

# Histogram of k Candidates from Top 50 Pairs 
k_freq = {}
for (a, b), c in sorted_counts:
    k = recover_k(a, b)
    if k is not None:
        k_freq[k] = k_freq.get(k, 0) + c

ks, cs = zip(*sorted(k_freq.items(), key=lambda kv: kv[1], reverse=True))
bars = plt.bar(ks, cs, color='crimson')
for bar, k in zip(bars, ks):
    if k == 7:
        bar.set_color('gold')

plt.title("Histogram of Recovered k Candidates (Top 50)")
plt.xlabel("Recovered k")
plt.ylabel("Aggregated Shot Count")
plt.xticks(range(0, ORDER))
plt.grid(True, linestyle='--', alpha=0.5)
plt.show()

# Heatmap of Invertibility and Signal Strength 
invert_map = np.zeros((ORDER, ORDER))
for (a, b), c in counts.items():
    invert_map[a, b] = c if gcd(b, ORDER) == 1 else 0

plt.figure(figsize=(8, 6))
plt.imshow(invert_map, origin='lower', cmap='plasma')
plt.colorbar(label='Shot Count (Invertible b only)')
plt.title("Signal Density from Invertible b Pairs Only")
plt.xlabel('b register')
plt.ylabel('a register')
plt.grid(False)
plt.show()

# Extract counts from bitstring results
counts = data['counts']
bit_counts = Counter(counts)

# Convert bitstrings to (a, b) integer tuples
pairs = []
pair_counts = []
for bitstring, count in bit_counts.items():
    bitstring = bitstring.replace(" ", "")
    a = int(bitstring[:4], 2)
    b = int(bitstring[4:], 2)
    pairs.append((a, b))
    pair_counts.append(count)

# Create a 16x16 matrix from the counts
counts_matrix = np.zeros((16, 16))
for (a, b), count in zip(pairs, pair_counts):
    counts_matrix[a, b] = count

# Define modular inverse function for Z_16
def modinv(x, mod=16):
    try:
        return pow(x, -1, mod)
    except ValueError:
        return None

# Entropy Heatmap: Information richness of each measurement 
entropy_map = np.zeros((16, 16))
total_shots = np.sum(counts_matrix)
prob_matrix = counts_matrix / total_shots
for a in range(16):
    for b in range(16):
        p = prob_matrix[a, b]
        entropy_map[a, b] = -p * np.log2(p) if p > 0 else 0

plt.figure()
plt.imshow(entropy_map, origin='lower', cmap='inferno')
plt.title('Local Entropy Map of Measurement Space')
plt.xlabel('b register')
plt.ylabel('a register')
plt.colorbar(label='Entropy Contribution')
plt.show()

# Invertibility Mask: which b values are algebraically usable 
invert_map = np.zeros((16, 16))
for a in range(16):
    for b in range(16):
        invert_map[a, b] = 1 if modinv(b) is not None else 0

plt.figure()
plt.imshow(invert_map, origin='lower', cmap='bone')
plt.title('Invertibility Mask for b Register')
plt.xlabel('b register')
plt.ylabel('a register')
plt.colorbar(label='Invertible (1 = yes)')
plt.show()

# Log-scaled signal intensity (invertible b only) 
success_score = np.zeros((16, 16))
for a in range(16):
    for b in range(16):
        count = counts_matrix[a, b]
        if count > 0 and modinv(b) is not None:
            success_score[a, b] = np.log(count)

plt.figure()
plt.imshow(success_score, origin='lower', cmap='viridis')
plt.title('Log Signal Strength (Invertible b Only)')
plt.xlabel('b register')
plt.ylabel('a register')
plt.colorbar(label='log(Shot Count)')
plt.show()

# Modular Aliasing Histogram: how often each k = -a * b⁻¹ mod 16 appears
k_map = np.zeros(16)
for (a, b), count in zip(pairs, pair_counts):
    b_inv = modinv(b)
    if b_inv is not None:
        k = (-a * b_inv) % 16
        k_map[k] += count

plt.figure()
plt.bar(range(16), k_map, color='teal')
plt.title('Modular Aliasing Frequency per k')
plt.xlabel('Recovered k')
plt.ylabel('Aggregated Count')
plt.xticks(range(16))
plt.show()

shot_dict = {bitstring: count for bitstring, count in raw_counts.items()}

# Prepare a, b value matrices
def parse_bitstring(s):
    s = s.replace(" ", "")[::-1]  # Endian reversal
    a = int(s[:4][::-1], 2)
    b = int(s[4:][::-1], 2)
    return a, b

a_vals, b_vals, counts = [], [], []
for bitstring, count in shot_dict.items():
    a, b = parse_bitstring(bitstring)
    a_vals.append(a)
    b_vals.append(b)
    counts.append(count)

A = np.array(a_vals)
B = np.array(b_vals)
C = np.array(counts)

# Phase-Collision Map that shows (a + k*b) mod 16 concentration with assumed k=7
plt.figure(figsize=(6, 5))
k_guess = 7
collision_bins = np.zeros((16, 16))
for a, b, c in zip(A, B, C):
    phase = (a + k_guess * b) % 16
    collision_bins[phase][b] += c
plt.imshow(collision_bins, cmap='viridis', origin='lower')
plt.colorbar(label='Interference Strength')
plt.title(f'Phase Collision Map assuming k={k_guess}')
plt.xlabel('b value')
plt.ylabel('(a + 7b) mod 16')
plt.show()

# Mutual Information Heatmap 
joint_dist = np.zeros((16, 16))
for a, b, c in zip(A, B, C):
    joint_dist[a][b] += c
joint_dist /= np.sum(joint_dist)
marg_a = np.sum(joint_dist, axis=1)
marg_b = np.sum(joint_dist, axis=0)
mi_map = np.zeros_like(joint_dist)
for i in range(16):
    for j in range(16):
        if joint_dist[i][j] > 0:
            mi_map[i][j] = joint_dist[i][j] * np.log2(joint_dist[i][j] / (marg_a[i] * marg_b[j]))
plt.figure(figsize=(6, 5))
plt.imshow(mi_map, cmap='magma', origin='lower')
plt.colorbar(label='Mutual Information Contribution')
plt.title('Mutual Information Heatmap (a ↔ b)')
plt.xlabel('b')
plt.ylabel('a')
plt.show()

# Qubit Parity Band Analysis 
parity_band = np.zeros(2)
for s, c in shot_dict.items():
    s_clean = s.replace(" ", "")
    parity = sum(int(bit) for bit in s_clean) % 2
    parity_band[parity] += c
plt.figure(figsize=(5, 4))
plt.bar(['Even Parity', 'Odd Parity'], parity_band, color=['steelblue', 'darkorange'])
plt.title("Parity Distribution Across All Qubits")
plt.ylabel("Shot Count")
plt.show()

# Error Surface - Expected vs Observed Ridge 
u_v_grid = np.zeros((16, 16))
for bitstring, count in shot_dict.items():
    s = bitstring.replace(" ", "")[::-1]
    a = int(s[:4][::-1], 2)
    b = int(s[4:][::-1], 2)
    k = 7
    if math.gcd(b, 16) == 1:
        inferred_u = a
        inferred_v = b
        delta = (inferred_u + k * inferred_v) % 16
        u_v_grid[inferred_u][inferred_v] += delta * count
plt.figure(figsize=(6, 5))
plt.imshow(u_v_grid, origin='lower', cmap='plasma')
plt.title("Deviation from Ideal u + kv ≡ 0 Surface (k = 7)")
plt.xlabel('v')
plt.ylabel('u')
plt.colorbar(label='Weighted Residue (mod 16)')
plt.show()

# End