Classical Steins algorithm

Author: fedimser

arithmetic/GCD using Stein's algorithm.ipynb

@@ -0,0 +1,113 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "e768812e",
+   "metadata": {},
+   "source": [
+    "Stein's algorithm for Greatest Common Divisor, adapted for implementation on a quantum computer.\n",
+    "\n",
+    "Source: [this paper](https://arxiv.org/abs/1304.7516)."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "id": "b7ab54fe",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Exhaustive tests for N=1 OK\n",
+      "Exhaustive tests for N=2 OK\n",
+      "Exhaustive tests for N=3 OK\n",
+      "Exhaustive tests for N=4 OK\n",
+      "Exhaustive tests for N=5 OK\n",
+      "Exhaustive tests for N=6 OK\n",
+      "Exhaustive tests for N=7 OK\n",
+      "Exhaustive tests for N=8 OK\n"
+     ]
+    },
+    {
+     "ename": "NameError",
+     "evalue": "name 'gcd_stein2' is not defined",
+     "output_type": "error",
+     "traceback": [
+      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
+      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
+      "Cell \u001b[0;32mIn[1], line 40\u001b[0m\n\u001b[1;32m     38\u001b[0m         b \u001b[38;5;241m=\u001b[39m random\u001b[38;5;241m.\u001b[39mrandint(\u001b[38;5;241m1\u001b[39m, \u001b[38;5;241m2\u001b[39m\u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39mN\u001b[38;5;241m-\u001b[39m\u001b[38;5;241m1\u001b[39m)\n\u001b[1;32m     39\u001b[0m         \u001b[38;5;28;01massert\u001b[39;00m gcd_stein(N,a,b) \u001b[38;5;241m==\u001b[39m math\u001b[38;5;241m.\u001b[39mgcd(a,b)\n\u001b[0;32m---> 40\u001b[0m         \u001b[38;5;28;01massert\u001b[39;00m \u001b[43mgcd_stein2\u001b[49m(a,b) \u001b[38;5;241m==\u001b[39m math\u001b[38;5;241m.\u001b[39mgcd(a,b)\n\u001b[1;32m     41\u001b[0m     \u001b[38;5;28mprint\u001b[39m(\u001b[38;5;124mf\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mRandom tests for N=\u001b[39m\u001b[38;5;132;01m{\u001b[39;00mN\u001b[38;5;132;01m}\u001b[39;00m\u001b[38;5;124m OK\u001b[39m\u001b[38;5;124m\"\u001b[39m)\n\u001b[1;32m     43\u001b[0m \u001b[38;5;28mprint\u001b[39m(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mOK\u001b[39m\u001b[38;5;124m\"\u001b[39m)\n",
+      "\u001b[0;31mNameError\u001b[0m: name 'gcd_stein2' is not defined"
+     ]
+    }
+   ],
+   "source": [
+    "def gcd_stein(N,a,b):\n",
+    "    assert 0 <= a < 2**N\n",
+    "    assert 0 < b < 2**N\n",
+    "    \n",
+    "    r = 1\n",
+    "    for i in range(2*N-3):\n",
+    "        a_is_not_zero = (a!=0)\n",
+    "        a_is_even = (a%2==0)\n",
+    "        b_is_even = (b%2==0)    \n",
+    "        if a_is_not_zero and b_is_even:\n",
+    "            b//=2\n",
+    "        if a_is_not_zero and a_is_even:\n",
+    "            a//=2\n",
+    "        a_less_b = (a<b)\n",
+    "        if a_is_not_zero and a_less_b:\n",
+    "            a,b=b,a\n",
+    "        assert(a==0 or a>=b)\n",
+    "        if (not a_is_even) and (not b_is_even):\n",
+    "            a -= b\n",
+    "            a //= 2\n",
+    "        if a_is_not_zero and a_is_even and b_is_even:\n",
+    "            r *= 2\n",
+    "    assert r<2**N\n",
+    "    return b * r\n",
+    "    \n",
+    "import math\n",
+    "import random\n",
+    "\n",
+    "for N in [1,2,3,4,5,6,7,8]:\n",
+    "    for a in range(0, 2**N):\n",
+    "        for b in range(1, 2**N):\n",
+    "            assert gcd_stein(N,a,b) == math.gcd(a,b)\n",
+    "    print(f\"Exhaustive tests for N={N} OK\")\n",
+    "        \n",
+    "for N in [9,10,16,32,64,128]:\n",
+    "    for _ in range(1000):\n",
+    "        a = random.randint(0, 2**N-1)\n",
+    "        b = random.randint(1, 2**N-1)\n",
+    "        assert gcd_stein(N,a,b) == math.gcd(a,b)\n",
+    "        assert gcd_stein2(a,b) == math.gcd(a,b)\n",
+    "    print(f\"Random tests for N={N} OK\")\n",
+    "    \n",
+    "print(\"OK\")"
+   ]
+  }
+ ],
+ "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.12.3"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}