1. Group Encoding
1. Group Encoding
Restrict attention to the order‑8 subgroup ⟨𝑃⟩ of an elliptic curve over 𝐹_11.
Map points to integers: 0𝑃 -> 0, 1𝑃 -> 1, …, 7𝑃 -> 7.
Group law becomes modular addition: (𝑥𝑃) + (𝑦𝑃) = ((𝑥 + 𝑦) mod 8))𝑃.
2. Quantum Registers
Register a: three qubits for the exponent 𝑎 ∈ {0, …, 7}.
Register b: three qubits for 𝑏 ∈ {0, …, 7}.
Register p: three qubits initialized to ∣0⟩ to hold the point index.
Classicals 𝑐_0 - 𝑐_5: record the measured values of a and b.
3. Superposition Preparation
Apply Hadamards to every qubit in a and b:
7
1/8 ∑ ∣a⟩_a ∣b⟩_b ∣0⟩_p
a, b=0
4. Oracle construction U_f
Goal: reversible map: ∣a⟩ ∣b⟩ ∣0⟩ -> ∣a⟩ ∣b⟩ ∣aP + bQ⟩.
Because indices add modulo 8, realize each controlled addition with a 3‑qubit permutation gate +cmod8.
Add aP: for each bit a_i (weight 2^𝑖), add (2^𝑖 mod 8)
Add bQ: compute 2^𝑖 𝑄 once classically ≡ (2^𝑖(7)) mod 8 and add it controlled on 𝑏_𝑖.
No gate ever references the secret k.
5. Global State after Oracle
1/8 ∑ ∣a⟩∣b⟩∣f(a, b)⟩, f(a, b) = a + kb (mod8).
a, b
6. Discard Point Register
The algorithm needs only the phase relation in 𝑎, 𝑏. A barrier isolates p.
7. Quantum Fourier Transform (QFT)
QFT 7
∣a⟩ -> 1/√8 ∑ e^(2πiau/8) ∣u⟩,
u=0
QFT 7
∣b⟩ -> 1/√8 ∑ e^(2πibv/8) ∣v⟩.
v=0
8. Interference Pattern
The joint amplitude for observing (u,v) is:
1/8√8 ∑ e^((2πi/8) (au + bv)) δ_(a + kb ≡ 0) = 1/√8 δ_(u + kv ≡ 0 (mod 8)),
which forms a diagonal “ridge” in the 8×8 outcome grid.
9. Measurement
Measure all six qubits. Outcomes concentrate on the eight pairs obeying:
u + kv ≡ 0 (mod 8).
10. Classical Post-Processing
Translate bitstrings (u,v) to integers (a,b) (due to endian reversal).
Keep rows with gcd(b, 8) = 1 so that b is invertible.
Recover the secret scalar via 𝑘 = (-𝑎)𝑏^-1 (mod 8)
With 8192 shots the most frequent valid row yields k = 7.
11. Verification and Storage
Print the result, looking for k = 7.
Save raw counts, qubit layout, and job ID to JSON for independent auditing.
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
# User settings
TOKEN = "IBMQ_KEY_O-`"
INSTANCE = "IBMQ_CRN"
BACKEND = "ibm_torino"
# Most recent calibration data
CAL_CSV = "/Users/steventippeconnic/Downloads/ibm_torino_calibrations_2025-05-02T21_47_05Z.csv"
SHOTS = 8192
# Toy‑curve parameters (order‑8 subgroup of E(F₁₁))
ORDER = 8 # |E(F_11)| = 8
P_IDX = 1 # generator P → index 1
Q_IDX = 7 # public point Q = k·P, here “7 mod 8”
# Only Q is needed to build the oracle
# 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 = 3 + 3 + 3 # a, b, point registers
PHYSICAL = best_qubits(CAL_CSV, N_Q_TOTAL)
# Constant‑adder modulo 8 as a reusable gate
def add_const_mod8_gate(c: int) -> UnitaryGate:
"""Return a 3‑qubit gate that maps |x⟩ ↦ |x+c (mod 8)⟩."""
mat = np.zeros((8, 8))
for x in range(8):
mat[(x + c) % 8, x] = 1
return UnitaryGate(mat, label=f"+{c}")
ADDERS = {c: add_const_mod8_gate(c) for c in range(1, 8)}
def controlled_add(qc: QuantumCircuit, ctrl_qubit, point_reg, constant):
"""Apply |x⟩ → |x+constant (mod 8)⟩ 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 8)
def ecdlp_oracle(qc, a_reg, b_reg, point_reg):
# Add a·P (P index = 1)
for i in range(3): # bit 2^i of a
constant = (P_IDX * (1 << i)) % ORDER
if constant:
controlled_add(qc, a_reg[i], point_reg, constant)
# Add b·Q (Q index known *publicly*)
for i in range(3):
constant = (Q_IDX * (1 << i)) % ORDER # 2^i Q mod 8
if constant:
controlled_add(qc, b_reg[i], point_reg, constant)
# Build the full Shor circuit
def shor_ecdlp_circuit() -> QuantumCircuit:
a = QuantumRegister(3, "a")
b = QuantumRegister(3, "b")
p = QuantumRegister(3, "p") # point index, init |0⟩
c = ClassicalRegister(6, "c") # will hold (a,b)
qc = QuantumCircuit(a, b, p, c, name="ECDLP_8pts")
# Prepare uniform superposition over (a,b)
qc.h(a); qc.h(b)
# Oracle U_f : adds aP + bQ
ecdlp_oracle(qc, a, b, p)
qc.barrier()
# Apply QFT on a and b
qc.append(QFT(3, do_swaps=False), a)
qc.append(QFT(3, do_swaps=False), b)
# Measure result registers
qc.measure(a, c[:3])
qc.measure(b, c[3:])
return qc
# IBMQ 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 of results
creg_name = trans.cregs[0].name
counts_raw = result[0].data.__getattribute__(creg_name).get_counts()
# Helper to interpret bitstrings (keep endian flip of Qiskit)
def bits_to_int(bs): return int(bs[::-1], 2)
counts = {(bits_to_int(k[3:]), bits_to_int(k[:3])): v
for k, v in counts_raw.items()}
top = sorted(counts.items(), key=lambda kv: kv[1], reverse=True)
k_found = None
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
log .info("Candidate k = %d from (a=%d, b=%d, count=%d)",
k_candidate, a_val, b_val, freq)
k_found = k_candidate
break
if k_found is not None:
print(f"\nSUCCESS — recovered k = {k_found}")
else:
print("\nNo invertible (a,b) pair found; try more shots.")
# Save raw data
out = {
"experiment": "ECDLP_8pts_Shors",
"backend": backend .name,
"physical_qubits": PHYSICAL,
"shots": SHOTS,
"counts": counts_raw
}
JSON_PATH = "/Users/steventippeconnic/Documents/ECDLP_8pts_Shors_Run_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, numpy as np, matplotlib.pyplot as plt
from pathlib import Path
from math import gcd, log2
# Ensure all legends/labels are fully visible
plt.rcParams.update({'figure.constrained_layout.use': True})
# Load & decode
F = Path('/Users/steventippeconnic/Documents/ECDLP_8pts_Shors_Run_0.json')
raw = json.loads(F.read_text())["counts"]
def bits_to_int(bs): return int(bs[::-1], 2) # endian flip
pairs = np.zeros((8, 8), dtype=int)
for bits, cnt in raw.items():
u = bits_to_int(bits[3:])
v = bits_to_int(bits[:3])
pairs[u, v] += cnt
total = pairs.sum()
probs = pairs / total
diag_id = np.eye(8) * 0.125
# Joint distribution heat‑map
plt.figure(figsize=(5,4))
plt.imshow(probs, cmap='viridis', vmin=0, vmax=probs.max())
plt.colorbar(label='P(u,v)')
plt.xticks(range(8)); plt.yticks(range(8))
plt.xlabel('v'); plt.ylabel('u')
plt.title('Joint QFT outcome $P(u,v)$')
plt.show()
# Histogram of (u−v) mod 8 (ridge test)
residue = np.zeros(8, dtype=int)
for u in range(8):
for v in range(8):
residue[(u-v) & 7] += pairs[u,v]
plt.figure(figsize=(6,3))
plt.bar(range(8), residue/total)
plt.xlabel('$(u-v)\\;\\text{mod}\\;8$')
plt.ylabel('Probability')
plt.title('Collapse onto $u\\equiv v\\;(\\text{mod}\\;8)$ ridge')
plt.show()
# Scatter of invertible rows & recovered k
us, vs, ks, ws = [], [], [], []
for u in range(8):
for v in range(8):
if gcd(v,8)==1 and pairs[u,v]:
ks.append((-u * pow(v,-1,8)) & 7)
us.append(u); vs.append(v); ws.append(pairs[u,v])
sizes = (np.array(ws)/max(ws))*600
plt.figure(figsize=(5,4))
sc = plt.scatter(us, vs, s=sizes, c=ks, cmap='twilight', edgecolors='k')
plt.colorbar(sc, label='$k$ candidate')
plt.xticks(range(8)); plt.yticks(range(8))
plt.xlabel('u'); plt.ylabel('v')
plt.title('Key candidates from invertible rows')
plt.show()
# Diagonal profile
diag_exp = np.array([pairs[d,d] for d in range(8)]) / total
diag_ideal = np.full(8, 1/8)
plt.figure(figsize=(6,3))
plt.plot(range(8), diag_exp, marker='o', label='experiment')
plt.plot(range(8), diag_ideal, linestyle='--', label='ideal')
plt.xlabel('d (where u=v=d)'); plt.ylabel('Diagonal probability')
plt.title('Diagonal strength & hardware bias')
plt.legend()
plt.show()
# Difference heat‑map
plt.figure(figsize=(5,4))
diff = probs - diag_id
mx = np.abs(diff).max()
plt.imshow(diff, cmap='coolwarm', vmin=-mx, vmax=mx)
plt.colorbar(label='$P_{\\text{obs}}-P_{\\text{ideal}}$')
plt.xticks(range(8)); plt.yticks(range(8))
plt.xlabel('v'); plt.ylabel('u')
plt.title('Deviation from ideal ridge')
plt.show()
# Marginal distributions of u and v
Pu = probs.sum(axis=1); Pv = probs.sum(axis=0)
x = np.arange(8)
plt.figure(figsize=(6,3))
plt.bar(x-0.15, Pu, width=0.3, label='$P(u)$')
plt.bar(x+0.15, Pv, width=0.3, label='$P(v)$')
plt.xticks(range(8))
plt.ylabel('Probability')
plt.title('Marginal distributions of $u$ and $v$')
plt.legend()
plt.show()
# Row‑wise signal‑to‑leak ratio
ratios = []
for d in range(8):
row_total = probs[d,:].sum()
off_diag_av = (row_total - probs[d,d]) / 7
ratios.append(probs[d,d] / off_diag_av if off_diag_av else 0)
plt.figure(figsize=(6,3))
plt.bar(range(8), ratios)
plt.xlabel('d (where u=v=d)')
plt.ylabel('Signal / row‑leak')
plt.title('Row‑wise ridge prominence')
plt.show()
# Histogram of recovered k values
k_counts = np.zeros(8, dtype=int)
for u in range(8):
for v in range(8):
if gcd(v,8)==1 and pairs[u,v]:
k = (-u * pow(v,-1,8)) & 7
k_counts[k] += pairs[u,v]
plt.figure(figsize=(6,3))
plt.bar(range(8), k_counts / k_counts.sum())
plt.xlabel('Recovered $k$ candidate')
plt.ylabel('Probability mass')
plt.title('Selectivity for the true key')
plt.show()
# Column‑normalised conditional P(u|v)
cond = np.where(pairs.sum(axis=0), probs / probs.sum(axis=0, keepdims=True), 0)
plt.figure(figsize=(5,4))
plt.imshow(cond, cmap='viridis')
plt.colorbar(label='$P(u\\mid v)$')
plt.xticks(range(8)); plt.yticks(range(8))
plt.xlabel('v'); plt.ylabel('u')
plt.title('Column‑normalised conditional $P(u\\mid v)$')
plt.show()
# Column‑wise ridge prominence
col_ratios = []
for d in range(8):
col_total = probs[:,d].sum()
off_diag_av = (col_total - probs[d,d]) / 7
col_ratios.append(probs[d,d] / off_diag_av if off_diag_av else 0)
plt.figure(figsize=(6,3))
plt.bar(range(8), col_ratios)
plt.xlabel('d (where u=v=d)')
plt.ylabel('Signal / column‑leak')
plt.title('Column‑wise ridge prominence')
plt.show()
# Probability vs |u-v| distance
dist = np.zeros(8, dtype=float)
for u in range(8):
for v in range(8):
dist[abs(u-v)] += probs[u,v]
plt.figure(figsize=(6,3))
plt.bar(range(8), dist)
plt.xlabel('$|u-v|$')
plt.ylabel('Probability')
plt.title('Leakage as function of distance from ridge')
plt.show()
# Row‑wise KL divergence from ideal uniform
kl = []
ideal_row = np.full(8, 1/8)
for u in range(8):
p_row = probs[u,:]
# Avoid log(0); treat zero cells as 0 contribution
kl_row = sum(p * log2(p/ideal_row[v]) for v,p in enumerate(p_row) if p)
kl.append(kl_row)
plt.figure(figsize=(6,3))
plt.bar(range(8), kl)
plt.xlabel('u (row index)')
plt.ylabel('$D_{\\mathrm{KL}}(P_{u}\\parallel U)$ [bits]')
plt.title('Row‑wise KL divergence from ideal')
plt.show()
# End
/////////////////////////////////////////////////////////////////////////////////
# Code to print top 50 pairs.
# Imports
import json
from pathlib import Path
FILE = Path("/Users/steventippeconnic/Documents/ECDLP_8pts_Shors_Run_0.json")
# Little‑endian helper
def bits_to_int(bitstr: str) -> int:
"""Convert a 3‑bit string (MSB left) to an int, reversing endianness."""
return int(bitstr[::-1], 2)
# Load and decode
with FILE.open() as fp:
counts_raw = json.load(fp)["counts"]
pairs = []
for bitstr, freq in counts_raw.items():
a = bits_to_int(bitstr[3:]) # last 3 bits -> a
b = bits_to_int(bitstr[:3]) # first 3 bits -> b
pairs.append(((a, b), freq))
# Sort by descending frequency and keep the top 50
top50 = sorted(pairs, key=lambda kv: kv[1], reverse=True)[:50]
# Print
print("rank\t(a,b)\tcount")
for idx, ((a, b), freq) in enumerate(top50, 1):
print(f"{idx:2d}\t({a},{b})\t{freq}")
# End