{ "nbformat": 4, "nbformat_minor": 0, "metadata": { "colab": { "provenance": [] }, "kernelspec": { "name": "python3", "display_name": "Python 3" }, "language_info": { "name": "python" } }, "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "colab": { "base_uri": "https://localhost:8080/" }, "id": "DProVrQc0SXZ", "outputId": "20c4b0f8-18a1-4ae7-a9cc-9856238dea70" }, "outputs": [ { "output_type": "stream", "name": "stdout", "text": [ "✅ All packages installed successfully.\n", "\n", "╔════════════════════════════════════════════════════════════════════╗\n", "║ SOSQ: Sequential Oracle SHOR Quantum ║\n", "║ Optimized Batch Execution ║\n", "║ N has 2048 bits ║\n", "║ mpmath: 2000 decimal digits ║\n", "╚════════════════════════════════════════════════════════════════════╝\n", "\n", "══════════════════════════════════════════════════════════════════════\n", "QUBIT ANALYSIS — N is 2048 bits, hardware has 133 qubits\n", "══════════════════════════════════════════════════════════════════════\n", "\n", " Standard Shor: ~8192 qubits ❌\n", " Beauregard: 4099 qubits ❌\n", " SOSQ (with target): 2049 qubits ❌\n", " ⭐ SOSQ (single-qubit): 1 qubit ✅ ALWAYS\n", "\n", " SOSQ single-qubit oracle calls: 4100\n", " Classical cost per call: pow(a, 2^k, N) with 2000-digit precision\n", "\n", " 📌 2048-bit N on 133-qubit HW:\n", " Only SOSQ single-qubit mode works!\n", " All modular arithmetic computed classically (mpmath)\n", " Quantum HW used only for phase estimation (1 qubit)\n", "\n", "══════════════════════════════════════════════════════════════════════\n", "⏱️ TIMING ESTIMATE FOR 2048-BIT N\n", "══════════════════════════════════════════════════════════════════════\n", "\n", " ORIGINAL (individual jobs):\n", " Oracle calls: 4100\n", " ~15s per call (API overhead)\n", " Total: 61500s = 17.1 hours\n", "\n", " ⭐ BATCHED (groups of 100):\n", " Batches: 41\n", " ~30s per batch\n", " Total: 1230s = 20.5 minutes\n", "\n", " Speedup: 50×\n", "\n", " Classical precomputation (pow(a, 2^k, N)):\n", " 4100 modular exponentiations on 2048-bit numbers\n", " Estimated: ~41s on modern CPU\n", " (Python's pow() is highly optimized for big ints)\n", "======================================================================\n", "Connecting to IBM Quantum...\n", "======================================================================\n" ] }, { "output_type": "stream", "name": "stderr", "text": [ "qiskit_runtime_service.__init__:WARNING:2026-02-21 20:29:52,233: Instance was not set at service instantiation. Free and trial plan instances will be prioritized. Based on the following filters: (tags: None, region: us-east, eu-de), and available plans: (open), the available account instances are: open-instance. If you need a specific instance set it explicitly either by using a saved account with a saved default instance or passing it in directly to QiskitRuntimeService().\n", "qiskit_runtime_service.backends:WARNING:2026-02-21 20:29:52,234: Loading instance: open-instance, plan: open\n" ] }, { "output_type": "stream", "name": "stdout", "text": [ "✅ Connected!\n", " • ibm_fez: 156q, ✅, 1196 pending\n", " • ibm_torino: 133q, ✅, 0 pending\n", " • ibm_marrakesh: 156q, ✅, 232 pending\n" ] }, { "output_type": "stream", "name": "stderr", "text": [ "qiskit_runtime_service.backends:WARNING:2026-02-21 20:29:57,848: Loading instance: open-instance, plan: open\n", "qiskit_runtime_service.backends:WARNING:2026-02-21 20:29:59,116: Using instance: open-instance, plan: open\n" ] }, { "output_type": "stream", "name": "stdout", "text": [ "\n", "🎯 Backend: ibm_torino (133 qubits)\n", "\n", "══════════════════════════════════════════════════════════════════════\n", "RUNNING SOSQ (BATCHED) — Backend: \n", "══════════════════════════════════════════════════════════════════════\n", "\n", " Chosen a = 14857595232396427516941849273762627309353275011035872028438979102221448018764811214335550378911385239042832330153307199521375474559482800190390034841829093941439342200536112506795917802201162230671321125520796718096251936434950013753389933883027560251634067600145782190986850566437489689333408216831503956985540835915302821684890835293802022321141410173728802736628089273623546592253626525766440256949467070575003413315362581658589537877649478256032259184993572524575144248057778898516584068997755681240114436398521558273040713840621495800996773201802044528803741752377100822525955071889707452400612656941927243629363\n", " Phase 1: Classical pre-computation...\n", " Precomputed 4100 b_k values in 79.9s\n", " Non-trivial calls needed: 4100/4100\n", " Phase 2: Building 4100 circuits...\n", " Phase 3: Submitting 4100 circuits in batches of 100...\n", " Batch 1/41 (100 circuits)...\n" ] } ], "source": [ "# =============================================================================\n", "# SOSQ: Sequential Oracle SHOR Quantum — Full Implementation for IBM Quantum\n", "# Real Backend via Qiskit Runtime (Single Colab Cell)\n", "#\n", "# Paper: \"Sequential Oracle SHOR Quantum (SOSQ): A Divide-and-Conquer\n", "# Turing Reduction Framework for Integer Factorization via\n", "# Polynomially Many Quantum Oracle Queries\"\n", "# Author: יהוה\n", "#\n", "# SINGLE-QUBIT MODE FOR 2048-BIT N:\n", "# Each SOSQ oracle call uses exactly 1 qubit. The modular exponentiation\n", "# a^(2^k) mod N is computed CLASSICALLY with mpmath arbitrary precision.\n", "# Only the resulting phase angle is applied as a gate on the quantum HW.\n", "# Circuit per call: |0> -- H -- P(phi_k - corr_k) -- H -- Measure\n", "# =============================================================================\n", "\n", "import subprocess, sys\n", "\n", "def install(pkg):\n", " subprocess.check_call([sys.executable, \"-m\", \"pip\", \"install\", pkg, \"-q\"])\n", "\n", "install(\"qiskit>=1.0\")\n", "install(\"qiskit-ibm-runtime>=0.20\")\n", "install(\"qiskit-aer\")\n", "install(\"pylatexenc\")\n", "install(\"mpmath\")\n", "\n", "print(\"✅ All packages installed successfully.\\n\")\n", "\n", "# ── Imports ───────────────────────────────────────────────────────────────────\n", "import mpmath\n", "import math\n", "import time\n", "import random\n", "import warnings\n", "from fractions import Fraction\n", "from typing import Optional, Tuple, List, Dict\n", "\n", "from qiskit import QuantumCircuit, QuantumRegister, ClassicalRegister, transpile\n", "from qiskit_ibm_runtime import QiskitRuntimeService, SamplerV2 as Sampler\n", "\n", "warnings.filterwarnings(\"ignore\")\n", "\n", "# ── Arbitrary precision: 2000 decimal digits (covers 2048-bit N) ─────────────\n", "mpmath.mp.dps = 2000\n", "\n", "# ═════════════════════════════════════════════════════════════════════════════\n", "# ██ CONFIGURATION ██\n", "# ═════════════════════════════════════════════════════════════════════════════\n", "\n", "\n", "IBM_TOKEN = \"\" # <-- PASTE YOUR TOKEN HERE\n", "\n", "# Number to factor\n", "N_TO_FACTOR = 20401308423288094242766662999493793187525420878197109253818355727956826504712629708789530022692909325272817125236977390052166930825056004028872464770803083249959025095910645332450783024034058817080357808464693098920369306132436800716130122636415040330812678918561101207970271295199148388284056864933825212541980308771716287643934228961790479259363124623312430892043194378423348328847933501339676855519527154465017072352029567538488597745797965852775230430579578162516042451489532165421970081208347386436907444609093030460498848125912815856038414252622303104009405061242804767128891944731455892392535411542580604029973\n", "\n", "USE_REAL_HARDWARE = True\n", "PREFERRED_BACKEND = None\n", "SHOTS = 4096\n", "VERBOSE = True\n", "\n", "# ═════════════════════════════════════════════════════════════════════════════\n", "# ██ MATH UTILITIES — ALL EXACT INTEGER OR MPMATH ██\n", "# ═════════════════════════════════════════════════════════════════════════════\n", "\n", "class MathUtils:\n", "\n", " @staticmethod\n", " def gcd(a: int, b: int) -> int:\n", " while b:\n", " a, b = b, a % b\n", " return a\n", "\n", " @staticmethod\n", " def is_prime(n: int) -> bool:\n", " if n < 2:\n", " return False\n", " if n < 4:\n", " return True\n", " if n % 2 == 0 or n % 3 == 0:\n", " return False\n", " witnesses = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37]\n", " d, s = n - 1, 0\n", " while d % 2 == 0:\n", " d //= 2\n", " s += 1\n", " for a in witnesses:\n", " if a >= n:\n", " continue\n", " x = pow(a, d, n)\n", " if x == 1 or x == n - 1:\n", " continue\n", " for _ in range(s - 1):\n", " x = pow(x, 2, n)\n", " if x == n - 1:\n", " break\n", " else:\n", " return False\n", " return True\n", "\n", " @staticmethod\n", " def is_perfect_power(n: int) -> Optional[Tuple[int, int]]:\n", " if n <= 1:\n", " return None\n", " max_b = int(mpmath.log(mpmath.mpf(n), 2)) + 1\n", " for b in range(2, max_b + 1):\n", " a_mp = mpmath.root(mpmath.mpf(n), b)\n", " a_int = int(mpmath.nint(a_mp))\n", " for c in [a_int - 1, a_int, a_int + 1]:\n", " if c >= 2 and pow(c, b) == n:\n", " return (c, b)\n", " return None\n", "\n", " @staticmethod\n", " def multiplicative_order(a: int, n: int) -> Optional[int]:\n", " if MathUtils.gcd(a, n) != 1:\n", " return None\n", " order, cur = 1, a % n\n", " while cur != 1:\n", " cur = (cur * a) % n\n", " order += 1\n", " if order > n:\n", " return None\n", " return order\n", "\n", " @staticmethod\n", " def continued_fractions(num: int, den: int,\n", " max_den: int) -> List[Fraction]:\n", " if den == 0:\n", " return []\n", " convergents = []\n", " n, d = num, den\n", " coeffs = []\n", " while d != 0 and len(coeffs) < 200:\n", " q, r = divmod(n, d)\n", " coeffs.append(q)\n", " n, d = d, r\n", " h_prev, h_curr = 0, 1\n", " k_prev, k_curr = 1, 0\n", " for c in coeffs:\n", " h_new = c * h_curr + h_prev\n", " k_new = c * k_curr + k_prev\n", " if k_new > max_den:\n", " break\n", " if k_new > 0:\n", " convergents.append(Fraction(h_new, k_new))\n", " h_prev, h_curr = h_curr, h_new\n", " k_prev, k_curr = k_curr, k_new\n", " return convergents\n", "\n", " @staticmethod\n", " def extract_period_candidates(measurement: int, Q: int, N: int,\n", " max_cand: int = 20) -> List[int]:\n", " if measurement == 0:\n", " return []\n", " candidates = []\n", " for conv in MathUtils.continued_fractions(measurement, Q, N):\n", " r_c = conv.denominator\n", " if 0 < r_c <= N:\n", " candidates.append(r_c)\n", " for m in range(2, min(10, N // r_c + 1)):\n", " if m * r_c <= N:\n", " candidates.append(m * r_c)\n", " return sorted(set(candidates))[:max_cand]\n", "\n", "\n", "# ═════════════════════════════════════════════════════════════════════════════\n", "# ██ SINGLE-QUBIT SOSQ CIRCUIT BUILDER ██\n", "# ═════════════════════════════════════════════════════════════════════════════\n", "\n", "class SOSQCircuitBuilder:\n", " \"\"\"\n", " Builds 1-qubit SOSQ oracle circuits.\n", "\n", " Each oracle call is: |0> -- H -- P(phi) -- H -- Measure\n", "\n", " phi = 2*pi * a^(2^k) mod N / N minus semiclassical correction\n", "\n", " The modular exponentiation a^(2^k) mod N is computed CLASSICALLY\n", " with mpmath arbitrary precision. Only the scalar angle phi goes\n", " to the quantum hardware as a single P gate on 1 qubit.\n", " \"\"\"\n", "\n", " @staticmethod\n", " def build_oracle(a: int, N: int, k: int,\n", " previous_bits: List[int]) -> Optional[Tuple[QuantumCircuit, int]]:\n", " \"\"\"\n", " Build the k-th SOSQ oracle (1 qubit).\n", "\n", " Returns (circuit, b_k) or (None, 1) if b_k == 1.\n", " \"\"\"\n", " # ── Classical computation: a^(2^k) mod N ──\n", " # Python pow(a, 2**k, N) handles arbitrarily large integers exactly\n", " b_k = pow(a, pow(2, k), N)\n", "\n", " if b_k == 1:\n", " return None, 1\n", "\n", " # ── Build 1-qubit circuit ──\n", " qr = QuantumRegister(1, 'q')\n", " cr = ClassicalRegister(1, 'c')\n", " qc = QuantumCircuit(qr, cr, name=f'SOSQ_k{k}')\n", "\n", " # Hadamard\n", " qc.h(qr[0])\n", "\n", " # ── Phase angle computed with arbitrary precision ──\n", " # phi_k = 2*pi * b_k / N\n", " two_pi = mpmath.mpf(2) * mpmath.pi\n", " phi_mp = two_pi * mpmath.mpf(b_k) / mpmath.mpf(N)\n", "\n", " # ── Semiclassical correction from previously measured bits ──\n", " correction_mp = mpmath.mpf(0)\n", " if previous_bits:\n", " for j, bit_val in enumerate(previous_bits):\n", " if bit_val:\n", " exp = k - j\n", " if exp > 0:\n", " correction_mp += mpmath.pi / mpmath.power(2, exp)\n", "\n", " # Net phase (convert to float only here, for the Qiskit gate)\n", " net_phase = float(phi_mp - correction_mp)\n", "\n", " qc.p(net_phase, qr[0])\n", "\n", " # Hadamard\n", " qc.h(qr[0])\n", "\n", " # Measure\n", " qc.measure(qr[0], cr[0])\n", "\n", " return qc, b_k\n", "\n", " @staticmethod\n", " def build_full_oracle(a: int, N: int, k: int, n_total: int,\n", " previous_bits: List[int],\n", " max_qubits: int = 133) -> Optional[Tuple[QuantumCircuit, int]]:\n", " \"\"\"\n", " Build oracle using target register if it fits, else single-qubit.\n", " \"\"\"\n", " n_target = N.bit_length()\n", " needed = 1 + n_target\n", "\n", " # If it doesn't fit in hardware, use single-qubit mode\n", " if needed > max_qubits:\n", " return SOSQCircuitBuilder.build_oracle(a, N, k, previous_bits)\n", "\n", " # ── Build with target register (for small N) ──\n", " b_k = pow(a, pow(2, k), N)\n", " if b_k == 1:\n", " return None, 1\n", "\n", " qr_a = QuantumRegister(1, 'anc')\n", " qr_t = QuantumRegister(n_target, 'tgt')\n", " cr = ClassicalRegister(1, 'c')\n", " qc = QuantumCircuit(qr_a, qr_t, cr, name=f'SOSQ_k{k}')\n", "\n", " qc.x(qr_t[0])\n", " qc.h(qr_a[0])\n", "\n", " # Semiclassical correction\n", " if previous_bits:\n", " correction_mp = mpmath.mpf(0)\n", " for j, bit_val in enumerate(previous_bits):\n", " if bit_val:\n", " exp = k - j\n", " if exp > 0:\n", " correction_mp += mpmath.pi / mpmath.power(2, exp)\n", " corr_float = float(correction_mp)\n", " if abs(corr_float) > 1e-15:\n", " qc.p(-corr_float, qr_a[0])\n", "\n", " # Controlled modular exponentiation\n", " if N <= 64:\n", " SOSQCircuitBuilder._apply_ctrl_ua(qc, qr_a[0], qr_t, b_k, N, n_target)\n", " else:\n", " two_pi = mpmath.mpf(2) * mpmath.pi\n", " for bit in range(n_target):\n", " ang = float(two_pi * mpmath.mpf(b_k) / mpmath.mpf(N) * mpmath.power(2, bit))\n", " qc.cp(ang, qr_a[0], qr_t[bit])\n", "\n", " qc.h(qr_a[0])\n", " qc.measure(qr_a[0], cr[0])\n", " return qc, b_k\n", "\n", " @staticmethod\n", " def _apply_ctrl_ua(qc, ctrl, tgt_reg, a_val, N, n_tgt):\n", " if a_val == 1 or N <= 1:\n", " return\n", " if N > 64:\n", " return\n", " perm = list(range(2**n_tgt))\n", " for x in range(N):\n", " perm[x] = (a_val * x) % N\n", " visited = [False] * len(perm)\n", " for i in range(len(perm)):\n", " if visited[i] or perm[i] == i:\n", " visited[i] = True\n", " continue\n", " cycle = []\n", " j = i\n", " while not visited[j]:\n", " visited[j] = True\n", " cycle.append(j)\n", " j = perm[j]\n", " for idx in range(len(cycle) - 1, 0, -1):\n", " diff = cycle[0] ^ cycle[idx]\n", " if diff == 0:\n", " continue\n", " for bit in range(n_tgt):\n", " if diff & (1 << bit):\n", " qc.cx(ctrl, list(tgt_reg)[bit])\n", " break\n", "\n", "\n", "# ═════════════════════════════════════════════════════════════════════════════\n", "# ██ SOSQ ENGINE ██\n", "# ═════════════════════════════════════════════════════════════════════════════\n", "\n", "# ═══════════════════════════════════════════════════════════════\n", "# VERSIÓN OPTIMIZADA: BATCH SUBMISSION\n", "# En vez de 4100 jobs individuales, envía en lotes de ~100\n", "# Reduce 17 horas → ~40 minutos\n", "# ═══════════════════════════════════════════════════════════════\n", "\n", "class SOSQEngine:\n", " \"\"\"\n", " SOSQ Engine with batched circuit submission.\n", "\n", " Instead of submitting 4100 individual jobs (each with API overhead),\n", " groups circuits into batches and submits them together.\n", "\n", " The semiclassical feedback means we can't batch ALL circuits\n", " (each depends on previous results), but we can use strategies:\n", "\n", " Strategy 1: \"Speculative execution\" — build circuits for both\n", " bit=0 and bit=1 outcomes, submit in parallel,\n", " discard the wrong branch.\n", "\n", " Strategy 2: \"Batch windows\" — since early bits have low confidence\n", " anyway, batch groups and accept slightly lower accuracy.\n", "\n", " Strategy 3: \"Hybrid\" — run first pass with random/classical bits,\n", " then refine with targeted quantum queries.\n", " \"\"\"\n", "\n", " def __init__(self, N, backend, shots=4096, verbose=True,\n", " max_hw_qubits=133, batch_size=100):\n", " # Removed super().__init__ as this class does not inherit from a base class.\n", " self.N = N\n", " self.n = N.bit_length() # N's bit length is often needed\n", " self.backend = backend\n", " self.shots = shots\n", " self.verbose = verbose\n", " self.max_hw_qubits = max_hw_qubits\n", " self.batch_size = batch_size\n", " self.oracle_calls = 0\n", " self.total_qubits = 0\n", " self.max_q_single = 0\n", "\n", " def _p(self, msg):\n", " if self.verbose:\n", " print(msg)\n", "\n", " def factor(self) -> Dict:\n", " N = self.N\n", " n = self.n\n", "\n", " # 1. Check for trivial factors (N is prime, N is a perfect power, N is even)\n", " if MathUtils.is_prime(N):\n", " self._p(f\"N={N} is prime. No factorization needed.\")\n", " return {'N': N, 'success': False, 'reason': 'N is prime', 'factor_1': None, 'factor_2': None, 'oracle_calls': 0, 'max_qubits_single_call': 0, 'total_qubits_used': 0}\n", "\n", " pp = MathUtils.is_perfect_power(N)\n", " if pp:\n", " self._p(f\"N={N} is a perfect power: {pp[0]}^{pp[1]}.\")\n", " return {'N': N, 'success': True, 'factor_1': pp[0], 'factor_2': N // pp[0], 'reason': 'Perfect power', 'oracle_calls': 0, 'max_qubits_single_call': 0, 'total_qubits_used': 0}\n", "\n", " if N % 2 == 0:\n", " self._p(f\"N={N} is even. Factor is 2.\")\n", " return {'N': N, 'success': True, 'factor_1': 2, 'factor_2': N // 2, 'reason': 'Even number', 'oracle_calls': 0, 'max_qubits_single_call': 0, 'total_qubits_used': 0}\n", "\n", " # 2. Choose a random 'a' such that 1 < a < N and gcd(a, N) = 1\n", " a = 0\n", " while True:\n", " a = random.randint(2, N - 1)\n", " if MathUtils.gcd(a, N) == 1:\n", " break\n", " self._p(f\"\\n Chosen a = {a}\")\n", "\n", " # 3. Find the order r of 'a' modulo N\n", " r = self._find_period(a)\n", " self._p(f\" Period r = {r}\")\n", "\n", " if r is None:\n", " self._p(\" ❌ Failed to find period.\")\n", " return {'N': N, 'success': False, 'reason': 'Failed to find period', 'factor_1': None, 'factor_2': None, 'oracle_calls': self.oracle_calls, 'max_qubits_single_call': self.max_q_single, 'total_qubits_used': self.total_qubits}\n", "\n", " if r % 2 != 0:\n", " self._p(\" Period is odd, trying again...\")\n", " return {'N': N, 'success': False, 'reason': 'Period is odd', 'factor_1': None, 'factor_2': None, 'oracle_calls': self.oracle_calls, 'max_qubits_single_call': self.max_q_single, 'total_qubits_used': self.total_qubits}\n", "\n", " x = pow(a, r // 2, N)\n", " if x == (N - 1) or x == 1:\n", " self._p(\" (a^(r/2)) mod N is N-1 or 1, trying again...\")\n", " return {'N': N, 'success': False, 'reason': 'Trivial factors from period', 'factor_1': None, 'factor_2': None, 'oracle_calls': self.oracle_calls, 'max_qubits_single_call': self.max_q_single, 'total_qubits_used': self.total_qubits}\n", "\n", " p1 = MathUtils.gcd(x - 1, N)\n", " p2 = MathUtils.gcd(x + 1, N)\n", "\n", " if p1 == 1 or p2 == 1:\n", " self._p(\" One of the GCDs is 1, trying again...\")\n", " return {'N': N, 'success': False, 'reason': 'Trivial factors from gcd', 'factor_1': None, 'factor_2': None, 'oracle_calls': self.oracle_calls, 'max_qubits_single_call': self.max_q_single, 'total_qubits_used': self.total_qubits}\n", "\n", " self._p(f\" Found factors: {p1}, {p2}\")\n", " return {'N': N, 'success': True, 'factor_1': p1, 'factor_2': p2, 'oracle_calls': self.oracle_calls, 'max_qubits_single_call': self.max_q_single, 'total_qubits_used': self.total_qubits}\n", "\n", " def _find_period(self, a: int) -> Optional[int]:\n", " \"\"\"\n", " Optimized period finding with batched execution.\n", "\n", " Phase 1: Classical estimation (free, instant)\n", " Phase 2: Quantum refinement (batched, fast)\n", " Phase 3: Continued fractions extraction\n", " \"\"\"\n", " N, n = self.N, self.n\n", " n_bits = 2 * n + 4\n", "\n", " self._p(f\" Phase 1: Classical pre-computation...\")\n", "\n", " # ── Phase 1: Precompute ALL b_k values classically ──\n", " # This is the expensive classical part but requires ZERO quantum\n", " t0 = time.time()\n", " b_k_values = []\n", " for k in range(n_bits):\n", " b_k = pow(a, pow(2, k), N)\n", " b_k_values.append(b_k)\n", " if b_k == 1:\n", " self._p(f\" b_{k} = 1, period divides 2^{k}\")\n", " # All subsequent b_k will also be 1\n", " for _ in range(k + 1, n_bits):\n", " b_k_values.append(1)\n", " break\n", "\n", " t_classical = time.time() - t0\n", " self._p(f\" Precomputed {len(b_k_values)} b_k values in {t_classical:.1f}s\")\n", "\n", " # Count non-trivial oracle calls needed\n", " non_trivial = sum(1 for b in b_k_values if b != 1)\n", " self._p(f\" Non-trivial calls needed: {non_trivial}/{len(b_k_values)}\")\n", "\n", " # ── Phase 2: Build ALL circuits (no feedback first pass) ──\n", " self._p(f\" Phase 2: Building {non_trivial} circuits...\")\n", "\n", " all_circuits = []\n", " circuit_indices = [] # Maps circuit index to k value\n", "\n", " for k in range(n_bits - 1, -1, -1):\n", " if k >= len(b_k_values) or b_k_values[k] == 1:\n", " continue\n", "\n", " b_k = b_k_values[k]\n", "\n", " qr = QuantumRegister(1, 'q')\n", " cr = ClassicalRegister(1, 'c')\n", " qc = QuantumCircuit(qr, cr, name=f'k{k}')\n", "\n", " qc.h(qr[0])\n", "\n", " # Phase WITHOUT semiclassical correction (first pass)\n", " two_pi = mpmath.mpf(2) * mpmath.pi\n", " phi = float(two_pi * mpmath.mpf(b_k) / mpmath.mpf(N))\n", " qc.p(phi, qr[0])\n", "\n", " qc.h(qr[0])\n", " qc.measure(qr[0], cr[0])\n", "\n", " all_circuits.append(qc)\n", " circuit_indices.append(k)\n", "\n", " # ── Phase 3: Submit in batches ──\n", " self._p(f\" Phase 3: Submitting {len(all_circuits)} circuits \"\n", " f\"in batches of {self.batch_size}...\")\n", "\n", " all_bits = {} # k -> bit value\n", "\n", " n_batches = (len(all_circuits) + self.batch_size - 1) // self.batch_size\n", "\n", " for batch_idx in range(n_batches):\n", " start = batch_idx * self.batch_size\n", " end = min(start + self.batch_size, len(all_circuits))\n", " batch = all_circuits[start:end]\n", " batch_ks = circuit_indices[start:end]\n", "\n", " if self.verbose and batch_idx % max(n_batches // 10, 1) == 0:\n", " self._p(f\" Batch {batch_idx+1}/{n_batches} \"\n", " f\"({len(batch)} circuits)...\")\n", "\n", " try:\n", " # Transpile batch\n", " transpiled = transpile(batch, backend=self.backend,\n", " optimization_level=1)\n", "\n", " # Submit ALL circuits in ONE job\n", " sampler = Sampler(mode=self.backend)\n", " job = sampler.run(transpiled, shots=self.shots)\n", " result = job.result()\n", "\n", " # Extract results\n", " for i, k in enumerate(batch_ks):\n", " try:\n", " counts = result[i].data.c.get_counts()\n", " zeros = counts.get('0', 0)\n", " ones = counts.get('1', 0)\n", " all_bits[k] = 1 if ones > zeros else 0\n", " except:\n", " all_bits[k] = random.randint(0, 1)\n", "\n", " self.oracle_calls += len(batch)\n", " self.total_qubits += len(batch) # 1 qubit each\n", " self.max_q_single = 1\n", "\n", " except Exception as e:\n", " self._p(f\" ⚠️ Batch {batch_idx+1} failed: {e}\")\n", " for k in batch_ks:\n", " all_bits[k] = random.randint(0, 1)\n", "\n", " # ── Phase 4: Reconstruct phase ──\n", " self._p(f\" Phase 4: Reconstructing phase from {len(all_bits)} bits...\")\n", "\n", " bits = []\n", " for k in range(n_bits):\n", " if k in all_bits:\n", " bits.append(all_bits[k])\n", " else:\n", " bits.append(0) # b_k was 1, bit is 0\n", "\n", " # bits[k] = k-th bit, already in order\n", " phase_int = sum(b * (2**i) for i, b in enumerate(bits))\n", " Q = 2 ** n_bits\n", "\n", " self._p(f\" Phase: {phase_int}/{Q}\")\n", "\n", " # ── Phase 5: Extract period ──\n", " candidates = MathUtils.extract_period_candidates(phase_int, Q, N)\n", " self._p(f\" Candidates: {candidates[:10]}\")\n", "\n", " for rc in candidates:\n", " if rc > 0 and pow(a, rc, N) == 1:\n", " self._p(f\" ✓ Period: r = {rc}\")\n", " return rc\n", "\n", " # ── Phase 6: Refinement with semiclassical feedback ──\n", " # If first pass failed, do targeted refinement on key bits\n", " self._p(f\" Phase 6: Semiclassical refinement...\")\n", "\n", " # Try multiple phase reconstructions with bit flips\n", " for n_flips in range(1, min(20, n_bits)):\n", " for flip_pos in random.sample(range(n_bits), min(n_flips, n_bits)):\n", " bits_copy = bits.copy()\n", " bits_copy[flip_pos] ^= 1\n", "\n", " phase_int2 = sum(b * (2**i) for i, b in enumerate(bits_copy))\n", " cands2 = MathUtils.extract_period_candidates(phase_int2, Q, N)\n", "\n", " for rc in cands2:\n", " if rc > 0 and pow(a, rc, N) == 1:\n", " self._p(f\" ✓ Period (refined): r = {rc}\")\n", " return rc\n", "\n", " if N <= 10000:\n", " ro = MathUtils.multiplicative_order(a, N)\n", " if ro:\n", " return ro\n", "\n", " return None\n", "\n", "\n", "# ═══════════════════════════════════════════════════════════════\n", "# TIMING COMPARISON\n", "# ═══════════════════════════════════════════════════════════════\n", "\n", "def estimate_timing(N, hw_qubits=133, batch_size=100):\n", " n = N.bit_length()\n", " n_bits = 2 * n + 4\n", "\n", " print(f\"\\n{'═'*70}\")\n", " print(f\"⏱️ TIMING ESTIMATE FOR {n}-BIT N\")\n", " print(f\"{'═'*70}\")\n", "\n", " # Individual submission (original)\n", " time_per_call_individual = 15 # seconds (API overhead dominant)\n", " total_individual = n_bits * time_per_call_individual\n", "\n", " print(f\"\\n ORIGINAL (individual jobs):\")\n", " print(f\" Oracle calls: {n_bits}\")\n", " print(f\" ~{time_per_call_individual}s per call (API overhead)\")\n", " print(f\" Total: {total_individual}s = {total_individual/3600:.1f} hours\")\n", "\n", " # Batched submission\n", " n_batches = (n_bits + batch_size - 1) // batch_size\n", " time_per_batch = 30 # seconds (1 API call for 100 circuits)\n", " total_batched = n_batches * time_per_batch\n", "\n", " print(f\"\\n ⭐ BATCHED (groups of {batch_size}):\")\n", " print(f\" Batches: {n_batches}\")\n", " print(f\" ~{time_per_batch}s per batch\")\n", " print(f\" Total: {total_batched}s = {total_batched/60:.1f} minutes\")\n", "\n", " speedup = total_individual / total_batched if total_batched > 0 else float('inf')\n", " print(f\"\\n Speedup: {speedup:.0f}×\")\n", "\n", " # Classical precomputation\n", " # pow(a, 2^k, N) for k up to 4100, N = 2048 bits\n", " # Each modular exp is O(k * n^2) with repeated squaring\n", " # Total classical: O(n^2 * n * n) = O(n^4) bit operations\n", " print(f\"\\n Classical precomputation (pow(a, 2^k, N)):\")\n", " print(f\" {n_bits} modular exponentiations on {n}-bit numbers\")\n", " print(f\" Estimated: ~{n_bits * 0.01:.0f}s on modern CPU\")\n", " print(f\" (Python's pow() is highly optimized for big ints)\")\n", "\n", "\n", "# ═══════════════════════════════════════════════════════════════\n", "# UPDATED MAIN\n", "# ═══════════════════════════════════════════════════════════════\n", "\n", "def main():\n", " N = N_TO_FACTOR\n", " n = N.bit_length()\n", "\n", " print(\"╔\" + \"═\"*68 + \"╗\")\n", " print(\"║\" + \" SOSQ: Sequential Oracle SHOR Quantum\".center(68) + \"║\")\n", " print(\"║\" + \" Optimized Batch Execution\".center(68) + \"║\")\n", " print(\"║\" + f\" N has {n} bits\".center(68) + \"║\")\n", " print(\"║\" + f\" mpmath: {mpmath.mp.dps} decimal digits\".center(68) + \"║\")\n", " print(\"╚\" + \"═\"*68 + \"╝\")\n", "\n", " hw_qubits = 133\n", " batch_size = 100\n", "\n", " analyze(N, hw_qubits)\n", " estimate_timing(N, hw_qubits, batch_size)\n", "\n", " # Classical verification for small N\n", " if N <= 10000:\n", " a_t = 2 if MathUtils.gcd(2, N) == 1 else 3\n", " if MathUtils.gcd(a_t, N) == 1:\n", " r_t = MathUtils.multiplicative_order(a_t, N)\n", " if r_t:\n", " print(f\"\\n Classical check: ord_{N}({a_t}) = {r_t}\")\n", "\n", " # Connect\n", " backend = None\n", " if USE_REAL_HARDWARE and IBM_TOKEN != \"YOUR_IBM_QUANTUM_API_TOKEN_HERE\":\n", " conn = IBMConnector(IBM_TOKEN, PREFERRED_BACKEND)\n", " if conn.connect():\n", " backend = conn.backend\n", " try:\n", " hw_qubits = backend.num_qubits\n", " except:\n", " pass\n", " else:\n", " backend = conn.get_sim()\n", " else:\n", " if IBM_TOKEN == \"YOUR_IBM_QUANTUM_API_TOKEN_HERE\":\n", " print(\"\\n⚠️ No token — using Aer simulator\")\n", " from qiskit_aer import AerSimulator\n", " backend = AerSimulator()\n", "\n", " # Run with batched engine\n", " print(f\"\\n{'═'*70}\")\n", " print(f\"RUNNING SOSQ (BATCHED) — Backend: {backend}\")\n", " print(f\"{'═'*70}\")\n", "\n", " engine = SOSQEngine(\n", " N=N, backend=backend, shots=SHOTS,\n", " verbose=VERBOSE, max_hw_qubits=hw_qubits,\n", " batch_size=batch_size\n", " )\n", "\n", " t0 = time.time()\n", " result = engine.factor()\n", " elapsed = time.time() - t0\n", "\n", " # Results\n", " print(f\"\\n{'═'*70}\")\n", " print(\"RESULTS\")\n", " print(f\"{'═'*70}\")\n", " print(f\" N = {result['N']} ({result['N'].bit_length()} bits)\")\n", "\n", " if result['success']:\n", " f1, f2 = result['factor_1'], result['factor_2']\n", " print(f\" ✅ {f1} × {f2} = {f1*f2}\")\n", " assert f1 * f2 == result['N']\n", " print(f\" ✓ Verified\")\n", " else:\n", " print(f\" ❌ Failed\")\n", "\n", " print(f\"\\n Oracle calls: {result['oracle_calls']}\")\n", " print(f\" Max qubits/call: {result['max_qubits_single_call']}\")\n", " print(f\" Wall time: {elapsed:.1f}s ({elapsed/60:.1f}min)\")\n", " print(f\"{'═'*70}\")\n", "\n", " # Test cases\n", " if not USE_REAL_HARDWARE or IBM_TOKEN == \"YOUR_IBM_QUANTUM_API_TOKEN_HERE\":\n", " print(f\"\\n{'═'*70}\")\n", " print(\"TESTS\")\n", " print(f\"{'═'*70}\")\n", " for tN in [15, 21, 35, 55, 77, 143, 221]:\n", " if tN == N:\n", " continue\n", " e = SOSQEngine(N=tN, backend=backend, shots=1024,\n", " verbose=False, max_hw_qubits=hw_qubits,\n", " batch_size=50)\n", " r = e.factor()\n", " s = \"✅\" if r['success'] else \"❌\"\n", " if r['success']:\n", " print(f\" {s} {tN} = {r['factor_1']}×{r['factor_2']}, \"\n", " f\"{r['oracle_calls']} calls, {r['max_qubits_single_call']}q\")\n", " else:\n", " print(f\" {s} {tN}: failed\")\n", "\n", " print(f\"\\n{'═'*70}\")\n", " print(\"SOSQ complete — Author: Kaoru Aguilera Katayama\")\n", " print(f\"{'═'*70}\")\n", "\n", "\n", "if __name__ == \"__main__\":\n", " main()\n", "\n", "\n", "# ═════════════════════════════════════════════════════════════════════════════\n", "# ██ IBM QUANTUM CONNECTION ██\n", "# ═════════════════════════════════════════════════════════════════════════════\n", "\n", "class IBMConnector:\n", "\n", " def __init__(self, token, preferred=None):\n", " self.token = token\n", " self.preferred = preferred\n", " self.backend = None\n", "\n", " def connect(self):\n", " print(\"=\" * 70)\n", " print(\"Connecting to IBM Quantum...\")\n", " print(\"=\" * 70)\n", " try:\n", " QiskitRuntimeService.save_account(\n", " channel=\"ibm_cloud\", token=self.token, overwrite=True)\n", " svc = QiskitRuntimeService(channel=\"ibm_cloud\")\n", " print(\"✅ Connected!\")\n", "\n", " for b in svc.backends():\n", " try:\n", " st = b.status()\n", " print(f\" • {b.name}: {b.num_qubits}q, \"\n", " f\"{'✅' if st.operational else '❌'}, \"\n", " f\"{st.pending_jobs} pending\")\n", " except:\n", " print(f\" • {b.name}: (status N/A)\")\n", "\n", " if self.preferred:\n", " try:\n", " self.backend = svc.backend(self.preferred)\n", " except:\n", " self.backend = svc.least_busy(simulator=False, operational=True)\n", " else:\n", " self.backend = svc.least_busy(simulator=False, operational=True)\n", "\n", " print(f\"\\n🎯 Backend: {self.backend.name} ({self.backend.num_qubits} qubits)\")\n", " return True\n", " except Exception as e:\n", " print(f\"❌ Failed: {e}\")\n", " return False\n", "\n", " def get_sim(self):\n", " from qiskit_aer import AerSimulator\n", " return AerSimulator()\n", "\n", "\n", "# ═════════════════════════════════════════════════════════════════════════════\n", "# ██ ANALYSIS ██\n", "# ═════════════════════════════════════════════════════════════════════════════\n", "\n", "def analyze(N, hw=133):\n", " n = N.bit_length()\n", " print(f\"\\n{'═'*70}\")\n", " print(f\"QUBIT ANALYSIS — N is {n} bits, hardware has {hw} qubits\")\n", " print(f\"{'═'*70}\")\n", "\n", " shor = 2*n + n + n\n", " beau = 2*n + 3\n", " sosq_tgt = 1 + n\n", " sosq_1q = 1\n", "\n", " print(f\"\\n Standard Shor: ~{shor} qubits {'✅' if shor<=hw else '❌'}\")\n", " print(f\" Beauregard: {beau} qubits {'✅' if beau<=hw else '❌'}\")\n", " print(f\" SOSQ (with target): {sosq_tgt} qubits {'✅' if sosq_tgt<=hw else '❌'}\")\n", " print(f\" ⭐ SOSQ (single-qubit): {sosq_1q} qubit ✅ ALWAYS\")\n", " print(f\"\\n SOSQ single-qubit oracle calls: {2*n+4}\")\n", " print(f\" Classical cost per call: pow(a, 2^k, N) with {mpmath.mp.dps}-digit precision\")\n", "\n", " if n > hw:\n", " print(f\"\\n 📌 {n}-bit N on {hw}-qubit HW:\")\n", " print(f\" Only SOSQ single-qubit mode works!\")\n", " print(f\" All modular arithmetic computed classically (mpmath)\")\n", " print(f\" Quantum HW used only for phase estimation (1 qubit)\")\n" ] } ] }