{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Factorisation des entiers"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Marc Lorenzi\n",
    "\n",
    "11 juillet 2016"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import random\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams['figure.figsize'] = (10, 6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1 L'algorithme naïf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.1 Trouver un diviseur d'un entier"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On recherche le plus petit diviseur de $n$ (il sera forcément premier) en essayant tous les entiers à partir de 2. On peut tout de même remarquer qu'un entier composé $n$ a un diviseur inférieur à $\\sqrt n$. En effet, si $a$ divise $n$ alors $\\frac n a$ aussi et l'un des deux entiers $a$ et $\\frac n a$ est inférieur à $\\sqrt n$ puisque leur produit vaut $n$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plus_petit_diviseur(n):\n",
    "    k = 2\n",
    "    while k * k <= n and n % k != 0:\n",
    "        k = k + 1\n",
    "    if k * k > n:\n",
    "        return n\n",
    "    else:\n",
    "        return k"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Cette fonction a clairement une complexité en pire cas en $\\mathcal O(\\sqrt n)$, atteinte lorsque $n$ est premier. Si l'on accepte de faire, disons, un million d'opérations, on pourra trouver des diviseurs de nombres ayant une douzaine de chiffres."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.2 Factoriser un entier"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pour factoriser l'entier $n$, on recherche son plus petit diviseur premier $p$. Ensuite, on factorise (récursivement) $\\frac n p$ et on ajuste ..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La fonction __factoriser__ prend un entier non nul $n$ en paramètre. Elle renvoie la liste $[(p_1,k_1),\\ldots,(p_s,k_s)]$ où les $p_i$ sont des nombres premiers distincts et $n=\\prod_{i=1}^s p_i^{k_i}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def factoriser(n):\n",
    "    if n == 1:\n",
    "        return []\n",
    "    else:\n",
    "        p = plus_petit_diviseur(n)\n",
    "        s = factoriser(n // p)\n",
    "        return ajuster_valuation(p, s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def ajuster_valuation(p, s):\n",
    "    if s == []:\n",
    "        return [(p, 1)]\n",
    "    else:\n",
    "        (q, k) = s[0]\n",
    "        if p == q:\n",
    "            return [(p, k + 1)] + s[1:]\n",
    "        else:\n",
    "            return [(p, 1)] + s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "factoriser(123456789)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Cet algorithme de factorisation permet de factoriser des entiers ayant une quinzaine de chiffres, mais guère plus. Il faut savoir qu'il existe des algorithmes très sophistiqués permettant de factoriser des nombres d'une cinquantaine de chiffres. Dans ce qui suit nous allons nous intéresser à un algorithme relativement simple, l'algorithme $\\rho$ de Pollard, qui permet la factorisation d'entiers d'une trentaine de chiffres. Mais avant cela nous allons dire quelques mots sur les suites récurrentes à valeurs dans un ensemble fini."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2 Détection d'un cycle dans une suite récurrente"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.1 L'algorithme de Floyd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Soit $f:E\\to E$ une fonction où $E$ est un ensemble fini. Soit $a\\in E$. La suite définie par $x_0=a$ et $x_{n+1}=f(x_n)$ est ultimement périodique : il existe un entier $n_0$ et un entier $T>0$ tels que pour tout $n\\ge n_0$ on ait $x_{n+T}=x_n$. Plus précisément, les valeurs $a,f(a),f(f(a))$, etc. sont distinctes, puis on tombe à l'indice $n_0$ sur un cycle. Si l'on relie graphiquement les valeurs successives de la suite, on ne forme pas un rond (la suite ne boucle pas sur son premier terme) mais un $\\rho$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "L'algorithme de Floyd permet de trouver une période ultime de $f$ (la longueur de la boucle d'un $\\rho$) à partir de la valeur $a$ en espace constant et dans un temps de l'ordre de cette période. On initialise une variable $x$ à la valeur $f(a)$ et une variable $y$ à la valeur $f(f(a))$. Puis, tant que $x\\ne y$, on remplace $x$ par $f(x)$ et $y$ par $f(f(y))$. Le point $y$ se déplace deux fois plus vite que $x$ et finira par \"rattraper\" $x$ le long d'un cycle. \n",
    "\n",
    "On montre facilement que le nombre d'itérations effectuées par l'algorithme est inférieur à la longueur du $\\rho$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La fonction __floyd__ ci-dessous prend en paramètres la fonction $f$ et un élément $a$ de l'ensemble de départ (et d'arrivée !) de $f$. Elle renvoie un triplet $(x,p,T)$ où $x$ est un point cyclique pour $f$ (le point où la boucle du $\\rho$ se rattache à sa jambe), $p$ est l'indice de $x$ dans la suite, et $T$ est une période ultime de $f$. Néanmoins, l'entier $T$ n'est pas nécessairement la plus petite période possible."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def floyd(f, a):\n",
    "    x = f(a)\n",
    "    p = 1\n",
    "    y = f(f(a))\n",
    "    q = 2\n",
    "    while x != y:\n",
    "        x = f(x)\n",
    "        p = p + 1\n",
    "        y = f(f(y))\n",
    "        q = q + 2\n",
    "    return (x, p, q - p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Juste pour savoir, peut-on obtenir la plus petite période ? Oui, c'est facile maintenant que l'on dispose d'un point cyclique ! La fonction __floyd_opt__ renvoie LA période du point cyclique $x$ en un temps au pire double du temps d'exécution de __floyd__."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def floyd_opt(f, a):\n",
    "    x, p, T = floyd(f, a)\n",
    "    y = f(x)\n",
    "    t = 1\n",
    "    while y != x:\n",
    "        y = f(y)\n",
    "        t = t + 1\n",
    "    return (x, p, t)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.2 Un exemple"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Prenons $E=\\mathbb Z/10403\\mathbb Z$. Prenons $f:E\\to E$ définie par $f(x)=x^2+1$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def f(x):\n",
    "    return (x * x + 1) % (10403)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La fonction ci-dessous renvoie la liste des $n$ premières valeurs après l'indice $n_0$ de la suite définie par $f$ et son premier terme $x_0$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def iterer (f, n0, n, x0):\n",
    "    for k in range(n0):\n",
    "        x0 = f(x0)\n",
    "    s = [x0]\n",
    "    for k in range(n):\n",
    "        s.append(f(s[-1]))\n",
    "    return s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "s = iterer(f, 9288, 300, 6)\n",
    "plt.grid()\n",
    "plt.plot(s, 'k')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pas facile de voir quoi que ce soit. Lançons l'algorithme de Floyd."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "floyd(f, 6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3 L'algorithme $\\rho$ de Pollard"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.1 L'implémentation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Soit $c$ un entier quelconque. Pour tout entier non nul $n$, soit $f_n:x\\mapsto x^2+c$, définie sur $\\mathbb Z/n\\mathbb Z$. Appelons $t_n$ le période ultime de $f_n$ obtenue à partir d'une valeur initiale $a$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Soit $N$ un entier. Supposons que $N$ est composé. Soit $p$ un facteur premier de $N$. Si l'on applique l'algorithme de détection de cycle à $f_p$ à partir d'un point $a$, on obtient un point cyclique pour $f_p$ de période $t_p$. Il y a de grandes chances pour que $t_p\\ne t_q$, ceci pour tous les autres facteurs premiers $q$ de $N$. Par exemple, $t_p<t_q$. Soit $a$ le point cyclique trouvé pour $f_p$. Soit $b=f_p^{t_p}(a)$. On a alors $a\\equiv b[p]$ mais $a\\not\\equiv b[q]$. Et donc le pgcd de $b-a$ et de $N$ est égal à $p$. \n",
    "\n",
    "En fait, on modifie légérement l'algorithme de Floyd pour trouver ce qui nous intéresse, à savoir non pas une période, mais un pgcd différent de 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def gcd(a, b):\n",
    "    while b != 0:\n",
    "        a, b = b, a % b\n",
    "    return a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def f(x, N, c):\n",
    "    return (x * x + c) % N"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def pollard(N):\n",
    "    x = random.randint(0, N - 1)\n",
    "    c = random.randint(0, N - 1)\n",
    "    y = f (x, N, c)\n",
    "    count = 0\n",
    "    while gcd(y - x, N) == 1:\n",
    "        x = f(x, N, c)\n",
    "        y = f(f(y, N, c), N, c)\n",
    "        count += 1\n",
    "    print(\"facteur trouvé en %d itérations\" % count)\n",
    "    return gcd(y - x, N)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.2 Tests"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 1234567907 * 10987654367\n",
    "N"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pollard(N)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Encore un essai, avec un nombre d'une trentaine de chiffres sans petits diviseurs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 2 ** 101 - 1\n",
    "N"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pollard(N)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.3 Complexité de l'algorithme"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Voici le graphe de la fonction $x\\mapsto(x^2+1)\\bmod 501$. Mise à part la symétrie évidente $f(501-x)=f(x)$ due au carré, le graphe apparaît très irrégulier. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "s = [(x ** 2 + 1) % 501 for x in range(501)]\n",
    "plt.plot(s, 'k')\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Nous n'allons pas rigoureusement prouver la complexité de l'algorithme de Pollard, mais juste donner quelques idées sur la question.\n",
    "\n",
    "Soit $N$ un entier composé. Soit $p$ le plus petit facteur premier de $N$. En admettant que la fonction $x\\mapsto x^2+c$ se comporte comme une fonction \"aléatoire\", on peut voir les termes de la suite $(x_n)$ de l'algorithme de Pollard comme des boules tirées au hasard avec remise dans une urne contenant $p$ boules. La question est de savoir au bout de combien de temps on aura tiré deux fois la même boule. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Proposition__ : On considère une urne contenant $p$ boules. On effectue des tirages avec remise dans l'urne. L'espérance du nombre de tirages avant que deux boules identiques aient été tirées est $\\mathcal O(\\sqrt p)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Preuve__ : Soit $s$ le nombre de tirages avant qu'une \"collision\" ne survienne. Pour $p\\ge j\\ge 2$, on a\n",
    "\n",
    "$$P(s\\ge j) = \\frac 1 {p^{j-1}} \\prod_{1\\le i < j} (p-(i-1)) = \\prod_{1\\le i < j} (1-\\frac{i-1}{p})$$\n",
    "\n",
    "C'est en effet la probabilité que les $j-1$ premières boules tirées soirent dictinctes : $p$ possibilités pour le tirage de la première boule, puis $p-1$ pour la deuxième boule, ainsi de suite. Puis on divise par $p^{j-1}$, le nombre de possibilités. Bien évidemment, si $j>p$, $P(s\\ge j)=0$.\n",
    "\n",
    "Utilisons l'inégalité $1-x\\le e^{-x}$, vraie pour tout réel $x$. Il vient\n",
    "\n",
    "$$P(s\\ge j) \\le \\prod_{1\\le i < j} e^{-\\frac{i-1}{p}}=e^{-(j-1)(j-2)/(2p)}\\le e^{-(j-2)^2/(2p)}$$\n",
    "\n",
    "L'espérance de $s$ vérifie donc\n",
    "\n",
    "$$E(s) = \\sum_{j\\ge 1} P(s\\ge j) \\le 1 +\\sum_{j\\ge 0} e^{-j^2/2p}\\le 2 + \\int_0^\\infty e^{-x^2/2p}dx$$\n",
    "\n",
    "la dernière inégalité provenant d'une comparaison série-intégrale. Le changement de variable $x=\\sqrt{2p}t$ dans l'intégrale nous donne \n",
    "\n",
    "$$E(s)\\le 2+\\sqrt{2p}\\int_0^\\infty e^{-t^2}dt=2+\\sqrt{\\frac{p\\pi}{2}}=\\mathcal O(\\sqrt p)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "_Ainsi, la complexité en moyenne de l'algorithme de Pollard est $\\mathcal O(\\sqrt p)$ où $p$ est le plus petit facteur premier de $N$. C'est à dire, au pire, en $\\mathcal O(N^{1/4})$, mais beaucoup mieux si $N$ a par chance un facteur premier pas trop grand._"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Cette complexité est à comparer avec le $\\mathcal O(N^{1/2})$ de l'algorithme naïf. Concrètement, cela signifie que l'on pourra avec l'agorithme $\\rho$ de Pollard factoriser des entiers deux fois plus longs qu'avec l'algorithme naïf. Si l'on accepte, mettons, un million d'opérations, l'algorithme naïf factorise des entiers d'une douzaine de chiffres et l'algorithme de Pollard factorise des entiers d'environ 25 chiffres. Évidemment, Avec un peu de chance on peut faire mieux. Tentons par exemple de trouver un facteur du 10ème nombre de Fermat, $F_{10}=2^{2^{10}}+1$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def fermat(n):\n",
    "    return 2 ** (2 ** n) + 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fermat(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pollard(fermat(10))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Voici le mieux que l'on puisse faire en quelques secondes. L'entier $N$ ci-dessous est un nombre de 80 bits, produit de deux nombres premiers de 40 bits chacun."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 740514396871 * 1069728598117\n",
    "N"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pollard(N)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}