{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Résidus quadratiques\n",
    "\n",
    "Marc Lorenzi - 22 juin 2018"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import random\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams['figure.figsize'] = (8, 8)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. De quoi allons-nous parler ?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Dans tout ce notebook, $p$ désigne un nombre premier impair. Rappelons que l'ensemble $\\mathbb F_p = \\mathbb Z/p\\mathbb Z$ est un corps. L'ensemble de ses éléments non nuls, $\\mathbb F_p^*$ est un groupe multiplicatif, que nous noterons $\\mathbb G$. On notera les entiers naturels et les éléments de $\\mathbb F_p$ de façon identique, le contexte permettant sans peine de décider si des calculs sont réalisés dans $\\mathbb N$ ou dans $\\mathbb F_p$.\n",
    "\n",
    "__Définition__ : Soit $x\\in\\mathbb G$. On dit que $x$ est un __résidu quadratique__ (modulo $p$) lorsqu'il existe $y\\in\\mathbb G$ tel que $x=y^2$.\n",
    "\n",
    "Voici les questions auxquelles nous allons répondre :\n",
    "\n",
    "- Combien y a-t-il de résidus quadratiques modulo $p$ ?\n",
    "- Étant donné $x\\in\\mathbb G$, comment déterminer si $x$ est un résidu quadratique ?\n",
    "- Si $x$ est un résidu quadratique, comment trouver $y\\in \\mathbb G$ tel que $y^2=x$ ? \n",
    "\n",
    "Pour faire un parallèle intéressant, remplacez $\\mathbb F_p$ par $\\mathbb R$. Les réponses à ces trois questions sont :\n",
    "\n",
    "- Une infinité. Philosophiquement parlant, on a envie de dire \"un sur deux\".\n",
    "- On regarde si le signe de $x$ est positif.\n",
    "- On appuie sur la touche $\\sqrt\\bullet$ de la calculatrice (il y a d'autres réponses plus satisfaisantes :-))\n",
    "\n",
    "Nous noterons $\\mathcal Q$ l'ensemble des résidus quadratiques. Nous dirons de temps en temps \"résidu\" au lieu de \"résidu quadratique\"."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Nombres premiers"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Nous aurons besoin de fabriquer des nombres premiers pour tester nos fonctions. Voici tout d'abord une fonction (inefficace) permettant de vérifier si un entier est premier."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def premier(n):\n",
    "    if n <= 1: return False\n",
    "    else:\n",
    "        k = 2\n",
    "        while k * k <= n and n % k != 0: k += 1\n",
    "        return k * k > n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print([n for n in range(100) if premier(n)])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Et voici une fonction `premier suivant` qui prend un entier $n$ en paramètre et renvoie le plus petit nombre premier strictement supérieur à $n$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def premier_suivant(n):\n",
    "    p = n + 1\n",
    "    while not premier(p): p = p +1\n",
    "    return p"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "premier_suivant(123456789)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Entrons maintenant dans le vif du sujet."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. La fonction $x\\mapsto x^2$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Nous allons étudier la fonction $\\varphi:\\mathbb G\\to\\mathbb G$ définie par $x\\mapsto x^2$. Commençons par l'évidence :\n",
    "\n",
    "__Proposition__ : $\\varphi$ est un endomorphisme du groupe $\\mathbb G$. On a $\\ker\\varphi=\\{-1,1\\}$ et $Im\\ \\varphi=\\mathcal Q$.\n",
    "\n",
    "__Démonstration__ : la propriété de morphisme est évidente : $(xy)^2=x^2y^2$. L'image de $\\varphi$ est $\\mathcal Q$ par définition même d'un résidu quadratique. Seul le noyau de $\\varphi$ mérite un peu d'attention. Soit $x\\in\\mathbb G$. On a $x\\in\\ker \\varphi$ si et seulement si $x^2-1=0$, ou encore $(x-1)(x+1)=0$. N'oublions pas que $\\mathbb G\\subset \\mathbb F_p$, qui est un corps. Or, dans un corps, un produit est nul si et seulement si l'un des facteurs est nul. D'où le résultat.\n",
    "\n",
    "__Corollaire__ : $\\mathcal Q$ est un groupe.\n",
    "\n",
    "__Démonstration__ : c'est l'image d'un morphisme de groupes."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Prenons un petit exemple, $p=17$. Traçons la fonction $\\varphi$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p = 17\n",
    "print('p = ', p)\n",
    "plt.plot(range(1, p), [(k * k) % p for k in range(1, p)], 'ok', markersize=5)\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " Regardez les ordonnées des points : les résidus modulo 17 sont : 1, 2, 4, 8, 9, 13, 15 et 16. Il y en a donc 8. Pourquoi 8 ? Et pourquoi ces nombres-ci et pas les autres ? Et pourquoi les racines carrées de 2 sont-elles 6 et 11 ??? Pour bien comprendre que la question n'est pas totalement évidente, prenons un plus gros exemple."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p = premier_suivant(6000)\n",
    "print('p = ', p)\n",
    "rg = range(1, p)\n",
    "plt.plot(rg, [(k * k) % p for k in rg], 'ok', markersize=1)\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "C'est très joli. Cela montre surtout que la fonction $x\\mapsto x^2$ a un comportement très compliqué. En y regardant de plus près, on voit des motifs apparaître, des \"bouts de paraboles\". Changez la valeur de `rg` ci-dessus pour faire des zooms sur des parties du graphe. Valeurs conseillées (rien ne vous empêche d'essayer autre chose) :\n",
    "\n",
    "- 0-100\n",
    "- 700-800\n",
    "- 2500-3500\n",
    "- etc."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Facile, évidemment, d'obtenir la liste des résidus modulo $p$. La fonction `liste_residus` fait le travail."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def liste_residus(p):\n",
    "    s = set()\n",
    "    for k in range(1, p):\n",
    "            s = s.union(set([(k * k) % p]))\n",
    "    return s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(liste_residus(17))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Combien de résidus ?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Proposition__ : Il y a exactement $\\frac{p-1} 2$ résidus quadratiques modulo $p$.\n",
    "\n",
    "__Démonstration__ : Reprenons notre fonction carré \"surjectivisée\", $\\varphi :\\mathbb G\\to\\mathcal Q$. On a $\\mathbb G=\\varphi^{-1}(\\mathcal Q)=\\varphi^{-1}(\\bigcup_{a\\in \\mathcal Q}\\{a\\})=\\bigcup_{a\\in \\mathcal Q}\\varphi^{-1}(\\{a\\})$. Les ensembles $\\varphi^{-1}(\\{a\\})$ sont disjoints, donc \n",
    "\n",
    "$$(1)\\quad |\\mathbb G|=\\sum_{a\\in \\mathcal Q}|\\varphi^{-1}(\\{a\\})|$$\n",
    "\n",
    "où les deux barres verticales désignent le cardinal de l'ensemble.\n",
    "\n",
    "Soit $a\\in\\mathcal Q$. Soient $x,y\\in\\mathbb G$. Supposons $\\varphi(x)=\\varphi(y)=a$. On a alors $x^2=y^2$, ou encore $(x-y)(x+y)=0$. Donc $y=\\pm x$ (corps !). Maintenant, par définition même d'un résidu quadratique, $a$ a un antécédent $x$ par $\\varphi$. Ainsi, $\\varphi^{-1}(\\{a\\})=\\{-x,x\\}$. Peut-on avoir $-x = x$ ? Non, car dans ce cas $2x=0$ dans $\\mathbb F_p$. Mais $2\\ne 0$ car $p$ est impair, et $x\\ne 0$ également. Et $\\mathbb F_p$ est un corps. Ainsi, $|\\varphi^{-1}(\\{a\\})|=2$. Réinjectons dans l'égalité $(1)$ pour obtenir $|\\mathbb G|=\\sum_{a\\in \\mathcal Q}2=2|\\mathcal Q|$. D'où le résultat souhaité, puisque $|\\mathbb G|=p-1$.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Creusons un peu : dans $\\mathbb F_p^*$, il y a $\\frac{p-1}{2}$ résidus, et donc $\\frac{p-1}{2}$ non résidus.\n",
    "\n",
    "- Si on multiplie entre-eux deux résidus, on obtient un résidu, puisque $\\mathbb G$ est un groupe.\n",
    "- Si on multiplie un résidu par un non-résidu, on obtient un non-résidu, toujours grâce à la structure de groupe de $\\mathbb G$ (supposez le contraire ...).\n",
    "\n",
    "Mais que se passe-t-il si on multiplie deux non-résidus ? Mais oui, rassurez-vous, cela fait un résidu, mais ce n'est pas totalement immédiat. Soit $\\alpha$ un non résidu quadratique. La fonction $f:x\\mapsto \\alpha x$ est une bijection de $\\mathbb F_p^*$ sur lui-même. qui envoie les résidus sur DES non-résidus, d'après ce que nous avons dit au point numéro 2. Mais il y a exactement le même nombre de résidus et de non-résidus, donc la fonction $f$ envoie les résidus sur LES non-résidus. Les non-résidus n'ont donc, par l'injectivité de $f$, plus le choix : leur image EST un résidu. Ainsi :\n",
    "\n",
    "- Si on multiplie un non-résidu par un non-résidu, on obtient un résidu."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5. Qui est résidu ? Qui ne l'est pas ?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5.1 La formule de Legendre"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Soit $x\\in\\mathbb G$. Comment savoir si $x$ est un carré ? On pourrait calculer tous les carrés possibles et voir si on trouve $x$ dans la liste, mais soyons lucides, ce serait totalement inefficace. Il nous faudrait faire $\\mathcal O(p)$ opérations. Peut-on faire mieux ? Oui. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Proposition__ : Soit $x\\in\\mathbb G$. Alors $x$ est un résidu quadratique si et seulement si $x^{\\frac{p-1}{2}}=1$.\n",
    "\n",
    "__Démonstration__ : Les deux implications sont intéressantes.\n",
    "\n",
    "- $(\\Rightarrow)$ Supposons que $x$ est un résidu. Il existe donc $y\\in\\mathbb G$ tel que $x=y^2$. Mais alors $x^{\\frac{p-1}{2}}=y^{p-1}=1$ par le petit théorème de Fermat.\n",
    "\n",
    "- $(\\Leftarrow)$ Considérons le polynôme $P=X^{\\frac{p-1}{2}}-1\\in\\mathbb F_p[X]$. Ce polynôme est de degré $\\frac{p-1} 2$, il a donc au plus $\\frac{p-1}2$ racines. Mais nous avons déjà, par l'implication directe montrée juste avant, trouvé $\\frac {p-1} 2$ racines de $P$ : les résidus quadratiques. Il ne peut donc pas y avoir d'autres racines, donc LES racines de $P$ sont LES résidus quadratiques. En prime, $P$ est scindé et toutes ses racines sont simples.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Et si $x$ n'est pas un résidu ? Que vaut $x^{\\frac{p-1}2}$ ?\n",
    "\n",
    "__Proposition__ : Soit $x\\in\\mathbb G$. Si $x$ n'est pas un résidu quadratique alors $x^{\\frac{p-1}{2}}=-1$.\n",
    "\n",
    "__Démonstration__ : Soit $y = x^{\\frac{p-1}{2}}$. On a $y^2=x^{p-1}=1$ par le petit théorème de Fermat. Donc $y=\\pm 1$ (corps !). Mais $y\\ne 1$ puisque $x$ n'est pas un résidu."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Nous voici donc en possession d'un test simple pour savoir si $x$ est un résidu. On élève $x$ à la puissance $\\frac {p-1} 2$. La réponse est OUI si et seulement si le résultat est 1 et NON si et seulement si le résultat n'est pas 1. Il est alors $-1$.\n",
    "\n",
    "Certes, bon, d'accord, mais cette élévation à une puissance va nous obliger à faire $\\mathcal O(p)$ multiplications ? Eh non, car l'exponentiation rapide est notre amie. Cet algorithme bien connu permet de calculer une puissance $x^n$ avec une nombre d'opérations de l'ordre de $\\log n$. Je ne détaillerai pas ici son fonctionnement, vous l'étudierez (ou l'avez étudié) en cours d'informatique. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La fonction `powermod` prend en paramètres trois entiers $x,n,p$. Elle renvoie $x^n$ modulo $p$ après avoir fait $\\mathcal O(\\log n)$ multiplications dans $\\mathbb F_p$. Bien entendu, si $p$ est grand, chaque multiplication est en soi coûteuse, mais je n'entrerai pas plus avant dans les détails. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def powermod(x, n, p):\n",
    "    if n == 0: return 1\n",
    "    else:\n",
    "        y = powermod(x, n // 2, p)\n",
    "        z = (y * y) % p\n",
    "        if n % 2 == 0: return z\n",
    "        else: return (x * z) % p"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Même si l'objet de ce notebook est l'étude de $x\\mapsto x^2$, ne résistons pas à la tentation de tracer la courbe de $x\\mapsto 2^x$, histoire de tester notre fonction `powermod`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p = premier_suivant(6000)\n",
    "plt.plot([powermod(2, k, p) for k in range(1, p)], 'ok', markersize=1)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On ne distingue pas grand chose, les points ont l'air vraiment répartis de façon très \"aléatoire\". Imaginons que l'on me donne $x\\in\\mathbb F_p^*$. Supposons que $x=2^k$ est une puissance de 2. Que vaut $k$ ? Ceci est une question diificile, on ne connaît aucun algorithme permettant d'y répondre en temps raisonnable (pour $p$ grand, cela va de soi). En fait, le problème du calcul de $k$ à partir de $x$ est appelé le problème du __logarithme discret__. Il est à la base d'un système cryptographique appelé le __cryptage ElGamal__."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Remarque__ : En fait on voit quelque chose. Il y a des répétions de 1000 en 1000, à vue de nez. Pourquoi ? Eh bien parce que l'ordre de 2 dans $\\mathbb F_p^*$ doit être de l'ordre de 1000. Quel est cet ordre ?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(p)\n",
    "c = 1\n",
    "x = 2\n",
    "while x != 1:\n",
    "    x = (x * 2) % p\n",
    "    c = c + 1\n",
    "print(c)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "L'ordre de 2 dans le groupe $\\mathbb F_{6007}^*$ est 1001. La fonction $x\\mapsto 2^x$ est donc de période 1001 pour $p=6007$. Remarquez que 1001 divise $p-1=6006$, ce qui est normal en vertu du théorème de Lagrange qui dit que l'ordre d'un élément divise l'ordre du groupe.\n",
    "\n",
    "__Expérience__ : Remplacez 2 par 3, vous ne verrez plus aucune périodicité parce que l'ordre de 3 est 6006 : 3 est un __générateur__ du groupe cyclique $\\mathbb F_{6007}^*$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Revenons à notre sujet. La fonction `est_residu` prend en paramètres deux entiers $x$ et $p$, où $p$ est premier et impair, et $x$ n'est pas un multiple de $p$. Elle renvoie `True` lorsque $x$ (sa classe modulo $p$) est un résidu quadratique, et `False` sinon."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def est_residu(x, p):\n",
    "    return x % p != 0 and powermod(x, (p - 1) // 2, p) == 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Petit test avec petit $p$ ?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "[x for x in range(17) if est_residu(x, 17)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Mise à l'épreuve avec un nombre premier un peu plus grand."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p = premier_suivant(100000)\n",
    "print(p)\n",
    "len([x for x in range(1, p) if est_residu(x, p)])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Le compte est bon, puisque $\\frac{100003-1}{2}=50001$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5.2 Peut-on faire mieux ?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "En utilisant des résultats mathématiques plus fins on peut écrire un algorithme qui s'exécute plus efficacement que notre algorithme utilisant l'exponentiation rapide. Le lecteur intéressé s'orientera vers les mots clés _symbole de Legendre_, _symbole de Jacobi_, _loi de réciprocité quadratique_."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5.3 Un petit quizz"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La fonction `chi(x,p)` renvoie 0 si $x$ est multiple de $p$, et $\\pm 1$ sinon, selon que est un résidu modulo $p$ ou pas."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def chi(x, p):\n",
    "    if x % p == 0: return 0\n",
    "    elif est_residu(x, p): return 1\n",
    "    else: return -1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On trace $\\chi(x, p)$ pour $p$ premier, $200\\le p<1300$ et $0\\le x < 200$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rg = [p for p in range(200, 1300) if premier(p)]\n",
    "A = [[chi(x, p) for x in range(0,200)] for p in rg]\n",
    "plt.imshow(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Question__ : Mais que sont ces bandes jaunes verticales ?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6. Calcul des racines carrées modulo $p$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "J'ai dit un peu plus haut que résoudre dans $\\mathbb F_p$ l'équation $2^y=x$ d'inconnue $y$ était un problème très difficile pour lequel on ne possède à l'heure actuelle aucune solution satisfaisante. Qu'en est-il de la résolution de l'équation $y^2=x$ ? Eh bien ce n'est pas évident, mais ON A une solution satisfaisante à ce problème.\n",
    "\n",
    "__Problème__ : $p$ est un nombre premier impair, $x$ est un résidu quadratique modulo $p$. Trouver une racine carrée de $x$.\n",
    "\n",
    "Nous allons décrire et écrire l'algorithme de __Tonelli-Shanks__, qui résout ce problème. Voici le principe de l'algorithme. La difficulté essentielle dans ce qui va suivre est le nombre de lettres que l'on doit introduire dans la démonstration. Prenez des notes si vous voulez comprendre !"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "L'entier $p$ est impair. Il s'écrit donc de façon unique $p=1+2^eq$ où $e\\ge 1$ et $q$ est un nombre impair. On a $x^{\\frac{p-1} 2}=(x^q)^{2^{e-1}}=1$ puisque $x$ est un résidu quadratique. Posons $b=x^q$. Comme l'ensemble $\\mathcal Q$ des résidus est un groupe, on a $b\\in\\mathcal Q$. De plus, $b^{2^{e-1}}=1$.\n",
    "\n",
    "Soit $g\\in\\mathbb G$ un non-résidu. Imaginons que l'on a fabriqué un entier pair $k$ tel que $bg^k=1$. Posons illico $r=x^{\\frac{q+1}2}g^{\\frac k 2}$. On a $r^2=x^{q+1}g^k=x x^qg^k=xbg^k=x$. Ainsi, $r$ est une racine carrée de $x$.\n",
    "\n",
    "Il nous reste évidemment deux problèmes non triviaux à résoudre.\n",
    "\n",
    "1. Trouver $g$.\n",
    "2. Trouver $k$.\n",
    "\n",
    "Commençons par $g$ ..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 6.1 Étape 1 : trouver un non-résidu"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Si je tire au hasard un entier entre 1 et $p-1$, la probabilité que ce soit un non-résidu est $\\frac 1 2$. Notre algorithme de recherche de non-résidu sera probabiliste : tirer des entiers au hasard jusqu'à ce qu'on tombe sur un non-résidu. Pourquoi cet algorithme va-t-il fonctionner ? Soit $n$ un entier $\\ge 1$. Considérons l'univers $\\Omega$ des $n$-uplets d'entiers entre 1 et $p-1$ muni de la probabilité uniforme. L'événement $A$ \"aucun des éléments du $n$-uplet n'est un résidu\" a une probabilité $(\\frac 1 2)^n$. Tirages avec remise ! Ainsi, par exemple, la probabilité de ne pas avoir obtenu un non-résidu après 20 tirages aléatoires est $\\frac 1 {2^{20}}\\simeq 10^{-6}$.\n",
    "\n",
    "__Théorème__ : En prenant $n$ suffisamment grand, la probabilité d'échec est inférieure à celle qu'un yack obtienne un doctorat de physique nucléaire.\n",
    "\n",
    "__Démonstration__ : La probabilité du yack est strictement positive (quoique petite), et $\\frac 1 {2^n}$ tend vers 0 lorsque $n$ tend vers l'infini.\n",
    "\n",
    "__Note aux amis des animaux__ : je suis un ami des yacks."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pour être encore plus précis, Soit $X$ le nombre de tirages nécessaires pour obtenir un non-résidu. Pour tout entier $n\\in\\mathbb N^*$, on a $P(X=n)=\\frac 1{2^n}$. La probabilité que l'algorithme ne renvoie JAMAIS de non-résidu est donc $1-\\sum_{n=1}^\\infty \\frac 1{2^n}=1-1=0$. L'événement \"la fonction ne termine pas\" est de probabilité NULLE. Cela ne veut pas dire qu'il est impossible, mais vous ne l'observerez jamais (avec probabilité 1 :-)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La fonction `trouver_non_residu` fait le travail. Elle prend également un paramètre optionnel `stats` qui permet d'afficher le nombre de tirages aléatoires avant succès. Essayez donc de dépasser 20 tirages !"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def trouver_non_residu(p, stats=False):\n",
    "    count = 0\n",
    "    while True:\n",
    "        r = random.randint(1, p - 1)\n",
    "        count=count + 1\n",
    "        if not est_residu(r, p):\n",
    "            if stats: print(count)\n",
    "            return r"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p = premier_suivant(10000000)\n",
    "trouver_non_residu(p, stats=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 6.2 Trouver $k$, calculer une racine carrée"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Là, ça se corse. On va s'intéresser à des couples $(e_1,e_2)$ tels que \n",
    "\n",
    "$$(1)\\quad x^{e_1}g^{e_2}=1$$\n",
    "\n",
    "Je rappelle que $x$ est le résidu dont on cherche la racine carrée, et que $g$ est un non-résidu. Remarquons que $e_2$ est pair. En effet, $x^{e_1}\\in\\mathcal Q$ et $1$ aussi. Comme $\\mathcal Q$ est un groupe il en résulte que $g^{e_2}$ est un résidu. Mais $g$ n'en est pas un, et une puissance impaire d'un non résidu est un non-résidu.\n",
    "\n",
    "On part de $e_1=\\frac{p-1}{2}$ et $e_2=0$. Tant que $e_1$ est pair, on le divise par 2 et on adapte $e_2$ pour maintenir la condition $(1)$. Lorsque la condition de parité sur $e_1$ n'est plus respectée, on a $e_1=q$, LE $q$ impair tel que $p=1+2^eq$. Et donc, $x^qg^{e_2}=1$, ou encore, en posant $b=x^q$ et $k=e_2$, $bg^k=1$.\n",
    "\n",
    "À chaque division de $e_1$ par 2, on cherche un $e_2$ tel que $(1)$ reste vraie. Remarquons que $(x^{\\frac {e_1} 2}g^{\\frac {e_2} 2})^2=x^{e_1}g^{e_2}=1$. Donc (corps !), $x^{\\frac {e_1} 2}g^{\\frac {e_2} 2}=\\pm 1$.\n",
    "\n",
    "- Si c'est 1, c'est le bonheur. En divisant $e_1$ et $e_2$ par 2, on obtient un nouveau couple qui satisfait $(1)$.\n",
    "- Si c'est $-1$, c'est pas grave. Parce que comme $g$ est un non résidu, on a $g^{\\frac{p-1}2}=-1$. On prend comme nouvelle valeur de $e_2$ l'entier $\\frac {e_2} 2 + \\frac {p-1}2$. Le nouveau couple ainsi créé vérifie encore $(1)$ (exercice !)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Voici donc notre fonction qui calcule une racine carrée du résidu $x$ modulo $p$. Elle reprend les notations de la démonstration.\n",
    "\n",
    "La présentation que j'ai choisie ici a l'avantage de la clarté (si, si, je vous assure). En revanche, le code ci-dessous effectue des exponentiations rapides dont on pourrait se dispenser, au prix d'un obscurcissement du code. La fonction `racine` garde tout de même une complexité très correcte : la boucle `while` s'exécute en pire cas $\\log p$ fois (et souvent (cf l'annexe !) beaucoup moins, juste une ou deux fois, voire pas du tout). Et chaque itération calcule deux puissances qui nécessitent $\\mathcal O(\\log p)$ opérations. Notre fonction est donc de complexité en pire cas $\\mathcal O((\\log p)^2)$. Des nombres premiers de plusieurs centaines de chiffres restent tout à fait envisageables.\n",
    "\n",
    "En fait, ce qui nous bride dans ce notebook est la fonction `premier` qui, elle, s'exécute en temps $\\mathcal O(\\sqrt p)$ et ne supportera pas d'être exécutée sur des nombres premiers de plus d'une quinzaine de chiffres. Mais là encore, il existe des algorithmes bien plus performants permettant de s'attaquer à des nombres premiers ayant des milliers de chiffres (voir le notebook sur l'algorithme de Miller-Rabin). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def racine(x, p):\n",
    "    g = trouver_non_residu(p)\n",
    "    q = (p - 1) // 2\n",
    "    e1, e2 = (q, 0)\n",
    "    while e1 % 2 == 0:\n",
    "        e1, e2 = (e1 // 2, e2 // 2)\n",
    "        z = (powermod(x, e1, p) * powermod(g, e2, p)) % p\n",
    "        if z == p - 1: e2 = e2 + q\n",
    "    return (powermod(x, (e1 + 1) // 2, p) * powermod(g, e2 // 2, p)) % p"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Prenons pour tester un nombre premier pas trop petit (mais pas trop grand :-))."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p = premier_suivant(1000000000)\n",
    "print(p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Prenons un résidu modulo $p$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "est_residu(123456789, p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Il y avait une chance sur deux que ça marche. Calculons une racine carrée de ce nombre."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "racine(123456789, p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Vérifions ?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(_ * _) % p"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Youpi."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 6.3 Le graphe !"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Traçons le graphe de la \"fonction\" racine carrée modulo $p$. Chaque résidu quadratique a deux racines, qui sont bien entendu opposées (dans $\\mathbb F_p$). En noir, l'une des deux. En rouge, l'autre. Ce n'est pas moi qui choisis les couleurs, c'est la fonction `racine`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def tracer_racine(p):\n",
    "    rg = [k for k in range(1, p) if est_residu(k, p)]\n",
    "    s1 = [racine(k, p) for k in rg]\n",
    "    s2 = [p - racine(k, p) for k in rg]\n",
    "    plt.plot(rg, s1, 'ok', markersize=1)\n",
    "    plt.plot(rg, s2, 'or', markersize=1)\n",
    "    plt.grid()\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p = premier_suivant(6000)\n",
    "print(p)\n",
    "tracer_racine(p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pour tout $x$ entre 1 et $p-1$, il y a\n",
    "\n",
    "- zéro point d'abscisse $x$ si $x$ n'est pas un résidu.\n",
    "- Un point rouge et un point noir d'abscisse $x$ sinon.\n",
    "\n",
    "Et bien entendu, le symétrique orthogonal de ce dessin par rapport à la droite d'équation $y=x$ est le dessin que nous avons fait à la section 3.\n",
    "\n",
    "\n",
    "Je vous laisse compter les points :-)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 6.4 Annexe : Le nombre moyen d'itérations dans l'algorithme de Tonelli-Shanks"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Dans ce qui va suivre, j'utiliserai un certain nombre de résultats de théorie des nombres absolument pas triviaux. Nous allons calculer le nombre moyen d'itérations dans la boucle `while` de la fonction `racine`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Posons, comme nous l'avons fait plus haut, $p=1 + 2^e q$ où $e\\ge 1$ et $q$ est un entier impair. Juste avant le début de la boucle `while`, on a $e_1 = 2^{e-1}q$. Et à chaque itération, $e_1$ est divisé par 2. Ainsi, le nombre d'itérations de la boucle est $e-1$.\n",
    "\n",
    "__Bilan : Si $p\\equiv 1[2^e]$ et $p\\not\\equiv 1[2^{e+1}]$, le nombre d'itérations est $e-1$.__\n",
    "\n",
    "Par exemple, il y a zéro itération si et seulement si $p\\equiv 1[2]$, ce qui est toujours vrai, et $p\\not\\equiv 1[4]$.\n",
    "\n",
    "Un nombre impair est congru soit à 1, soit à 3 modulo 4. Mais combien y a-t-il de nombres premiers congrus à 3 modulo 4 ? Eh bien une infinité. Ceci peut se montrer de façon élémentaire par l'absurde. Et les nombres premiers congrus à un modulo 4 ? Il y en a une infinité aussi. Cela peut aussi se montrer, mais un peu moins facilement que pour les nombres congrus à 3 modulo 4.\n",
    "\n",
    "__Posons la question autrement__ : parmi tous les nombres premiers, quelle est la \"proportion\" des nombres premiers congurs à 3 modulo 4 ? L'idéal serait de pouvoir définir une probabilité uniforme sur l'ensemble $\\mathcal P$ des nombres premiers, mais ceci est impossible. Il est en revanche possible, je ne le ferai pas, de définir une \"quasi-probabilité\" (une mesure finiment additive) qui réponde à nos besoins. Nous la noterons $\\mu$. Un résultat très général, dû à Dirichlet et de la Vallée Poussin, dit :\n",
    "\n",
    "__Théorème__ : Soient $a\\ge 2$ et $1\\le b\\le a$ tels que $pgcd(a,b)=1$. Notons $\\mathcal P_{a,b}$ l'ensemble des nombres premiers $p$ tels que $p\\equiv b [a]$. Alors, $\\mu(\\mathcal P_{a,b})= \\frac 1{\\varphi(a)}$ où $\\varphi$ est la fonction indicatrice d'Euler.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Exemple__ : Pour $a=4$ et $b=1$ ou 3, on a $\\mu(\\mathcal P_{4,1})=\\mu(\\mathcal P_{4,3})=\\frac 1 2$.\n",
    "\n",
    "__Notons__ $X:\\mathcal P\\to \\mathbb N$ la fonction qui à tout nombre premier $p$ (impair) associe le nombre d'itérations de la boucle `while`. $X$ n'est pas une variable aléatoire mais il s'en faut de peu. En utilisant les notations standard des probabilités, nous avons donc $\\mu(X=0)=\\frac 1 2$. Pour la \"moitié\" des nombres premiers, la boucle ne s'exécute pas !"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Prenons maintenant $p\\equiv 1[4]$. Il y a deux possibilités : $p\\equiv 1[8]$ ou $p\\equiv 5[8]$. En termes d'ensembles, cela signifie que $p\\in\\mathcal P_{8,1}$ ou $p\\in\\mathcal P_{8, 5}$. $\\frac {p-1}4$ est impair si et seulement si on est dans le second cas.\n",
    "\n",
    "Ainsi, $\\mu(X=1)=\\mu(\\mathcal P_{8,5})=\\frac 1 {\\varphi(8)}=\\frac 1 4$. Je ne détaillerai pas le reste, mais on peut montrer assez facilement grâce au théorème de Dirichlet que pour tout entier $n$, $\\mu(X=n)=\\frac 1{2^{n+1}}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Quelle est l'\"espérance\" de $X$ ? Eh bien $E(X)=\\sum_{n=0}^\\infty n\\mu(X=n)=\\sum_{n=0}^\\infty\\frac{n}{2^{n+1}}$. Il nous reste à calculer la somme de cette série. Posons $f(x)=\\frac 1 2 \\sum_{n=0}^\\infty (\\frac x 2)^n$. Cette série géométrique est convergente lorsque $|x|<2$ et sa somme est $f(x)=\\frac 1{2-x}$. Admettons ici que $f$ est dérivable sur $]-2,2[$ et que sa dérivée est obtenue en dérivant la série terme à terme (résultat bien connu sur les séries entières). On a alors pour tout $x\\in]-2,2[$, $f'(x)=\\frac 1 2\\sum_{n=1}^\\infty \\frac {nx^{n-1}}{2^n}$ et on a donc $f'(1)=E(X)$. Or, pour tout $x\\in]-2,2[$, $f'(x)=\\frac 1{(2-x)^2}$. Donc :\n",
    "\n",
    "$$E(X)=1$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Théorème : le nombre moyen d'itérations de la boucle `while` est ... 1.__"
   ]
  },
  {
   "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": 2
}
