Breaking a 6-Bit Elliptic Curve Key using IBM’s 133-Qubit Quantum Computer

Conclusion
In the end, this experiment successfully broke a 6‑bit elliptic‑curve key with a Shor‑style quantum attack on IBM’s 133‑qubit processor. A 18‑qubit circuit (12 logical + 6 ancilla) encoded the oracle over ℤ₆₄, never inserting the secret k and instead weaving it into the phase of the superposition. Despite a 340k‑gate, 3.4 × 10⁵‑layer depth the device generated a discernible interference valley only ± 3 steps wide around the ideal line a + 42b ≡ 0. Classical post‑processing, mod‑inverse filtering and top‑100 enumeration, recovered k = 42 three times within the 100 highest‑count ridge cells, confirming success after 16,384 shots. 2‑D and 3‑D visualizations exposed the signed‑distance trench, and hot‑spot control registers b, confirming quantum coherence amplifies the correct modular relation. This shows that Shor’s algorithm continues to scale under deeper circuit regimes and that dictionary-based key recovery strategies (top 100 enumeration) remain viable as bit-length increases. All code, visualizations, calibration data, and raw backend results are available at https://github.com/SteveTipp/Qwork.github.io or via the project website www.qubits.work.

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
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-07-28T01_40_46Z.csv"
SHOTS    = 16384

# Toy‑curve parameters (order‑64 subgroup of E(F_p))
ORDER = 64  # |E(F_p)| = 64
P_IDX = 11  # Generator P  -> index 11
Q_IDX = 14  # Public point Q = kP, here “462 mod 64” for k = 42

# 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 = 6
N_Q_TOTAL = N_Q * 3  # a, b, point
PHYSICAL = best_qubits(CAL_CSV, N_Q_TOTAL)

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

ADDERS = {c: add_const_mod64_gate(c) for c in range(1, 64)}

def controlled_add(qc: QuantumCircuit, ctrl_qubit, point_reg, constant):
    """Apply |x⟩ → |x+constant (mod 64)⟩ 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 64)
def ecdlp_oracle(qc, a_reg, b_reg, point_reg):
    for i in range(N_Q):
        constant = (P_IDX * (1 << i)) % ORDER
        if constant:
            controlled_add(qc, a_reg[i], point_reg, constant)

    for i in range(N_Q):
        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(N_Q, "a")
    b = QuantumRegister(N_Q, "b")
    p = QuantumRegister(N_Q, "p")
    c = ClassicalRegister(N_Q * 2, "c")
    qc = QuantumCircuit(a, b, p, c, name="ECDLP_64pts")

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

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

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

    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[N_Q:]), bits_to_int(k[:N_Q])): v
          for k, v in counts_raw.items()}
top = sorted(counts.items(), key=lambda kv: kv[1], reverse=True)

# Success criteria. Check top 100 invertible rows for k = 42
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) == 100:
        break

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

if found_k:
    print("\nSUCCESS — k = 42 found in top 100 results\n")
else:
    print("\nWARNING — k = 42 NOT found in top 100 results\n")

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

# Save raw data
out = {
    "experiment": "ECDLP_64pts_Shors",
    "backend": backend.name,
    "physical_qubits": PHYSICAL,
    "shots": SHOTS,
    "counts": counts_raw
}
JSON_PATH = "/Users/steventippeconnic/Documents/QC/Shors_ECC_6_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, math, numpy as np, matplotlib.pyplot as plt
from collections import Counter, defaultdict
from mpl_toolkits.mplot3d import Axes3D  

PATH = '/Users/steventippeconnic/Documents/QC/Shors_ECC_6_Bit_Key_0.json'
ORDER = 64                                                         
N_Q   = 6                                                         
K     = 42 

# Helpers
def bits_to_int(bs: str) -> int:          
    return int(bs[::-1], 2)

with open(PATH) as fp:
    raw = json.load(fp)['counts']

# Parse raw quantum counts into heatmap, key histogram, invertibility mask, and ridge hits
def parse_counts(path):
    with open(path) as fp:
        raw = json.load(fp)['counts']
    heat   = np.zeros((ORDER, ORDER), dtype=int)
    k_hist = Counter()
    b_bar  = defaultdict(int)
    ridge  = []                            

    for bitstr, cnt in raw.items():
        a = bits_to_int(bitstr[N_Q:])     
        b = bits_to_int(bitstr[:N_Q])      
        heat[a, b] += cnt
        if math.gcd(b, ORDER) == 1:
            k = (-a * pow(b, -1, ORDER)) % ORDER
            k_hist[k] += cnt
            if k == 42:
                ridge.append((a, b, cnt))
        else:
            b_bar['non‑invertible'] += cnt
            continue
        b_bar['invertible'] += cnt
    return heat, k_hist, b_bar, ridge

heat, k_hist, b_bar, ridge = parse_counts(PATH)

# Parse once into arrays
all_counts   = []
invert_heat  = np.zeros((ORDER, ORDER), dtype=int)   
b_counts_inv = np.zeros(ORDER, dtype=int)           
invmap_heat  = np.zeros((ORDER, ORDER), dtype=int)   
zipf_vals    = []

for bitstring, cnt in raw.items():
    a = bits_to_int(bitstring[N_Q:])
    b = bits_to_int(bitstring[:N_Q])
    all_counts.append(cnt)

    if math.gcd(b, ORDER) == 1:
        invert_heat[a, b] += cnt                 
        b_counts_inv[b]   += cnt                 
        inv_b             = pow(b, -1, ORDER)
        invmap_heat[b, inv_b] += cnt       

# Collect ridge stats, distance distributions, and top (a,b) cells for 3D scatter plotting
ridge_counts   = np.zeros(ORDER, dtype=int)              
row_totals     = np.zeros(ORDER, dtype=int)               
dist_tot       = np.zeros(ORDER//2 + 1, dtype=int)        
dist_cells     = np.zeros_like(dist_tot, dtype=int)      
top_cells      = []                                       

for bitstr, cnt in raw.items():
    a = bits_to_int(bitstr[N_Q:])
    b = bits_to_int(bitstr[:N_Q])
    row_totals[b] += cnt

    d = (a + K*b) % ORDER;  d = min(d, ORDER - d)
    dist_tot[d]   += cnt
    dist_cells[d] += 1

    if d == 0:
        ridge_counts[b] += cnt

    # collect for 3-D scatter
    if math.gcd(b, ORDER) == 1:
        k_val = (-a * pow(b, -1, ORDER)) % ORDER
        top_cells.append((cnt, a, b, k_val == K))

# keep only the 200 loudest for readability
top_cells.sort(reverse=True)
top_cells = top_cells[:200]      

residue_hist = Counter()                        
polar_theta  = []                              
polar_size   = []

for bitstr, cnt in raw.items():
    a = bits_to_int(bitstr[N_Q:])
    b = bits_to_int(bitstr[:N_Q])
    heat[a, b] += cnt

    # ridge bookkeeping
    if a == (-K * b) % ORDER:
        ridge_counts[b] += cnt

    # residue distribution
    residue_hist[(a + K * b) % ORDER] += cnt

    # polar points only for recovered k==42
    if math.gcd(b, ORDER) == 1:
        if (-a * pow(b, -1, ORDER)) % ORDER == K:
            theta = np.arctan2(b, -a)           
            polar_theta.append(theta)
            polar_size.append(cnt)

# Parse raw bitstring data into heatmap, ridge, and distance histograms
dmap = np.zeros((ORDER, ORDER),    dtype=int)   
ridge_pairs = []                                
dist_hist = np.zeros(ORDER//2 + 1, dtype=int)   

for bitstr, c in raw.items():
    a = bits_to_int(bitstr[N_Q:])
    b = bits_to_int(bitstr[:N_Q])
    heat[a, b] += c

    d = (a + K * b) % ORDER
    d = min(d, ORDER - d)                      
    dmap[d, b] += c
    dist_hist[d] += c

    if d == 0:
        ridge_pairs.append((c, a, b))

ridge_pairs.sort(reverse=True)                  
cum = np.cumsum([c for c, _, _ in ridge_pairs])

# Compute signed distance surface and modular key estimates
dist_surf   = np.zeros((ORDER, ORDER), dtype=int)  
k_b_counts  = np.zeros((ORDER, ORDER), dtype=int)
ridge_line  = np.zeros(ORDER, dtype=int)
sc_a, sc_b, sc_c = [], [], []                       

for bits, cnt in raw.items():
    a, b = bits_to_int(bits[N_Q:]), bits_to_int(bits[:N_Q])

    d = (a + K*b) % ORDER           
    if d >= ORDER//2:               
        d -= ORDER
    d_idx = d + ORDER//2            
    dist_surf[d_idx, b] += cnt

    if math.gcd(b, ORDER) == 1:
        k_val = (-a * pow(b, -1, ORDER)) % ORDER
        k_b_counts[k_val, b] += cnt
        if k_val == K:
            sc_a.append(a); sc_b.append(b); sc_c.append(cnt)

    if d == 0:
        ridge_line[b] += cnt

# 64x64 heatmap of (a,b) counts 
plt.figure(figsize=(6,6))
plt.title('Heatmap of raw (a, b) counts')
plt.imshow(heat, cmap='viridis', origin='lower', interpolation='nearest')
plt.colorbar(label='counts')
plt.xlabel('b  (0-63)');  plt.ylabel('a  (0-63)')
plt.show()

# Histogram of recovered k values 
plt.figure()
plt.title('Recovered k distribution')
plt.bar(list(k_hist.keys()), list(k_hist.values()))
plt.xticks(range(0, 64, 8))
plt.xlabel('k'); plt.ylabel('total counts')
plt.show()

# Invertible vs non‑invertible b counts
plt.figure()
plt.title('Counts by b invertibility')
plt.bar(b_bar.keys(), b_bar.values(), color=['tab:green', 'tab:red'])
plt.ylabel('total counts')
plt.show()

# Scatter of (a, b) pairs decoding to k = 42 
if ridge:
    a_vals, b_vals, sizes = zip(*ridge)
    sizes = np.array(sizes) * 4           
    plt.figure(figsize=(6,6))
    plt.title('(a, b) pairs yielding k = 42')
    plt.scatter(b_vals, a_vals, s=sizes, alpha=0.7, edgecolor='k')
    plt.xlabel('b'); plt.ylabel('a')
    plt.xlim(-1, 64); plt.ylim(-1, 64)
    plt.grid(True, lw=0.3)
    plt.show()
else:
    print('No (a, b) pair produced k = 42.')

# Zipf tail (rank‑vs‑frequency)
zipf_vals = sorted(all_counts, reverse=True)
ranks     = np.arange(1, len(zipf_vals)+1)

# Heatmap over invertible b (gcd(b,64)=1)
plt.figure(figsize=(6,6))
plt.title('Heatmap over invertible b (gcd(b,64)=1)')
plt.imshow(invert_heat, cmap='plasma', origin='lower', interpolation='nearest')
plt.colorbar(label='counts')
plt.xlabel('b  (invertible only)'); plt.ylabel('a')
plt.show()

# Counts per invertible b
plt.figure()
plt.title('Counts per invertible b')
plt.bar(np.arange(ORDER), b_counts_inv, color='tab:purple')
plt.xticks(range(0,64,8))
plt.xlabel('b'); plt.ylabel('total counts')
plt.show()

# Mod‑inverse map:  b  versus  b⁻¹  (counts as intensity)
plt.figure(figsize=(6,6))
plt.title('Mod‑inverse map:  b  versus  b⁻¹  (counts as intensity)')
plt.imshow(invmap_heat, cmap='viridis', origin='lower', interpolation='nearest')
plt.colorbar(label='counts')
plt.xlabel('b'); plt.ylabel('b⁻¹   (mod 64)')
plt.show()

# Zipf‑like rank plot of bitstring counts
plt.figure()
plt.title('Zipf‑like rank plot of bitstring counts')
plt.loglog(ranks, zipf_vals, marker='.', linestyle='none')
plt.xlabel('rank'); plt.ylabel('count')
plt.grid(True, which='both', ls='--', lw=0.3)
plt.show()

# 3‑D surface of full (a,b) counts  
fig = plt.figure(figsize=(7,5))
ax  = fig.add_subplot(111, projection='3d')
A, B = np.meshgrid(np.arange(ORDER), np.arange(ORDER), indexing='ij')
ax.plot_surface(B, A, heat, cmap='viridis', linewidth=0, antialiased=False)
ax.set_title('3‑D surface of raw counts')
ax.set_xlabel('b'); ax.set_ylabel('a'); ax.set_zlabel('counts')
plt.tight_layout(); plt.show()

# 3‑D ridge profile  a = −Kb (mod 64)  
fig = plt.figure()
ax  = fig.add_subplot(111, projection='3d')
b_vals = np.arange(ORDER)
a_vals = (-K * b_vals) % ORDER
ax.bar3d(b_vals, a_vals, np.zeros(ORDER), 0.7, 0.7, ridge_counts,
         shade=True, color='crimson')
ax.set_title('Amplitude along the k = 42 ridge')
ax.set_xlabel('b'); ax.set_ylabel('a = −42·b (mod 64)'); ax.set_zlabel('counts')
plt.tight_layout(); plt.show()

# Residue histogram  r = (a + K b) mod 64  
plt.figure()
r_keys = list(range(ORDER))
r_vals = [residue_hist[r] for r in r_keys]
plt.bar(r_keys, r_vals, color='tab:orange')
plt.title('Residue distribution  (a + 42 b) mod 64')
plt.xticks(range(0,64,8)); plt.xlabel('residue r'); plt.ylabel('counts')
plt.tight_layout(); plt.show()

# Polar scatter of vectors (−a, b) decoding to k = 42 
plt.figure(figsize=(6,6))
ax = plt.subplot(projection='polar')
sizes = np.array(polar_size) * 3                 
ax.scatter(polar_theta, np.ones_like(polar_theta),
           s=sizes, alpha=0.7, c='teal')
ax.set_rticks([]); ax.set_title('Phase‑ridge angles for k = 42')
plt.show()

# Heatmap with ridge overlay 
plt.figure(figsize=(6,6))
plt.title('Heatmap with ridge cells highlighted')
plt.imshow(heat, cmap='inferno', origin='lower', interpolation='nearest')
ra, rb = zip(*[(a,b) for _, a, b in ridge_pairs])
plt.scatter(rb, ra, marker='x', c='white', s=30, lw=0.8)
plt.colorbar(label='counts'); plt.xlabel('b'); plt.ylabel('a')
plt.show()

# Histogram of |distance| to ridge   
plt.figure()
plt.title('|a + 42 b|₍₆₄₎ distance histogram')
plt.bar(range(len(dist_hist)), dist_hist, color='tab:cyan')
plt.yscale('log'); plt.xlabel('distance d'); plt.ylabel('counts (log‐scale)')
plt.show()

# Cumulative energy captured by top‑n ridge points 
top_n   = np.arange(1, len(cum)+1)
plt.figure()
plt.title('Cumulative ridge power vs  n  (sorted by count)')
plt.plot(top_n, cum / cum[-1] * 100, marker='.')
plt.xscale('log'); plt.xlabel('top‑n ridge pairs'); plt.ylabel('% of ridge counts')
plt.grid(True, which='both', ls='--', lw=0.3)
plt.show()

# Ridge amplitude vs b 
plt.figure()
plt.title('Ridge amplitude  (a = −42·b mod 64)')
plt.plot(range(ORDER), ridge_counts, marker='o', lw=1)
plt.xticks(range(0,64,8)); plt.xlabel('b'); plt.ylabel('counts along ridge')
plt.tight_layout(); plt.show()

# Signed‑distance surface
fig = plt.figure(figsize=(7,5)); ax = fig.add_subplot(111, projection='3d')
signed_range = np.arange(-ORDER//2, ORDER//2)       
B, D = np.meshgrid(np.arange(ORDER), signed_range, indexing='xy')
ax.plot_surface(B, D, dist_surf, cmap='cividis', linewidth=0)
ax.set_title('Signed‑distance surface to ridge'); ax.set_xlabel('b')
ax.set_ylabel('signed distance'); ax.set_zlabel('counts')
plt.tight_layout(); plt.show()

# Recovered‑k vs b surface
fig = plt.figure(figsize=(7,5)); ax = fig.add_subplot(111, projection='3d')
KX, BX = np.meshgrid(np.arange(ORDER), np.arange(ORDER), indexing='ij')
ax.plot_surface(BX, KX, k_b_counts, cmap='plasma', linewidth=0)
ax.set_title('Counts by recovered k and b'); ax.set_xlabel('b')
ax.set_ylabel('recovered k'); ax.set_zlabel('counts')
plt.tight_layout(); plt.show()

# Ridge wire
fig = plt.figure(); ax = fig.add_subplot(111, projection='3d')
a_vals = (-K * np.arange(ORDER)) % ORDER
ax.plot(np.arange(ORDER), a_vals, ridge_line, color='gold', lw=2)
ax.set_title('Ridge wire  (a = −42·b mod 64)')
ax.set_xlabel('b'); ax.set_ylabel('a'); ax.set_zlabel('counts')
plt.tight_layout(); plt.show()

# Scatter of cells decoding to k = 42
fig = plt.figure(figsize=(7,5)); ax = fig.add_subplot(111, projection='3d')
ax.scatter(sc_b, sc_a, sc_c, c='lime', s=40, alpha=0.8)
ax.set_title('All cells yielding k = 42'); ax.set_xlabel('b')
ax.set_ylabel('a'); ax.set_zlabel('counts')
plt.tight_layout(); plt.show()

# End

# Code to compute recovery rate for the correct key k = 42 
# Imports
import json
import math
from collections import Counter

# Experiment parameters
ORDER = 64
N_Q = 6  # Number of qubits per register
K = 42   # Target scalar key

def bits_to_int(bits):
    return int(bits, 2)

# Compute K recovery rate
def compute_k_recovery_rate(path):
    with open(path) as f:
        data = json.load(f)

    counts = data["counts"]
    total_invertible = 0
    total_k42 = 0

    for bitstr, cnt in counts.items():
        a = bits_to_int(bitstr[N_Q:])   # rightmost bits = a
        b = bits_to_int(bitstr[:N_Q])   # leftmost bits = b

        if math.gcd(b, ORDER) == 1:
            total_invertible += cnt
            k = (-a * pow(b, -1, ORDER)) % ORDER
            if k == K:
                total_k42 += cnt

    recovery_rate = total_k42 / total_invertible if total_invertible else 0

    print(f"Correct key (k = {K}) total count: {total_k42}")
    print(f"Total invertible (a, b) counts: {total_invertible}")
    print(f"Recovery rate: {recovery_rate:.4f} ({recovery_rate*100:.2f}%)")

# Run the analysis on 6-bit ECC experiment data
compute_k_recovery_rate('/Users/steventippeconnic/Documents/QC/Shors_ECC_6_Bit_Key_0.json')

# End

// Code for Three-js visual from experiment JSON
// Load JSON data - 6-bit backend result
async function loadQuantumData() {
  const response = await fetch('Shors_ECC_6_Bit_Key_0.json');
  const data     = await response.json();
  return data.counts;
}

// Three.js setup
const scene    = new THREE.Scene();
const camera   = new THREE.PerspectiveCamera(75, window.innerWidth / window.innerHeight, 0.1, 1000);
const renderer = new THREE.WebGLRenderer({ antialias: true });
renderer.setSize(window.innerWidth, window.innerHeight);
document.body.appendChild(renderer.domElement);

// Orbit controls
const controls = new THREE.OrbitControls(camera, renderer.domElement);

// Lighting
const light = new THREE.DirectionalLight(0xffffff, 1);
light.position.set(10, 10, 10);
scene.add(light);

// Camera
camera.position.set(0, 40, 75);
camera.lookAt(0, 0, 0);

// Globals
let mesh;
let redDots          = [];
let highlightSpheres = [];
let time = 0;

// Helpers 
// All ridge hits measured
function getAllCorrectSet(counts, halfBits, K) {
  const mod = 1 << halfBits;
  const set = new Set();
  for (const bit of Object.keys(counts)) {
    const u = parseInt(bit.slice(halfBits), 2);
    const v = parseInt(bit.slice(0, halfBits), 2);
    if ((u + K * v) % mod === 0) set.add(`${u},${v}`);
  }
  return set;
}

// Top‑3 highest‑count ridge hits
function getTop3CorrectSet(counts, halfBits, K) {
  const mod = 1 << halfBits;
  return new Set(
    Object.entries(counts)
      .filter(([bit]) => {
        const u = parseInt(bit.slice(halfBits), 2);
        const v = parseInt(bit.slice(0, halfBits), 2);
        return (u + K * v) % mod === 0;
      })
      .sort((a, b) => b[1] - a[1])
      .slice(0, 3)
      .map(([bit]) => {
        const u = parseInt(bit.slice(halfBits), 2);
        const v = parseInt(bit.slice(0, halfBits), 2);
        return `${u},${v}`;
      })
  );
}

// Parse counts -> matrix
function parseCountsToGrid(counts) {
  const sampleKey = Object.keys(counts)[0];
  const halfBits  = sampleKey.length / 2;        
  const gridSize  = 1 << halfBits;               
  const matrix    = Array.from({ length: gridSize }, () => Array(gridSize).fill(0));
  let   maxCount  = 0;

  for (const [bit, cnt] of Object.entries(counts)) {
    const u = parseInt(bit.slice(halfBits), 2);
    const v = parseInt(bit.slice(0, halfBits), 2);
    matrix[u][v] = cnt;
    if (cnt > maxCount) maxCount = cnt;
  }
  return { matrix, maxCount, gridSize, halfBits };
}

// Build surface mesh + dots + highlight spheres
function createWaveSurface(matrix, maxCount, gridSize, correctSet, top3Set) {
  const WIDTH    = gridSize * 2;     
  const SEGMENTS = WIDTH;            
  const STEP     = WIDTH / (gridSize - 1);   

  const geometry = new THREE.PlaneGeometry(WIDTH, WIDTH, SEGMENTS, SEGMENTS);
  const colors   = [];
  const gridX    = SEGMENTS + 1;
  const gridY    = SEGMENTS + 1;

  for (let y = 0; y < gridY; y++) {
    for (let x = 0; x < gridX; x++) {
      const u = Math.floor(x / (gridX / gridSize));
      const v = Math.floor(y / (gridY / gridSize));
      const key = `${u},${v}`;

      // default green
      let color = new THREE.Color(0x00ff00);
      // ridge hits -> blue
      if (correctSet.has(key)) color = new THREE.Color(0x0000ff);

      colors.push(color.r, color.g, color.b);
    }
  }

  geometry.setAttribute('color', new THREE.Float32BufferAttribute(colors, 3));
  const material = new THREE.MeshBasicMaterial({ vertexColors: true, wireframe: true });

  mesh = new THREE.Mesh(geometry, material);
  mesh.rotation.x = -Math.PI / 2;
  mesh.userData   = { matrix, maxCount, gridSize, STEP };
  scene.add(mesh);

  // red dots every (u,v)
  const dotMat = new THREE.MeshBasicMaterial({ color: 0xff0000 });
  for (let v = 0; v < gridSize; v++) {
    for (let u = 0; u < gridSize; u++) {
      const dot = new THREE.Mesh(new THREE.SphereGeometry(0.15, 8, 8), dotMat);
      dot.position.set(
        (u - (gridSize - 1) / 2) * STEP,
        0,
        (v - (gridSize - 1) / 2) * STEP
      );
      redDots.push(dot);
      scene.add(dot);
    }
  }

  // yellow spheres on top‑3
  const hMat = new THREE.MeshBasicMaterial({ color: 0xffff00 });
  top3Set.forEach(key => {
    const [uStr, vStr] = key.split(',');
    const u = parseInt(uStr, 10);
    const v = parseInt(vStr, 10);
    const amp = matrix[u][v] / maxCount;
    const sph = new THREE.Mesh(new THREE.SphereGeometry(0.8, 20, 20), hMat);
    sph.position.set(
      (u - (gridSize - 1) / 2) * STEP,
      amp * 25 + 1.2,
      (v - (gridSize - 1) / 2) * STEP
    );
    sph.userData = { u, v };
    highlightSpheres.push(sph);
    scene.add(sph);
  });
}

// Animate
function animate() {
  requestAnimationFrame(animate);
  time += 0.2;
  if (!mesh) return;

  const { matrix, maxCount, gridSize, STEP } = mesh.userData;
  const pos   = mesh.geometry.attributes.position;
  const gridX = mesh.geometry.parameters.widthSegments + 1;
  const gridY = mesh.geometry.parameters.heightSegments + 1;

  for (let i = 0; i < pos.count; i++) {
    const x   = i % gridX;
    const y   = Math.floor(i / gridX);
    const u   = Math.floor(x / (gridX / gridSize));
    const v   = Math.floor(y / (gridY / gridSize));
    const amp = matrix[u]?.[v] ? matrix[u][v] / maxCount : 0;
    const wave = Math.sin((x + time) * 0.4) * Math.cos((y + time) * 0.4);
    pos.setZ(i, amp * wave * 25);
  }

  redDots.forEach((dot, idx) => {
    const u = idx % gridSize;
    const v = Math.floor(idx / gridSize);
    const amp = matrix[u]?.[v] ? matrix[u][v] / maxCount : 0;
    dot.position.y = amp * 25 + 0.6;
  });

  highlightSpheres.forEach(sph => {
    const { u, v } = sph.userData;
    const amp = matrix[u]?.[v] ? matrix[u][v] / maxCount : 0;
    sph.position.y = amp * 25 + 1.2;
  });

  pos.needsUpdate = true;
  controls.update();
  renderer.render(scene, camera);
}

// Main
loadQuantumData().then(counts => {
  const K = 42;
  const { matrix, maxCount, gridSize, halfBits } = parseCountsToGrid(counts);

  const correctSet = getAllCorrectSet(counts, halfBits, K);
  const top3Set    = getTop3CorrectSet(counts, halfBits, K);

  createWaveSurface(matrix, maxCount, gridSize, correctSet, top3Set);
  animate();
});


<!-- Index.html -->
<!DOCTYPE html>
<html lang="en">
<head>
  <meta charset="UTF-8">
  <meta name="viewport" content="width=device-width, initial-scale=1.0">
  <title>Quantum Wavefunction Visualization</title>

  <!-- Three.js -->
  <script src="https://cdnjs.cloudflare.com/ajax/libs/three.js/r128/three.min.js"></script>

  <!-- OrbitControls -->
  <script src="https://cdn.jsdelivr.net/npm/three@0.128.0/examples/js/controls/OrbitControls.js"></script>

  <style>
    body { margin: 0; overflow: hidden; }
    canvas { display: block; }
  </style>
</head>
<body>
  <script src="6_bit_ECC_break_visual.js"></script>
</body>
</html>