{
"cells": [
{
"cell_type": "markdown",
"id": "8c28da72",
"metadata": {},
"source": [
"# Factoring and Shor's Algorithm\n",
"\n",
"Quantum computers can factor large composite numbers in polynomial time using an algorithm initially introduced in [S1994] and improved over the years. Notably, [GE2019] optimized and compiled a variant of Shor's algorithm for the surface code. We'll first investigate how modular exponentiation can be used to back out factors of composite numbers then show the reference decomposition for modular exponentiation.\n",
"\n",
"\n",
"### References:\n",
"\n",
" - [S1994] Algorithms for quantum computation: discrete logarithms and factoring.\n",
"[10.1109/SFCS.1994.365700](https://dx.doi.org/10.1109/SFCS.1994.365700). Shor. 1994.\n",
"\n",
" - [GE2019] How to factor 2048 bit RSA integers in 8 hours using 20 million noisy qubits.\n",
"[arxiv:1905.09749](https://arxiv.org/abs/1905.09749). Gidney and Ekerå. 2019."
]
},
{
"cell_type": "markdown",
"id": "74e6cb0c",
"metadata": {},
"source": [
"## Walk through the factoring of 13 * 17\n",
"\n",
"We can factor a large number by finding the period of modular exponentiation. We set our input variables with demo values that we can manage classically.\n",
"\n",
" - `N` is the composite number to factor. Here, we set it to $13*17$\n",
" - `n` is its bitsize"
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "c182f5d8",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(221, 8)"
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"import numpy as np\n",
"N = 13*17\n",
"n = int(np.ceil(np.log2(N)))\n",
"N, n"
]
},
{
"cell_type": "markdown",
"id": "d54206cf",
"metadata": {},
"source": [
"Our modular exponentiation uses the composite number as the modulus and we will try a variety of different exponent values (to find the period). The base of the exponent doesn't matter that much and we choose a random base $g$."
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "47fb7255",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"8"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"g = 8\n",
"g"
]
},
{
"cell_type": "markdown",
"id": "7c542b60",
"metadata": {},
"source": [
"### Period Finding\n",
"\n",
"We need to find the period of the exponentiation -- namely how long it takes to cycle back to 1 under modular arithmetic. Below, we do this with a classical `for` loop for demonstration. On the quantum computer, we will execute modular exponentiation using a superposition of exponent values."
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "50a027d8",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" 0 1 *\n",
" 1 8\n",
" 2 64\n",
" 3 70\n",
" 4 118\n",
" 5 60\n",
" 6 38\n",
" 7 83\n",
" 8 1 *\n",
" 9 8\n",
" 10 64\n",
" 11 70\n",
" 12 118\n",
" 13 60\n",
" 14 38\n",
" 15 83\n",
" 16 1 *\n",
" 17 8\n",
" 18 64\n",
" 19 70\n"
]
}
],
"source": [
"for e in range(20):\n",
" f = (g ** e) % N\n",
" \n",
" star = ' *' if f == 1 else ''\n",
" print(f'{e:5d} {f:5d}{star}')"
]
},
{
"cell_type": "markdown",
"id": "7cb006f7",
"metadata": {},
"source": [
"The period is indeed a consistent value of 8."
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "b55445ff",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(8, 8)"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"16-8, 8-0"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "d7c8cf3f",
"metadata": {},
"outputs": [],
"source": [
"period = 8"
]
},
{
"cell_type": "markdown",
"id": "367ca473",
"metadata": {},
"source": [
"### Use the period to find factors\n",
"\n",
"We can use some numerical tricks find our two factors from this period. Consult the references if you'd like an explanation for why this works. Note that we make some assertions about $g$ and the period. If these are violated, you must re-run the period finding with a different choice for $g$."
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "9aaa0618",
"metadata": {},
"outputs": [],
"source": [
"assert period % 2 == 0\n",
"assert g**(period//2) != -1\n",
"\n",
"half_period = g**(period//2)\n",
"p1 = half_period + 1\n",
"m1 = half_period - 1\n",
"\n",
"assert (p1*m1) % N == 0"
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "caea53b2",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"gcd(119, 221), gcd(117, 221)\n"
]
},
{
"data": {
"text/plain": [
"(17, 13)"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"print(f'gcd{p1%N, N}, gcd{m1%N, N}')\n",
"import math\n",
"math.gcd(p1%N, N), math.gcd(m1%N, N)"
]
},
{
"cell_type": "markdown",
"id": "a07f7aa8",
"metadata": {},
"source": [
"We've recovered the two numbers and the factoring has been successfull."
]
},
{
"cell_type": "markdown",
"id": "79cb559b",
"metadata": {},
"source": [
"## Using the `ModExp` Bloq\n",
"\n",
"We can use the classical simulation capabilities of the library to show that our `ModExp` bloq performs the same arithmetic. We can test that its decomposition into sub-bloqs is correct by classically simulating the decomposition. When executed on a quantum computer, the `exponent` register will be in superposition leading to exponentially faster period finding."
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "4ce6ef64",
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
""
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"from qualtran.bloqs.cryptography.rsa.rsa_mod_exp import ModExp\n",
"from qualtran.drawing import show_bloq\n",
"\n",
"mod_exp = ModExp(base=g, mod=N, exp_bitsize=32, x_bitsize=32)\n",
"show_bloq(mod_exp)"
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "3d366054",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"{CModMulK(dtype=QUInt(bitsize=32), k=1, mod=221): 29,\n",
" Join(dtype=QUInt(bitsize=32)): 1,\n",
" CModMulK(dtype=QUInt(bitsize=32), k=118, mod=221): 1,\n",
" Split(dtype=QUInt(bitsize=32)): 1,\n",
" CModMulK(dtype=QUInt(bitsize=32), k=64, mod=221): 1,\n",
" CModMulK(dtype=QUInt(bitsize=32), k=8, mod=221): 1,\n",
" IntState(val=1, bitsize=32): 1}"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"mod_exp_decomp = mod_exp.decompose_bloq()\n",
"mod_exp_decomp.bloq_counts()"
]
},
{
"cell_type": "code",
"execution_count": 10,
"id": "680d8a30",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" 0 1 1 1 *\n",
" 1 8 8 8 \n",
" 2 64 64 64 \n",
" 3 70 70 70 \n",
" 4 118 118 118 \n",
" 5 60 60 60 \n",
" 6 38 38 38 \n",
" 7 83 83 83 \n",
" 8 1 1 1 *\n",
" 9 8 8 8 \n",
" 10 64 64 64 \n",
" 11 70 70 70 \n",
" 12 118 118 118 \n",
" 13 60 60 60 \n",
" 14 38 38 38 \n",
" 15 83 83 83 \n",
" 16 1 1 1 *\n",
" 17 8 8 8 \n",
" 18 64 64 64 \n",
" 19 70 70 70 \n"
]
}
],
"source": [
"from qualtran import QUInt\n",
"for e in range(20):\n",
" ref = (g ** e) % N\n",
" _, bloq_eval = mod_exp.call_classically(exponent=e)\n",
" _, decomp_eval = mod_exp_decomp.call_classically(exponent=e)\n",
" \n",
" star = ' *' if ref == 1 else ''\n",
" print(f'{e:5d} {ref:5d} {bloq_eval:5d} {decomp_eval:5d} {star}')"
]
}
],
"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.11.8"
}
},
"nbformat": 4,
"nbformat_minor": 5
}