Causal Quantum Timekeeping with Chained Bloch Clocks while Breaking a 3-Bit Elliptic Curve Key via Shor’s Algorithm on a 133-Qubit Quantum Computer

Code:

# Imports
import logging, json
from math import gcd
import numpy as np
import pandas as pd
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, plot_bloch_multivector
from qiskit.quantum_info import Statevector
import matplotlib.pyplot as plt

# Settings
TOKEN = "YOUR_IBMQ_API_KEY"
INSTANCE = "YOUR_IBMQ_CRN"
BACKEND  = "ibm_torino"
CAL_CSV = "/Users/steventippeconnic/Downloads/ibm_torino_calibrations_2025-07-01T20_15_50Z.csv"
SHOTS = 8192
ORDER = 8
P_IDX = 1
Q_IDX = 7
TARGET_K = 7

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

# Calibration Based Qubit Selection
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
PHYSICAL = best_qubits(CAL_CSV, N_Q_TOTAL)

# Adders
def add_const_mod8_gate(c: int) -> UnitaryGate:
    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, ctrl_qubit, point_reg, constant):
    qc.append(ADDERS[constant].control(), [ctrl_qubit, *point_reg])

def ecdlp_oracle(qc, a_reg, b_reg, point_reg):
    for i in range(3):
        constant = (P_IDX * (1 << i)) % ORDER
        if constant:
            controlled_add(qc, a_reg[i], point_reg, constant)
    for i in range(3):
        constant = (Q_IDX * (1 << i)) % ORDER
        if constant:
            controlled_add(qc, b_reg[i], point_reg, constant)

# Circuit
def shor_ecdlp_with_chained_clocks():
    clk2 = QuantumRegister(1, "clk2")
    clk1 = QuantumRegister(1, "clk1")
    a = QuantumRegister(3, "a")
    b = QuantumRegister(3, "b")
    p = QuantumRegister(3, "p")
    c = ClassicalRegister(6, "c")
    clk1_c = ClassicalRegister(1, "c_clk1")
    clk2_c = ClassicalRegister(1, "c_clk2")
    qc = QuantumCircuit(clk2, clk1, a, b, p, c, clk1_c, clk2_c)

    # Clock 1 evolution
    qc.h(clk1[0])
    qc.rx(np.pi / 4, clk1[0])
    qc.ry(np.pi / 4, clk1[0])
    qc.rz(np.pi / 4, clk1[0])

    # ECC logic
    qc.h(a)
    qc.h(b)
    ecdlp_oracle(qc, a, b, p)
    qc.barrier()

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

    # ECC measurements
    qc.measure(a, c[:3])
    qc.measure(b, c[3:])
    qc.measure(clk1[0], clk1_c[0])

    # Clock 2 evolution
    qc.h(clk2[0])
    qc.rx(np.pi / 4, clk2[0])
    qc.ry(np.pi / 4, clk2[0])
    qc.rz(np.pi / 4, clk2[0])

    # Entangle clocks: If clk1 = 1, flip clk2
    qc .cx(clk1[0], clk2[0])
    qc.measure(clk2[0], clk2_c[0])
    return qc

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

qc = shor_ecdlp_with_chained_clocks()

# Select physical qubits for clk1 and clk2 that do not overlap with ECC qubits
PHYSICAL_CLOCKS = [q for q in range(127) if q not in PHYSICAL][:2]
FULL_LAYOUT = PHYSICAL + PHYSICAL_CLOCKS

trans = transpile(qc, backend=backend, initial_layout=FULL_LAYOUT, 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()

# Extract counts
cregs = {creg .name: creg.size for creg in trans.cregs}
clk1_name = [name for name in cregs if "clk1" in name][0]
clk2_name = [name for name in cregs if "clk2" in name][0]
main_name = [name for name in cregs if name == "c"][0]

clk1_counts = result[0].data.__getattribute__(clk1_name).get_counts()
clk2_counts = result[0].data.__getattribute__(clk2_name).get_counts()
main_counts_raw = result[0].data.__getattribute__(main_name).get_counts()

# ECC Recovery
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 main_counts_raw.items()}
top = sorted(counts.items(), key=lambda kv: kv[1], reverse=True)

k_found = None
for i, ((a_val, b_val), freq) in enumerate(top[:10]):
    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)
    if k_candidate == TARGET_K:
        k_found = k_candidate
        print(f"\nSUCCESS — recovered correct k = {k_candidate} from top-{i+1} (a={a_val}, b={b_val}, count={freq})")
        break

if k_found is None:
    print("\nNo correct key found in top 10; try more shots or examine further.")

# Save Results
out = {
    "experiment": "ECDLP_8pts_Shors_BlochClock_CausalChain",
    "backend": backend .name,
    "physical_qubits": FULL_LAYOUT,
    "shots": SHOTS,
    "counts": main_counts_raw,
    "clock1_counts": clk1_counts,
    "clock2_counts": clk2_counts
}
JSON_PATH = "/Users/steventippeconnic/Documents/QC/Shors_ECC_3_Bit_Key_Two_Bloch_Clocks_0.json"
with open(JSON_PATH, "w") as fp:
    json.dump(out, fp, indent=4)
log .info("Results saved → %s", JSON_PATH)

# Visualize ideal state of Clock 2
qc_bloch2 = QuantumCircuit(1)
qc_bloch2.h(0)
qc_bloch2.rx(np.pi / 4, 0)
qc_bloch2.ry(np.pi / 4, 0)
qc_bloch2.rz(np.pi / 4, 0)
state2 = Statevector.from_instruction(qc_bloch2)
plot_bloch_multivector(state2)
plt .show()

# End 




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

          
# Code for all visuals from experiment JSON

# Imports
import json
import matplotlib.pyplot as plt
import numpy as np
from collections import defaultdict
from collections import Counter
from itertools import product
from mpl_toolkits.mplot3d import Axes3D
from collections import defaultdict
from itertools import product

# Load results
with open('/Users/steventippeconnic/Documents/QC/Shors_ECC_3_Bit_Key_Two_Bloch_Clocks_0.json') as f:
    data = json.load(f)

counts = data["counts"]
shots = sum(counts.values())
clk1_counts = data["clock1_counts"]
clk2_counts = data["clock2_counts"]

# Clock Collapse Comparison
plt.figure()
plt.bar(["Clock 1 = 0", "Clock 1 = 1"], [clk1_counts["0"], clk1_counts["1"]], color="blue", alpha=0.6, label="Clock 1")
plt.bar(["Clock 2 = 0", "Clock 2 = 1"], [clk2_counts["0"], clk2_counts["1"]], color="green", alpha=0.6, label="Clock 2")
plt.title("Collapse Probabilities of Clock 1 and Clock 2")
plt.ylabel("Counts")
plt.legend()
plt.show()

# ECC Bitstring Histogram
plt.figure()
sorted_items = sorted(counts.items(), key=lambda x: x[1], reverse=True)
bitstrings = [k for k, _ in sorted_items]
frequencies = [v for _, v in sorted_items]
plt.bar(range(len(bitstrings)), frequencies)
plt.xticks(range(len(bitstrings)), bitstrings, rotation='vertical', fontsize=6)
plt.title("ECC Bitstring Frequencies (a ⊕ b)")
plt.ylabel("Counts")
plt.tight_layout()
plt.show()

# (a, b) Heatmap
ab_matrix = np.zeros((8, 8))
for bitstring, count in counts.items():
    a = int(bitstring[3:], 2)
    b = int(bitstring[:3], 2)
    ab_matrix[a, b] += count
plt.figure()
plt.imshow(ab_matrix, cmap='viridis', interpolation='nearest')
plt.colorbar(label="Counts")
plt.title("(a, b) Frequency Heatmap")
plt.xlabel("b")
plt.ylabel("a")
plt.xticks(range(8))
plt.yticks(range(8))
plt.show()

# Ridge Alignment for k = 7 (u + 7v ≡ 0 mod 8)
ridge_matrix = np.zeros((8, 8))
for bitstring, count in counts.items():
    a = int(bitstring[3:], 2)
    b = int(bitstring[:3], 2)
    try:
        inv_b = pow(b, -1, 8)
        k = (-a * inv_b) % 8
        if k == 7:
            ridge_matrix[a, b] += count
    except ValueError:
        continue
plt.figure()
plt.imshow(ridge_matrix, cmap='plasma', interpolation='nearest')
plt.colorbar(label="Ridge-Aligned Counts for k = 7")
plt.title("Interference Ridge Map (k = 7)")
plt.xlabel("b")
plt.ylabel("a")
plt.xticks(range(8))
plt.yticks(range(8))
plt.show()

# Clock Mutual Collapse Matrix
joint = np.zeros((2, 2))  # rows clk1, cols clk2
for bitstring, count in counts.items():
    clk1_bit = bitstring[-1]  # last two bits are clock measurements
    clk2_bit = bitstring[-2]
    i = int(clk1_bit)
    j = int(clk2_bit)
    joint[i, j] += count
plt.figure()
plt.imshow(joint, cmap="cividis")
plt.colorbar(label="Counts")
plt.xticks([0, 1], ["clk2 = 0", "clk2 = 1"])
plt.yticks([0, 1], ["clk1 = 0", "clk1 = 1"])
plt.title("Joint Clock Collapse Correlation")
plt.xlabel("Clock 2")
plt.ylabel("Clock 1")
plt.show()

# Bitwise 1-Bias in ECC Register
bitwise = [0] * 6
total = [0] * 6
for bitstring, count in counts.items():
    for i in range(6):
        if bitstring[5 - i] == '1':
            bitwise[i] += count
        total[i] += count
bias = [bitwise[i] / total[i] for i in range(6)]
plt.figure()
plt.bar(range(6), bias)
plt.title("Bitwise 1-Bias in ECC Register")
plt.xlabel("Bit Index (0 = LSB of b)")
plt.ylabel("P(bit = 1)")
plt.ylim(0.4, 0.6)
plt.xticks(range(6), [f"b{i}" if i < 3 else f"a{i-3}" for i in range(6)])
plt.show()

# Number of Unique (a, b) Pairs Supporting Each Candidate k
from math import gcd
support = defaultdict(set)
for bitstring in counts:
    a = int(bitstring[3:], 2)
    b = int(bitstring[:3], 2)
    if gcd(b, 8) != 1:
        continue
    inv_b = pow(b, -1, 8)
    k = (-a * inv_b) % 8
    support[k].add((a, b))
ks = sorted(support)
num_supports = [len(support[k]) for k in ks]
plt.figure()
plt.bar(ks, num_supports)
plt.title("Support Count per Candidate k")
plt.xlabel("Candidate k")
plt.ylabel("Distinct (a, b) Pairs")
plt.xticks(range(8))
plt.show()

# ECC Output Histogram by Clock 1 = 1 and Clock 2 = 1
filtered_counts = defaultdict(int)
for bitstring, count in counts.items():
    if bitstring[-1] == '1' and bitstring[-2] == '1':
        filtered_counts[bitstring[:6]] += count
top = sorted(filtered_counts.items(), key=lambda x: x[1], reverse=True)
keys = [x[0] for x in top]
vals = [x[1] for x in top]
plt.figure()
plt.bar(range(len(keys)), vals)
plt.xticks(range(len(keys)), keys, rotation='vertical', fontsize=6)
plt.title("ECC Bitstring Histogram | clk1 = 1 and clk2 = 1")
plt.ylabel("Counts")
plt.tight_layout()
plt.show()

# Helper to decode bitstrings
def decode_bitstring(bitstring):
    b = int(bitstring[:3], 2)
    a = int(bitstring[3:6], 2)
    return a, b

decoded = [decode_bitstring(k) for k in counts.keys()]
frequencies = [v for v in counts.values()]
ab_counts = dict(zip(decoded, frequencies))

ridge_k = 7

# 3D Ridge Surface Plot 

X, Y = np.meshgrid(range(8), range(8))
Z = np.zeros_like(X, dtype=float)

for (a, b), count in ab_counts.items():
    Z[a, b] = count

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(X, Y, Z, cmap='viridis', edgecolor='k', alpha=0.9)
ax.set_title("3D Surface of (a, b) Frequencies")
ax.set_xlabel("b")
ax.set_ylabel("a")
ax.set_zlabel("Counts")
plt.tight_layout()
plt.show()

# Ridge Residue Intensity Map (Phase Residue) 

residue_matrix = np.zeros((8, 8))

for (a, b), count in ab_counts.items():
    residue = (a + ridge_k * b) % 8
    residue_matrix[a, b] = residue

plt.figure()
plt.imshow(residue_matrix, cmap='twilight', interpolation='nearest')
plt.title("Phase Residue Map (a + 7b mod 8)")
plt.xlabel("b")
plt.ylabel("a")
plt.colorbar(label="Residue ∈ [0,7]")
plt.show()

# Complex Projection on the Circle (Phasor Plot) 

phasor_sum = defaultdict(complex)
for (a, b), count in ab_counts.items():
    theta = 2 * np.pi * ((a + ridge_k * b) % 8) / 8
    phasor_sum[(a, b)] = count * np.exp(1j * theta)

phasor_total = sum(phasor_sum.values())

plt.figure(figsize=(6,6))
plt.polar([np.angle(v) for v in phasor_sum.values()],
          [np.abs(v) for v in phasor_sum.values()],
          'o', alpha=0.7)
plt.title("Complex Projection of (a + 7b) mod 8 onto Unit Circle")
plt.annotate("TPD",
             xy=(np.angle(phasor_total), np.abs(phasor_total)),
             xytext=(np.angle(phasor_total), np.abs(phasor_total)+20),
             arrowprops=dict(arrowstyle="->", color='red'))
plt.show()

# 3D Ridge Fidelity Hill 

fidelity_surface = np.zeros((8, 8))
for (a, b), count in ab_counts.items():
    theta = 2 * np.pi * ((a + ridge_k * b) % 8) / 8
    fidelity = np.cos(theta) ** 2
    fidelity_surface[a, b] = fidelity * count

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(X, Y, fidelity_surface, cmap='plasma', edgecolor='k')
ax.set_title("3D Ridge Fidelity Landscape (cos² weighted)")
ax.set_xlabel("b")
ax.set_ylabel("a")
ax.set_zlabel("Weighted Fidelity")
plt.tight_layout()
plt.show()

# End.


//////////////////////////////////////////////////////////
Code to analyze Bloch Clock 
# Imports
import base64, io, json, pathlib, zlib
from collections import Counter
import numpy as np
import matplotlib.pyplot as plt
from qiskit import QuantumCircuit
from qiskit.quantum_info import Statevector

# File location
RESULT_PATH = pathlib.Path("/Users/steventippeconnic/Documents/QC/job-d1i47se5no2c739tvkeg-result.json")

# Helpers
def decode_field(field_dict) -> np.ndarray:
    b64 = field_dict["__value__"]["array"]["__value__"]
    raw = zlib.decompress(base64.b64decode(b64))
    return np.load(io.BytesIO(raw), allow_pickle=False)

def unpack_bits(packed: np.ndarray) -> np.ndarray:
    bits = ((packed[:, None] >> np.arange(6)) & 1).astype(np.uint8)
    return bits

def ab_from_bits(bits: np.ndarray) -> tuple[int,int]:
    a = bits[0] | (bits[1] << 1) | (bits[2] << 2)
    b = bits[3] | (bits[4] << 1) | (bits[5] << 2)
    return a, b

# Load result JSON
with RESULT_PATH.open() as f:
    res_json = json.load(f)

fields = res_json["__value__"]["pub_results"][0]["__value__"]["data"]["__value__"]["fields"]

c_arr     = decode_field(fields["c"])
clk1_arr  = decode_field(fields["c_clk1"])
clk2_arr  = decode_field(fields["c_clk2"])

if c_arr.ndim == 2 and c_arr.shape[1] == 1:
    c_bits = unpack_bits(c_arr[:,0])
else:
    c_bits = c_arr

clk1_bits = clk1_arr.reshape(-1).astype(np.uint8)
clk2_bits = clk2_arr.reshape(-1).astype(np.uint8)
shots = c_bits.shape[0]

assert clk1_bits.size == shots and clk2_bits.size == shots, "Clock register size mismatch"

# Build conditionals
cnt_clk1_0, cnt_clk1_1 = Counter(), Counter()
cnt_clk2_0, cnt_clk2_1 = Counter(), Counter()

for bits, c1, c2 in zip(c_bits, clk1_bits, clk2_bits):
    a, b = ab_from_bits(bits)
    (cnt_clk1_1 if c1 else cnt_clk1_0)[(a,b)] += 1
    (cnt_clk2_1 if c2 else cnt_clk2_0)[(a,b)] += 1

tot1_0, tot1_1 = sum(cnt_clk1_0.values()), sum(cnt_clk1_1.values())
tot2_0, tot2_1 = sum(cnt_clk2_0.values()), sum(cnt_clk2_1.values())

# Report
def top(counter, tot, name):
    print(f"\n{name}  (shots = {tot})")
    for (a,b), c in counter.most_common(10):
        print(f"  (a={a}, b={b})   {c:>4}   {100*c/tot:5.2f}%")

print(f"\nTotal shots       : {shots}")
print(f"Clock1 |1⟩        : {tot1_1}  ({100*tot1_1/shots:.1f} %)")
print(f"Clock1 |0⟩        : {tot1_0}  ({100*tot1_0/shots:.1f} %)")
print(f"Clock2 |1⟩        : {tot2_1}  ({100*tot2_1/shots:.1f} %)")
print(f"Clock2 |0⟩        : {tot2_0}  ({100*tot2_0/shots:.1f} %)")

top(cnt_clk1_1, tot1_1, "Top-10  |  clk1 = 1")
top(cnt_clk1_0, tot1_0, "Top-10  |  clk1 = 0")
top(cnt_clk2_1, tot2_1, "Top-10  |  clk2 = 1")
top(cnt_clk2_0, tot2_0, "Top-10  |  clk2 = 0")

# Heatmap visualization
def heat(counter, title):
    grid = np.zeros((8,8), dtype=int)
    for (a,b), n in counter.items():
        grid[a,b] = n
    plt.figure(figsize=(4,4))
    plt.title(title)
    plt.imshow(grid, origin="lower")
    plt.xlabel("b")
    plt.ylabel("a")
    plt.colorbar(label="counts")
    plt.tight_layout()

heat(cnt_clk1_1, "Heat-map   clk1 = 1")
heat(cnt_clk1_0, "Heat-map   clk1 = 0")
heat(cnt_clk2_1, "Heat-map   clk2 = 1")
heat(cnt_clk2_0, "Heat-map   clk2 = 0")
plt.show()

# Ideal Bloch vector
qc = QuantumCircuit(1)
qc.h(0)
qc.rx(np.pi/4, 0)
qc.ry(np.pi/4, 0)
qc.rz(np.pi/4, 0)
bloch_sv = Statevector.from_instruction(qc)
bloch_xyz = (
    bloch_sv.expectation_value([[0,1],[1,0]]),
    bloch_sv.expectation_value([[0,-1j],[1j,0]]),
    bloch_sv.expectation_value([[1,0],[0,-1]])
)
print(f"\nIdeal Bloch-clock  ⟨Z⟩  =  {np.real(bloch_xyz[2]): .4f}")

# End