{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Le nième nombre premier\n",
    "\n",
    "Marc Lorenzi - 5 juin 2018"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "from math import log, sqrt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams['figure.figsize'] = (12, 6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pour $n\\ge 1$, notons $p_n$ le $n$-ième nombre premier. On a $p_1=2$, $p_2=3$, $p_3=5$. Mais que vaut $p_n$ pour $n$ quelconque ? Peut-on avoir un équivalent de $p_n$ ? Un développement asymptotique de $p_n$ ?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Nombres premiers"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Les fonctions ci-dessous ne sont pas efficaces mais elles suffiront pour faire quelques simulations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La fonction `smallest_divisor` prend un entier $n\\ge 2$ en paramètre. Elle renvoie le plus petit diviseur premier de $n$. Sa complexité est $O(\\sqrt(n))$, ce qui ne permet de l'utiliser que pour des entiers d'au plus une dizaine de chiffres."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def smallest_divisor(n):\n",
    "    k = 2\n",
    "    while n % k != 0 and k * k <= n: k += 1\n",
    "    if k * k > n: return n\n",
    "    else: return k"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "smallest_divisor(1234567)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Un entier naturel $n$ est premier si et seulement si $n\\ge 2$ et $n$ n'admet que 1 et $n$ comme diviseurs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def is_prime(n):\n",
    "    if n >= 2: return smallest_divisor(n) == n\n",
    "    else: return False"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(is_prime(1234567891))\n",
    "print(is_prime(1234567893))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print([n for n in range(100) if is_prime(n)])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La fonction `primes_from_to` prend deux entiers $m$ et $n$ en paramètres et renvoie la liste des nombres premiers dans l'intervalle $[m,n]$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def primes_from_to(m, n):\n",
    "    s = []\n",
    "    for k in range(m, n + 1):\n",
    "        if is_prime(k):\n",
    "            s.append(k)\n",
    "    return s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "primes = primes_from_to(1, 100)\n",
    "print(primes)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Le théorème des nombres premiers"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pour tout $x\\ge 1$, notons $\\pi(x)$ le nombre de nombres premiers inférieurs ou égaux à $x$. On a $\\pi(2)=1,\\pi(3)=2,\\pi(4)=2,\\pi(5)=3$, etc.\n",
    "\n",
    "__Théorème__ : On a\n",
    "\n",
    "$$\\pi(x)\\sim\\frac x {\\ln x}$$\n",
    "\n",
    "Cet équivalent a tout d'abord été conjecturé par Adrien Marie Legendre à la fin du XVIIIe siècle. L'histoire de sa démonstration est plutôt longue, mais s'est achevée de façon heureuse : le théorème a été démontré indépendamment en 1896 par Jacques Hadamard et Charles Jean de la Vallée Poussin. Sa démonstration originale fait intervenir des résultats sur les fonctions de variables complexe. Une démonstration \"élémentaire\" (c'est à dire ne passant pas par le plan complexe) en a été faite en 1948 par, entre-autres, Paul Erdös. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Voyons un peu ce que donne concrètement tout cela. La fonction`pi` prend un entier $n$ en paramètre et renvoie la liste des $\\pi(k)$ pour $k=2,\\ldots,n$. Un compteur est mis au départ à 0. Chaque fois qu'un nombre premier est rencontré, on incrémente le compteur. Évidemment, cette fonction ne fonctionne que pour des petites valeurs de $n$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def pi(n):\n",
    "    s = []\n",
    "    c = 0\n",
    "    for k in range(2, n + 1):\n",
    "        if is_prime(k): c += 1\n",
    "        s.append(c)\n",
    "    return s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(pi(100))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Traçons sur un même graphe $\\pi(x)$ et $\\frac x {\\ln x}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 10000\n",
    "ks = range(2, N + 1)\n",
    "plt.plot(ks, pi(N), 'k')\n",
    "plt.plot(ks, [k / log(k) for k in ks], 'r')\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ce n'est pas mal mais ce n'est pas si bien que cela. En fait, un « meilleur équivalent » de $\\pi(x)$ est \n",
    "\n",
    "$$Li(x)=\\int_2^x\\frac{dt}{\\ln t}$$\n",
    "\n",
    "Mais c'est quoi un meilleur équivalent ? Après tout, tous les équivalents sont ... équivalents ? Cela signifie que \n",
    "\n",
    "$$\\pi(x)=Li(x)+o(Li(x))$$\n",
    "\n",
    "où le petit o est un petit o du petit o qui apparaît dans\n",
    "\n",
    "$$\\pi(x)=\\frac x {\\ln x}+o(\\frac x{\\ln x})$$\n",
    "\n",
    "Ma phrase est rigolote mais incompréhensible. Pour être plus clairs :\n",
    "\n",
    "$$\\pi(x)-Li(x)=o\\left(\\pi(x)-\\frac x {\\ln x}\\right)$$\n",
    "\n",
    "Lowell Schoenfeld a démontré en 1976 que si une certaine conjecture appelée __l'hypothèse de Riemann__ est vérifiée (voir plus loin !), alors, pour tout $x\\ge 2657$ on a\n",
    "$$|\\pi(x)-Li(x)|\\le \\frac{\\sqrt x\\ln x}{8\\pi}$$\n",
    "\n",
    "\n",
    "Oui, bon, mais cette intégrale ne peut pas s'évaluer à l'aide des fonctions usuelles ? Eh bien non, mais il s'avère qu'on a un développement asymptotique de $Li(x)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Proposition__ : On a\n",
    "\n",
    "$$Li(x)=\\sum_{k=1}^n(k-1)!\\frac{x}{\\ln^k x} + o\\left(\\frac{x}{\\ln^n x}\\right)$$\n",
    "\n",
    "__Démonstration__ : Intégrer par parties."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Remarquons que \n",
    "\n",
    "$$\\frac{\\sqrt x\\ln x}{\\frac{x}{\\ln^n x}}=\\frac{\\ln^{n-1}x}{\\sqrt x}\\to 0$$\n",
    "\n",
    "orsque $x$ tend vers l'infini. En d'autres termes, le majorant de Schoenfeld est meilleur que toute approximation que nous pourrions obtenir à l'aide du développement asymptotique de $Li(x)$. Difficile de faire mieux ..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La fonction `approx_Li` prend un réel $x$ et un entier $n$ en paramètre et renvoie la somme des $n$ premiers termes du développement asymptotique ci-dessus."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def approx_Li(x, n):\n",
    "    s = 0\n",
    "    l = log(x)\n",
    "    p = 1\n",
    "    for k in range(1, n + 1):\n",
    "        s = s + p * x / l ** k\n",
    "        p = p * k\n",
    "    return s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 10000\n",
    "plt.plot(pi(N), 'k')\n",
    "plt.plot([0,0] + [approx_Li(k, 4) for k in range(2, N)], 'r')\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "C'est remarquable. Reprenez les valeurs dans la cellule ci-dessus et faites vos propres tests. C'est tout l'intérêt des notebooks ...\n",
    "\n",
    "Un graphique encore plus parlant est celui de l'erreur commise en approchant $\\pi(x)$ par $Li(x)$. Le voici."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 100000\n",
    "s1 = pi(N)\n",
    "s2 = [approx_Li(k, 5) for k in range(2, N + 1)]\n",
    "s = [s1[k] - s2[k] for k in range(3, len(s1))]\n",
    "plt.plot(s, 'k')\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Remarquez la petitesse de l'erreur : au plus une quarantaine. La majoration de Schoenfeld donne quant à elle :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sqrt(1e5) * log(1e5) / (8 * 3.1415926536)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Je vous laisse essayer avec $N=10^6$. Remarquons que l'erreur commise est plus ou moins croissante. C'est normal : regardez le petit $o$ dans le développement asymptotique, c'est quelque chose de négligeable devant un truc qui tend vers l'infini ! Si vous espériez obtenir une somme \"convergente\", vous serez déçus. Les nombres premiers ne se laissent pas attraper aussi facilement."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Le $n$ième nombre premier"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.1 Un équivalent de $p_n$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Nous y voilà. Que vaut $p_n$, le $n$-ième nombre premier ?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Théorème__ : On a\n",
    "\n",
    "$$p_n\\sim n\\ln n$$\n",
    "\n",
    "__Démonstration__ : on a \n",
    "\n",
    "$$\\pi(p_n) = n\\sim\\frac{p_n}{\\ln p_n}$$\n",
    "\n",
    "Ainsi, $p_n\\sim n \\ln p_n$. Mais on peut prendre des logarithmes dans des équivalents qui tendent vers l'infini. Donc \n",
    "\n",
    "$$\\ln n\\sim\\ln\\left(\\frac{p_n}{\\ln p_n}\\right)=\\ln p_n-\\ln\\ln p_n\\sim\\ln p_n$$\n",
    "\n",
    "d'où le résultat. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Une petite simulation ? En noir, le $n$ième nombre premier. En rouge, $n\\ln n$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 100000\n",
    "primes = primes_from_to(2, N)\n",
    "s = [t * log(t) for t in range(2, len(primes))]\n",
    "plt.plot(primes, 'k')\n",
    "plt.plot([0,0] + s, 'r')\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Oui, bon, c'est un peu comme $\\pi(n)$. Pas mal mais peut mieux faire. Vu que $\\pi(p_n)=n\\sim Li(p_n)$ est une meilleure approximation, on se doute que si on calcule la réciproque de la fonction $Li$ on sera sur la bonne voie. Cette fonction est appelée $ali$ par les experts. On a donc $p_n\\sim ali(n)$. Mais comment calculer la réciproque d'une fonction que l'on ne sait pas elle-même calculer ?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.2 Quelques résultats"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Voici sans preuve (on s'en doute, ce n'est pas trivial), un développement asymptotique de $p_n$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Théorème__ : On a\n",
    "\n",
    "$$p_n = n\\ln n \\left(1+\\sum_{k=1}^N\\frac{P_{k-1}(\\ln\\ln n)}{\\ln^kn}+O\\left(n\\left(\\frac{\\ln\\ln n}{\\ln n}\\right)^N\\right)\\right)$$\n",
    "\n",
    "où les $P_k$ sont des polynômes, de degré $\\le k$ si $k\\ge 1$. On a en particulier :\n",
    "\n",
    "- $P_0=X-1$\n",
    "- $P_1=X-2$\n",
    "- $P_2=\\frac 1 2(X^2-6X+11)$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Testons."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def f(n):\n",
    "    if n <= 2: return 0\n",
    "    else:\n",
    "        u = log(n)\n",
    "        x = log(u)\n",
    "        P0 = x - 1\n",
    "        P1 = x - 2\n",
    "        P2 = -(x ** 2 - 6 * x + 11) / 2\n",
    "        return n * (u + P0 + P1 / u  + P2 / (2 * u ** 2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Par exemple, le $2\\ 10^{17}$ème nombre premier est (paraît-il :-)) 8512677386048191063."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f(2e17)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Une erreur relative de l'ordre de $10^{-4}$ ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 10000\n",
    "primes = primes_from_to(2, N)\n",
    "s = [f(k + 1) for k in range(len(primes))]\n",
    "plt.plot(primes, 'k')\n",
    "plt.plot(s, 'r')\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Les deux courbes se superposent quasiment. Enfin, traçon le graphe de l'erreur commise."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 100000\n",
    "primes = primes_from_to(2, N)\n",
    "s = [f(t + 1) for t in range(len(primes))]\n",
    "s1 = [abs(primes[k] - s[k]) for k in range(2, len(primes))]\n",
    "plt.plot(s1, 'k')\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "L'erreur ne dépasse pas ici quelques centaines."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.3 Une majoration théorique de l'erreur"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "L'__hypothèse de Riemann__ affirme que les zéros complexes non triviaux d'une certaine fonction (la fonction zeta) sont tous situés sur une certaine droite. Cette conjecture a été énoncée par B. Riemann en 1859, elle est à ce jour non démontrée. Cela dit, beaucoup d'énoncés d'arithmétique commencent par \"Si l\"hypothèse de Riemann est vérifiée alors\" ... Conclusion :\n",
    "\n",
    "- Si quelqu'un démontre un jour que l'hypothèse de Riemann est vraie nous aurons plein de théorèmes démontrés d'un coup.\n",
    "- Si quelqu'un démontre un jour que l'hypothèse de Riemann est fausse nous aurons plein de théorèmes à mettre à la poubelle, ou, plutôt, plein de théorèmes qu'il nous faudra essayer de prouver autrement."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Théorème__ : Si l'Hypothèse de Riemann est vérifiée, alors pour tout $n\\ge 11$\n",
    "\n",
    "$$|p_n-ali(n)|\\le \\frac 1 \\pi \\sqrt n (\\ln n)^{\\frac 5 2}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Conclusion"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Un notebook sur la fonction $\\zeta$ semble tout indiqué."
   ]
  }
 ],
 "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}