{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Pohlig-Hellman algorithm\n",
    "\n",
    "We will perform the Pohlig-Hellman algorithm to compute $\\log_{a} b$ in $\\mathbf{Z}_{p}^*$, where $a = 11$, $b = 2020$, and $p = 15121$ is a prime number."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from algorithms.factorization import totalFactorization\n",
    "from algorithms.discreteLogarithm import babyStepGiantStep\n",
    "from algorithms.modular import crt\n",
    "\n",
    "a = 11\n",
    "b = 2020\n",
    "p = 15121"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We notice that $p-1$ can be efficiently factorized. This allows the Pohlig-Hellman algorithm to break the problem into multiple smaller problems which can be solved, say, using the giant-step/baby-step algorithm."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{3: 3, 2: 4, 5: 1, 7: 1}"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "n = p - 1\n",
    "f = totalFactorization(n)\n",
    "f"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now perform the main loop of the algorithm. The idea is to compute discrete logarithms with base $(p-1)/q$ for each prime power factor $q^e$ of $p-1$ which give a $q$-ary representation of the sought result modulo $q^e$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "prime power factor: 3^3\n",
      "a^(15120/3) = 7292\n",
      "c_0 = c_-1 * a^-k_-1 = 1\n",
      "b_0 = (b_-1*c_0)^(15120/3^1) = 7828\n",
      "z = log_7292(7828) mod 15121, order = 3\n",
      "k_0 = z*q^0 = 2\n",
      "c_1 = c_0 * a^-k_0 = 11372\n",
      "b_1 = (b_0*c_1)^(15120/3^2) = 1\n",
      "z = log_7292(1) mod 15121, order = 3\n",
      "k_1 = z*q^1 = 0\n",
      "c_2 = c_1 * a^-k_1 = 11372\n",
      "b_2 = (b_1*c_2)^(15120/3^3) = 1\n",
      "z = log_7292(1) mod 15121, order = 3\n",
      "k_2 = z*q^2 = 0\n",
      "prime power factor: 2^4\n",
      "a^(15120/2) = 15120\n",
      "c_0 = c_-1 * a^-k_-1 = 1\n",
      "b_0 = (b_-1*c_0)^(15120/2^1) = 15120\n",
      "z = log_15120(15120) mod 15121, order = 2\n",
      "k_0 = z*q^0 = 1\n",
      "c_1 = c_0 * a^-k_0 = 4124\n",
      "b_1 = (b_0*c_1)^(15120/2^2) = 1\n",
      "z = log_15120(1) mod 15121, order = 2\n",
      "k_1 = z*q^1 = 0\n",
      "c_2 = c_1 * a^-k_1 = 4124\n",
      "b_2 = (b_1*c_2)^(15120/2^3) = 15120\n",
      "z = log_15120(15120) mod 15121, order = 2\n",
      "k_2 = z*q^2 = 4\n",
      "c_3 = c_2 * a^-k_2 = 8938\n",
      "b_3 = (b_2*c_3)^(15120/2^4) = 15120\n",
      "z = log_15120(15120) mod 15121, order = 2\n",
      "k_3 = z*q^3 = 8\n",
      "prime power factor: 5^1\n",
      "a^(15120/5) = 11389\n",
      "c_0 = c_-1 * a^-k_-1 = 1\n",
      "b_0 = (b_-1*c_0)^(15120/5^1) = 1383\n",
      "z = log_11389(1383) mod 15121, order = 5\n",
      "k_0 = z*q^0 = 2\n",
      "prime power factor: 7^1\n",
      "a^(15120/7) = 10187\n",
      "c_0 = c_-1 * a^-k_-1 = 1\n",
      "b_0 = (b_-1*c_0)^(15120/7^1) = 7205\n",
      "z = log_10187(7205) mod 15121, order = 7\n",
      "k_0 = z*q^0 = 6\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "{27: 2, 16: 13, 5: 2, 7: 6}"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ai = a^-1 % p\n",
    "v = {}\n",
    "for q, e in f.items():\n",
    "    print(\"prime power factor: %d^%d\" % (q, e))\n",
    "    nq = n // q\n",
    "    ap = pow(a, nq, p)\n",
    "    print(\"a^(%d/%d) = %d\" % (n, q, ap))\n",
    "    bp, qp = b, 1\n",
    "    c, k, x = 1, 0, 0\n",
    "    for i in range(e):\n",
    "        c = c * pow(ai, k, p) % p\n",
    "        bp = pow(b*c, nq, p)\n",
    "        print(\"c_%d = c_%d * a^-k_%d = %d\" % (i, i-1, i-1, c))\n",
    "        print(\"b_%d = (b_%d*c_%d)^(%d/%d^%d) = %d\" %\n",
    "                                            (i, i-1, i, n, q, i+1, bp))\n",
    "        print(\"z = log_%d(%d) mod %d, order = %d\" % (ap, bp, p, q))\n",
    "        k = (babyStepGiantStep(ap, bp, p, order = q) % q) * qp\n",
    "        x += k\n",
    "        print(\"k_%d = z*q^%d = %d\" % (i, i, k))\n",
    "        qp *= q\n",
    "        nq //= q\n",
    "    v[qp] = x\n",
    "v"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This gives a system of congruences with pairwise coprime moduli, which can then be solved using the Chinese remainder theorem."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "12557"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = crt(v)\n",
    "x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us check that we get the correct answer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2020"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pow(a, x, p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Running time comparison\n",
    "\n",
    "We will now use the function `pohligHellman` to compute some more discrete logarithms and measure the evaluation times."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "from algorithms.discreteLogarithm import pohligHellman"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 30.6 ms, sys: 0 ns, total: 30.6 ms\n",
      "Wall time: 28.8 ms\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "6935101882"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%time pohligHellman(47, 191, 100000000003)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 29 ms, sys: 0 ns, total: 29 ms\n",
      "Wall time: 28.5 ms\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "6780065714244411514"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%time pohligHellman(47, 191, 10000000000000000051)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 33.9 ms, sys: 0 ns, total: 33.9 ms\n",
      "Wall time: 32.6 ms\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "354067223054486232245"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%time pohligHellman(47, 191, 1000000000000000000117)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%time pohligHellman(47, 191, 10000000000000000000000343, trace = 2)"
   ]
  }
 ],
 "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
}