{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Integer factorization\n",
    "\n",
    "We will see three algorithms for integer factorization: the Pollard $p-1$ algorithm, the Pollard $\\rho$ algorithm, and Dixon's random squares algorithm."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from algorithms.factorization import pollardP1, pollardRho, totalFactorization, factorizeByBase"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Pollard $p-1$ algorithm\n",
    "\n",
    "Let us try to determine the bound $B$ needed to factor $262063$ and $9420457$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "521"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pollardP1(262063, 13)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "521.is_prime()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2^3 * 5 * 13"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "520.factor()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "503"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "262063 / 521"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "503.is_prime()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2 * 251"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "502.factor()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2351"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pollardP1(9420457, 47)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "2351.is_prime()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2 * 5^2 * 47"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "2350.factor()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Pollard $\\rho$ algorithm\n",
    "\n",
    "Let us try to factor $221$ using the function $f(x) = (x^2 + 1) \\bmod{221}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x = 194\n",
      "f(x) = 67\n",
      "gcd(x-f(x), n) = 1\n",
      "f^1(x) = 67\n",
      "f^3(x) = 39\n",
      "gcd(f^1(x)-f^3(x), n) = 1\n",
      "f^2(x) = 70\n",
      "f^5(x) = 184\n",
      "gcd(f^2(x)-f^5(x), n) = 1\n",
      "f^3(x) = 39\n",
      "f^7(x) = 169\n",
      "gcd(f^3(x)-f^7(x), n) = 13\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "13"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "n = 221\n",
    "f = lambda x: (x^2 + 1) % n\n",
    "pollardRho(n, f, trace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "17"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "221 / 13"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Dixon's random squares algorithm\n",
    "\n",
    "We will factorize $n = 15770708441$ using the random integers $14056852462$, $9004436975$, $5286195500$, $11335959904$ and $7052612564$. Let us factorize the squares of the random integers modulo $n$. We notice that they can be expressed as products of powers of the smallest seven primes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[[1, 3, 0, 0, 0, 3, 2],\n",
       " [0, 2, 3, 8, 0, 0, 0],\n",
       " [2, 3, 0, 3, 0, 2, 0],\n",
       " [4, 0, 4, 0, 1, 0, 4],\n",
       " [0, 8, 1, 0, 1, 2, 0]]"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "n = 15770708441\n",
    "r = [14056852462, 9004436975, 5286195500, 11335959904, 7052612564]\n",
    "b = [2, 3, 5, 7, 11, 13, 17]\n",
    "fac = [factorizeByBase(x^2 % n, b, n) for x in r]\n",
    "fac"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us gather the exponents into a matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1 3 0 0 0 3 2]\n",
       "[0 2 3 8 0 0 0]\n",
       "[2 3 0 3 0 2 0]\n",
       "[4 0 4 0 1 0 4]\n",
       "[0 8 1 0 1 2 0]"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "M = Matrix(fac)\n",
    "M"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(4, 10, 8, 8, 2, 2, 4)"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sum(M[[1, 3, 4], :])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We notice that summing the second and the last two rows gives a vector consisting entirely of even numbers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[(0), (0), (0), (0), (0), (0), (0)]"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rows = (1, 3, 4)\n",
    "[sum(M[rows, i]) % 2 for i in range(len(b))]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use this to find a factor of $n$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(135979, 115979)"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pr = prod(r[i] for i in rows) % n\n",
    "pb = prod(p^e for p, (e, ) in zip(b, (sum(M[rows, i])/2 for i in range(len(b))))) % n\n",
    "g = gcd(abs(pr - pb), n)\n",
    "(g, n/g)"
   ]
  }
 ],
 "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}