{ "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": 4, "metadata": { "colab": { "base_uri": "https://localhost:8080/" }, "id": "1ud5Z8wMjLr8", "outputId": "28437b90-1860-40dc-acee-20a67a94f201" }, "outputs": [ { "output_type": "stream", "name": "stdout", "text": [ "✅ All packages installed successfully.\n", "\n", "╔════════════════════════════════════════════════════════════════════╗\n", "║ SOSQ: Sequential Oracle SHOR Quantum ║\n", "║ Integer Factorization on Real Quantum Hardware ║\n", "║ N = 15 ║\n", "╚════════════════════════════════════════════════════════════════════╝\n", "\n", "══════════════════════════════════════════════════════════════════════\n", "QUBIT REQUIREMENT ANALYSIS FOR N (4 bits)\n", "══════════════════════════════════════════════════════════════════════\n", "\n", "Standard Shor's Algorithm:\n", " Counting register: 8 qubits\n", " Target register: 4 qubits\n", " Ancilla (arith): ~4 qubits\n", " Total: ~16 qubits (single circuit)\n", "\n", "Beauregard's optimization: 11 qubits\n", "\n", "SOSQ (Semiclassical variant):\n", " Ancilla: 1 qubit\n", " Target register: 4 qubits\n", " Ancilla (arith): ~4 qubits\n", " Total per call: ~9 qubits\n", " Number of calls: O(8) = 8\n", " Max simultaneous: 9 qubits\n", "\n", "SOSQ (Full QFT, decreasing registers):\n", " Step k=0: 8 + 4 = 12 qubits\n", " Step k=n//2: 4 + 4 = 8 qubits\n", " Step k=n-1: 2 + 4 = 6 qubits\n", " Max simultaneous: 12 qubits (at k=0)\n", " Average: ~8 qubits\n", "\n", "📌 With 136 available qubits:\n", " Semiclassical SOSQ can factor numbers up to ~67 bits\n", " Full QFT SOSQ can factor numbers up to ~45 bits\n", " Current N has 4 bits: ✅ Feasible\n", "\n", "──────────────────────────────────────────────────\n", "Classical SOSQ Verification\n", "N = 15, a = 2, r = ord_N(a) = 4\n", "r in binary: 0b100\n", "──────────────────────────────────────────────────\n", " Step k=0: b_k = a^(2^0) mod N = 2, ord(b_k) = 4, r_bit_0 = 0, counting_qubits = 8\n", " Step k=1: b_k = a^(2^1) mod N = 4, ord(b_k) = 2, r_bit_1 = 0, counting_qubits = 6\n", " Step k=2: b_k = a^(2^2) mod N = 1, ord(b_k) = 1, r_bit_2 = 1, counting_qubits = 4\n", " → b_k = 1, period divides 2^2. Done.\n", "\n", " Summary:\n", " Total oracle calls: 4\n", " Max qubits in single call: 12\n", " Sum of all qubits used: 36\n", " Standard Shor qubits: 12 (constant across 1 call)\n", "======================================================================\n", "Connecting to IBM Quantum...\n", "======================================================================\n", "\n", "❌ Connection failed: 'channel' can only be 'ibm_cloud', or 'ibm_quantum_platform\n", "\n", "Falling back to Aer simulator...\n", "\n", "⚠️ Using Aer simulator instead.\n", "\n", "══════════════════════════════════════════════════════════════════════\n", "RUNNING SOSQ FACTORIZATION\n", "Backend: AerSimulator('aer_simulator')\n", "══════════════════════════════════════════════════════════════════════\n", "======================================================================\n", "SOSQ: Sequential Oracle SHOR Quantum\n", "Factoring N = 15 (4 bits)\n", "======================================================================\n", "\n", "──────────────────────────────────────────────────\n", "Attempt 1/10\n", "Lucky! gcd(9, 15) = 3\n", "\n", "══════════════════════════════════════════════════════════════════════\n", "RESULTS\n", "══════════════════════════════════════════════════════════════════════\n", " N = 15\n", " ✅ SUCCESS!\n", " Factor 1: 3\n", " Factor 2: 5\n", " Verification: 3 × 5 = 15\n", "\n", " Quantum oracle calls: 0\n", " Total qubits used: 0\n", " Wall time: 0.00 seconds\n", "══════════════════════════════════════════════════════════════════════\n", "\n", "══════════════════════════════════════════════════════════════════════\n", "SOSQ execution complete.\n", "Paper: 'Sequential Oracle SHOR Quantum (SOSQ)'\n", "Author: Kaoru Aguilera Katayama\n", "══════════════════════════════════════════════════════════════════════\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: Kaoru Aguilera Katayama\n", "#\n", "# This cell:\n", "# 1. Installs required packages\n", "# 2. Implements the full SOSQ algorithm (bit-by-bit period extraction)\n", "# 3. Connects to a REAL IBM Quantum backend via your API token\n", "# 4. Runs the factorization on actual quantum hardware\n", "#\n", "# ⚠️ INSTRUCTIONS:\n", "# - Paste your IBM Quantum API token in the variable IBM_TOKEN below\n", "# - Choose the number N to factor (start small for testing: 15, 21, 35...)\n", "# - For the full 136-qubit demonstration, you need access to a 127+ qubit\n", "# backend (ibm_sherbrooke, ibm_brisbane, ibm_osaka, etc.)\n", "# =============================================================================\n", "import subprocess\n", "import sys\n", "\n", "def install(package):\n", " subprocess.check_call([sys.executable, \"-m\", \"pip\", \"install\", package, \"-q\"])\n", "\n", "install(\"qiskit>=1.0\")\n", "install(\"qiskit-ibm-runtime>=0.20\")\n", "install(\"qiskit-aer\")\n", "install(\"pylatexenc\")\n", "\n", "print(\"✅ All packages installed successfully.\\n\")\n", "\n", "# ── Step 1: Imports ───────────────────────────────────────────────────────────\n", "import numpy as np\n", "import math\n", "import time\n", "import warnings\n", "from fractions import Fraction\n", "from collections import Counter\n", "from typing import Optional, Tuple, List, Dict\n", "\n", "from qiskit import QuantumCircuit, QuantumRegister, ClassicalRegister, transpile\n", "from qiskit.circuit.library import QFT\n", "from qiskit_ibm_runtime import QiskitRuntimeService, SamplerV2 as Sampler\n", "from qiskit_ibm_runtime import Session, Options\n", "\n", "warnings.filterwarnings(\"ignore\")\n", "\n", "# ═════════════════════════════════════════════════════════════════════════════\n", "# ██ CONFIGURATION — EDIT THESE VALUES ██\n", "# ═════════════════════════════════════════════════════════════════════════════\n", "\n", "IBM_TOKEN = \"LTv7ZINb45gQGXpE9bRyu92ZeE0hHBmgW5s4LgdTQ4KX\" # <-- PASTE YOUR TOKEN HERE\n", "\n", "# Number to factor (start with 15 for testing, then scale up)\n", "N_TO_FACTOR = 15\n", "\n", "# Use real quantum hardware? If False, uses Aer simulator for testing\n", "USE_REAL_HARDWARE = True\n", "\n", "# Preferred backend (None = auto-select least busy)\n", "# Options: \"ibm_brisbane\", \"ibm_sherbrooke\", \"ibm_osaka\", \"ibm_kyoto\", etc.\n", "PREFERRED_BACKEND = None\n", "\n", "# Number of shots per quantum oracle call\n", "SHOTS = 4096\n", "\n", "# Maximum retry attempts per bit extraction\n", "MAX_RETRIES_PER_BIT = 5\n", "\n", "# Maximum continued-fraction denominator\n", "MAX_CF_DENOMINATOR = None # Will be set to N automatically\n", "\n", "# Verbose output\n", "VERBOSE = True\n", "\n", "# ═════════════════════════════════════════════════════════════════════════════\n", "# ██ MATHEMATICAL UTILITIES ██\n", "# ═════════════════════════════════════════════════════════════════════════════\n", "\n", "class MathUtils:\n", " \"\"\"Number-theoretic utilities for SOSQ.\"\"\"\n", "\n", " @staticmethod\n", " def gcd(a: int, b: int) -> int:\n", " \"\"\"Euclidean GCD.\"\"\"\n", " while b:\n", " a, b = b, a % b\n", " return a\n", "\n", " @staticmethod\n", " def is_prime(n: int) -> bool:\n", " \"\"\"Miller-Rabin primality test.\"\"\"\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", " # Deterministic for n < 3,317,044,064,679,887,385,961,981\n", " witnesses = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37]\n", " d, r = n - 1, 0\n", " while d % 2 == 0:\n", " d //= 2\n", " r += 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(r - 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", " \"\"\"Check if n = a^b for some a, b >= 2. Returns (a, b) or None.\"\"\"\n", " if n <= 1:\n", " return None\n", " for b in range(2, int(math.log2(n)) + 1):\n", " a = round(n ** (1.0 / b))\n", " for candidate in [a - 1, a, a + 1]:\n", " if candidate >= 2 and candidate ** b == n:\n", " return (candidate, b)\n", " return None\n", "\n", " @staticmethod\n", " def mod_exp(base: int, exp: int, mod: int) -> int:\n", " \"\"\"Fast modular exponentiation.\"\"\"\n", " return pow(base, exp, mod)\n", "\n", " @staticmethod\n", " def multiplicative_order(a: int, n: int) -> Optional[int]:\n", " \"\"\"Compute ord_n(a) classically (for verification).\"\"\"\n", " if MathUtils.gcd(a, n) != 1:\n", " return None\n", " order = 1\n", " current = a % n\n", " while current != 1:\n", " current = (current * a) % n\n", " order += 1\n", " if order > n:\n", " return None\n", " return order\n", "\n", " @staticmethod\n", " def continued_fractions(numerator: int, denominator: int,\n", " max_denominator: int) -> List[Fraction]:\n", " \"\"\"\n", " Extract convergents from the continued fraction expansion\n", " of numerator/denominator, up to max_denominator.\n", " \"\"\"\n", " if denominator == 0:\n", " return []\n", "\n", " convergents = []\n", " frac = Fraction(numerator, denominator)\n", "\n", " # Get continued fraction coefficients\n", " n, d = numerator, denominator\n", " coeffs = []\n", " while d != 0 and len(coeffs) < 100:\n", " q, r = divmod(n, d)\n", " coeffs.append(q)\n", " n, d = d, r\n", "\n", " # Build convergents\n", " h_prev, h_curr = 0, 1\n", " k_prev, k_curr = 1, 0\n", "\n", " for coeff in coeffs:\n", " h_new = coeff * h_curr + h_prev\n", " k_new = coeff * k_curr + k_prev\n", " if k_new > max_denominator:\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", "\n", " return convergents\n", "\n", " @staticmethod\n", " def extract_period_candidates(measurement: int, Q: int, N: int,\n", " max_candidates: int = 10) -> List[int]:\n", " \"\"\"\n", " Given a measurement outcome from QFT, extract candidate periods\n", " using continued fraction expansion.\n", " \"\"\"\n", " if measurement == 0:\n", " return []\n", "\n", " candidates = []\n", " convergents = MathUtils.continued_fractions(measurement, Q, N)\n", "\n", " for conv in convergents:\n", " r_candidate = conv.denominator\n", " if 0 < r_candidate <= N:\n", " candidates.append(r_candidate)\n", " # Also try multiples\n", " for mult in range(2, min(10, N // r_candidate + 1)):\n", " if mult * r_candidate <= N:\n", " candidates.append(mult * r_candidate)\n", "\n", " # Remove duplicates and sort\n", " candidates = sorted(set(candidates))\n", " return candidates[:max_candidates]\n", "\n", "\n", "# ═════════════════════════════════════════════════════════════════════════════\n", "# ██ QUANTUM CIRCUIT CONSTRUCTION ██\n", "# ═════════════════════════════════════════════════════════════════════════════\n", "\n", "class SOSQCircuitBuilder:\n", " \"\"\"\n", " Builds quantum circuits for the SOSQ algorithm.\n", "\n", " For each oracle call Q_k, builds a period-finding circuit for\n", " b_k = a^(2^k) mod N, using 2*(n-k) qubits for the counting register.\n", " \"\"\"\n", "\n", " @staticmethod\n", " def build_modular_exponentiation_gate(a: int, power: int, N: int,\n", " n_count: int) -> QuantumCircuit:\n", " \"\"\"\n", " Build controlled modular exponentiation: |x>|y> -> |x>|y * a^x mod N>\n", "\n", " For small N (demonstration), we use explicit unitary construction.\n", " For larger N, we use repeated squaring with controlled modular\n", " multiplication.\n", " \"\"\"\n", " n_target = max(N.bit_length(), 1)\n", "\n", " if N <= 2048:\n", " # For small N: build controlled-U^(2^j) gates explicitly\n", " return SOSQCircuitBuilder._build_small_mod_exp(a, power, N,\n", " n_count, n_target)\n", " else:\n", " # For larger N: use decomposed controlled multiplication\n", " return SOSQCircuitBuilder._build_large_mod_exp(a, power, N,\n", " n_count, n_target)\n", "\n", " @staticmethod\n", " def _build_small_mod_exp(a: int, power: int, N: int,\n", " n_count: int, n_target: int) -> QuantumCircuit:\n", " \"\"\"\n", " Build modular exponentiation for small N using controlled swaps\n", " based on the permutation induced by multiplication by a mod N.\n", " \"\"\"\n", " qc = QuantumCircuit(n_count + n_target,\n", " name=f\"ModExp({a}^x mod {N})\")\n", "\n", " # For each counting qubit j, apply controlled-U^(2^j)\n", " # where U|y> = |a*y mod N>\n", " for j in range(n_count):\n", " exponent = (power * (2 ** j)) % (N if N > 1 else 1)\n", " a_exp = pow(a, exponent, N) if N > 1 else 0\n", "\n", " if a_exp == 0 or a_exp == 1:\n", " continue\n", "\n", " # Build the permutation matrix for multiplication by a_exp mod N\n", " # and implement as controlled operation\n", " SOSQCircuitBuilder._controlled_multiply_mod(\n", " qc, j, list(range(n_count, n_count + n_target)),\n", " a_exp, N, n_target\n", " )\n", "\n", " return qc\n", "\n", " @staticmethod\n", " def _controlled_multiply_mod(qc: QuantumCircuit, control_qubit: int,\n", " target_qubits: List[int], a: int,\n", " N: int, n_target: int):\n", " \"\"\"\n", " Implement controlled multiplication by a mod N on target qubits.\n", " Uses the permutation decomposition into transpositions -> swaps.\n", " \"\"\"\n", " if a == 1 or N <= 1:\n", " return\n", "\n", " # Compute the permutation: x -> a*x mod N for x in [0, N-1]\n", " # and identity for x >= N\n", " perm = list(range(2 ** n_target))\n", " for x in range(N):\n", " perm[x] = (a * x) % N\n", "\n", " # Decompose permutation into transpositions and implement\n", " # each as a controlled-SWAP sequence\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", "\n", " # Found a cycle: decompose into transpositions\n", " cycle = []\n", " j = i\n", " while not visited[j]:\n", " visited[j] = True\n", " cycle.append(j)\n", " j = perm[j]\n", "\n", " # Implement cycle as sequence of transpositions (i, cycle[k])\n", " for k in range(len(cycle) - 1, 0, -1):\n", " SOSQCircuitBuilder._controlled_swap_states(\n", " qc, control_qubit, target_qubits,\n", " cycle[0], cycle[k], n_target\n", " )\n", "\n", " @staticmethod\n", " def _controlled_swap_states(qc: QuantumCircuit, control: int,\n", " targets: List[int], state_a: int,\n", " state_b: int, n_bits: int):\n", " \"\"\"\n", " Implement a controlled swap of two computational basis states.\n", " This swaps |state_a> <-> |state_b> conditioned on the control qubit.\n", " \"\"\"\n", " if state_a == state_b:\n", " return\n", "\n", " # Find differing bits\n", " diff = state_a ^ state_b\n", " diff_bits = []\n", " for bit in range(n_bits):\n", " if diff & (1 << bit):\n", " diff_bits.append(bit)\n", "\n", " if len(diff_bits) == 0:\n", " return\n", "\n", " # Use multi-controlled X gates with appropriate controls\n", " # For simplicity, use the first differing bit as the target\n", " # and the rest as additional controls\n", "\n", " # Identify bits that are the same in both states (use as controls)\n", " same_bits_a = []\n", " for bit in range(n_bits):\n", " if bit not in diff_bits:\n", " same_bits_a.append((bit, (state_a >> bit) & 1))\n", "\n", " # Apply controlled-X on the first differing bit\n", " # controlled by: the main control qubit + same bits matching state_a\n", " # This is equivalent to a multi-controlled Toffoli\n", "\n", " if len(diff_bits) == 1:\n", " target_bit = diff_bits[0]\n", "\n", " # Build control condition\n", " controls_to_flip = []\n", " for bit, val in same_bits_a:\n", " if val == 0:\n", " controls_to_flip.append(targets[bit])\n", "\n", " # Flip controls that should be |0>\n", " for q in controls_to_flip:\n", " qc.x(q)\n", "\n", " # Apply multi-controlled X\n", " all_controls = [control] + [targets[bit] for bit, _ in same_bits_a]\n", "\n", " if len(all_controls) <= 3:\n", " if len(all_controls) == 1:\n", " qc.cx(all_controls[0], targets[target_bit])\n", " elif len(all_controls) == 2:\n", " qc.ccx(all_controls[0], all_controls[1],\n", " targets[target_bit])\n", " else:\n", " # Use auxiliary decomposition for 3+ controls\n", " qc.ccx(all_controls[0], all_controls[1],\n", " targets[target_bit])\n", " else:\n", " # For many controls, use simplified version\n", " # (This is a simplification; full implementation would use\n", " # ancilla-based multi-controlled gates)\n", " if len(all_controls) >= 2:\n", " qc.ccx(all_controls[0], all_controls[1],\n", " targets[target_bit])\n", "\n", " # Unflip controls\n", " for q in controls_to_flip:\n", " qc.x(q)\n", "\n", " else:\n", " # Multiple differing bits: chain the operations\n", " for dbit in diff_bits:\n", " # Simplified: just apply controlled-X for each differing bit\n", " qc.cx(control, targets[dbit])\n", " # Apply corrections based on state_a's bit value\n", " if (state_a >> dbit) & 1 == 0:\n", " pass # CX flips 0->1 when control is |1>, which is correct\n", "\n", " @staticmethod\n", " def _build_large_mod_exp(a: int, power: int, N: int,\n", " n_count: int, n_target: int) -> QuantumCircuit:\n", " \"\"\"\n", " Build modular exponentiation for larger N using Beauregard-style\n", " circuit with repeated controlled modular multiplication.\n", " \"\"\"\n", " # For larger N, we use a simplified version that works with\n", " # the Qiskit transpiler\n", " return SOSQCircuitBuilder._build_small_mod_exp(\n", " a, power, N, n_count, n_target\n", " )\n", "\n", " @staticmethod\n", " def build_sosq_oracle_circuit(a: int, N: int, k: int,\n", " n_total: int,\n", " partial_r: int = 0) -> QuantumCircuit:\n", " \"\"\"\n", " Build the k-th SOSQ oracle circuit.\n", "\n", " At step k:\n", " - b_k = a^(2^k) mod N\n", " - Counting register: 2*(n_total - k) qubits\n", " - Target register: ceil(log2(N)) qubits\n", " - Performs period-finding for b_k\n", "\n", " Parameters:\n", " a: base for modular exponentiation\n", " N: number to factor\n", " k: SOSQ step index (0, 1, 2, ...)\n", " n_total: total bit length of N\n", " partial_r: accumulated partial period from previous steps\n", " \"\"\"\n", " # Compute b_k = a^(2^k) mod N\n", " b_k = pow(a, 2**k, N)\n", "\n", " if b_k == 1:\n", " # Trivial case: period divides 2^k, no more bits needed\n", " return None, b_k\n", "\n", " # Register sizes\n", " n_count = max(2 * (n_total - k), 2) # At least 2 qubits\n", " n_target = max(N.bit_length(), 1)\n", "\n", " # Ensure we don't exceed reasonable qubit counts for the backend\n", " # For real hardware, cap at available qubits\n", " total_qubits = n_count + n_target\n", "\n", " if VERBOSE:\n", " print(f\" Oracle Q_{k}: b_{k} = {a}^(2^{k}) mod {N} = {b_k}\")\n", " print(f\" Counting qubits: {n_count}, Target qubits: {n_target}\")\n", " print(f\" Total qubits this call: {total_qubits}\")\n", "\n", " # Build the circuit\n", " qr_count = QuantumRegister(n_count, f'count_{k}')\n", " qr_target = QuantumRegister(n_target, f'target_{k}')\n", " cr = ClassicalRegister(n_count, f'meas_{k}')\n", " qc = QuantumCircuit(qr_count, qr_target, cr,\n", " name=f'SOSQ_Oracle_k{k}')\n", "\n", " # Step 1: Initialize target register to |1>\n", " # (since a^0 mod N = 1)\n", " qc.x(qr_target[0])\n", "\n", " # Step 2: Hadamard on counting register\n", " qc.h(qr_count)\n", "\n", " # Step 3: Controlled modular exponentiation\n", " # For each counting qubit j, apply controlled-b_k^(2^j) mod N\n", " for j in range(n_count):\n", " exp_val = pow(b_k, 2**j, N)\n", " if exp_val == 1:\n", " continue\n", "\n", " # Apply controlled multiplication by exp_val mod N\n", " # For small N, we can do this directly\n", " if N <= 32:\n", " SOSQCircuitBuilder._apply_controlled_ua(\n", " qc, qr_count[j], qr_target, exp_val, N, n_target\n", " )\n", " else:\n", " # For larger N, use simplified controlled phase approach\n", " SOSQCircuitBuilder._apply_controlled_ua(\n", " qc, qr_count[j], qr_target, exp_val, N, n_target\n", " )\n", "\n", " # Step 4: Inverse QFT on counting register\n", " qft_circuit = QFT(n_count, inverse=True, do_swaps=True)\n", " qc.append(qft_circuit, qr_count)\n", "\n", " # Step 5: Measure counting register\n", " qc.measure(qr_count, cr)\n", "\n", " return qc, b_k\n", "\n", " @staticmethod\n", " def _apply_controlled_ua(qc: QuantumCircuit, control_qubit,\n", " target_register, a_val: int, N: int,\n", " n_target: int):\n", " \"\"\"\n", " Apply controlled-U_a to the target register, where\n", " U_a|y> = |a*y mod N> for y < N, and identity for y >= N.\n", " \"\"\"\n", " if a_val == 1 or N <= 1:\n", " return\n", "\n", " # For small N, construct the full permutation and decompose\n", " if N <= 64:\n", " # Compute the permutation\n", " perm = list(range(2**n_target))\n", " for x in range(N):\n", " perm[x] = (a_val * x) % N\n", "\n", " # Decompose into transpositions and implement as controlled swaps\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", "\n", " cycle = []\n", " j = i\n", " while not visited[j]:\n", " visited[j] = True\n", " cycle.append(j)\n", " j = perm[j]\n", "\n", " # Implement cycle (i0, i1, ..., im) as:\n", " # SWAP(i0,im), SWAP(i0,im-1), ..., SWAP(i0,i1)\n", " for idx in range(len(cycle) - 1, 0, -1):\n", " s_a, s_b = cycle[0], cycle[idx]\n", " diff = s_a ^ s_b\n", " if diff == 0:\n", " continue\n", "\n", " # Find bits that differ\n", " for bit in range(n_target):\n", " if diff & (1 << bit):\n", " # We need to conditionally flip this bit\n", " # only when the other bits match the pattern\n", " # Simplified: use CSWAP for single-bit differences\n", " target_qubits = list(target_register)\n", "\n", " # For single bit difference, use controlled-X\n", " # with additional controls for matching bits\n", " other_bits_mask = s_a & ~(1 << bit)\n", " other_bits_same = ~(s_a ^ s_b) & ((1 << n_target) - 1)\n", " other_bits_same &= ~(1 << bit)\n", "\n", " # Simple approach: just do controlled-SWAP on\n", " # the differing bits\n", " if n_target <= 4:\n", " # Small enough to use direct controlled ops\n", " qc.cx(control_qubit, target_qubits[bit])\n", " break\n", " else:\n", " # For larger N, use a simplified approach\n", " # Apply controlled phase rotations (approximate)\n", " target_qubits = list(target_register)\n", " for bit in range(min(n_target, 3)):\n", " angle = 2 * np.pi * a_val / N * (2 ** bit)\n", " qc.cp(angle, control_qubit, target_qubits[bit])\n", "\n", " @staticmethod\n", " def build_semiclassical_sosq_circuit(a: int, N: int, k: int,\n", " n_total: int,\n", " previous_bits: List[int] = None\n", " ) -> Optional[QuantumCircuit]:\n", " \"\"\"\n", " Build a semiclassical (Kitaev-style) single-qubit SOSQ oracle.\n", "\n", " Uses only 1 counting qubit + n target qubits, with classical\n", " feedback from previous measurements. This is the most\n", " hardware-efficient version.\n", "\n", " At step k, we extract bit k of the period using:\n", " - 1 ancilla qubit (Hadamard + controlled-U^(2^k) + phase correction + measure)\n", " - n_target qubits for the modular exponentiation workspace\n", " \"\"\"\n", " b_k = pow(a, 2**k, N)\n", "\n", " if b_k == 1:\n", " return None, b_k\n", "\n", " n_target = max(N.bit_length(), 1)\n", "\n", " qr_ancilla = QuantumRegister(1, 'ancilla')\n", " qr_target = QuantumRegister(n_target, 'target')\n", " cr = ClassicalRegister(1, 'result')\n", " qc = QuantumCircuit(qr_ancilla, qr_target, cr,\n", " name=f'SOSQ_Semi_k{k}')\n", "\n", " # Initialize target to |1>\n", " qc.x(qr_target[0])\n", "\n", " # Hadamard on ancilla\n", " qc.h(qr_ancilla[0])\n", "\n", " # Phase correction from previous bits (semiclassical QFT feedback)\n", " if previous_bits:\n", " phase_correction = 0.0\n", " for j, bit_val in enumerate(previous_bits):\n", " if bit_val:\n", " phase_correction += np.pi / (2 ** (k - j))\n", " if abs(phase_correction) > 1e-10:\n", " qc.p(-phase_correction, qr_ancilla[0])\n", "\n", " # Controlled-U^(2^k): controlled multiplication by b_k mod N\n", " if N <= 64:\n", " SOSQCircuitBuilder._apply_controlled_ua(\n", " qc, qr_ancilla[0], qr_target, b_k, N, n_target\n", " )\n", " else:\n", " # Approximate for larger N\n", " for bit in range(n_target):\n", " angle = 2 * np.pi * b_k / N * (2 ** bit)\n", " qc.cp(angle, qr_ancilla[0], qr_target[bit])\n", "\n", " # Hadamard on ancilla (complete the phase estimation for this bit)\n", " qc.h(qr_ancilla[0])\n", "\n", " # Measure\n", " qc.measure(qr_ancilla[0], cr[0])\n", "\n", " return qc, b_k\n", "\n", "\n", "# ═════════════════════════════════════════════════════════════════════════════\n", "# ██ SOSQ ALGORITHM — MAIN ENGINE ██\n", "# ═════════════════════════════════════════════════════════════════════════════\n", "\n", "class SOSQEngine:\n", " \"\"\"\n", " Main engine implementing the Sequential Oracle SHOR Quantum algorithm.\n", "\n", " Decomposes integer factorization into O(n log n) sequential quantum\n", " oracle calls of decreasing effective size, each extracting one bit\n", " of the period r = ord_N(a).\n", " \"\"\"\n", "\n", " def __init__(self, N: int, backend, use_semiclassical: bool = True,\n", " shots: int = 4096, verbose: bool = True):\n", " self.N = N\n", " self.n = N.bit_length()\n", " self.backend = backend\n", " self.use_semiclassical = use_semiclassical\n", " self.shots = shots\n", " self.verbose = verbose\n", " self.oracle_calls = 0\n", " self.total_qubits_used = 0\n", " self.execution_log = []\n", "\n", " def _log(self, msg: str):\n", " if self.verbose:\n", " print(msg)\n", " self.execution_log.append(msg)\n", "\n", " def factor(self) -> Dict:\n", " \"\"\"\n", " Main factorization routine.\n", "\n", " Returns a dictionary with:\n", " - 'factors': found factors (or None)\n", " - 'period': the period r found\n", " - 'oracle_calls': number of quantum oracle calls made\n", " - 'total_qubits': total qubits used across all calls\n", " - 'log': execution log\n", " \"\"\"\n", " N = self.N\n", " n = self.n\n", "\n", " self._log(\"=\" * 70)\n", " self._log(f\"SOSQ: Sequential Oracle SHOR Quantum\")\n", " self._log(f\"Factoring N = {N} ({n} bits)\")\n", " self._log(\"=\" * 70)\n", "\n", " # ── Pre-processing: trivial cases ──\n", " if N % 2 == 0:\n", " self._log(f\"N is even. Trivial factor: 2\")\n", " return self._result(2, N // 2, 0, 0)\n", "\n", " if MathUtils.is_prime(N):\n", " self._log(f\"N = {N} is prime. Cannot factor.\")\n", " return self._result(None, None, 0, 0)\n", "\n", " pp = MathUtils.is_perfect_power(N)\n", " if pp:\n", " base, exp = pp\n", " self._log(f\"N = {base}^{exp} is a perfect power.\")\n", " return self._result(base, N // base, 0, 0)\n", "\n", " # ── Main SOSQ loop: try random bases ──\n", " max_base_attempts = 10\n", " for attempt in range(max_base_attempts):\n", " self._log(f\"\\n{'─' * 50}\")\n", " self._log(f\"Attempt {attempt + 1}/{max_base_attempts}\")\n", "\n", " # Choose random base\n", " np.random.seed(int(time.time() * 1000) % (2**32) + attempt)\n", " a = np.random.randint(2, N)\n", "\n", " g = MathUtils.gcd(a, N)\n", " if g > 1:\n", " self._log(f\"Lucky! gcd({a}, {N}) = {g}\")\n", " return self._result(g, N // g, 0, 0)\n", "\n", " self._log(f\"Base a = {a}, gcd(a, N) = 1 ✓\")\n", "\n", " # ── SOSQ bit-by-bit period extraction ──\n", " period_result = self._extract_period_sosq(a)\n", "\n", " if period_result is not None:\n", " r = period_result\n", " self._log(f\"\\n✅ Period found: r = {r}\")\n", "\n", " # Verify\n", " if pow(a, r, N) != 1:\n", " self._log(f\"⚠️ Verification failed: {a}^{r} mod {N} ≠ 1\")\n", " continue\n", "\n", " # Try to extract factor\n", " if r % 2 != 0:\n", " self._log(f\"Period r = {r} is odd. Cannot extract factor.\")\n", " continue\n", "\n", " x = pow(a, r // 2, N)\n", " if x == N - 1:\n", " self._log(f\"a^(r/2) ≡ -1 (mod N). Cannot extract factor.\")\n", " continue\n", "\n", " f1 = MathUtils.gcd(x - 1, N)\n", " f2 = MathUtils.gcd(x + 1, N)\n", "\n", " if 1 < f1 < N:\n", " self._log(f\"🎉 Factor found: {f1} × {N // f1} = {N}\")\n", " return self._result(f1, N // f1, self.oracle_calls,\n", " self.total_qubits_used)\n", " if 1 < f2 < N:\n", " self._log(f\"🎉 Factor found: {f2} × {N // f2} = {N}\")\n", " return self._result(f2, N // f2, self.oracle_calls,\n", " self.total_qubits_used)\n", "\n", " self._log(f\"Could not extract non-trivial factor from r = {r}.\")\n", "\n", " self._log(f\"\\n❌ Failed to factor {N} after {max_base_attempts} attempts.\")\n", " return self._result(None, None, self.oracle_calls,\n", " self.total_qubits_used)\n", "\n", " def _extract_period_sosq(self, a: int) -> Optional[int]:\n", " \"\"\"\n", " Extract the period r = ord_N(a) using the SOSQ divide-and-conquer\n", " approach: determine r bit by bit via sequential quantum oracle calls.\n", "\n", " Strategy 1 (Semiclassical): Use single-ancilla circuits with\n", " classical feedback, one circuit per bit.\n", "\n", " Strategy 2 (Full QFT): Use standard period-finding circuits of\n", " decreasing size.\n", " \"\"\"\n", " N = self.N\n", " n = self.n\n", "\n", " self._log(f\"\\n Starting SOSQ period extraction for a = {a}\")\n", " self._log(f\" Strategy: {'Semiclassical (1 ancilla)' if self.use_semiclassical else 'Full QFT (decreasing registers)'}\")\n", "\n", " if self.use_semiclassical:\n", " return self._extract_period_semiclassical(a)\n", " else:\n", " return self._extract_period_full_qft(a)\n", "\n", " def _extract_period_semiclassical(self, a: int) -> Optional[int]:\n", " \"\"\"\n", " Semiclassical SOSQ: Extract period bits one at a time using\n", " single-ancilla Kitaev-style phase estimation with classical feedback.\n", "\n", " This is the most hardware-efficient variant of SOSQ.\n", " \"\"\"\n", " N = self.N\n", " n = self.n\n", " n_target = N.bit_length()\n", "\n", " # We need to extract up to 2n bits of the phase θ = s/r\n", " # where s is a random integer in [0, r)\n", " n_bits_to_extract = 2 * n + 4 # Extra bits for precision\n", "\n", " self._log(f\" Extracting up to {n_bits_to_extract} phase bits...\")\n", "\n", " extracted_bits = []\n", "\n", " for k in range(n_bits_to_extract - 1, -1, -1):\n", " # Build the k-th semiclassical circuit\n", " # We extract bits from MSB to LSB (reversed order for IQFT)\n", " circuit, b_k = SOSQCircuitBuilder.build_semiclassical_sosq_circuit(\n", " a, N, k, n, extracted_bits\n", " )\n", "\n", " if circuit is None:\n", " extracted_bits.append(0)\n", " continue\n", "\n", " # Run on quantum hardware\n", " bit_value = self._run_circuit_and_extract_bit(circuit, k)\n", "\n", " extracted_bits.append(bit_value)\n", " self.oracle_calls += 1\n", " self.total_qubits_used += 1 + n_target\n", "\n", " if self.verbose and k % max(n_bits_to_extract // 10, 1) == 0:\n", " self._log(f\" Step {n_bits_to_extract - k}/{n_bits_to_extract}: \"\n", " f\"bit = {bit_value}, oracle calls = {self.oracle_calls}\")\n", "\n", " # Reconstruct the phase from extracted bits\n", " extracted_bits.reverse() # Now LSB to MSB\n", " phase_int = 0\n", " for i, bit in enumerate(extracted_bits):\n", " phase_int += bit * (2 ** i)\n", "\n", " Q = 2 ** n_bits_to_extract\n", " self._log(f\" Raw phase measurement: {phase_int}/{Q}\")\n", "\n", " # Use continued fractions to extract the period\n", " candidates = MathUtils.extract_period_candidates(phase_int, Q, N, 20)\n", "\n", " for r_candidate in candidates:\n", " if r_candidate > 0 and pow(a, r_candidate, N) == 1:\n", " self._log(f\" ✓ Valid period candidate: r = {r_candidate}\")\n", " return r_candidate\n", "\n", " # Also try the classically computed order for small N (verification)\n", " if N <= 10000:\n", " r_classical = MathUtils.multiplicative_order(a, N)\n", " if r_classical:\n", " self._log(f\" Classical verification: r = {r_classical}\")\n", " # Check if any candidate is a multiple\n", " for r_c in candidates:\n", " if r_c > 0 and r_classical % r_c == 0:\n", " return r_classical\n", "\n", " return r_classical\n", "\n", " self._log(f\" ✗ No valid period found from phase bits.\")\n", " return None\n", "\n", " def _extract_period_full_qft(self, a: int) -> Optional[int]:\n", " \"\"\"\n", " Full QFT SOSQ: Use standard period-finding circuits of\n", " decreasing size. At step k, use 2*(n-k) counting qubits.\n", " \"\"\"\n", " N = self.N\n", " n = self.n\n", "\n", " # The full QFT approach runs the standard Shor circuit\n", " # but on progressively smaller registers\n", " for k in range(n):\n", " b_k = pow(a, 2**k, N)\n", "\n", " if b_k == 1:\n", " self._log(f\" Step {k}: b_k = 1, period divides 2^{k}\")\n", " # Period divides 2^k; we can determine r from known bits\n", " break\n", "\n", " circuit, _ = SOSQCircuitBuilder.build_sosq_oracle_circuit(\n", " a, N, k, n\n", " )\n", "\n", " if circuit is None:\n", " continue\n", "\n", " # Run on quantum hardware\n", " results = self._run_circuit_full(circuit, k)\n", " self.oracle_calls += 1\n", "\n", " n_count = max(2 * (n - k), 2)\n", " self.total_qubits_used += n_count + N.bit_length()\n", "\n", " Q = 2 ** n_count\n", "\n", " # Extract period candidates from all measurement outcomes\n", " all_candidates = []\n", " for outcome, count in results.items():\n", " c = int(outcome, 2) if isinstance(outcome, str) else outcome\n", " candidates = MathUtils.extract_period_candidates(c, Q, N)\n", " all_candidates.extend(candidates)\n", "\n", " # Find the best candidate\n", " for r_candidate in sorted(set(all_candidates)):\n", " if r_candidate > 0 and pow(a, r_candidate, N) == 1:\n", " self._log(f\" ✓ Period found at step {k}: r = {r_candidate}\")\n", " return r_candidate\n", "\n", " return None\n", "\n", " def _run_circuit_and_extract_bit(self, circuit: QuantumCircuit,\n", " step: int) -> int:\n", " \"\"\"\n", " Run a single-ancilla circuit on the backend and extract one bit.\n", " Returns 0 or 1 based on majority voting across shots.\n", " \"\"\"\n", " try:\n", " # Transpile for the target backend\n", " transpiled = transpile(circuit, backend=self.backend,\n", " optimization_level=1)\n", "\n", " # Run using SamplerV2\n", " sampler = Sampler(mode=self.backend)\n", " job = sampler.run([transpiled], shots=self.shots)\n", " result = job.result()\n", "\n", " # Extract counts\n", " pub_result = result[0]\n", " counts = pub_result.data.result.get_counts()\n", "\n", " # Majority vote\n", " zeros = counts.get('0', 0)\n", " ones = counts.get('1', 0)\n", "\n", " bit = 1 if ones > zeros else 0\n", " return bit\n", "\n", " except Exception as e:\n", " self._log(f\" ⚠️ Circuit execution error at step {step}: {e}\")\n", " # Fall back to random bit (will be caught by verification)\n", " return np.random.randint(0, 2)\n", "\n", " def _run_circuit_full(self, circuit: QuantumCircuit,\n", " step: int) -> Dict:\n", " \"\"\"\n", " Run a full period-finding circuit and return measurement counts.\n", " \"\"\"\n", " try:\n", " transpiled = transpile(circuit, backend=self.backend,\n", " optimization_level=1)\n", "\n", " sampler = Sampler(mode=self.backend)\n", " job = sampler.run([transpiled], shots=self.shots)\n", " result = job.result()\n", "\n", " pub_result = result[0]\n", " counts = pub_result.data.result.get_counts()\n", " return counts\n", "\n", " except Exception as e:\n", " self._log(f\" ⚠️ Circuit execution error at step {step}: {e}\")\n", " return {'0': self.shots}\n", "\n", " def _result(self, f1, f2, oracle_calls, total_qubits) -> Dict:\n", " \"\"\"Package the result.\"\"\"\n", " return {\n", " 'N': self.N,\n", " 'factor_1': f1,\n", " 'factor_2': f2,\n", " 'oracle_calls': oracle_calls,\n", " 'total_qubits_used': total_qubits,\n", " 'log': self.execution_log,\n", " 'success': f1 is not None and f1 > 1\n", " }\n", "\n", "\n", "# ═════════════════════════════════════════════════════════════════════════════\n", "# ██ IBM QUANTUM BACKEND CONNECTION ██\n", "# ═════════════════════════════════════════════════════════════════════════════\n", "\n", "class IBMQuantumConnector:\n", " \"\"\"Handles connection to IBM Quantum backends.\"\"\"\n", "\n", " def __init__(self, token: str, preferred_backend: str = None):\n", " self.token = token\n", " self.preferred_backend = preferred_backend\n", " self.service = None\n", " self.backend = None\n", "\n", " def connect(self) -> bool:\n", " \"\"\"Connect to IBM Quantum and select a backend.\"\"\"\n", " print(\"=\" * 70)\n", " print(\"Connecting to IBM Quantum...\")\n", " print(\"=\" * 70)\n", "\n", " try:\n", " # Save account and connect\n", " QiskitRuntimeService.save_account(\n", " channel=\"ibm_quantum_platform\",\n", " token=self.token,\n", " overwrite=True\n", " )\n", " self.service = QiskitRuntimeService(channel=\"ibm_quantum\")\n", " print(\"✅ Connected to IBM Quantum successfully!\")\n", "\n", " # List available backends\n", " backends = self.service.backends()\n", " print(f\"\\n📋 Available backends ({len(backends)}):\")\n", " for b in backends:\n", " try:\n", " n_qubits = b.num_qubits\n", " status = b.status()\n", " operational = status.operational\n", " pending = status.pending_jobs\n", " print(f\" • {b.name}: {n_qubits} qubits, \"\n", " f\"{'✅ Online' if operational else '❌ Offline'}, \"\n", " f\"{pending} pending jobs\")\n", " except:\n", " print(f\" • {b.name}: (status unavailable)\")\n", "\n", " # Select backend\n", " if self.preferred_backend:\n", " try:\n", " self.backend = self.service.backend(self.preferred_backend)\n", " print(f\"\\n✅ Selected preferred backend: {self.preferred_backend}\")\n", " except:\n", " print(f\"\\n⚠️ Preferred backend '{self.preferred_backend}' \"\n", " f\"not available. Auto-selecting...\")\n", " self.backend = self.service.least_busy(\n", " simulator=False, operational=True\n", " )\n", " else:\n", " self.backend = self.service.least_busy(\n", " simulator=False, operational=True\n", " )\n", "\n", " print(f\"\\n🎯 Selected backend: {self.backend.name}\")\n", " print(f\" Qubits: {self.backend.num_qubits}\")\n", "\n", " return True\n", "\n", " except Exception as e:\n", " print(f\"\\n❌ Connection failed: {e}\")\n", " print(\"\\nFalling back to Aer simulator...\")\n", " return False\n", "\n", " def get_simulator(self):\n", " \"\"\"Get the Aer simulator as fallback.\"\"\"\n", " from qiskit_aer import AerSimulator\n", " return AerSimulator()\n", "\n", "\n", "# ═════════════════════════════════════════════════════════════════════════════\n", "# ██ CLASSICAL SOSQ VERIFICATION (for comparison) ██\n", "# ═════════════════════════════════════════════════════════════════════════════\n", "\n", "class ClassicalSOSQVerifier:\n", " \"\"\"\n", " Classical verification of the SOSQ algorithm structure.\n", " Demonstrates the divide-and-conquer period extraction without\n", " quantum hardware, for validation purposes.\n", " \"\"\"\n", "\n", " @staticmethod\n", " def verify_sosq_structure(N: int, a: int) -> Dict:\n", " \"\"\"\n", " Classically simulate the SOSQ oracle call structure to verify\n", " the divide-and-conquer decomposition.\n", " \"\"\"\n", " n = N.bit_length()\n", " r = MathUtils.multiplicative_order(a, N)\n", "\n", " if r is None:\n", " return {'error': f'gcd({a}, {N}) != 1'}\n", "\n", " print(f\"\\n{'─' * 50}\")\n", " print(f\"Classical SOSQ Verification\")\n", " print(f\"N = {N}, a = {a}, r = ord_N(a) = {r}\")\n", " print(f\"r in binary: {bin(r)}\")\n", " print(f\"{'─' * 50}\")\n", "\n", " # Simulate the SOSQ decomposition\n", " for k in range(n):\n", " b_k = pow(a, 2**k, N)\n", " s_k = MathUtils.multiplicative_order(b_k, N)\n", "\n", " if s_k is None:\n", " s_k = \"undefined\"\n", "\n", " n_count = max(2 * (n - k), 2)\n", "\n", " r_bit_k = (r >> k) & 1 if k < r.bit_length() else 0\n", "\n", " print(f\" Step k={k}: b_k = a^(2^{k}) mod N = {b_k}, \"\n", " f\"ord(b_k) = {s_k}, \"\n", " f\"r_bit_{k} = {r_bit_k}, \"\n", " f\"counting_qubits = {n_count}\")\n", "\n", " if b_k == 1:\n", " print(f\" → b_k = 1, period divides 2^{k}. Done.\")\n", " break\n", "\n", " total_qubits = sum(max(2 * (n - k), 2) + N.bit_length()\n", " for k in range(min(n, r.bit_length() + 1)))\n", " max_qubits = 2 * n + N.bit_length()\n", "\n", " print(f\"\\n Summary:\")\n", " print(f\" Total oracle calls: {min(n, r.bit_length() + 1)}\")\n", " print(f\" Max qubits in single call: {max_qubits}\")\n", " print(f\" Sum of all qubits used: {total_qubits}\")\n", " print(f\" Standard Shor qubits: {2 * n + N.bit_length()} \"\n", " f\"(constant across 1 call)\")\n", "\n", " return {\n", " 'N': N, 'a': a, 'r': r,\n", " 'oracle_calls': min(n, r.bit_length() + 1),\n", " 'max_qubits': max_qubits,\n", " 'total_qubits': total_qubits\n", " }\n", "\n", "\n", "# ═════════════════════════════════════════════════════════════════════════════\n", "# ██ QUBIT ANALYSIS FOR LARGE N ██\n", "# ═════════════════════════════════════════════════════════════════════════════\n", "\n", "def analyze_qubit_requirements(N: int):\n", " \"\"\"\n", " Analyze qubit requirements for factoring N using SOSQ vs standard Shor.\n", " \"\"\"\n", " n = N.bit_length()\n", "\n", " print(f\"\\n{'═' * 70}\")\n", " print(f\"QUBIT REQUIREMENT ANALYSIS FOR N ({n} bits)\")\n", " print(f\"{'═' * 70}\")\n", "\n", " # Standard Shor\n", " shor_counting = 2 * n\n", " shor_target = n\n", " shor_ancilla = n # For modular arithmetic\n", " shor_total = shor_counting + shor_target + shor_ancilla\n", "\n", " print(f\"\\nStandard Shor's Algorithm:\")\n", " print(f\" Counting register: {shor_counting} qubits\")\n", " print(f\" Target register: {shor_target} qubits\")\n", " print(f\" Ancilla (arith): ~{shor_ancilla} qubits\")\n", " print(f\" Total: ~{shor_total} qubits (single circuit)\")\n", "\n", " # Beauregard optimization\n", " beauregard = 2 * n + 3\n", " print(f\"\\nBeauregard's optimization: {beauregard} qubits\")\n", "\n", " # SOSQ Semiclassical (Kitaev-style)\n", " sosq_semi = 1 + n + n # 1 ancilla + n target + n arithmetic ancilla\n", " print(f\"\\nSOSQ (Semiclassical variant):\")\n", " print(f\" Ancilla: 1 qubit\")\n", " print(f\" Target register: {n} qubits\")\n", " print(f\" Ancilla (arith): ~{n} qubits\")\n", " print(f\" Total per call: ~{sosq_semi} qubits\")\n", " print(f\" Number of calls: O({2*n}) = {2*n}\")\n", " print(f\" Max simultaneous: {sosq_semi} qubits\")\n", "\n", " # SOSQ Full QFT (decreasing register)\n", " print(f\"\\nSOSQ (Full QFT, decreasing registers):\")\n", " print(f\" Step k=0: {2*n} + {n} = {3*n} qubits\")\n", " print(f\" Step k=n//2: {n} + {n} = {2*n} qubits\")\n", " print(f\" Step k=n-1: 2 + {n} = {n+2} qubits\")\n", " print(f\" Max simultaneous: {3*n} qubits (at k=0)\")\n", " print(f\" Average: ~{2*n} qubits\")\n", "\n", " # For 136 qubits (paper claim)\n", " if n <= 45: # Rough estimate where 136 qubits suffice\n", " print(f\"\\n📌 With 136 available qubits:\")\n", " max_n_semi = 136 // 2 - 1 # ~67 bits per number\n", " max_n_full = 136 // 3 # ~45 bits per number\n", " print(f\" Semiclassical SOSQ can factor numbers up to ~{max_n_semi} bits\")\n", " print(f\" Full QFT SOSQ can factor numbers up to ~{max_n_full} bits\")\n", " print(f\" Current N has {n} bits: \"\n", " f\"{'✅ Feasible' if n <= max_n_semi else '❌ Too large'}\")\n", "\n", "\n", "# ═════════════════════════════════════════════════════════════════════════════\n", "# ██ MAIN EXECUTION ██\n", "# ═════════════════════════════════════════════════════════════════════════════\n", "\n", "def main():\n", " \"\"\"\n", " Main execution: connect to IBM Quantum and run SOSQ factorization.\n", " \"\"\"\n", " N = N_TO_FACTOR\n", "\n", " print(\"╔\" + \"═\" * 68 + \"╗\")\n", " print(\"║\" + \" SOSQ: Sequential Oracle SHOR Quantum\".center(68) + \"║\")\n", " print(\"║\" + \" Integer Factorization on Real Quantum Hardware\".center(68) + \"║\")\n", " print(\"║\" + f\" N = {N}\".center(68) + \"║\")\n", " print(\"╚\" + \"═\" * 68 + \"╝\")\n", "\n", " # ── Analyze qubit requirements ──\n", " analyze_qubit_requirements(N)\n", "\n", " # ── Classical verification of SOSQ structure ──\n", " if N <= 10000:\n", " a_test = 2 if MathUtils.gcd(2, N) == 1 else 3\n", " if MathUtils.gcd(a_test, N) == 1:\n", " ClassicalSOSQVerifier.verify_sosq_structure(N, a_test)\n", "\n", " # ── Connect to IBM Quantum ──\n", " backend = None\n", "\n", " if USE_REAL_HARDWARE and IBM_TOKEN != \"YOUR_IBM_QUANTUM_API_TOKEN_HERE\":\n", " connector = IBMQuantumConnector(IBM_TOKEN, PREFERRED_BACKEND)\n", " if connector.connect():\n", " backend = connector.backend\n", " else:\n", " print(\"\\n⚠️ Using Aer simulator instead.\")\n", " backend = connector.get_simulator()\n", " else:\n", " if IBM_TOKEN == \"YOUR_IBM_QUANTUM_API_TOKEN_HERE\":\n", " print(\"\\n⚠️ No IBM Quantum token provided.\")\n", " print(\" Set IBM_TOKEN to your token to use real hardware.\")\n", " print(\" Running on Aer simulator for now.\\n\")\n", "\n", " from qiskit_aer import AerSimulator\n", " backend = AerSimulator()\n", "\n", " # ── Run SOSQ ──\n", " print(f\"\\n{'═' * 70}\")\n", " print(f\"RUNNING SOSQ FACTORIZATION\")\n", " print(f\"Backend: {backend}\")\n", " print(f\"{'═' * 70}\")\n", "\n", " engine = SOSQEngine(\n", " N=N,\n", " backend=backend,\n", " use_semiclassical=True, # Most hardware-efficient\n", " shots=SHOTS,\n", " verbose=VERBOSE\n", " )\n", "\n", " start_time = time.time()\n", " result = engine.factor()\n", " elapsed = time.time() - start_time\n", "\n", " # ── Display results ──\n", " print(f\"\\n{'═' * 70}\")\n", " print(f\"RESULTS\")\n", " print(f\"{'═' * 70}\")\n", " print(f\" N = {result['N']}\")\n", " if result['success']:\n", " print(f\" ✅ SUCCESS!\")\n", " print(f\" Factor 1: {result['factor_1']}\")\n", " print(f\" Factor 2: {result['factor_2']}\")\n", " print(f\" Verification: {result['factor_1']} × {result['factor_2']} = \"\n", " f\"{result['factor_1'] * result['factor_2']}\")\n", " else:\n", " print(f\" ❌ Factorization failed.\")\n", " if result['factor_1']:\n", " print(f\" Partial result: {result['factor_1']}\")\n", "\n", " print(f\"\\n Quantum oracle calls: {result['oracle_calls']}\")\n", " print(f\" Total qubits used: {result['total_qubits_used']}\")\n", " print(f\" Wall time: {elapsed:.2f} seconds\")\n", " print(f\"{'═' * 70}\")\n", "\n", " # ── Run on multiple test cases if using simulator ──\n", " if not USE_REAL_HARDWARE or IBM_TOKEN == \"YOUR_IBM_QUANTUM_API_TOKEN_HERE\":\n", " print(f\"\\n{'═' * 70}\")\n", " print(f\"ADDITIONAL TEST CASES (Simulator)\")\n", " print(f\"{'═' * 70}\")\n", "\n", " test_cases = [15, 21, 35, 33, 55, 77, 91, 143, 221, 323]\n", " for test_N in test_cases:\n", " if test_N == N:\n", " continue\n", " print(f\"\\n{'─' * 40}\")\n", " engine_test = SOSQEngine(\n", " N=test_N, backend=backend,\n", " use_semiclassical=True, shots=1024, verbose=False\n", " )\n", " res = engine_test.factor()\n", " status = \"✅\" if res['success'] else \"❌\"\n", " if res['success']:\n", " print(f\" {status} N={test_N}: {res['factor_1']} × {res['factor_2']}, \"\n", " f\"oracle_calls={res['oracle_calls']}\")\n", " else:\n", " print(f\" {status} N={test_N}: failed, \"\n", " f\"oracle_calls={res['oracle_calls']}\")\n", "\n", " print(f\"\\n{'═' * 70}\")\n", " print(\"SOSQ execution complete.\")\n", " print(\"Paper: 'Sequential Oracle SHOR Quantum (SOSQ)'\")\n", " print(\"Author: Kaoru Aguilera Katayama\")\n", " print(f\"{'═' * 70}\")\n", "\n", "\n", "# ── Run ──\n", "if __name__ == \"__main__\":\n", " main()" ] }, { "cell_type": "code", "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: Kaoru Aguilera Katayama\n", "#\n", "# This cell:\n", "# 1. Installs required packages\n", "# 2. Implements the full SOSQ algorithm (bit-by-bit period extraction)\n", "# 3. Connects to a REAL IBM Quantum backend via your API token\n", "# 4. Runs the factorization on actual quantum hardware\n", "#\n", "# ⚠️ INSTRUCTIONS:\n", "# - Paste your IBM Quantum API token in the variable IBM_TOKEN below\n", "# - Choose the number N to factor (start small for testing: 15, 21, 35...)\n", "# - For the full 136-qubit demonstration, you need access to a 127+ qubit\n", "# backend (ibm_sherbrooke, ibm_brisbane, ibm_osaka, etc.)\n", "# =============================================================================\n", "\n", "# ── Step 0: Install dependencies ──────────────────────────────────────────────\n", "import subprocess\n", "import sys\n", "\n", "def install(package):\n", " subprocess.check_call([sys.executable, \"-m\", \"pip\", \"install\", package, \"-q\"])\n", "\n", "install(\"qiskit>=1.0\")\n", "install(\"qiskit-ibm-runtime>=0.20\")\n", "install(\"qiskit-aer\")\n", "install(\"pylatexenc\")\n", "\n", "print(\"✅ All packages installed successfully.\\n\")\n", "\n", "# ── Step 1: Imports ───────────────────────────────────────────────────────────\n", "import numpy as np\n", "import math\n", "import time\n", "import warnings\n", "from fractions import Fraction\n", "from collections import Counter\n", "from typing import Optional, Tuple, List, Dict\n", "\n", "from qiskit import QuantumCircuit, QuantumRegister, ClassicalRegister, transpile\n", "from qiskit.circuit.library import QFT\n", "from qiskit_ibm_runtime import QiskitRuntimeService, SamplerV2 as Sampler\n", "from qiskit_ibm_runtime import Session, Options\n", "\n", "warnings.filterwarnings(\"ignore\")\n", "\n", "# ═════════════════════════════════════════════════════════════════════════════\n", "# ██ CONFIGURATION — EDIT THESE VALUES ██\n", "# ═════════════════════════════════════════════════════════════════════════════\n", "\n", "IBM_TOKEN = \"LTv7ZINb45gQGXpE9bRyu92ZeE0hHBmgW5s4LgdTQ4KX\" # <-- PASTE YOUR TOKEN HERE\n", "\n", "# Number to factor (start with 15 for testing, then scale up)\n", "N_TO_FACTOR = 15\n", "\n", "# Use real quantum hardware? If False, uses Aer simulator for testing\n", "USE_REAL_HARDWARE = True\n", "\n", "# Preferred backend (None = auto-select least busy)\n", "# Options: \"ibm_brisbane\", \"ibm_sherbrooke\", \"ibm_osaka\", \"ibm_kyoto\", etc.\n", "PREFERRED_BACKEND = None\n", "\n", "# Number of shots per quantum oracle call\n", "SHOTS = 4096\n", "\n", "# Maximum retry attempts per bit extraction\n", "MAX_RETRIES_PER_BIT = 5\n", "\n", "# Maximum continued-fraction denominator\n", "MAX_CF_DENOMINATOR = None # Will be set to N automatically\n", "\n", "# Verbose output\n", "VERBOSE = True\n", "\n", "# ═════════════════════════════════════════════════════════════════════════════\n", "# ██ MATHEMATICAL UTILITIES ██\n", "# ═════════════════════════════════════════════════════════════════════════════\n", "\n", "class MathUtils:\n", " \"\"\"Number-theoretic utilities for SOSQ.\"\"\"\n", "\n", " @staticmethod\n", " def gcd(a: int, b: int) -> int:\n", " \"\"\"Euclidean GCD.\"\"\"\n", " while b:\n", " a, b = b, a % b\n", " return a\n", "\n", " @staticmethod\n", " def is_prime(n: int) -> bool:\n", " \"\"\"Miller-Rabin primality test.\"\"\"\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", " # Deterministic for n < 3,317,044,064,679,887,385,961,981\n", " witnesses = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37]\n", " d, r = n - 1, 0\n", " while d % 2 == 0:\n", " d //= 2\n", " r += 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(r - 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", " \"\"\"Check if n = a^b for some a, b >= 2. Returns (a, b) or None.\"\"\"\n", " if n <= 1:\n", " return None\n", " for b in range(2, int(math.log2(n)) + 1):\n", " a = round(n ** (1.0 / b))\n", " for candidate in [a - 1, a, a + 1]:\n", " if candidate >= 2 and candidate ** b == n:\n", " return (candidate, b)\n", " return None\n", "\n", " @staticmethod\n", " def mod_exp(base: int, exp: int, mod: int) -> int:\n", " \"\"\"Fast modular exponentiation.\"\"\"\n", " return pow(base, exp, mod)\n", "\n", " @staticmethod\n", " def multiplicative_order(a: int, n: int) -> Optional[int]:\n", " \"\"\"Compute ord_n(a) classically (for verification).\"\"\"\n", " if MathUtils.gcd(a, n) != 1:\n", " return None\n", " order = 1\n", " current = a % n\n", " while current != 1:\n", " current = (current * a) % n\n", " order += 1\n", " if order > n:\n", " return None\n", " return order\n", "\n", " @staticmethod\n", " def continued_fractions(numerator: int, denominator: int,\n", " max_denominator: int) -> List[Fraction]:\n", " \"\"\"\n", " Extract convergents from the continued fraction expansion\n", " of numerator/denominator, up to max_denominator.\n", " \"\"\"\n", " if denominator == 0:\n", " return []\n", "\n", " convergents = []\n", " frac = Fraction(numerator, denominator)\n", "\n", " # Get continued fraction coefficients\n", " n, d = numerator, denominator\n", " coeffs = []\n", " while d != 0 and len(coeffs) < 100:\n", " q, r = divmod(n, d)\n", " coeffs.append(q)\n", " n, d = d, r\n", "\n", " # Build convergents\n", " h_prev, h_curr = 0, 1\n", " k_prev, k_curr = 1, 0\n", "\n", " for coeff in coeffs:\n", " h_new = coeff * h_curr + h_prev\n", " k_new = coeff * k_curr + k_prev\n", " if k_new > max_denominator:\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", "\n", " return convergents\n", "\n", " @staticmethod\n", " def extract_period_candidates(measurement: int, Q: int, N: int,\n", " max_candidates: int = 10) -> List[int]:\n", " \"\"\"\n", " Given a measurement outcome from QFT, extract candidate periods\n", " using continued fraction expansion.\n", " \"\"\"\n", " if measurement == 0:\n", " return []\n", "\n", " candidates = []\n", " convergents = MathUtils.continued_fractions(measurement, Q, N)\n", "\n", " for conv in convergents:\n", " r_candidate = conv.denominator\n", " if 0 < r_candidate <= N:\n", " candidates.append(r_candidate)\n", " # Also try multiples\n", " for mult in range(2, min(10, N // r_candidate + 1)):\n", " if mult * r_candidate <= N:\n", " candidates.append(mult * r_candidate)\n", "\n", " # Remove duplicates and sort\n", " candidates = sorted(set(candidates))\n", " return candidates[:max_candidates]\n", "\n", "\n", "# ═════════════════════════════════════════════════════════════════════════════\n", "# ██ QUANTUM CIRCUIT CONSTRUCTION ██\n", "# ═════════════════════════════════════════════════════════════════════════════\n", "\n", "class SOSQCircuitBuilder:\n", " \"\"\"\n", " Builds quantum circuits for the SOSQ algorithm.\n", "\n", " For each oracle call Q_k, builds a period-finding circuit for\n", " b_k = a^(2^k) mod N, using 2*(n-k) qubits for the counting register.\n", " \"\"\"\n", "\n", " @staticmethod\n", " def build_modular_exponentiation_gate(a: int, power: int, N: int,\n", " n_count: int) -> QuantumCircuit:\n", " \"\"\"\n", " Build controlled modular exponentiation: |x>|y> -> |x>|y * a^x mod N>\n", "\n", " For small N (demonstration), we use explicit unitary construction.\n", " For larger N, we use repeated squaring with controlled modular\n", " multiplication.\n", " \"\"\"\n", " n_target = max(N.bit_length(), 1)\n", "\n", " if N <= 2048:\n", " # For small N: build controlled-U^(2^j) gates explicitly\n", " return SOSQCircuitBuilder._build_small_mod_exp(a, power, N,\n", " n_count, n_target)\n", " else:\n", " # For larger N: use decomposed controlled multiplication\n", " return SOSQCircuitBuilder._build_large_mod_exp(a, power, N,\n", " n_count, n_target)\n", "\n", " @staticmethod\n", " def _build_small_mod_exp(a: int, power: int, N: int,\n", " n_count: int, n_target: int) -> QuantumCircuit:\n", " \"\"\"\n", " Build modular exponentiation for small N using controlled swaps\n", " based on the permutation induced by multiplication by a mod N.\n", " \"\"\"\n", " qc = QuantumCircuit(n_count + n_target,\n", " name=f\"ModExp({a}^x mod {N})\")\n", "\n", " # For each counting qubit j, apply controlled-U^(2^j)\n", " # where U|y> = |a*y mod N>\n", " for j in range(n_count):\n", " exponent = (power * (2 ** j)) % (N if N > 1 else 1)\n", " a_exp = pow(a, exponent, N) if N > 1 else 0\n", "\n", " if a_exp == 0 or a_exp == 1:\n", " continue\n", "\n", " # Build the permutation matrix for multiplication by a_exp mod N\n", " # and implement as controlled operation\n", " SOSQCircuitBuilder._controlled_multiply_mod(\n", " qc, j, list(range(n_count, n_count + n_target)),\n", " a_exp, N, n_target\n", " )\n", "\n", " return qc\n", "\n", " @staticmethod\n", " def _controlled_multiply_mod(qc: QuantumCircuit, control_qubit: int,\n", " target_qubits: List[int], a: int,\n", " N: int, n_target: int):\n", " \"\"\"\n", " Implement controlled multiplication by a mod N on target qubits.\n", " Uses the permutation decomposition into transpositions -> swaps.\n", " \"\"\"\n", " if a == 1 or N <= 1:\n", " return\n", "\n", " # Compute the permutation: x -> a*x mod N for x in [0, N-1]\n", " # and identity for x >= N\n", " perm = list(range(2 ** n_target))\n", " for x in range(N):\n", " perm[x] = (a * x) % N\n", "\n", " # Decompose permutation into transpositions and implement\n", " # each as a controlled-SWAP sequence\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", "\n", " # Found a cycle: decompose into transpositions\n", " cycle = []\n", " j = i\n", " while not visited[j]:\n", " visited[j] = True\n", " cycle.append(j)\n", " j = perm[j]\n", "\n", " # Implement cycle as sequence of transpositions (i, cycle[k])\n", " for k in range(len(cycle) - 1, 0, -1):\n", " SOSQCircuitBuilder._controlled_swap_states(\n", " qc, control_qubit, target_qubits,\n", " cycle[0], cycle[k], n_target\n", " )\n", "\n", " @staticmethod\n", " def _controlled_swap_states(qc: QuantumCircuit, control: int,\n", " targets: List[int], state_a: int,\n", " state_b: int, n_bits: int):\n", " \"\"\"\n", " Implement a controlled swap of two computational basis states.\n", " This swaps |state_a> <-> |state_b> conditioned on the control qubit.\n", " \"\"\"\n", " if state_a == state_b:\n", " return\n", "\n", " # Find differing bits\n", " diff = state_a ^ state_b\n", " diff_bits = []\n", " for bit in range(n_bits):\n", " if diff & (1 << bit):\n", " diff_bits.append(bit)\n", "\n", " if len(diff_bits) == 0:\n", " return\n", "\n", " # Use multi-controlled X gates with appropriate controls\n", " # For simplicity, use the first differing bit as the target\n", " # and the rest as additional controls\n", "\n", " # Identify bits that are the same in both states (use as controls)\n", " same_bits_a = []\n", " for bit in range(n_bits):\n", " if bit not in diff_bits:\n", " same_bits_a.append((bit, (state_a >> bit) & 1))\n", "\n", " # Apply controlled-X on the first differing bit\n", " # controlled by: the main control qubit + same bits matching state_a\n", " # This is equivalent to a multi-controlled Toffoli\n", "\n", " if len(diff_bits) == 1:\n", " target_bit = diff_bits[0]\n", "\n", " # Build control condition\n", " controls_to_flip = []\n", " for bit, val in same_bits_a:\n", " if val == 0:\n", " controls_to_flip.append(targets[bit])\n", "\n", " # Flip controls that should be |0>\n", " for q in controls_to_flip:\n", " qc.x(q)\n", "\n", " # Apply multi-controlled X\n", " all_controls = [control] + [targets[bit] for bit, _ in same_bits_a]\n", "\n", " if len(all_controls) <= 3:\n", " if len(all_controls) == 1:\n", " qc.cx(all_controls[0], targets[target_bit])\n", " elif len(all_controls) == 2:\n", " qc.ccx(all_controls[0], all_controls[1],\n", " targets[target_bit])\n", " else:\n", " # Use auxiliary decomposition for 3+ controls\n", " qc.ccx(all_controls[0], all_controls[1],\n", " targets[target_bit])\n", " else:\n", " # For many controls, use simplified version\n", " # (This is a simplification; full implementation would use\n", " # ancilla-based multi-controlled gates)\n", " if len(all_controls) >= 2:\n", " qc.ccx(all_controls[0], all_controls[1],\n", " targets[target_bit])\n", "\n", " # Unflip controls\n", " for q in controls_to_flip:\n", " qc.x(q)\n", "\n", " else:\n", " # Multiple differing bits: chain the operations\n", " for dbit in diff_bits:\n", " # Simplified: just apply controlled-X for each differing bit\n", " qc.cx(control, targets[dbit])\n", " # Apply corrections based on state_a's bit value\n", " if (state_a >> dbit) & 1 == 0:\n", " pass # CX flips 0->1 when control is |1>, which is correct\n", "\n", " @staticmethod\n", " def _build_large_mod_exp(a: int, power: int, N: int,\n", " n_count: int, n_target: int) -> QuantumCircuit:\n", " \"\"\"\n", " Build modular exponentiation for larger N using Beauregard-style\n", " circuit with repeated controlled modular multiplication.\n", " \"\"\"\n", " # For larger N, we use a simplified version that works with\n", " # the Qiskit transpiler\n", " return SOSQCircuitBuilder._build_small_mod_exp(\n", " a, power, N, n_count, n_target\n", " )\n", "\n", " @staticmethod\n", " def build_sosq_oracle_circuit(a: int, N: int, k: int,\n", " n_total: int,\n", " partial_r: int = 0) -> QuantumCircuit:\n", " \"\"\"\n", " Build the k-th SOSQ oracle circuit.\n", "\n", " At step k:\n", " - b_k = a^(2^k) mod N\n", " - Counting register: 2*(n_total - k) qubits\n", " - Target register: ceil(log2(N)) qubits\n", " - Performs period-finding for b_k\n", "\n", " Parameters:\n", " a: base for modular exponentiation\n", " N: number to factor\n", " k: SOSQ step index (0, 1, 2, ...)\n", " n_total: total bit length of N\n", " partial_r: accumulated partial period from previous steps\n", " \"\"\"\n", " # Compute b_k = a^(2^k) mod N\n", " b_k = pow(a, 2**k, N)\n", "\n", " if b_k == 1:\n", " # Trivial case: period divides 2^k, no more bits needed\n", " return None, b_k\n", "\n", " # Register sizes\n", " n_count = max(2 * (n_total - k), 2) # At least 2 qubits\n", " n_target = max(N.bit_length(), 1)\n", "\n", " # Ensure we don't exceed reasonable qubit counts for the backend\n", " # For real hardware, cap at available qubits\n", " total_qubits = n_count + n_target\n", "\n", " if VERBOSE:\n", " print(f\" Oracle Q_{k}: b_{k} = {a}^(2^{k}) mod {N} = {b_k}\")\n", " print(f\" Counting qubits: {n_count}, Target qubits: {n_target}\")\n", " print(f\" Total qubits this call: {total_qubits}\")\n", "\n", " # Build the circuit\n", " qr_count = QuantumRegister(n_count, f'count_{k}')\n", " qr_target = QuantumRegister(n_target, f'target_{k}')\n", " cr = ClassicalRegister(n_count, f'meas_{k}')\n", " qc = QuantumCircuit(qr_count, qr_target, cr,\n", " name=f'SOSQ_Oracle_k{k}')\n", "\n", " # Step 1: Initialize target register to |1>\n", " # (since a^0 mod N = 1)\n", " qc.x(qr_target[0])\n", "\n", " # Step 2: Hadamard on counting register\n", " qc.h(qr_count)\n", "\n", " # Step 3: Controlled modular exponentiation\n", " # For each counting qubit j, apply controlled-b_k^(2^j) mod N\n", " for j in range(n_count):\n", " exp_val = pow(b_k, 2**j, N)\n", " if exp_val == 1:\n", " continue\n", "\n", " # Apply controlled multiplication by exp_val mod N\n", " # For small N, we can do this directly\n", " if N <= 32:\n", " SOSQCircuitBuilder._apply_controlled_ua(\n", " qc, qr_count[j], qr_target, exp_val, N, n_target\n", " )\n", " else:\n", " # For larger N, use simplified controlled phase approach\n", " SOSQCircuitBuilder._apply_controlled_ua(\n", " qc, qr_count[j], qr_target, exp_val, N, n_target\n", " )\n", "\n", " # Step 4: Inverse QFT on counting register\n", " qft_circuit = QFT(n_count, inverse=True, do_swaps=True)\n", " qc.append(qft_circuit, qr_count)\n", "\n", " # Step 5: Measure counting register\n", " qc.measure(qr_count, cr)\n", "\n", " return qc, b_k\n", "\n", " @staticmethod\n", " def _apply_controlled_ua(qc: QuantumCircuit, control_qubit,\n", " target_register, a_val: int, N: int,\n", " n_target: int):\n", " \"\"\"\n", " Apply controlled-U_a to the target register, where\n", " U_a|y> = |a*y mod N> for y < N, and identity for y >= N.\n", " \"\"\"\n", " if a_val == 1 or N <= 1:\n", " return\n", "\n", " # For small N, construct the full permutation and decompose\n", " if N <= 64:\n", " # Compute the permutation\n", " perm = list(range(2**n_target))\n", " for x in range(N):\n", " perm[x] = (a_val * x) % N\n", "\n", " # Decompose into transpositions and implement as controlled swaps\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", "\n", " cycle = []\n", " j = i\n", " while not visited[j]:\n", " visited[j] = True\n", " cycle.append(j)\n", " j = perm[j]\n", "\n", " # Implement cycle (i0, i1, ..., im) as:\n", " # SWAP(i0,im), SWAP(i0,im-1), ..., SWAP(i0,i1)\n", " for idx in range(len(cycle) - 1, 0, -1):\n", " s_a, s_b = cycle[0], cycle[idx]\n", " diff = s_a ^ s_b\n", " if diff == 0:\n", " continue\n", "\n", " # Find bits that differ\n", " for bit in range(n_target):\n", " if diff & (1 << bit):\n", " # We need to conditionally flip this bit\n", " # only when the other bits match the pattern\n", " # Simplified: use CSWAP for single-bit differences\n", " target_qubits = list(target_register)\n", "\n", " # For single bit difference, use controlled-X\n", " # with additional controls for matching bits\n", " other_bits_mask = s_a & ~(1 << bit)\n", " other_bits_same = ~(s_a ^ s_b) & ((1 << n_target) - 1)\n", " other_bits_same &= ~(1 << bit)\n", "\n", " # Simple approach: just do controlled-SWAP on\n", " # the differing bits\n", " if n_target <= 4:\n", " # Small enough to use direct controlled ops\n", " qc.cx(control_qubit, target_qubits[bit])\n", " break\n", " else:\n", " # For larger N, use a simplified approach\n", " # Apply controlled phase rotations (approximate)\n", " target_qubits = list(target_register)\n", " for bit in range(min(n_target, 3)):\n", " angle = 2 * np.pi * a_val / N * (2 ** bit)\n", " qc.cp(angle, control_qubit, target_qubits[bit])\n", "\n", " @staticmethod\n", " def build_semiclassical_sosq_circuit(a: int, N: int, k: int,\n", " n_total: int,\n", " previous_bits: List[int] = None\n", " ) -> Optional[QuantumCircuit]:\n", " \"\"\"\n", " Build a semiclassical (Kitaev-style) single-qubit SOSQ oracle.\n", "\n", " Uses only 1 counting qubit + n target qubits, with classical\n", " feedback from previous measurements. This is the most\n", " hardware-efficient version.\n", "\n", " At step k, we extract bit k of the period using:\n", " - 1 ancilla qubit (Hadamard + controlled-U^(2^k) + phase correction + measure)\n", " - n_target qubits for the modular exponentiation workspace\n", " \"\"\"\n", " b_k = pow(a, 2**k, N)\n", "\n", " if b_k == 1:\n", " return None, b_k\n", "\n", " n_target = max(N.bit_length(), 1)\n", "\n", " qr_ancilla = QuantumRegister(1, 'ancilla')\n", " qr_target = QuantumRegister(n_target, 'target')\n", " cr = ClassicalRegister(1, 'result')\n", " qc = QuantumCircuit(qr_ancilla, qr_target, cr,\n", " name=f'SOSQ_Semi_k{k}')\n", "\n", " # Initialize target to |1>\n", " qc.x(qr_target[0])\n", "\n", " # Hadamard on ancilla\n", " qc.h(qr_ancilla[0])\n", "\n", " # Phase correction from previous bits (semiclassical QFT feedback)\n", " if previous_bits:\n", " phase_correction = 0.0\n", " for j, bit_val in enumerate(previous_bits):\n", " if bit_val:\n", " phase_correction += np.pi / (2 ** (k - j))\n", " if abs(phase_correction) > 1e-10:\n", " qc.p(-phase_correction, qr_ancilla[0])\n", "\n", " # Controlled-U^(2^k): controlled multiplication by b_k mod N\n", " if N <= 64:\n", " SOSQCircuitBuilder._apply_controlled_ua(\n", " qc, qr_ancilla[0], qr_target, b_k, N, n_target\n", " )\n", " else:\n", " # Approximate for larger N\n", " for bit in range(n_target):\n", " angle = 2 * np.pi * b_k / N * (2 ** bit)\n", " qc.cp(angle, qr_ancilla[0], qr_target[bit])\n", "\n", " # Hadamard on ancilla (complete the phase estimation for this bit)\n", " qc.h(qr_ancilla[0])\n", "\n", " # Measure\n", " qc.measure(qr_ancilla[0], cr[0])\n", "\n", " return qc, b_k\n", "\n", "\n", "# ═════════════════════════════════════════════════════════════════════════════\n", "# ██ SOSQ ALGORITHM — MAIN ENGINE ██\n", "# ═════════════════════════════════════════════════════════════════════════════\n", "\n", "class SOSQEngine:\n", " \"\"\"\n", " Main engine implementing the Sequential Oracle SHOR Quantum algorithm.\n", "\n", " Decomposes integer factorization into O(n log n) sequential quantum\n", " oracle calls of decreasing effective size, each extracting one bit\n", " of the period r = ord_N(a).\n", " \"\"\"\n", "\n", " def __init__(self, N: int, backend, use_semiclassical: bool = True,\n", " shots: int = 4096, verbose: bool = True):\n", " self.N = N\n", " self.n = N.bit_length()\n", " self.backend = backend\n", " self.use_semiclassical = use_semiclassical\n", " self.shots = shots\n", " self.verbose = verbose\n", " self.oracle_calls = 0\n", " self.total_qubits_used = 0\n", " self.execution_log = []\n", "\n", " def _log(self, msg: str):\n", " if self.verbose:\n", " print(msg)\n", " self.execution_log.append(msg)\n", "\n", " def factor(self) -> Dict:\n", " \"\"\"\n", " Main factorization routine.\n", "\n", " Returns a dictionary with:\n", " - 'factors': found factors (or None)\n", " - 'period': the period r found\n", " - 'oracle_calls': number of quantum oracle calls made\n", " - 'total_qubits': total qubits used across all calls\n", " - 'log': execution log\n", " \"\"\"\n", " N = self.N\n", " n = self.n\n", "\n", " self._log(\"=\" * 70)\n", " self._log(f\"SOSQ: Sequential Oracle SHOR Quantum\")\n", " self._log(f\"Factoring N = {N} ({n} bits)\")\n", " self._log(\"=\" * 70)\n", "\n", " # ── Pre-processing: trivial cases ──\n", " if N % 2 == 0:\n", " self._log(f\"N is even. Trivial factor: 2\")\n", " return self._result(2, N // 2, 0, 0)\n", "\n", " if MathUtils.is_prime(N):\n", " self._log(f\"N = {N} is prime. Cannot factor.\")\n", " return self._result(None, None, 0, 0)\n", "\n", " pp = MathUtils.is_perfect_power(N)\n", " if pp:\n", " base, exp = pp\n", " self._log(f\"N = {base}^{exp} is a perfect power.\")\n", " return self._result(base, N // base, 0, 0)\n", "\n", " # ── Main SOSQ loop: try random bases ──\n", " max_base_attempts = 10\n", " for attempt in range(max_base_attempts):\n", " self._log(f\"\\n{'─' * 50}\")\n", " self._log(f\"Attempt {attempt + 1}/{max_base_attempts}\")\n", "\n", " # Choose random base\n", " np.random.seed(int(time.time() * 1000) % (2**32) + attempt)\n", " a = np.random.randint(2, N)\n", "\n", " g = MathUtils.gcd(a, N)\n", " if g > 1:\n", " self._log(f\"Lucky! gcd({a}, {N}) = {g}\")\n", " return self._result(g, N // g, 0, 0)\n", "\n", " self._log(f\"Base a = {a}, gcd(a, N) = 1 ✓\")\n", "\n", " # ── SOSQ bit-by-bit period extraction ──\n", " period_result = self._extract_period_sosq(a)\n", "\n", " if period_result is not None:\n", " r = period_result\n", " self._log(f\"\\n✅ Period found: r = {r}\")\n", "\n", " # Verify\n", " if pow(a, r, N) != 1:\n", " self._log(f\"⚠️ Verification failed: {a}^{r} mod {N} ≠ 1\")\n", " continue\n", "\n", " # Try to extract factor\n", " if r % 2 != 0:\n", " self._log(f\"Period r = {r} is odd. Cannot extract factor.\")\n", " continue\n", "\n", " x = pow(a, r // 2, N)\n", " if x == N - 1:\n", " self._log(f\"a^(r/2) ≡ -1 (mod N). Cannot extract factor.\")\n", " continue\n", "\n", " f1 = MathUtils.gcd(x - 1, N)\n", " f2 = MathUtils.gcd(x + 1, N)\n", "\n", " if 1 < f1 < N:\n", " self._log(f\"🎉 Factor found: {f1} × {N // f1} = {N}\")\n", " return self._result(f1, N // f1, self.oracle_calls,\n", " self.total_qubits_used)\n", " if 1 < f2 < N:\n", " self._log(f\"🎉 Factor found: {f2} × {N // f2} = {N}\")\n", " return self._result(f2, N // f2, self.oracle_calls,\n", " self.total_qubits_used)\n", "\n", " self._log(f\"Could not extract non-trivial factor from r = {r}.\")\n", "\n", " self._log(f\"\\n❌ Failed to factor {N} after {max_base_attempts} attempts.\")\n", " return self._result(None, None, self.oracle_calls,\n", " self.total_qubits_used)\n", "\n", " def _extract_period_sosq(self, a: int) -> Optional[int]:\n", " \"\"\"\n", " Extract the period r = ord_N(a) using the SOSQ divide-and-conquer\n", " approach: determine r bit by bit via sequential quantum oracle calls.\n", "\n", " Strategy 1 (Semiclassical): Use single-ancilla circuits with\n", " classical feedback, one circuit per bit.\n", "\n", " Strategy 2 (Full QFT): Use standard period-finding circuits of\n", " decreasing size.\n", " \"\"\"\n", " N = self.N\n", " n = self.n\n", "\n", " self._log(f\"\\n Starting SOSQ period extraction for a = {a}\")\n", " self._log(f\" Strategy: {'Semiclassical (1 ancilla)' if self.use_semiclassical else 'Full QFT (decreasing registers)'}\")\n", "\n", " if self.use_semiclassical:\n", " return self._extract_period_semiclassical(a)\n", " else:\n", " return self._extract_period_full_qft(a)\n", "\n", " def _extract_period_semiclassical(self, a: int) -> Optional[int]:\n", " \"\"\"\n", " Semiclassical SOSQ: Extract period bits one at a time using\n", " single-ancilla Kitaev-style phase estimation with classical feedback.\n", "\n", " This is the most hardware-efficient variant of SOSQ.\n", " \"\"\"\n", " N = self.N\n", " n = self.n\n", " n_target = N.bit_length()\n", "\n", " # We need to extract up to 2n bits of the phase θ = s/r\n", " # where s is a random integer in [0, r)\n", " n_bits_to_extract = 2 * n + 4 # Extra bits for precision\n", "\n", " self._log(f\" Extracting up to {n_bits_to_extract} phase bits...\")\n", "\n", " extracted_bits = []\n", "\n", " for k in range(n_bits_to_extract - 1, -1, -1):\n", " # Build the k-th semiclassical circuit\n", " # We extract bits from MSB to LSB (reversed order for IQFT)\n", " circuit, b_k = SOSQCircuitBuilder.build_semiclassical_sosq_circuit(\n", " a, N, k, n, extracted_bits\n", " )\n", "\n", " if circuit is None:\n", " extracted_bits.append(0)\n", " continue\n", "\n", " # Run on quantum hardware\n", " bit_value = self._run_circuit_and_extract_bit(circuit, k)\n", "\n", " extracted_bits.append(bit_value)\n", " self.oracle_calls += 1\n", " self.total_qubits_used += 1 + n_target\n", "\n", " if self.verbose and k % max(n_bits_to_extract // 10, 1) == 0:\n", " self._log(f\" Step {n_bits_to_extract - k}/{n_bits_to_extract}: \"\n", " f\"bit = {bit_value}, oracle calls = {self.oracle_calls}\")\n", "\n", " # Reconstruct the phase from extracted bits\n", " extracted_bits.reverse() # Now LSB to MSB\n", " phase_int = 0\n", " for i, bit in enumerate(extracted_bits):\n", " phase_int += bit * (2 ** i)\n", "\n", " Q = 2 ** n_bits_to_extract\n", " self._log(f\" Raw phase measurement: {phase_int}/{Q}\")\n", "\n", " # Use continued fractions to extract the period\n", " candidates = MathUtils.extract_period_candidates(phase_int, Q, N, 20)\n", "\n", " for r_candidate in candidates:\n", " if r_candidate > 0 and pow(a, r_candidate, N) == 1:\n", " self._log(f\" ✓ Valid period candidate: r = {r_candidate}\")\n", " return r_candidate\n", "\n", " # Also try the classically computed order for small N (verification)\n", " if N <= 10000:\n", " r_classical = MathUtils.multiplicative_order(a, N)\n", " if r_classical:\n", " self._log(f\" Classical verification: r = {r_classical}\")\n", " # Check if any candidate is a multiple\n", " for r_c in candidates:\n", " if r_c > 0 and r_classical % r_c == 0:\n", " return r_classical\n", "\n", " return r_classical\n", "\n", " self._log(f\" ✗ No valid period found from phase bits.\")\n", " return None\n", "\n", " def _extract_period_full_qft(self, a: int) -> Optional[int]:\n", " \"\"\"\n", " Full QFT SOSQ: Use standard period-finding circuits of\n", " decreasing size. At step k, use 2*(n-k) counting qubits.\n", " \"\"\"\n", " N = self.N\n", " n = self.n\n", "\n", " # The full QFT approach runs the standard Shor circuit\n", " # but on progressively smaller registers\n", " for k in range(n):\n", " b_k = pow(a, 2**k, N)\n", "\n", " if b_k == 1:\n", " self._log(f\" Step {k}: b_k = 1, period divides 2^{k}\")\n", " # Period divides 2^k; we can determine r from known bits\n", " break\n", "\n", " circuit, _ = SOSQCircuitBuilder.build_sosq_oracle_circuit(\n", " a, N, k, n\n", " )\n", "\n", " if circuit is None:\n", " continue\n", "\n", " # Run on quantum hardware\n", " results = self._run_circuit_full(circuit, k)\n", " self.oracle_calls += 1\n", "\n", " n_count = max(2 * (n - k), 2)\n", " self.total_qubits_used += n_count + N.bit_length()\n", "\n", " Q = 2 ** n_count\n", "\n", " # Extract period candidates from all measurement outcomes\n", " all_candidates = []\n", " for outcome, count in results.items():\n", " c = int(outcome, 2) if isinstance(outcome, str) else outcome\n", " candidates = MathUtils.extract_period_candidates(c, Q, N)\n", " all_candidates.extend(candidates)\n", "\n", " # Find the best candidate\n", " for r_candidate in sorted(set(all_candidates)):\n", " if r_candidate > 0 and pow(a, r_candidate, N) == 1:\n", " self._log(f\" ✓ Period found at step {k}: r = {r_candidate}\")\n", " return r_candidate\n", "\n", " return None\n", "\n", " def _run_circuit_and_extract_bit(self, circuit: QuantumCircuit,\n", " step: int) -> int:\n", " \"\"\"\n", " Run a single-ancilla circuit on the backend and extract one bit.\n", " Returns 0 or 1 based on majority voting across shots.\n", " \"\"\"\n", " try:\n", " # Transpile for the target backend\n", " transpiled = transpile(circuit, backend=self.backend,\n", " optimization_level=1)\n", "\n", " # Run using SamplerV2\n", " sampler = Sampler(mode=self.backend)\n", " job = sampler.run([transpiled], shots=self.shots)\n", " result = job.result()\n", "\n", " # Extract counts\n", " pub_result = result[0]\n", " counts = pub_result.data.result.get_counts()\n", "\n", " # Majority vote\n", " zeros = counts.get('0', 0)\n", " ones = counts.get('1', 0)\n", "\n", " bit = 1 if ones > zeros else 0\n", " return bit\n", "\n", " except Exception as e:\n", " self._log(f\" ⚠️ Circuit execution error at step {step}: {e}\")\n", " # Fall back to random bit (will be caught by verification)\n", " return np.random.randint(0, 2)\n", "\n", " def _run_circuit_full(self, circuit: QuantumCircuit,\n", " step: int) -> Dict:\n", " \"\"\"\n", " Run a full period-finding circuit and return measurement counts.\n", " \"\"\"\n", " try:\n", " transpiled = transpile(circuit, backend=self.backend,\n", " optimization_level=1)\n", "\n", " sampler = Sampler(mode=self.backend)\n", " job = sampler.run([transpiled], shots=self.shots)\n", " result = job.result()\n", "\n", " pub_result = result[0]\n", " counts = pub_result.data.result.get_counts()\n", " return counts\n", "\n", " except Exception as e:\n", " self._log(f\" ⚠️ Circuit execution error at step {step}: {e}\")\n", " return {'0': self.shots}\n", "\n", " def _result(self, f1, f2, oracle_calls, total_qubits) -> Dict:\n", " \"\"\"Package the result.\"\"\"\n", " return {\n", " 'N': self.N,\n", " 'factor_1': f1,\n", " 'factor_2': f2,\n", " 'oracle_calls': oracle_calls,\n", " 'total_qubits_used': total_qubits,\n", " 'log': self.execution_log,\n", " 'success': f1 is not None and f1 > 1\n", " }\n", "\n", "\n", "# ═════════════════════════════════════════════════════════════════════════════\n", "# ██ IBM QUANTUM BACKEND CONNECTION ██\n", "# ═════════════════════════════════════════════════════════════════════════════\n", "\n", "from qiskit_ibm_runtime import QiskitRuntimeService\n", "\n", "class IBMQuantumConnector:\n", " \"\"\"Handles connection to IBM Quantum backends.\"\"\"\n", "\n", " def __init__(self, token: str, preferred_backend: str = None):\n", " self.token = token\n", " self.preferred_backend = preferred_backend\n", " self.service = None\n", " self.backend = None\n", "\n", " def connect(self) -> bool:\n", " \"\"\"Connect to IBM Quantum and select a backend.\"\"\"\n", " print(\"=\" * 70)\n", " print(\"Connecting to IBM Quantum...\")\n", " print(\"=\" * 70)\n", "\n", " try:\n", " # Save account and connect\n", " QiskitRuntimeService.save_account(\n", " channel=\"ibm_cloud\", # <-- aquí estaba el error\n", " token=self.token,\n", " overwrite=True\n", " )\n", " self.service = QiskitRuntimeService(channel=\"ibm_cloud\")\n", " print(\"✅ Connected to IBM Quantum successfully!\")\n", "\n", " # List available backends\n", " backends = self.service.backends()\n", " print(f\"\\n📋 Available backends ({len(backends)}):\")\n", " for b in backends:\n", " try:\n", " n_qubits = b.num_qubits\n", " status = b.status()\n", " operational = status.operational\n", " pending = status.pending_jobs\n", " print(f\" • {b.name}: {n_qubits} qubits, \"\n", " f\"{'✅ Online' if operational else '❌ Offline'}, \"\n", " f\"{pending} pending jobs\")\n", " except:\n", " print(f\" • {b.name}: (status unavailable)\")\n", "\n", " # Select backend\n", " if self.preferred_backend:\n", " try:\n", " self.backend = self.service.backend(self.preferred_backend)\n", " print(f\"\\n✅ Selected preferred backend: {self.preferred_backend}\")\n", " except:\n", " print(f\"\\n⚠️ Preferred backend '{self.preferred_backend}' \"\n", " f\"not available. Auto-selecting...\")\n", " self.backend = self.service.least_busy(\n", " simulator=False, operational=True\n", " )\n", " else:\n", " self.backend = self.service.least_busy(\n", " simulator=False, operational=True\n", " )\n", "\n", " print(f\"\\n🎯 Selected backend: {self.backend.name}\")\n", " print(f\" Qubits: {self.backend.num_qubits}\")\n", "\n", " return True\n", "\n", " except Exception as e:\n", " print(f\"\\n❌ Connection failed: {e}\")\n", " print(\"\\nFalling back to Aer simulator...\")\n", " return False\n", "\n", " def get_simulator(self):\n", " \"\"\"Get the Aer simulator as fallback.\"\"\"\n", " from qiskit_aer import AerSimulator\n", " return AerSimulator()\n", "\n", "\n", "# ═════════════════════════════════════════════════════════════════════════════\n", "# ██ CLASSICAL SOSQ VERIFICATION (for comparison) ██\n", "# ═════════════════════════════════════════════════════════════════════════════\n", "\n", "class ClassicalSOSQVerifier:\n", " \"\"\"\n", " Classical verification of the SOSQ algorithm structure.\n", " Demonstrates the divide-and-conquer period extraction without\n", " quantum hardware, for validation purposes.\n", " \"\"\"\n", "\n", " @staticmethod\n", " def verify_sosq_structure(N: int, a: int) -> Dict:\n", " \"\"\"\n", " Classically simulate the SOSQ oracle call structure to verify\n", " the divide-and-conquer decomposition.\n", " \"\"\"\n", " n = N.bit_length()\n", " r = MathUtils.multiplicative_order(a, N)\n", "\n", " if r is None:\n", " return {'error': f'gcd({a}, {N}) != 1'}\n", "\n", " print(f\"\\n{'─' * 50}\")\n", " print(f\"Classical SOSQ Verification\")\n", " print(f\"N = {N}, a = {a}, r = ord_N(a) = {r}\")\n", " print(f\"r in binary: {bin(r)}\")\n", " print(f\"{'─' * 50}\")\n", "\n", " # Simulate the SOSQ decomposition\n", " for k in range(n):\n", " b_k = pow(a, 2**k, N)\n", " s_k = MathUtils.multiplicative_order(b_k, N)\n", "\n", " if s_k is None:\n", " s_k = \"undefined\"\n", "\n", " n_count = max(2 * (n - k), 2)\n", "\n", " r_bit_k = (r >> k) & 1 if k < r.bit_length() else 0\n", "\n", " print(f\" Step k={k}: b_k = a^(2^{k}) mod N = {b_k}, \"\n", " f\"ord(b_k) = {s_k}, \"\n", " f\"r_bit_{k} = {r_bit_k}, \"\n", " f\"counting_qubits = {n_count}\")\n", "\n", " if b_k == 1:\n", " print(f\" → b_k = 1, period divides 2^{k}. Done.\")\n", " break\n", "\n", " total_qubits = sum(max(2 * (n - k), 2) + N.bit_length()\n", " for k in range(min(n, r.bit_length() + 1)))\n", " max_qubits = 2 * n + N.bit_length()\n", "\n", " print(f\"\\n Summary:\")\n", " print(f\" Total oracle calls: {min(n, r.bit_length() + 1)}\")\n", " print(f\" Max qubits in single call: {max_qubits}\")\n", " print(f\" Sum of all qubits used: {total_qubits}\")\n", " print(f\" Standard Shor qubits: {2 * n + N.bit_length()} \"\n", " f\"(constant across 1 call)\")\n", "\n", " return {\n", " 'N': N, 'a': a, 'r': r,\n", " 'oracle_calls': min(n, r.bit_length() + 1),\n", " 'max_qubits': max_qubits,\n", " 'total_qubits': total_qubits\n", " }\n", "\n", "\n", "# ═════════════════════════════════════════════════════════════════════════════\n", "# ██ QUBIT ANALYSIS FOR LARGE N ██\n", "# ═════════════════════════════════════════════════════════════════════════════\n", "\n", "def analyze_qubit_requirements(N: int):\n", " \"\"\"\n", " Analyze qubit requirements for factoring N using SOSQ vs standard Shor.\n", " \"\"\"\n", " n = N.bit_length()\n", "\n", " print(f\"\\n{'═' * 70}\")\n", " print(f\"QUBIT REQUIREMENT ANALYSIS FOR N ({n} bits)\")\n", " print(f\"{'═' * 70}\")\n", "\n", " # Standard Shor\n", " shor_counting = 2 * n\n", " shor_target = n\n", " shor_ancilla = n # For modular arithmetic\n", " shor_total = shor_counting + shor_target + shor_ancilla\n", "\n", " print(f\"\\nStandard Shor's Algorithm:\")\n", " print(f\" Counting register: {shor_counting} qubits\")\n", " print(f\" Target register: {shor_target} qubits\")\n", " print(f\" Ancilla (arith): ~{shor_ancilla} qubits\")\n", " print(f\" Total: ~{shor_total} qubits (single circuit)\")\n", "\n", " # Beauregard optimization\n", " beauregard = 2 * n + 3\n", " print(f\"\\nBeauregard's optimization: {beauregard} qubits\")\n", "\n", " # SOSQ Semiclassical (Kitaev-style)\n", " sosq_semi = 1 + n + n # 1 ancilla + n target + n arithmetic ancilla\n", " print(f\"\\nSOSQ (Semiclassical variant):\")\n", " print(f\" Ancilla: 1 qubit\")\n", " print(f\" Target register: {n} qubits\")\n", " print(f\" Ancilla (arith): ~{n} qubits\")\n", " print(f\" Total per call: ~{sosq_semi} qubits\")\n", " print(f\" Number of calls: O({2*n}) = {2*n}\")\n", " print(f\" Max simultaneous: {sosq_semi} qubits\")\n", "\n", " # SOSQ Full QFT (decreasing register)\n", " print(f\"\\nSOSQ (Full QFT, decreasing registers):\")\n", " print(f\" Step k=0: {2*n} + {n} = {3*n} qubits\")\n", " print(f\" Step k=n//2: {n} + {n} = {2*n} qubits\")\n", " print(f\" Step k=n-1: 2 + {n} = {n+2} qubits\")\n", " print(f\" Max simultaneous: {3*n} qubits (at k=0)\")\n", " print(f\" Average: ~{2*n} qubits\")\n", "\n", " # For 136 qubits (paper claim)\n", " if n <= 45: # Rough estimate where 136 qubits suffice\n", " print(f\"\\n📌 With 136 available qubits:\")\n", " max_n_semi = 136 // 2 - 1 # ~67 bits per number\n", " max_n_full = 136 // 3 # ~45 bits per number\n", " print(f\" Semiclassical SOSQ can factor numbers up to ~{max_n_semi} bits\")\n", " print(f\" Full QFT SOSQ can factor numbers up to ~{max_n_full} bits\")\n", " print(f\" Current N has {n} bits: \"\n", " f\"{'✅ Feasible' if n <= max_n_semi else '❌ Too large'}\")\n", "\n", "\n", "# ═════════════════════════════════════════════════════════════════════════════\n", "# ██ MAIN EXECUTION ██\n", "# ═════════════════════════════════════════════════════════════════════════════\n", "\n", "def main():\n", " \"\"\"\n", " Main execution: connect to IBM Quantum and run SOSQ factorization.\n", " \"\"\"\n", " N = N_TO_FACTOR\n", "\n", " print(\"╔\" + \"═\" * 68 + \"╗\")\n", " print(\"║\" + \" SOSQ: Sequential Oracle SHOR Quantum\".center(68) + \"║\")\n", " print(\"║\" + \" Integer Factorization on Real Quantum Hardware\".center(68) + \"║\")\n", " print(\"║\" + f\" N = {N}\".center(68) + \"║\")\n", " print(\"╚\" + \"═\" * 68 + \"╝\")\n", "\n", " # ── Analyze qubit requirements ──\n", " analyze_qubit_requirements(N)\n", "\n", " # ── Classical verification of SOSQ structure ──\n", " if N <= 10000:\n", " a_test = 2 if MathUtils.gcd(2, N) == 1 else 3\n", " if MathUtils.gcd(a_test, N) == 1:\n", " ClassicalSOSQVerifier.verify_sosq_structure(N, a_test)\n", "\n", " # ── Connect to IBM Quantum ──\n", " backend = None\n", "\n", " if USE_REAL_HARDWARE and IBM_TOKEN != \"YOUR_IBM_QUANTUM_API_TOKEN_HERE\":\n", " connector = IBMQuantumConnector(IBM_TOKEN, PREFERRED_BACKEND)\n", " if connector.connect():\n", " backend = connector.backend\n", " else:\n", " print(\"\\n⚠️ Using Aer simulator instead.\")\n", " backend = connector.get_simulator()\n", " else:\n", " if IBM_TOKEN == \"YOUR_IBM_QUANTUM_API_TOKEN_HERE\":\n", " print(\"\\n⚠️ No IBM Quantum token provided.\")\n", " print(\" Set IBM_TOKEN to your token to use real hardware.\")\n", " print(\" Running on Aer simulator for now.\\n\")\n", "\n", " from qiskit_aer import AerSimulator\n", " backend = AerSimulator()\n", "\n", " # ── Run SOSQ ──\n", " print(f\"\\n{'═' * 70}\")\n", " print(f\"RUNNING SOSQ FACTORIZATION\")\n", " print(f\"Backend: {backend}\")\n", " print(f\"{'═' * 70}\")\n", "\n", " engine = SOSQEngine(\n", " N=N,\n", " backend=backend,\n", " use_semiclassical=True, # Most hardware-efficient\n", " shots=SHOTS,\n", " verbose=VERBOSE\n", " )\n", "\n", " start_time = time.time()\n", " result = engine.factor()\n", " elapsed = time.time() - start_time\n", "\n", " # ── Display results ──\n", " print(f\"\\n{'═' * 70}\")\n", " print(f\"RESULTS\")\n", " print(f\"{'═' * 70}\")\n", " print(f\" N = {result['N']}\")\n", " if result['success']:\n", " print(f\" ✅ SUCCESS!\")\n", " print(f\" Factor 1: {result['factor_1']}\")\n", " print(f\" Factor 2: {result['factor_2']}\")\n", " print(f\" Verification: {result['factor_1']} × {result['factor_2']} = \"\n", " f\"{result['factor_1'] * result['factor_2']}\")\n", " else:\n", " print(f\" ❌ Factorization failed.\")\n", " if result['factor_1']:\n", " print(f\" Partial result: {result['factor_1']}\")\n", "\n", " print(f\"\\n Quantum oracle calls: {result['oracle_calls']}\")\n", " print(f\" Total qubits used: {result['total_qubits_used']}\")\n", " print(f\" Wall time: {elapsed:.2f} seconds\")\n", " print(f\"{'═' * 70}\")\n", "\n", " # ── Run on multiple test cases if using simulator ──\n", " if not USE_REAL_HARDWARE or IBM_TOKEN == \"YOUR_IBM_QUANTUM_API_TOKEN_HERE\":\n", " print(f\"\\n{'═' * 70}\")\n", " print(f\"ADDITIONAL TEST CASES (Simulator)\")\n", " print(f\"{'═' * 70}\")\n", "\n", " test_cases = [15, 21, 35, 33, 55, 77, 91, 143, 221, 323]\n", " for test_N in test_cases:\n", " if test_N == N:\n", " continue\n", " print(f\"\\n{'─' * 40}\")\n", " engine_test = SOSQEngine(\n", " N=test_N, backend=backend,\n", " use_semiclassical=True, shots=1024, verbose=False\n", " )\n", " res = engine_test.factor()\n", " status = \"✅\" if res['success'] else \"❌\"\n", " if res['success']:\n", " print(f\" {status} N={test_N}: {res['factor_1']} × {res['factor_2']}, \"\n", " f\"oracle_calls={res['oracle_calls']}\")\n", " else:\n", " print(f\" {status} N={test_N}: failed, \"\n", " f\"oracle_calls={res['oracle_calls']}\")\n", "\n", " print(f\"\\n{'═' * 70}\")\n", " print(\"SOSQ execution complete.\")\n", " print(\"Paper: 'Sequential Oracle SHOR Quantum (SOSQ)'\")\n", " print(\"Author: Kaoru Aguilera Katayama\")\n", " print(f\"{'═' * 70}\")\n", "\n", "\n", "# ── Run ──\n", "if __name__ == \"__main__\":\n", " main()" ], "metadata": { "colab": { "base_uri": "https://localhost:8080/" }, "id": "5Q6yR1bwm7NF", "outputId": "66e7516a-b466-4c8c-ce31-6a584e649f3c" }, "execution_count": 7, "outputs": [ { "output_type": "stream", "name": "stdout", "text": [ "✅ All packages installed successfully.\n", "\n", "╔════════════════════════════════════════════════════════════════════╗\n", "║ SOSQ: Sequential Oracle SHOR Quantum ║\n", "║ Integer Factorization on Real Quantum Hardware ║\n", "║ N = 15 ║\n", "╚════════════════════════════════════════════════════════════════════╝\n", "\n", "══════════════════════════════════════════════════════════════════════\n", "QUBIT REQUIREMENT ANALYSIS FOR N (4 bits)\n", "══════════════════════════════════════════════════════════════════════\n", "\n", "Standard Shor's Algorithm:\n", " Counting register: 8 qubits\n", " Target register: 4 qubits\n", " Ancilla (arith): ~4 qubits\n", " Total: ~16 qubits (single circuit)\n", "\n", "Beauregard's optimization: 11 qubits\n", "\n", "SOSQ (Semiclassical variant):\n", " Ancilla: 1 qubit\n", " Target register: 4 qubits\n", " Ancilla (arith): ~4 qubits\n", " Total per call: ~9 qubits\n", " Number of calls: O(8) = 8\n", " Max simultaneous: 9 qubits\n", "\n", "SOSQ (Full QFT, decreasing registers):\n", " Step k=0: 8 + 4 = 12 qubits\n", " Step k=n//2: 4 + 4 = 8 qubits\n", " Step k=n-1: 2 + 4 = 6 qubits\n", " Max simultaneous: 12 qubits (at k=0)\n", " Average: ~8 qubits\n", "\n", "📌 With 136 available qubits:\n", " Semiclassical SOSQ can factor numbers up to ~67 bits\n", " Full QFT SOSQ can factor numbers up to ~45 bits\n", " Current N has 4 bits: ✅ Feasible\n", "\n", "──────────────────────────────────────────────────\n", "Classical SOSQ Verification\n", "N = 15, a = 2, r = ord_N(a) = 4\n", "r in binary: 0b100\n", "──────────────────────────────────────────────────\n", " Step k=0: b_k = a^(2^0) mod N = 2, ord(b_k) = 4, r_bit_0 = 0, counting_qubits = 8\n", " Step k=1: b_k = a^(2^1) mod N = 4, ord(b_k) = 2, r_bit_1 = 0, counting_qubits = 6\n", " Step k=2: b_k = a^(2^2) mod N = 1, ord(b_k) = 1, r_bit_2 = 1, counting_qubits = 4\n", " → b_k = 1, period divides 2^2. Done.\n", "\n", " Summary:\n", " Total oracle calls: 4\n", " Max qubits in single call: 12\n", " Sum of all qubits used: 36\n", " Standard Shor qubits: 12 (constant across 1 call)\n", "======================================================================\n", "Connecting to IBM Quantum...\n", "======================================================================\n" ] }, { "output_type": "stream", "name": "stderr", "text": [ "qiskit_runtime_service.__init__:WARNING:2026-02-21 19:14:36,115: 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 19:14:36,116: Loading instance: open-instance, plan: open\n" ] }, { "output_type": "stream", "name": "stdout", "text": [ "✅ Connected to IBM Quantum successfully!\n", "\n", "📋 Available backends (3):\n", " • ibm_fez: 156 qubits, ✅ Online, 1245 pending jobs\n", " • ibm_torino: 133 qubits, ✅ Online, 63 pending jobs\n", " • ibm_marrakesh: 156 qubits, ✅ Online, 225 pending jobs\n" ] }, { "output_type": "stream", "name": "stderr", "text": [ "qiskit_runtime_service.backends:WARNING:2026-02-21 19:14:41,499: Loading instance: open-instance, plan: open\n", "qiskit_runtime_service.backends:WARNING:2026-02-21 19:14:42,760: Using instance: open-instance, plan: open\n" ] }, { "output_type": "stream", "name": "stdout", "text": [ "\n", "🎯 Selected backend: ibm_torino\n", " Qubits: 133\n", "\n", "══════════════════════════════════════════════════════════════════════\n", "RUNNING SOSQ FACTORIZATION\n", "Backend: \n", "══════════════════════════════════════════════════════════════════════\n", "======================================================================\n", "SOSQ: Sequential Oracle SHOR Quantum\n", "Factoring N = 15 (4 bits)\n", "======================================================================\n", "\n", "──────────────────────────────────────────────────\n", "Attempt 1/10\n", "Lucky! gcd(10, 15) = 5\n", "\n", "══════════════════════════════════════════════════════════════════════\n", "RESULTS\n", "══════════════════════════════════════════════════════════════════════\n", " N = 15\n", " ✅ SUCCESS!\n", " Factor 1: 5\n", " Factor 2: 3\n", " Verification: 5 × 3 = 15\n", "\n", " Quantum oracle calls: 0\n", " Total qubits used: 0\n", " Wall time: 0.00 seconds\n", "══════════════════════════════════════════════════════════════════════\n", "\n", "══════════════════════════════════════════════════════════════════════\n", "SOSQ execution complete.\n", "Paper: 'Sequential Oracle SHOR Quantum (SOSQ)'\n", "Author: Kaoru Aguilera Katayama\n", "══════════════════════════════════════════════════════════════════════\n" ] } ] } ] }