Recovering a 5-Bit Scalar in an Order-32 Elliptic Curve Subgroup via a Phase-Only Shor-Style Construction on IBM’s 156-Qubit ibm_marrakesh Quantum Processor

Code:

# Main circuit
# Imports
import logging, json
import numpy as np
from qiskit import QuantumCircuit, QuantumRegister, ClassicalRegister, transpile
from qiskit.circuit.library import QFT
from qiskit_ibm_runtime import QiskitRuntimeService, SamplerV2
import pandas as pd

# IBMQ
TOKEN = "YOUR_IBMQ_API_KEY"
INSTANCE = "YOUR_IBMQ_CRN"
BACKEND  = "ibm_marrakesh"
CAL_CSV  = "MOST_RECENT_ibm_marrakesh_calibrations.csv"
SHOTS    = 16384

# Toy parameters (order-32 subgroup)
ORDER = 32
K_SECRET = 7  # k = 7

# 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 = 5
N_Q_TOTAL = N_Q * 2
PHYSICAL = best_qubits(CAL_CSV, N_Q_TOTAL)

# Phase Only Oracle. Implements: |a⟩|b⟩ → exp(2πi(a + k b)/32) |a⟩|b⟩
def phase_oracle(qc, a_reg, b_reg):

    # Phase contribution from 'a'
    for i in range(N_Q):
        weight = (1 << i)
        angle = 2 * np.pi * weight / ORDER
        qc.rz(angle, a_reg[i])

    # Phase contribution from 'k * b'
    for i in range(N_Q):
        weight = (K_SECRET * (1 << i)) % ORDER
        angle = 2 * np.pi * weight / ORDER
        qc.rz(angle, b_reg[i])

# Build the Shor-style circuit
def shor_ecdlp_phase_circuit() -> QuantumCircuit:

    a = QuantumRegister(N_Q, "a")
    b = QuantumRegister(N_Q, "b")
    c = ClassicalRegister(N_Q * 2, "c")

    qc = QuantumCircuit(a, b, c, name="ECDLP_32pts_PhaseOnly")

    # Uniform superposition
    qc.h(a)
    qc.h(b)

    # Phase encoding of a + k b
    phase_oracle(qc, a, b)
    qc.barrier()

    # QFT on both registers
    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_phase_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)

# Phase-only success criteria: for ideal phase-only: a should lock to 31 and b should lock to (-k) mod 32.
A_TARGET = ORDER - 1           # 31
B_TARGET = (-K_SECRET) % ORDER # 25

a_counts = {a_val: 0 for a_val in range(ORDER)}
for (a_val, b_val), freq in counts.items():
    a_counts[a_val] += freq

p_a_target = a_counts[A_TARGET] / SHOTS

cond_total = 0
cond_b_target = 0
k_spectrum = {k_val: 0 for k_val in range(ORDER)}

for (a_val, b_val), freq in counts.items():
    if a_val != A_TARGET:
        continue
    cond_total += freq
    if b_val == B_TARGET:
        cond_b_target += freq
    k_hat = (-b_val) % ORDER
    k_spectrum[k_hat] += freq

p_b_target_given_a = (cond_b_target / cond_total) if cond_total > 0 else 0.0

k_sorted = sorted(k_spectrum.items(), key=lambda kv: kv[1], reverse=True)
k_best, k_best_count = k_sorted[0]

found_k = (k_best == K_SECRET)

print(f"\nShots: {SHOTS}")
top_ab, top_ct = top[0]
print(f"Top (a,b): {top_ab} count: {top_ct}")
print(f"P(a={A_TARGET}) = {p_a_target:.4f}")
print(f"P(b={B_TARGET} | a={A_TARGET}) = {p_b_target_given_a:.4f}")
print(f"Best k̂ by peak: {k_best} count: {k_best_count}")

if found_k:
    print(f"\nSUCCESS — k = {K_SECRET} is the dominant phase-only key peak\n")
else:
    print(f"\nWARNING — k = {K_SECRET} is NOT the dominant phase-only key peak\n")

print("Top 10 k̂ values from phase-only key spectrum (conditioned on a=31):")
for k_val, ct in k_sorted[:10]:
    tag = " <<<" if k_val == K_SECRET else ""
    print(f"  k̂ = {k_val:2}   (count = {ct}){tag}")

# Save raw data
out = {
    "experiment": "ECDLP_32pts_PhaseOnly",
    "backend": backend.name,
    "physical_qubits": PHYSICAL,
    "shots": SHOTS,
    "counts": counts_raw
}

JSON_PATH = "FILE_PATH_TO_SAVE_BACKEND_RESULT_JSON.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 numpy as np
import matplotlib.pyplot as plt

FILE_PATH = "FILE_PATH_TO_IMPORT_BACKEND_RESULT_JSON.json"

ORDER = 32
N_Q = 5

A_DOM    = 31
A_TARGET = 31
B_TARGET = 25
K_TRUE   = 7

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

# Load and build C[a,b]
with open(FILE_PATH, "r") as f:
    data = json.load(f)

counts_raw = data["counts"]

C = np.zeros((ORDER, ORDER), dtype=np.int64)
for bitstring, v in counts_raw.items():
    a_val = bits_to_int(bitstring[N_Q:2*N_Q])
    b_val = bits_to_int(bitstring[:N_Q])
    C[a_val, b_val] += int(v)

shots = int(C.sum())
print(f"Shots: {shots}")

Pab = C / shots
Pa = Pab.sum(axis=1, keepdims=True)   
Pb = Pab.sum(axis=0, keepdims=True)   

# Heatmap
plt.figure()
plt.imshow(np.log10(C + 1), origin="lower")
plt.title("Phase-Only Output Heatmap: log10(C[a,b]+1)")
plt.xlabel("b")
plt.ylabel("a")
plt.colorbar(label="log10(count+1)")
plt.show()

# Marginal P(a)
marg_a = C.sum(axis=1)
plt.figure()
plt.bar(range(ORDER), marg_a)
plt.title("Marginal P(a) from Phase-Only Run")
plt.xlabel("a")
plt.ylabel("counts")
plt.show()

# Conditional slice at a=31
row = C[A_DOM, :]
row_total = int(row.sum())

if row_total > 0:
    plt.figure()
    plt.bar(range(ORDER), row)
    plt.title(f"Conditional Slice: P(b | a={A_DOM}) (counts within a={A_DOM})")
    plt.xlabel("b")
    plt.ylabel(f"counts (within a={A_DOM})")
    plt.show()

# Key spectrum (phase-only)
k_spec = np.zeros(ORDER, dtype=np.int64)
if row_total > 0:
    for b in range(ORDER):
        k_hat = (-b) % ORDER
        k_spec[k_hat] += int(row[b])

plt.figure()
plt.bar(range(ORDER), k_spec)
plt.title(f"Key Spectrum (Phase-Only): k̂ = (-b) mod 32, conditioned on a={A_DOM}")
plt.xlabel("k̂")
plt.ylabel(f"counts (within a={A_DOM})")
plt.show()

# Data Prints 
top_idx = np.unravel_index(np.argmax(C), C.shape)
a_top, b_top = top_idx
c_top = C[a_top, b_top]

print(f"Top (a,b): ({a_top}, {b_top}) count: {c_top}")

p_a_dom = marg_a[A_DOM] / shots
p_b_top_given_a = C[A_DOM, b_top] / row_total

print(f"P(a={A_DOM}) = {p_a_dom:.4f}")
print(f"P(b={b_top} | a={A_DOM}) = {p_b_top_given_a:.4f}")

best_k = int(np.argmax(k_spec))
print(f"Best k̂ by peak (phase-only rule) = {best_k}  count: {int(k_spec[best_k])}")

# Lorenz and Gini
p_flat = Pab.flatten()
p_sorted = np.sort(p_flat)
cum_p = np.cumsum(p_sorted)
n = len(p_sorted)
x = np.arange(1, n + 1) / n

lorenz_x = np.concatenate(([0.0], x))
lorenz_y = np.concatenate(([0.0], cum_p))

area = np.trapezoid(lorenz_y, lorenz_x)
gini = 1.0 - 2.0 * area

plt.figure()
plt.plot(lorenz_x, lorenz_y)
plt.plot([0, 1], [0, 1])
plt.title(f"Lorenz Curve of P(a,b) (Gini ≈ {gini:.4f})")
plt.xlabel("Fraction of outcomes (sorted)")
plt.ylabel("Cumulative probability mass")
plt.show()

print(f"Gini concentration index ≈ {gini:.6f}  (0=uniform, 1=all mass in one cell)")

# PMI
eps = 1e-15
den = np.maximum(Pa @ Pb, eps)
ratio = np.where(Pab > 0, Pab / den, np.nan)
PMI = np.log2(ratio)

plt.figure()
plt.imshow(PMI, origin="lower", aspect="auto")
plt.title("PMI Heatmap: log2(P(a,b)/(P(a)P(b)))")
plt.xlabel("b")
plt.ylabel("a")
plt.colorbar(label="PMI (bits)")
plt.show()

mask = Pab > 0
MI = float(np.sum(Pab[mask] * np.log2(np.maximum(Pab[mask] / den[mask], eps))))
print(f"Mutual Information MI(A;B) ≈ {MI:.6f} bits")

# Key Projection Map
Kmat = np.zeros((ORDER, ORDER), dtype=np.int64)
for a in range(ORDER):
    for b in range(ORDER):
        k_hat = (-b) % ORDER
        Kmat[a, k_hat] += C[a, b]

plt.figure()
plt.imshow(np.log10(Kmat + 1), origin="lower", aspect="auto")
plt.title("Key-Projection Map: log10(Kmat[a,k̂]+1), where k̂=(-b) mod 32")
plt.xlabel("k̂")
plt.ylabel("a")
plt.colorbar(label="log10(count+1)")
plt.show()

kcol = Kmat[:, K_TRUE]
print(f"Total mass at k̂={K_TRUE}: {int(kcol.sum())}  ({kcol.sum()/shots:.4f} of shots)")
print(f"Mass at (a=31, k̂=7): {int(Kmat[A_TARGET, K_TRUE])}")

# Toroidal Blur (ring mass)
def torus_dist(a, b, a0, b0, mod=32):
    da = min((a - a0) % mod, (a0 - a) % mod)
    db = min((b - b0) % mod, (b0 - b) % mod)
    return float(np.sqrt(da*da + db*db))

D = np.zeros((ORDER, ORDER))
for a in range(ORDER):
    for b in range(ORDER):
        D[a, b] = torus_dist(a, b, A_TARGET, B_TARGET, ORDER)

dvals = D.flatten()
pvals = Pab.flatten()

dmax = int(np.ceil(dvals.max()))
mass_by_r = np.zeros(dmax + 1)

for d, p in zip(dvals, pvals):
    r = int(np.floor(d + 1e-12))
    mass_by_r[r] += p

plt.figure()
plt.bar(np.arange(len(mass_by_r)), mass_by_r)
plt.title("Toroidal Blur Spectrum Around Ideal Peak (31,25)")
plt.xlabel("Radius bucket r = floor(distance)")
plt.ylabel("Probability mass in ring")
plt.show()

print(f"P(ideal peak a=31,b=25) = {C[A_TARGET,B_TARGET]/shots:.6f}  (count={int(C[A_TARGET,B_TARGET])})")
print(f"Mass within r<=1: {mass_by_r[:2].sum():.6f}")
print(f"Mass within r<=2: {mass_by_r[:3].sum():.6f}")

# Per-row Key Fidelity
fidelity = np.zeros(ORDER)
row_mass = C.sum(axis=1)

for a in range(ORDER):
    if row_mass[a] > 0:
        fidelity[a] = C[a, (-K_TRUE) % ORDER] / row_mass[a]

plt.figure()
plt.bar(range(ORDER), fidelity)
plt.title("Per-row Key Fidelity: P(k̂=7 | a) = C[a,25] / sum_b C[a,b]")
plt.xlabel("a")
plt.ylabel("P(k̂=7 | a)")
plt.show()

print(f"P(k̂=7 | a=31) = {fidelity[A_TARGET]:.6f}")

# SNR curve
dmax = int(np.ceil(max(
    torus_dist(a, b, A_TARGET, B_TARGET, ORDER)
    for a in range(ORDER) for b in range(ORDER)
)))

signal = np.zeros(dmax + 1, dtype=np.float64)

for r in range(dmax + 1):
    s = 0.0
    for a in range(ORDER):
        for b in range(ORDER):
            if torus_dist(a, b, A_TARGET, B_TARGET, ORDER) <= r + 1e-12:
                s += Pab[a, b]
    signal[r] = s

noise = 1.0 - signal
snr = np.where(noise > 0, signal / noise, np.nan)

plt.figure()
plt.plot(range(dmax + 1), snr, marker="o")
plt.title("Peak-vs-Halo SNR on Torus: SNR(r) = mass(d<=r) / mass(d>r)")
plt.xlabel("radius r")
plt.ylabel("SNR(r)")
plt.show()

print(f"SNR(r=0) = {snr[0]:.6f}  (pure peak vs rest)")
print(f"SNR(r=1) = {snr[1]:.6f}")
print(f"SNR(r=2) = {snr[2]:.6f}")

# Edge Map
def wrap(i, mod=32):
    return i % mod

G = np.zeros_like(C)
for a in range(ORDER):
    for b in range(ORDER):
        ap = wrap(a + 1)
        am = wrap(a - 1)
        bp = wrap(b + 1)
        bm = wrap(b - 1)
        G[a, b] = abs(int(C[ap, b]) - int(C[am, b])) + \
                  abs(int(C[a, bp]) - int(C[a, bm]))

plt.figure()
plt.imshow(np.log10(G + 1), origin="lower")
plt.title("Interference Edge Map: log10(|Δa|+|Δb| + 1) (noise-sharpness proxy)")
plt.xlabel("b")
plt.ylabel("a")
plt.colorbar(label="log10(edge+1)")
plt.show()

# Hamming Diffusion
a_star, b_star = np.unravel_index(np.argmax(C), C.shape)
print(f"Dominant peak (a*,b*) = ({a_star},{b_star}) count={int(C[a_star,b_star])}")

def popcount(x):
    return int(bin(x).count("1"))

hd_mass = np.zeros(2 * N_Q + 1)

for a in range(ORDER):
    for b in range(ORDER):
        hd = popcount(a ^ a_star) + popcount(b ^ b_star)
        hd_mass[hd] += Pab[a, b]

plt.figure()
plt.bar(range(len(hd_mass)), hd_mass)
plt.title("Hamming-Diffusion Spectrum around Dominant Peak (a*,b*)")
plt.xlabel("Hamming distance hd = wt(a⊕a*) + wt(b⊕b*)")
plt.ylabel("Probability mass")
plt.show()

print(f"Mass at hd=0 (exact peak) = {hd_mass[0]:.6f}")
print(f"Mass at hd<=1 = {hd_mass[:2].sum():.6f}")
print(f"Mass at hd<=2 = {hd_mass[:3].sum():.6f}")

# End