shor.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "f0f43643",
   "metadata": {},
   "outputs": [],
   "source": [
    "import fractions\n",
    "import math\n",
    "import random\n",
    "\n",
    "import numpy as np\n",
    "import sympy\n",
    "from typing import Callable, List, Optional, Sequence, Union\n",
    "\n",
    "import cirq"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "02a734e4",
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"Function to compute the elements of Z_n.\"\"\"\n",
    "def multiplicative_group(n: int) -> List[int]:\n",
    "    \"\"\"Returns the multiplicative group modulo n.\n",
    "\n",
    "    Args:\n",
    "        n: Modulus of the multiplicative group.\n",
    "    \"\"\"\n",
    "    assert n > 1\n",
    "    group = [1]\n",
    "    for x in range(2, n):\n",
    "        if math.gcd(x, n) == 1:\n",
    "            group.append(x)\n",
    "    return group"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "53ded347",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The multiplicative group modulo n = 15 is:\n",
      "[1, 2, 4, 7, 8, 11, 13, 14]\n"
     ]
    }
   ],
   "source": [
    "# multiplicative group modulo n = 15\n",
    "\"\"\"Example of a multiplicative group.\"\"\"\n",
    "n = 15\n",
    "print(f\"The multiplicative group modulo n = {n} is:\")\n",
    "print(multiplicative_group(n))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "b143006c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# classical\n",
    "\"\"\"Function for classically computing the order of an element of Z_n.\"\"\"\n",
    "def classical_order_finder(x: int, n: int) -> Optional[int]:\n",
    "    \"\"\"Computes smallest positive r such that x**r mod n == 1.\n",
    "\n",
    "    Args:\n",
    "        x: Integer whose order is to be computed, must be greater than one\n",
    "           and belong to the multiplicative group of integers modulo n (which\n",
    "           consists of positive integers relatively prime to n),\n",
    "        n: Modulus of the multiplicative group.\n",
    "\n",
    "    Returns:\n",
    "        Smallest positive integer r such that x**r == 1 mod n.\n",
    "        Always succeeds (and hence never returns None).\n",
    "\n",
    "    Raises:\n",
    "        ValueError when x is 1 or not an element of the multiplicative\n",
    "        group of integers modulo n.\n",
    "    \"\"\"\n",
    "    # Make sure x is both valid and in Z_n.\n",
    "    if x < 2 or x >= n or math.gcd(x, n) > 1:\n",
    "        raise ValueError(f\"Invalid x={x} for modulus n={n}.\")\n",
    "\n",
    "    # Determine the order.\n",
    "    r, y = 1, x\n",
    "    while y != 1:\n",
    "        y = (x * y) % n\n",
    "        r += 1\n",
    "    return r"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "8d2fae9a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x^r mod n = 8^4 mod 15 = 1\n"
     ]
    }
   ],
   "source": [
    "\"\"\"Example of (classically) computing the order of an element.\"\"\"\n",
    "n = 15  # The multiplicative group is [1, 2, 4, 7, 8, 11, 13, 14].\n",
    "x = 8\n",
    "r = classical_order_finder(x, n)\n",
    "\n",
    "# Check that the order is indeed correct.\n",
    "print(f\"x^r mod n = {x}^{r} mod {n} = {x**r % n}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "40c466d3",
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"Example of defining an arithmetic (quantum) operation in Cirq.\"\"\"\n",
    "class Adder(cirq.ArithmeticOperation):\n",
    "    \"\"\"Quantum addition.\"\"\"\n",
    "    def __init__(self, target_register, input_register):\n",
    "        self.input_register = input_register\n",
    "        self.target_register = target_register\n",
    "\n",
    "    def registers(self):\n",
    "        return self.target_register, self.input_register\n",
    "\n",
    "    def with_registers(self, *new_registers):\n",
    "        return Adder(*new_registers)\n",
    "\n",
    "    def apply(self, target_value, input_value):\n",
    "        return target_value + input_value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "3c72a24d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Circuit:\n",
      "\n",
      "0: ───X───#3──────────────────────────────────────────M───\n",
      "          │\n",
      "1: ───────#4──────────────────────────────────────────M───\n",
      "          │\n",
      "2: ───────<__main__.Adder object at 0x7f7859d223a0>───M───\n",
      "          │\n",
      "3: ───X───#2──────────────────────────────────────────M───\n",
      "\n",
      "\n",
      "Measurement outcomes:\n",
      "\n",
      "   0  1  2  3\n",
      "0  1  0  1  1\n",
      "1  1  0  1  1\n",
      "2  1  0  1  1\n",
      "3  1  0  1  1\n",
      "4  1  0  1  1\n"
     ]
    }
   ],
   "source": [
    "\"\"\"Example of using an Adder in a circuit.\"\"\"\n",
    "# Two qubit registers.\n",
    "qreg1 = cirq.LineQubit.range(2)\n",
    "qreg2 = cirq.LineQubit.range(2, 4)\n",
    "\n",
    "# Define the circuit.\n",
    "circ = cirq.Circuit(\n",
    "    cirq.ops.X.on(qreg1[0]),\n",
    "    cirq.ops.X.on(qreg2[1]),\n",
    "    Adder(input_register=qreg1, target_register=qreg2),\n",
    "    cirq.measure_each(*qreg1),\n",
    "    cirq.measure_each(*qreg2)\n",
    ")\n",
    "\n",
    "# Display it.\n",
    "print(\"Circuit:\\n\")\n",
    "print(circ)\n",
    "\n",
    "# Print the measurement outcomes.\n",
    "print(\"\\n\\nMeasurement outcomes:\\n\")\n",
    "print(cirq.sample(circ, repetitions=5).data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "d465b00f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0., 0., 0., 1.],\n",
       "       [1., 0., 0., 0.],\n",
       "       [0., 1., 0., 0.],\n",
       "       [0., 0., 1., 0.]])"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\"\"\"Example of the unitary of an Adder operation.\"\"\"\n",
    "cirq.unitary(\n",
    "    Adder(target_register=cirq.LineQubit.range(2),\n",
    "          input_register=1)\n",
    ").real"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "1a738d5c",
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"Defines the modular exponential operation used in Shor's algorithm.\"\"\"\n",
    "class ModularExp(cirq.ArithmeticOperation):\n",
    "    \"\"\"Quantum modular exponentiation.\n",
    "\n",
    "    This class represents the unitary which multiplies base raised to exponent\n",
    "    into the target modulo the given modulus. More precisely, it represents the\n",
    "    unitary V which computes modular exponentiation x**e mod n:\n",
    "\n",
    "        V|y⟩|e⟩ = |y * x**e mod n⟩ |e⟩     0 <= y < n\n",
    "        V|y⟩|e⟩ = |y⟩ |e⟩                  n <= y\n",
    "\n",
    "    where y is the target register, e is the exponent register, x is the base\n",
    "    and n is the modulus. Consequently,\n",
    "\n",
    "        V|y⟩|e⟩ = (U**e|y)|e⟩\n",
    "\n",
    "    where U is the unitary defined as\n",
    "\n",
    "        U|y⟩ = |y * x mod n⟩      0 <= y < n\n",
    "        U|y⟩ = |y⟩                n <= y\n",
    "    \"\"\"\n",
    "    def __init__(\n",
    "        self, \n",
    "        target: Sequence[cirq.Qid],\n",
    "        exponent: Union[int, Sequence[cirq.Qid]], \n",
    "        base: int,\n",
    "        modulus: int\n",
    "    ) -> None:\n",
    "        if len(target) < modulus.bit_length():\n",
    "            raise ValueError(f'Register with {len(target)} qubits is too small '\n",
    "                             f'for modulus {modulus}')\n",
    "        self.target = target\n",
    "        self.exponent = exponent\n",
    "        self.base = base\n",
    "        self.modulus = modulus\n",
    "\n",
    "    def registers(self) -> Sequence[Union[int, Sequence[cirq.Qid]]]:\n",
    "        return self.target, self.exponent, self.base, self.modulus\n",
    "\n",
    "    def with_registers(\n",
    "            self,\n",
    "            *new_registers: Union[int, Sequence['cirq.Qid']],\n",
    "    ) -> cirq.ArithmeticOperation:\n",
    "        if len(new_registers) != 4:\n",
    "            raise ValueError(f'Expected 4 registers (target, exponent, base, '\n",
    "                             f'modulus), but got {len(new_registers)}')\n",
    "        target, exponent, base, modulus = new_registers\n",
    "        if not isinstance(target, Sequence):\n",
    "            raise ValueError(\n",
    "                f'Target must be a qubit register, got {type(target)}')\n",
    "        if not isinstance(base, int):\n",
    "            raise ValueError(\n",
    "                f'Base must be a classical constant, got {type(base)}')\n",
    "        if not isinstance(modulus, int):\n",
    "            raise ValueError(\n",
    "                f'Modulus must be a classical constant, got {type(modulus)}')\n",
    "        return ModularExp(target, exponent, base, modulus)\n",
    "\n",
    "    def apply(self, *register_values: int) -> int:\n",
    "        assert len(register_values) == 4\n",
    "        target, exponent, base, modulus = register_values\n",
    "        if target >= modulus:\n",
    "            return target\n",
    "        return (target * base**exponent) % modulus\n",
    "\n",
    "    def _circuit_diagram_info_(\n",
    "            self,\n",
    "            args: cirq.CircuitDiagramInfoArgs,\n",
    "    ) -> cirq.CircuitDiagramInfo:\n",
    "        assert args.known_qubits is not None\n",
    "        wire_symbols: List[str] = []\n",
    "        t, e = 0, 0\n",
    "        for qubit in args.known_qubits:\n",
    "            if qubit in self.target:\n",
    "                if t == 0:\n",
    "                    if isinstance(self.exponent, Sequence):\n",
    "                        e_str = 'e'\n",
    "                    else:\n",
    "                        e_str = str(self.exponent)\n",
    "                    wire_symbols.append(\n",
    "                        f'ModularExp(t*{self.base}**{e_str} % {self.modulus})')\n",
    "                else:\n",
    "                    wire_symbols.append('t' + str(t))\n",
    "                t += 1\n",
    "            if isinstance(self.exponent, Sequence) and qubit in self.exponent:\n",
    "                wire_symbols.append('e' + str(e))\n",
    "                e += 1\n",
    "        return cirq.CircuitDiagramInfo(wire_symbols=tuple(wire_symbols))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "4009a307",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "To factor n = 15 which has L = 4 bits, we need 3L + 3 = 15 qubits.\n"
     ]
    }
   ],
   "source": [
    "\"\"\"Create the target and exponent registers for phase estimation,\n",
    "and see the number of qubits needed for Shor's algorithm.\n",
    "\"\"\"\n",
    "n = 15\n",
    "L = n.bit_length()\n",
    "\n",
    "# The target register has L qubits.\n",
    "target = cirq.LineQubit.range(L)\n",
    "\n",
    "# The exponent register has 2L + 3 qubits.\n",
    "exponent = cirq.LineQubit.range(L, 3 * L + 3)\n",
    "\n",
    "# Display the total number of qubits to factor this n.\n",
    "print(f\"To factor n = {n} which has L = {L} bits, we need 3L + 3 = {3 * L + 3} qubits.\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "abd46fea",
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"See (part of) the unitary for a modular exponential operation.\"\"\"\n",
    "# Pick some element of the multiplicative group modulo n.\n",
    "x = 5\n",
    "\n",
    "# Display (part of) the unitary. Uncomment if n is small enough.\n",
    "# cirq.unitary(ModularExp(target, exponent, x, n))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "10fb28db",
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"Function to make the quantum circuit for order finding.\"\"\"\n",
    "def make_order_finding_circuit(x: int, n: int) -> cirq.Circuit:\n",
    "    \"\"\"Returns quantum circuit which computes the order of x modulo n.\n",
    "\n",
    "    The circuit uses Quantum Phase Estimation to compute an eigenvalue of\n",
    "    the unitary\n",
    "\n",
    "        U|y⟩ = |y * x mod n⟩      0 <= y < n\n",
    "        U|y⟩ = |y⟩                n <= y\n",
    "\n",
    "    Args:\n",
    "        x: positive integer whose order modulo n is to be found\n",
    "        n: modulus relative to which the order of x is to be found\n",
    "\n",
    "    Returns:\n",
    "        Quantum circuit for finding the order of x modulo n\n",
    "    \"\"\"\n",
    "    L = n.bit_length()\n",
    "    target = cirq.LineQubit.range(L)\n",
    "    exponent = cirq.LineQubit.range(L, 3 * L + 3)\n",
    "    return cirq.Circuit(\n",
    "        cirq.X(target[L - 1]),\n",
    "        cirq.H.on_each(*exponent),\n",
    "        ModularExp(target, exponent, x, n),\n",
    "        cirq.qft(*exponent, inverse=True),\n",
    "        cirq.measure(*exponent, key='exponent'),\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "c9463756",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0: ────────ModularExp(t*7**e % 15)────────────────────────────\n",
      "           │\n",
      "1: ────────t1─────────────────────────────────────────────────\n",
      "           │\n",
      "2: ────────t2─────────────────────────────────────────────────\n",
      "           │\n",
      "3: ────X───t3─────────────────────────────────────────────────\n",
      "           │\n",
      "4: ────H───e0────────────────────────qft^-1───M('exponent')───\n",
      "           │                         │        │\n",
      "5: ────H───e1────────────────────────#2───────M───────────────\n",
      "           │                         │        │\n",
      "6: ────H───e2────────────────────────#3───────M───────────────\n",
      "           │                         │        │\n",
      "7: ────H───e3────────────────────────#4───────M───────────────\n",
      "           │                         │        │\n",
      "8: ────H───e4────────────────────────#5───────M───────────────\n",
      "           │                         │        │\n",
      "9: ────H───e5────────────────────────#6───────M───────────────\n",
      "           │                         │        │\n",
      "10: ───H───e6────────────────────────#7───────M───────────────\n",
      "           │                         │        │\n",
      "11: ───H───e7────────────────────────#8───────M───────────────\n",
      "           │                         │        │\n",
      "12: ───H───e8────────────────────────#9───────M───────────────\n",
      "           │                         │        │\n",
      "13: ───H───e9────────────────────────#10──────M───────────────\n",
      "           │                         │        │\n",
      "14: ───H───e10───────────────────────#11──────M───────────────\n"
     ]
    }
   ],
   "source": [
    "\"\"\"Example of the quantum circuit for period finding.\"\"\"\n",
    "n = 15\n",
    "x = 7\n",
    "circuit = make_order_finding_circuit(x, n)\n",
    "print(circuit)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "788fa76d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# example of measurement outcomes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "9bbd5511",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Raw measurements:\n",
      "exponent=01110111, 00000000, 00000000, 00000000, 00000000, 00000000, 00000000, 00000000, 00000000\n",
      "\n",
      "Integer in exponent register:\n",
      "   exponent\n",
      "0         0\n",
      "1       256\n",
      "2       256\n",
      "3       256\n",
      "4         0\n",
      "5       256\n",
      "6       256\n",
      "7       256\n"
     ]
    }
   ],
   "source": [
    "\"\"\"Measuring Shor's period finding circuit.\"\"\"\n",
    "circuit = make_order_finding_circuit(x=5, n=6)\n",
    "res = cirq.sample(circuit, repetitions=8)\n",
    "\n",
    "print(\"Raw measurements:\")\n",
    "print(res)\n",
    "\n",
    "print(\"\\nInteger in exponent register:\")\n",
    "print(res.data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "45aa4b66",
   "metadata": {},
   "outputs": [],
   "source": [
    "# classical post processing\n",
    "def process_measurement(result: cirq.Result, x: int, n: int) -> Optional[int]:\n",
    "    \"\"\"Interprets the output of the order finding circuit.\n",
    "\n",
    "    Specifically, it determines s/r such that exp(2πis/r) is an eigenvalue\n",
    "    of the unitary\n",
    "\n",
    "        U|y⟩ = |xy mod n⟩  0 <= y < n\n",
    "        U|y⟩ = |y⟩         n <= y\n",
    "\n",
    "    then computes r (by continued fractions) if possible, and returns it.\n",
    "\n",
    "    Args:\n",
    "        result: result obtained by sampling the output of the\n",
    "            circuit built by make_order_finding_circuit\n",
    "\n",
    "    Returns:\n",
    "        r, the order of x modulo n or None.\n",
    "    \"\"\"\n",
    "    # Read the output integer of the exponent register.\n",
    "    exponent_as_integer = result.data[\"exponent\"][0]\n",
    "    exponent_num_bits = result.measurements[\"exponent\"].shape[1]\n",
    "    eigenphase = float(exponent_as_integer / 2**exponent_num_bits)\n",
    "\n",
    "    # Run the continued fractions algorithm to determine f = s / r.\n",
    "    f = fractions.Fraction.from_float(eigenphase).limit_denominator(n)\n",
    "\n",
    "    # If the numerator is zero, the order finder failed.\n",
    "    if f.numerator == 0:\n",
    "        return None\n",
    "\n",
    "    # Else, return the denominator if it is valid.\n",
    "    r = f.denominator\n",
    "    if x**r % n != 1:\n",
    "        return None\n",
    "    return r"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "bb3a31b6",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Finding the order of x = 5 modulo n = 6\n",
      "\n",
      "Raw measurements:\n",
      "exponent=1, 0, 0, 0, 0, 0, 0, 0, 0\n",
      "\n",
      "Integer in exponent register:\n",
      "   exponent\n",
      "0       256\n",
      "\n",
      "Order r = 2\n",
      "x^r mod n = 5^2 mod 6 = 1\n"
     ]
    }
   ],
   "source": [
    "\"\"\"Example of the classical post-processing.\"\"\"\n",
    "# Set n and x here\n",
    "n = 6\n",
    "x = 5\n",
    "\n",
    "print(f\"Finding the order of x = {x} modulo n = {n}\\n\")\n",
    "measurement = cirq.sample(circuit, repetitions=1)\n",
    "print(\"Raw measurements:\")\n",
    "print(measurement)\n",
    "\n",
    "print(\"\\nInteger in exponent register:\")\n",
    "print(measurement.data)\n",
    "\n",
    "r = process_measurement(measurement, x, n)\n",
    "print(\"\\nOrder r =\", r)\n",
    "if r is not None:\n",
    "    print(f\"x^r mod n = {x}^{r} mod {n} = {x**r % n}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "a907a7fb",
   "metadata": {},
   "outputs": [],
   "source": [
    "# quantum order finder\n",
    "def quantum_order_finder(x: int, n: int) -> Optional[int]:\n",
    "    \"\"\"Computes smallest positive r such that x**r mod n == 1.\n",
    "\n",
    "    Args:\n",
    "        x: integer whose order is to be computed, must be greater than one\n",
    "           and belong to the multiplicative group of integers modulo n (which\n",
    "           consists of positive integers relatively prime to n),\n",
    "        n: modulus of the multiplicative group.\n",
    "    \"\"\"\n",
    "    # Check that the integer x is a valid element of the multiplicative group\n",
    "    # modulo n.\n",
    "    if x < 2 or n <= x or math.gcd(x, n) > 1:\n",
    "        raise ValueError(f'Invalid x={x} for modulus n={n}.')\n",
    "\n",
    "    # Create the order finding circuit.\n",
    "    circuit = make_order_finding_circuit(x, n)\n",
    "\n",
    "    # Sample from the order finding circuit.\n",
    "    measurement = cirq.sample(circuit)\n",
    "\n",
    "    # Return the processed measurement result.\n",
    "    return process_measurement(measurement, x, n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "e4683269",
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"Functions for factoring from start to finish.\"\"\"\n",
    "def find_factor_of_prime_power(n: int) -> Optional[int]:\n",
    "    \"\"\"Returns non-trivial factor of n if n is a prime power, else None.\"\"\"\n",
    "    for k in range(2, math.floor(math.log2(n)) + 1):\n",
    "        c = math.pow(n, 1 / k)\n",
    "        c1 = math.floor(c)\n",
    "        if c1**k == n:\n",
    "            return c1\n",
    "        c2 = math.ceil(c)\n",
    "        if c2**k == n:\n",
    "            return c2\n",
    "    return None\n",
    "\n",
    "\n",
    "def find_factor(\n",
    "    n: int,\n",
    "    order_finder: Callable[[int, int], Optional[int]] = quantum_order_finder,\n",
    "    max_attempts: int = 30\n",
    ") -> Optional[int]:\n",
    "    \"\"\"Returns a non-trivial factor of composite integer n.\n",
    "\n",
    "    Args:\n",
    "        n: Integer to factor.\n",
    "        order_finder: Function for finding the order of elements of the\n",
    "            multiplicative group of integers modulo n.\n",
    "        max_attempts: number of random x's to try, also an upper limit\n",
    "            on the number of order_finder invocations.\n",
    "\n",
    "    Returns:\n",
    "        Non-trivial factor of n or None if no such factor was found.\n",
    "        Factor k of n is trivial if it is 1 or n.\n",
    "    \"\"\"\n",
    "    # If the number is prime, there are no non-trivial factors.\n",
    "    if sympy.isprime(n):\n",
    "        print(\"n is prime!\")\n",
    "        return None\n",
    "\n",
    "    # If the number is even, two is a non-trivial factor.\n",
    "    if n % 2 == 0:\n",
    "        return 2\n",
    "\n",
    "    # If n is a prime power, we can find a non-trivial factor efficiently.\n",
    "    c = find_factor_of_prime_power(n)\n",
    "    if c is not None:\n",
    "        return c\n",
    "\n",
    "    for _ in range(max_attempts):\n",
    "        # Choose a random number between 2 and n - 1.\n",
    "        x = random.randint(2, n - 1)\n",
    "\n",
    "        # Most likely x and n will be relatively prime.\n",
    "        c = math.gcd(x, n)\n",
    "\n",
    "        # If x and n are not relatively prime, we got lucky and found\n",
    "        # a non-trivial factor.\n",
    "        if 1 < c < n:\n",
    "            return c\n",
    "\n",
    "        # Compute the order r of x modulo n using the order finder.\n",
    "        r = order_finder(x, n)\n",
    "\n",
    "        # If the order finder failed, try again.\n",
    "        if r is None:\n",
    "            continue\n",
    "\n",
    "        # If the order r is even, try again.\n",
    "        if r % 2 != 0:\n",
    "            continue\n",
    "\n",
    "        # Compute the non-trivial factor.\n",
    "        y = x**(r // 2) % n\n",
    "        assert 1 < y < n\n",
    "        c = math.gcd(y - 1, n)\n",
    "        if 1 < c < n:\n",
    "            return c\n",
    "\n",
    "    print(f\"Failed to find a non-trivial factor in {max_attempts} attempts.\")\n",
    "    return None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "97a1b28e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Factoring n = pq = 184573\n",
      "p = 487\n",
      "q = 379\n"
     ]
    }
   ],
   "source": [
    "\"\"\"Example of factoring via Shor's algorithm (order finding).\"\"\"\n",
    "# Number to factor\n",
    "n = 184573\n",
    "\n",
    "# Attempt to find a factor\n",
    "p = find_factor(n, order_finder=classical_order_finder)\n",
    "q = n // p\n",
    "\n",
    "print(\"Factoring n = pq =\", n)\n",
    "print(\"p =\", p)\n",
    "print(\"q =\", q)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "f69291ec",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\"\"\"Check the answer is correct.\"\"\"\n",
    "p * q == n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "81439a17",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}