{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Index calculus algorithm\n",
    "\n",
    "We will perform the Index calculus algorithm to compute $\\log_{a} b_i$ ($i = 1, 2$) in $\\mathbf{Z}_{p}^*$, where $a = 5$, $b_1 = 4389733$, $b_2 = 1234567$, and $p = 9330887$ is a prime number."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from algorithms.euclidean import gcd\n",
    "from algorithms.factorization import factorizeByBase\n",
    "\n",
    "p = 9330887\n",
    "a = 5\n",
    "b1 = 4389733\n",
    "b2 = 1234567"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will use a base consisting of $-1$ and all primes less than $50$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[-1, 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47]"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "base = [-1] + prime_range(50)\n",
    "base"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us now try to find the logarithms table. We construct a matrix by trying random powers of $a$ and factoring them with the numbers in our base. We stop when the matrix has full rank."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1 0 2 1 0 0 1 0 0 0 0 0 0 0 2 0]\n",
      "[1 0 1 0 0 0 0 0 0 0 0 1 3 0 0 0]\n",
      "[0 3 8 0 1 0 0 0 0 1 0 0 0 0 0 0]\n",
      "[0 1 2 1 2 0 0 0 0 0 1 0 1 0 0 0]\n",
      "[1 0 0 0 1 3 1 0 0 0 0 0 0 1 0 0]\n",
      "[1 4 0 0 1 0 0 1 0 0 1 0 0 0 0 0]\n",
      "[1 2 0 3 3 0 0 0 0 1 0 0 0 0 0 0]\n",
      "[1 2 2 0 0 1 1 1 0 1 0 0 0 0 0 0]\n",
      "[0 8 1 0 0 2 0 0 0 0 0 0 0 0 1 0]\n",
      "[0 3 1 0 0 1 0 0 1 1 0 0 0 0 0 0]\n",
      "[0 5 0 1 0 0 0 1 0 0 1 0 1 0 0 0]\n",
      "[0 0 0 1 1 0 0 0 1 0 0 0 0 2 0 0]\n",
      "[0 0 0 0 2 0 0 0 1 0 1 0 0 0 0 1]\n",
      "[0 1 3 0 4 0 1 0 0 0 0 0 0 0 0 0]\n",
      "[1 4 1 5 0 0 0 0 0 1 0 0 0 0 0 0]\n",
      "[1 1 1 2 0 2 1 0 1 0 0 0 0 0 0 0]\n",
      "[1342485, 669118, 7237936, 7922551, 9149713, 4288717, 563606, 4482702, 7606125, 6167680, 7916700, 6825736, 4712055, 7120183, 7905940, 4584527]\n"
     ]
    }
   ],
   "source": [
    "set_random_seed(0)\n",
    "\n",
    "v = []\n",
    "s = set()\n",
    "r, l = 0, len(base)\n",
    "M = Matrix(nrows=0, ncols=l)\n",
    "while r < l:\n",
    "    i = randint(1, p-1)\n",
    "    if i in s:\n",
    "        continue\n",
    "    s.add(i)\n",
    "    x = pow(a, i, p)\n",
    "    f = factorizeByBase(Integer(x), base, p)\n",
    "    if not f:\n",
    "        continue\n",
    "    MM = Matrix(list(M) + [f])\n",
    "    rr = MM.rank()\n",
    "    if rr > r:\n",
    "        M = MM\n",
    "        r = rr\n",
    "        v.append(i)\n",
    "print(M)\n",
    "print(v)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now solve the system of equations we have obtained. The solution represents our table of logarithms and can be used to find logarithms of multiple numbers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[4665443]\n",
       "[6670912]\n",
       "[2030334]\n",
       "[      1]\n",
       "[5786904]\n",
       "[1078197]\n",
       "[8534197]\n",
       "[7606749]\n",
       "[8519903]\n",
       "[2519168]\n",
       "[6200403]\n",
       "[9068634]\n",
       "[7409417]\n",
       "[5590350]\n",
       "[6037417]\n",
       "[6410599]"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "t = M^-1 * Matrix(zip(v)) % (p-1)\n",
    "t"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us verify that the table is correct."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[9330886, 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47]"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "[pow(a, i, p) for i, in t]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now find a power of $b_i$ that can be factorized with our base. We will use the following function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def findPower(b, p, base):\n",
    "    while True:\n",
    "        i = randint(1, p-1)\n",
    "        x = pow(b, i, p)\n",
    "        f = factorizeByBase(Integer(x), base, p)\n",
    "        if not f:\n",
    "            continue\n",
    "        return (i, f)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us find such a power of $b_1$ and its factorization."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(3157805, [0, 0, 3, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 0])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "i1, f1 = findPower(b1, p, base)\n",
    "i1, f1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now use the logarithm table to obtain a congruence in $\\log_a b_1$. However, there may be multiple solutions, and we need to check which one is our answer. The number of solutions is given by the greatest common divisor of $p-1$ and the obtained exponent."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "g1 = gcd(i1, p-1)\n",
    "g1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can thus obtain a solution modulo $p-1$ divided by this GCD."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5753305"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m1 = ((p-1)/g1)\n",
    "(x1,), = Matrix([f1])*t / i1 % m1\n",
    "x1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us now put this into a function which will also verify the potential solutions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "5^5753305 mod 9330887 = 4389733\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "5753305"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def checkSolutions(a, b, p, t, i, f):\n",
    "    g = gcd(i, p-1)\n",
    "    m = (p-1)/g\n",
    "    (x,), = Matrix([f])*t / i % m\n",
    "    for j in range(g):\n",
    "        y = pow(a, x, p)\n",
    "        print(\"%d^%d mod %d = %d\" % (a, x, p, y))\n",
    "        if y == b:\n",
    "            return x\n",
    "        x += m\n",
    "\n",
    "checkSolutions(a, b1, p, t, i1, f1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have thus computed $\\log_a b_1$. Let us now compute $\\log_a b_2$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(5353502, [0, 0, 2, 1, 0, 0, 2, 1, 0, 0, 0, 0, 0, 1, 0, 0])"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "i2, f2 = findPower(b2, p, base)\n",
    "i2, f2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "5^2608331 mod 9330887 = 8096320\n",
      "5^7273774 mod 9330887 = 1234567\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "7273774"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "checkSolutions(a, b2, p, t, i2, f2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Running time comparison\n",
    "\n",
    "We will now use the function `logarithmTable` to compute tables of logarithms, which will then be used to compute some more discrete logarithms with the function `indexCalculus`. We will measure the evaluation times."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "from algorithms.discreteLogarithm import logarithmTable, indexCalculus\n",
    "\n",
    "aa = 47\n",
    "bb = 191"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "found factorization 47^8920294660 = -1^1 * 2^2 * 3^2 * 7^1 * 11^1 * 13^1 * 19^1 * 41^1 * 43^1 * 47^1 (mod 100000000003)\n",
      "found factorization 47^86731650251 = -1^1 * 2^2 * 3^2 * 7^1 * 11^3 * 13^2 * 19^1 * 31^1 (mod 100000000003)\n",
      "found factorization 47^39752784926 = 2^7 * 3^1 * 5^1 * 11^3 * 19^2 * 23^1 (mod 100000000003)\n",
      "found factorization 47^69419634802 = 2^1 * 3^4 * 5^7 * 19^1 * 37^1 (mod 100000000003)\n",
      "found factorization 47^75320357659 = -1^1 * 3^1 * 5^3 * 7^1 * 11^2 * 13^2 * 31^1 * 47^1 (mod 100000000003)\n",
      "found factorization 47^40532049132 = -1^1 * 2^4 * 5^1 * 11^1 * 19^1 * 31^1 * 37^1 * 41^1 * 47^1 (mod 100000000003)\n",
      "found factorization 47^32937591819 = 3^3 * 5^1 * 11^4 * 13^1 * 31^2 (mod 100000000003)\n",
      "found factorization 47^6602367151 = -1^1 * 2^1 * 7^2 * 17^1 * 29^4 * 43^1 (mod 100000000003)\n",
      "found factorization 47^70383598779 = -1^1 * 2^1 * 11^1 * 17^1 * 31^2 * 37^1 * 41^1 * 43^1 (mod 100000000003)\n",
      "found factorization 47^82842048017 = 5^1 * 13^3 * 17^1 * 29^1 (mod 100000000003)\n",
      "found factorization 47^93452946324 = 2^3 * 3^1 * 5^1 * 11^1 * 19^1 * 23^1 * 29^1 * 31^1 * 41^1 (mod 100000000003)\n",
      "found factorization 47^85301660587 = 3^1 * 11^2 * 19^4 * 23^1 * 47^1 (mod 100000000003)\n",
      "found factorization 47^95880313549 = -1^1 * 2^1 * 3^4 * 13^4 * 23^3 (mod 100000000003)\n",
      "found factorization 47^32201305943 = 2^9 * 11^3 * 17^1 * 31^1 * 41^1 (mod 100000000003)\n",
      "found factorization 47^26135740991 = -1^1 * 2^2 * 3^8 * 5^1 * 7^1 * 11^2 * 13^1 * 37^1 (mod 100000000003)\n",
      "found factorization 47^47177975199 = -1^1 * 2^5 * 7^2 * 11^1 * 13^1 * 23^1 * 31^1 * 47^1 (mod 100000000003)\n",
      "CPU times: user 2.47 s, sys: 10.5 ms, total: 2.48 s\n",
      "Wall time: 2.51 s\n",
      "[50000000001, 88948660783, 74403602789, 76906212201, 2995382539, 11520765452, 33861633403, 95303711677, 78308239391, 74623569337, 9047223934, 26437938024, 66204838085, 75359412859, 30170333882, 1]\n",
      "found factorization 191^92538805852 = -1^1 * 2^2 * 3^3 * 23^2 * 29^1 * 47^2 (mod 100000000003)\n",
      "checking 2 solutions for 47^x = 191, x = 6935101882 (mod 50000000001)\n",
      "47^6935101882 mod 100000000003 = 191\n",
      "CPU times: user 86.4 ms, sys: 0 ns, total: 86.4 ms\n",
      "Wall time: 89.4 ms\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "6935101882"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pp = 100000000003\n",
    "table = %time logarithmTable(aa, pp, base, trace = True)\n",
    "print(table)\n",
    "%time indexCalculus(aa, bb, pp, base, table, trace = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pp = 10000000000000000051\n",
    "base = [-1] + prime_range(5000)\n",
    "table = %time logarithmTable(aa, pp, base, trace = True)\n",
    "%time indexCalculus(aa, bb, pp, base, table, trace = True)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "SageMath 9.2.rc2",
   "language": "sage",
   "name": "sagemath"
  },
  "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": 2
}