{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Tests de primalité\n",
    "\n",
    "Marc Lorenzi\n",
    "\n",
    "Juillet 2016"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## 1 Le test naïf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La première idée qui vient à l'esprit pour savoir si un nombre $p$ est premier est de chercher le plus petit diviseur de $p$ (autre que $1$, évidemment). En remarquant que si $a$ divise $p$ alors $\\frac p a$ divise $p$, on voit que si un entier $p$ est composé alors $p$ possède un diviseur inférieur à $\\sqrt p$.\n",
    "\n",
    "La fonction ci-dessous renvoie le plus petit diviseur d'un entier $p$. Elle effectue dans le pire cas, en particulier lorsque $p$ est premier, $\\mathcal O(\\sqrt p)$ opérations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plus_petit_diviseur(p):\n",
    "    k = 2\n",
    "    while k * k <= p and p % k != 0: k += 1\n",
    "    if k * k > p: return p\n",
    "    else: return k"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plus_petit_diviseur(91)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On en déduit immédiatement une fonction prenant en paramètre un entier $p$ et renvoyant True si et seulement si $p$ est premier."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def is_prime(p):\n",
    "    if p <= 1: return False\n",
    "    else: return plus_petit_diviseur(p) == p"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print([p for p in range(100) if is_prime(p)])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Quels sont les plus grands nombres $p$ dont on puisse décider la primalité avec la fonction __est_premier__ en un temps raisonnable ? En acceptant de faire, disons, de l'ordre de un milliard d'opérations, on peut tester des nombres $p$ inférieurs à $10^{18}$ soit des nombres d'une vingtaine de chiffres.\n",
    "\n",
    "Comment faire, alors, pour savoir si un nombre de 100, 1000, 1000000 de chiffres est premier ? Certainement pas en utilisant cet algorithme. Il va valloir essayer autre chose."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2 Première tentative ... et échec"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.1 Le petit théorème de Fermat"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La petit théorème de Fermat affirme que si $p$ est un nombre premier alors pour tout entier $a$ on a $a^p\\equiv a [p]$.\n",
    "\n",
    "Soit $p$ un entier naturel. Soit $a$ un entier. Nous dirons que $a$ est un témoin de non-primalité de $p$ lorsque $a^p\\not\\equiv a [p]$. Si un entier $p$ possède un témoin de non-primalité alors, d'après le petit théorème de Fermat, il n'est pas premier.\n",
    "\n",
    "La réciproque est-elle vraie ? Cela vaut la peine de tenter quelques expériences."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.2 Le test de Fermat"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Définissons tout d'abord une fonction d'exponentiation rapide modulo $p$. La fonction ci-dessous prend en paramètres trois entiers $x, n, p$. Elle renvoie $x^n$ modulo $p$ en effectuant un nombre de multiplications de l'ordre de $\\log p$. D'où son appellation de \"rapide\". C'est la clé de tout ce qui va suivre : on peut élever un nombre à des puissances énormes quasi-instantanément."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def power_mod(x, n, p):\n",
    "    z = 1\n",
    "    m = n\n",
    "    y = x\n",
    "    while m != 0:\n",
    "        if m % 2 == 1: z = (z * y) % p\n",
    "        m = m // 2\n",
    "        y = (y * y) % p\n",
    "    return z"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "power_mod(2, 3 ** 1000, 561)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Définissons ensuite la fonction __temoin_fermat__. Cette fonction prend un entier $a$ en paramètre. Elle renvoie True si l'entier $a$ est un témoin de non-primalité de $p$, et False sinon. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def temoin_fermat(a, p):\n",
    "    x = power_mod(a, p, p)\n",
    "    return x != a"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Voici enfin la fonction __test_fermat__. Cette fonction prend un entier $p$ en paramètre. Elle tire 20 nombres au hasard entre 2 et $p-1$. Si l'un d'entre eux est un témoin de non primalité de $p$, on sait que $p$ est composé. Si aucun n'est un témoin de primalité de $p$ ... eh bien on ne sait pas trop."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import random"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def test_fermat(p):\n",
    "    c = 20\n",
    "    while c > 0:\n",
    "        a = random.randint(2, p - 1)\n",
    "        c = c - 1\n",
    "        if temoin_fermat(a, p): return False\n",
    "    return True"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "test_fermat(13)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def is_fermat_prime(p):\n",
    "    if p <= 1: return False \n",
    "    elif p == 2: return True\n",
    "    else: return test_fermat(p) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print([p for p in range(100) if is_fermat_prime(p)])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.3 Les nombres de Carmichael"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Tout a l'air d'aller pour le mieux. Essayons jusqu'à 10000 ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print([p for p in range(10000) if is_fermat_prime(p) and not is_prime(p)])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Malheur ! Il y a 7 nombres entre 2 et 10000 qui ne sont pas premiers et qui, pourtant, passent le test. Nous n'y pouvons pas grand chose : ils font partie de la famille des nombres de Carmichaël. Bien que \"rares\" il y en a quand même une infinité : ce résultat a été démontré en 1994 par Alford, Granville et Pomerance. Pour les courageux, voici [l'article original](paper95.pdf)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Tout est-il perdu pour autant ? Non, car en adaptant ce qui vient d'être fait on obtient quelque chose qui fonctionne : l'algorithme de Miller-Rabin. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3 Miller-Rabin"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.1 Le principe"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Soit $p$ un nombre impair, $p\\ge 3$. On écrit tout d'abord $p=1 + q 2^t$ où $q$ est un nombre impair et $t$ un entier naturel, $t\\ge 1$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def factor2(p):\n",
    "    q = p - 1\n",
    "    t = 0\n",
    "    while q % 2 == 0:\n",
    "        q = q // 2\n",
    "        t = t + 1\n",
    "    return (q, t)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p = 10013\n",
    "print(is_prime(p))\n",
    "q, t = factor2(p)\n",
    "print(q, t)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Supposons que $p=1+q2^t$ est un nombre premier. Soit $a$ un entier entre 2 et $p-1$. D'après le petit théorème de Fermat, on a $a^{p-1}\\equiv 1 [p]$, ou encore $a^{q2^t}\\equiv 1[p]$. Soit $x=a^{q2^{t-1}}$. On a $x^2\\equiv 1$ mais comme $p$ est premier, l'anneau $\\mathbb Z / p\\mathbb Z$ est un corps. Donc, de $x^2-1=(x-1)(x+1)=0$ on déduit $x=\\pm 1$. Si $x\\equiv 1$, soit $y=a^{q2^{t-2}}$. On a $y^2\\equiv 1$ donc, encore une fois, $y\\equiv\\pm 1$. Ainsi de suite ... Ainsi, si $p$ est premier, alors\n",
    "- $a^q\\equiv 1[p]$, ou alors\n",
    "- il existe $e\\in[0,t-1]$ tel que $a^{q2^e}\\equiv -1[p]$.\n",
    "\n",
    "Nous définissons donc une fonction __temoin_miller__ qui prend en paramètres les entiers $a,q,t,p$, et qui renvoie True si aucune des deux conditions ci-dessus n'est réalisée, et False sinon. Si, pour un certain entier $a$, la fonction renvoie False, alors $a$ est dit témoin de Miller (de non primalité) pour $p$ : l'existence d'un tel $a$ prouve que $p$ n'est pas premier."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def temoin_miller(a, q, t, p):\n",
    "    b = power_mod(a, q, p)\n",
    "    if b == 1: return False\n",
    "    e = 0\n",
    "    while b != 1 and b != p - 1 and e <= t - 2:\n",
    "        b = (b * b) % p\n",
    "        e = e + 1\n",
    "    return b != p - 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.2 Cette fois-ci, cela fonctionne"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La fonction __temoin_miller__ est un raffinement de la fonction __temoin_fermat__. Et alors ? Eh bien nous avons le résultat suivant (que nous admettrons) :\n",
    "\n",
    "Soit $p\\ge 3$ un nombre impair non premier. Soit $L$ l'ensemble des entiers entre 2 et $p-1$ _qui ne sont pas_ des témoins de Miller de non primalité pour $p$. Alors, le cardinal de $L$ vérifie $|L|\\le \\frac{p-1}{4}$.\n",
    "\n",
    "Ainsi,\n",
    "- Si $p$ est premier, la fonction __temoin_miller__ renvoie toujours False\n",
    "- Si $p$ est composé, la fonction __temoin_miller(a)__, avec en argument un entier $a$ aléatoire entre 2 et $p-1$, renvoie False avec une probabilité inférieure à $1/4$.\n",
    "\n",
    "En admettant que les non-témoins de Miller sont uniformément répartis entre 2 et $p-1$, si l'on exécute la fonction __temoin_miller__ en prenant $c$ nombres $a$ tirés au hasard, la probabilité que l'on obtienne tout le temps False alors que $p$ est composé est donc inférieure à $(\\frac 1 4)^c$. Pour $c=20$, par exemple, cette probabilité est inférieure à : "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "0.25 ** 20"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "À titre de comparaison, c'est environ la probabilité de gagner cent mille fois de suite au loto."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p = 8191 * 23\n",
    "print(p)\n",
    "print(is_prime(p))\n",
    "q, t = factor2(p)\n",
    "s = [a for a in range(2, p) if not temoin_miller(a, q, t, p)]\n",
    "print(len(s) / float(p))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Regroupons tout cela dans une unique fonction."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def test_miller(p):\n",
    "    q, t = factor2(p)\n",
    "    c = 20\n",
    "    while c > 0:\n",
    "        a = random.randint(2, p - 1)\n",
    "        c = c - 1\n",
    "        if temoin_miller(a, q, t, p): return False\n",
    "    return True"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def is_probably_prime(p):\n",
    "    if p <= 1: return False \n",
    "    elif p == 2: return True\n",
    "    elif p % 2 == 0: return False\n",
    "    else: return test_miller(p) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "s1 = [p for p in range(10000) if is_probably_prime(p)]\n",
    "s2 = [p for p in range(10000) if is_prime(p)]\n",
    "print(s1 == s2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p = 1195068768795265792518361315725116351898245581\n",
    "print(is_probably_prime(p))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = 24444516448431392447461\n",
    "b = 48889032896862784894921\n",
    "print(is_probably_prime(a))\n",
    "print(is_probably_prime(b))\n",
    "print(a * b - p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.3 Test sur de \"grands\" nombres"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Posons, pour tout entier $n$, $M_n=2^n-1$. Il n'est pas difficile de montrer que si $M_n$ est premier, alors $n$ est aussi un nombre premier."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def mersenne(p):\n",
    "    return 2 ** p - 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for p in range(103):\n",
    "    if is_prime(p):\n",
    "        print(p, mersenne(p))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Les nombres de Mersenne deviennent vite grands lorsque $n$ augmente. Par exemple, le 1000 ème nombre de Mersenne possède 1000 chiffres en base 2, soit environ 300 chiffres décimaux."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from math import log\n",
    "1000 * log(2) / log(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(mersenne(1000))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Quels sont, les nombres de Mersenne $M_p$ probablement premiers, pour $p$ premier, $p\\le 2000$ ?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "s = [p for p in range(2000) if is_prime(p)]\n",
    "print(s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "s1 = [(p, mersenne(p)) for p in s]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "s2 = [p for (p, n) in s1 if is_probably_prime(n)]\n",
    "print(s2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(mersenne(1279))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Un dernier essai, avec $M_{4253}$ qui est le premier nombre de Mersenne premier ayant plus de 1000 chiffres."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p = mersenne(4253)\n",
    "print(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "is_probably_prime(p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "À titre de comparaison, l'algorithme naïf sur une machine effectuant $10^{10}$ opérations par secondes demanderait un temps de ... $10^{3000}$ millénaires."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## 4 Fabriquer de grands nombres premiers"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.1 Créer un nombre premier de taille donnée"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Comment fabriquer un nombre premier ayant une certaine taille, disons $L$ bits ? Un algorithme naïf consite à tirer un entier $n$ au hasard entre $2^{L-1}$ et $2^L$. Puis, tant que $n$ n'est pas premier, à incrémenter $n$. Quoique naïf, cet algorithme fonctionne plutôt bien jusqu'à des valeurs de $L$ de l'ordre du millier. Je n'entrerai pas plus avant dans les détails."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def random_prime(L):\n",
    "    a = 2 ** (L - 1)\n",
    "    b = 2 * a - 1\n",
    "    n = random.randint(a, b)\n",
    "    if n % 2 == 0: n = n + 1\n",
    "    while not is_probably_prime(n):\n",
    "        n = n + 2\n",
    "    return n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p = random_prime(30)\n",
    "print(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(is_probably_prime(p))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.2 Créer des clés RSA"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Le chiffrement crypographique RSA nécessite la création de clés. Créer une clé de taille $N$ (bits) consiste à créer deux nombres premiers de taille $N/2$. Vous désirez créer votre clé RSA-2048 ? rien de plus facile ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p = random_prime(1024)\n",
    "q = random_prime(1024)\n",
    "n = p * q\n",
    "print(n)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La connaissance de $n$ permet le cryptage de messages. La connaissance de $p$ et $q$ permet le décryptage. Avis aux amateurs : factorisez le nombre ci-DESSOUS. Parce que ci-dessus c'est facile : print(p, q) :-)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = random_prime(1024) * random_prime(1024)\n",
    "print(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
}