{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# La millionième décimale de $\\pi$\n",
    "\n",
    "Marc Lorenzi\n",
    "\n",
    "15 janvier 2019"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "_Il fit la cuve en métal fondu. Elle faisait 5 mètres d'un bord à l'autre et était entièrement ronde. Elle faisait 2 mètres et demi de haut et on pouvait mesurer sa circonférence avec un cordon de 15 mètres (Ancien Testament, 1 Rois, 7:23)_\n",
    "\n",
    "_Exemple pour faire un champ circulaire de 9 khet. Combien est-ce en surface ? Enlève lui $\\frac 1 9$, c'est à dire 1. Le reste est 8. Fais la multiplication : $8\\times 8$, cela devient $64$. Sa quantité c'est, en surface, 64 setat (Papyrus Rhind, problème 50)_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sympy import *\n",
    "import matplotlib.pyplot as plt\n",
    "import math\n",
    "import gc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 0. Remarques"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Nous allons dans ce notebook effectuer des calculs intensifs sur des nombres ayant un grand nombre de décimales. Lors de l'installation du module Python `sympy`, le module `mpmath` a été aussi installé. Ce module permet de calculer en virgule flottante à des précisions quelconques. Une remarque s'impose cependant.\n",
    "\n",
    "`mpmath` utilise par défaut les entiers longs de Python pour faire ses calculs. Il n'est donc pas extrêmement efficace. Pour ceux d'entre-vous qui voudraient calculer très vite, il convient d'installer `gmpy`. Cette installation se fait en deux temps :\n",
    "\n",
    "1. Installer la bibliothèque C `gmp` (_GNU MultiPrecision Library_). Plusieurs possibilités :\n",
    "\n",
    "    - Votre ordinateur tourne sous Windows. Je ne sais pas comment il faut faire.\n",
    "    - Votre ordinateur tourne sous Linux. Votre gestionnaire de packages installera `gmp` sans difficulté.\n",
    "    - Votre ordinateur tourne sous MacOS. Vous devez installer un gestionnaire de packages. Personnellement j'utilise MacPorts. Une fois MacPorts installé, ouvrir un terminal et exécuter `sudo port install gmp`.\n",
    "    \n",
    "2. Installer `gmpy`. Ça c'est facile. Ouvrir un terminal et exécuter `pip3 install gmpy`.\n",
    "\n",
    "Une fois `gmpy` installé, le module `mpmath` s'en aperçoit et utilise la bibliothèque `gmp` pour les calculs en virgule flottante. Le gain en vitesse est énorme, en gros un facteur 100. Conclusion :\n",
    "\n",
    "Si vous n'installez pas `gmpy`, d'ici la fin de ce notebook vous pourrez calculer $\\pi$ avec 100000 décimales en 3 secondes et avec 1 million de décimales en 3 minutes. Si vous installez `gmpy` vous pourrez calculer en 3 secondes 1 million de décimales de $\\pi$ et en 40 secondes 10 millions de décimales. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Les valeurs antiques de $\\pi$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Les premières valeurs approchées historiquement attestées du nombre $\\pi$ se trouvent sur des tablettes babyloniennes. Différentes valeurs y sont présentes, en particulier $3+\\frac 1 8\\simeq 3.125$. \n",
    "\n",
    "Les Égyptiens, utilisaient une valeur approchée de $\\pi$ intéressante. Le papyrus Rhind, qui date du XVIème siècle avant notre ère, en fournit la preuve (problèmes 41, 50, etc.). Lors de calculs de surfaces ou de volumes, nos amis du Nil utilisaient la valeur $4\\left(\\frac 8 9\\right)^2\\simeq 3.16$.\n",
    "\n",
    "On trouve également dans la Bible un passage où le nombre $\\pi$ apparaît de façon implicite. Il s'agit du premier livre des Rois, 7:23, parlant de la construction du palais de Salomon. On y trouve la valeur approchée $\\pi\\simeq \\frac{15} 5 = 3$. Cela ne veut pas dire que les Hébreux avaient une approximation de $\\pi$ moins précise que les Égyptiens, mais tout simplement que la valeur qu'ils utilisaient n'était pas d'une importance capitale dans le passage en question de l'ancien testament !\n",
    "\n",
    "Si l'on peut se prêter à sourire de ces approximations, il faut toutefois souligner que les besoins de précision des Égyptiens pour la construction de greniers à blé n'étaient peut-être pas aussi importants que les nôtres. Une approximation à $2\\times 10^{-2}$ près était sans aucun doute suffisante pour prévoir les capacités de stockage des récoltes. Par ailleurs, il n'est, ni chez les Hébreux ni chez les Égyptiens, fait mention du nombre $\\pi$ en tant que tel. Nous trouvons dans les textes antiques des \"recettes\" pour calculer l'aire d'un disque, sans que le coefficient multiplicateur y soit clairement identifié comme un nombre remarquable.\n",
    "\n",
    "D'excellents livres sur l'histoire de $\\pi$ existent, et ce n'est pas le but de ce notebook. Nous allons dans ce qui suit étudier quelques algorithmes liés au calcul de valeurs approchées \"modernes\" de $\\pi$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pi_babylone = 3 + 1 / 8\n",
    "pi_babylone"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pi_egypte = 4 * (8 / 9) ** 2\n",
    "pi_egypte"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Passons maintenant à un algorithme qui permet de calculer quelques centaines, voire milliers de décimales de $\\pi$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. La formule de Machin"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Euh, c'est quoi ce truc ? John Machin est un mathématicien anglais qui a calculé en 1706 les 100 premières décimales du nombre $\\pi$. Nous allons voir comment, et calculer 1000 décimales de $\\pi$ avec sa formule."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.1 Approximons la fonction arc tangente"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Je m'excuse à l'avance des quelques préliminaires mathématiques qui vont suivre. Si vous ne voyez pas le rapport entre $\\arctan$ et $\\pi$, soyez patients."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Soit $n$ un entier non nul. Soit $x\\in]0, 1]$.On a \n",
    "\n",
    "$$\\int_0^x\\sum_{k=0}^{n-1}(-1)^kt^{2k}dt=\\sum_{k=0}^{n-1}(-1)^{k}\\frac{x^{2k+1}}{2k+1}$$\n",
    "\n",
    "Par ailleurs, pour tout $t\\in [0, x]$, \n",
    "\n",
    "$$\\sum_{k=0}^{n-1}(-1)^kt^{2k}=\\frac{1-(-1)^{n}t^{2n}}{1+t^2}$$\n",
    "\n",
    "Ainsi,\n",
    "\n",
    "$$\\sum_{k=0}^{n-1}(-1)^{k}\\frac{x^{2k+1}}{2k+1}=\\int_0^x\\frac{dt}{1+t^2}-(-1)^{n}\\int_0^x \\frac{t^{2n}dt}{1+t^2}=\\arctan x-(-1)^{n}\\int_0^x \\frac{t^{2n}dt}{1+t^2}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "L'intégrale du membre de droite est facile à majorer, en remarquant que le dénominateur de la fraction est au moins égal à 1 :\n",
    "\n",
    "$$\\left|\\int_0^x \\frac{t^{2n}dt}{1+t^2}\\right|\\le \\int_0^x t^{2n}dt=\\frac{x^{2n+1}}{2n+1}$$\n",
    "\n",
    "On en déduit la \n",
    "\n",
    "__Proposition__ : \n",
    "\n",
    "$$\\sum_{k=0}^{n-1}(-1)^{k}\\frac{x^{2k+1}}{2k+1}=\\arctan x + \\varepsilon_n$$\n",
    "\n",
    "où\n",
    "\n",
    "$$|\\varepsilon_n|\\le \\frac{x^{2n+1}}{2n+1}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.2 Qu'obtient-on pour $x=1$ ?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "En prenant $x=1$ dans la proposition ci-dessus, on obtient \n",
    "\n",
    "$$\\sum_{k=0}^{n-1}\\frac{(-1)^{k}}{2k+1}=\\frac\\pi 4 + \\varepsilon_n$$\n",
    "\n",
    "où\n",
    "\n",
    "$$|\\varepsilon_n|\\le \\frac{1}{2n+1}$$\n",
    "\n",
    "Pour obtenir 10 décimales de $\\pi$ il faudrait (suffirait, pour être exact) donc additionner $n$ termes, où $n$ vérifie $\\frac 4 {2n+1}\\le 10^{-10}$, c'est à dire environs 20 milliards de termes. Contentons-nous d'une approximation à $10^{-5}$ près en ajoutant 200000 termes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def approx_mauvaise(n):\n",
    "    s = 0\n",
    "    for k in range(n):\n",
    "        s = s + (-1) ** k / (2 * k + 1)\n",
    "    return 4 * s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "approx_mauvaise(200000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "approx_mauvaise(200001)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.3 Passons à des choses plus sérieuses"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Dorénavant, prenons $0<x<1$. \n",
    "\n",
    "Soit $p\\in\\mathbb N$. Combien suffit-il d'additionner de termes dans la formule vue plus haut pour obtenir une approximation de $\\arctan x$ à $10^{-p}$ près ? Il nous suffit évidemment de choisir $n$ tel que $\\frac{x^{2n+1}}{2n+1}\\le 10^{-p}$, ce qui s'écrit encore \n",
    "\n",
    "$$(2n+1)\\log \\frac 1 x + \\log(2n+1)\\ge p$$\n",
    "\n",
    "où $\\log$ désigne le logarithme en base 10.\n",
    "\n",
    "Cette inégalité ne peut pas être résolue de façon exacte. Une majoration plus crue nous dit qu'il suffit de choisir $n$ tel que \n",
    "\n",
    "$$(2n+1)\\log \\frac 1 x\\ge p$$\n",
    "\n",
    "ou encore \n",
    "\n",
    "$$n \\ge \\frac 1 2\\left(\\frac{p}{\\log \\frac 1 x}-1\\right)$$\n",
    "\n",
    "Par exemple, pour $p=1000$ et $x=\\frac 1 5$, cela donne"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(1000 / math.log(5, 10) - 1) / 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Le facteur $2n+1$ est-il si important que cela ? Écrivons une fonction qui renvoie le plus petit entier $n$ tel que $(2n+1)\\log \\frac 1 x + \\log(2n+1)\\ge p$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def erreur(x, p):\n",
    "    n = 1\n",
    "    while (2 * n + 1) * math.log(1 / x, 10) + math.log(2 * n + 1, 10) <= p + 1: n = n + 1\n",
    "    return n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "erreur(1 / 5, 1000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Moralité : le facteur $2n+1$ n'est pas important du tout :-).\n",
    "\n",
    "Bon, au travail ! La fonction `approx_atan` prend en paramètres un flottant $x$ et un entier $p$. Elle renvoie une approximation de $\\arctan x$ à $10^{-p}$ près."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def approx_atan(x, p):\n",
    "    x = N(x, p + 10)\n",
    "    n = erreur(x, p)\n",
    "    s = 0\n",
    "    for k in range(n):\n",
    "        s = s + (-1) ** k * x ** (2 * k + 1) / (2 * k + 1)\n",
    "    return s"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Testons $\\arctan \\frac 1 5$ avec 1000 chiffres après la virgule."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x15 = S(1) / 5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "approx_atan(x15, 1000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Et vérifions avec `sympy`!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "atan(x15).evalf(1000) - approx_atan(x15, 1000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Tout va bien. Passons à la formule de Machin."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Proposition (Machin)__ : $$\\frac\\pi 4 = 4 \\arctan\\frac 1 5 - \\arctan \\frac 1 {239}$$\n",
    "\n",
    "__Démonstration__ :  Voici tout d'abord un exercice laissé au lecteur : pour $a,b\\in \\mathbb R_+$ tels que $ab<1$, on a \n",
    "\n",
    "$$\\arctan a + \\arctan b = \\arctan \\frac{a+b}{1-ab}$$\n",
    "\n",
    "On en déduit que $2\\arctan \\frac 1 5 = \\arctan \\frac{\\frac 2 5}{1-\\frac 2 5 \\frac 2 5}=\\arctan \\frac 5 {12}$.\n",
    "\n",
    "De là, $4\\arctan \\frac 1 5 = 2 \\arctan \\frac 5 {12} = \\arctan\\frac{\\frac{10}{12}}{1-\\frac{25}{144}}$.\n",
    "\n",
    "Au secours, Python !"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Rational(10, 12) / (1 - Rational(25, 144))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Bon, super. $4\\arctan \\frac 1 5=\\arctan \\frac{120}{119}$. Encore un petit effort !\n",
    "\n",
    "$\\arctan\\frac {120}{119}-\\arctan \\frac 1 {239}=\\arctan\\frac{\\frac {120}{119}-\\frac 1 {239}}{1+\\frac {120}{119}\\frac 1 {239}}$.\n",
    "\n",
    "Encore une fois, Python est notre ami."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = Rational(120, 119)\n",
    "y = Rational(1, 239)\n",
    "(x - y) / (1  + x * y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Bref, $4 \\arctan \\frac 1 5 - \\arctan \\frac 1 {239} = \\arctan 1 =\\frac \\pi 4$.\n",
    "\n",
    "Cette formule nous suggère une méthode pour calculer une approximation de $\\pi$. Combien doit-on additionner de termes dans notre formule fétiche ? Majorons grossièrement (pas trop quand même) l'erreur $\\varepsilon_n$ commise en additionnant $n$ termes pour approcher $\\arctan \\frac 1 5$ et $\\arctan \\frac 1 {239}$, puis en multipliant par 4 :"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\\varepsilon_n\\le \\frac 4 {2n+1}\\left(\\frac 4 {5^{2n+1}} + \\frac 1 {239^{2n+1}}\\right)\\le\\frac 4 {2n+1}\\frac 1 {5^{2n}}$$ \n",
    "\n",
    "car $239 > 5$ et $4+1=5$. Pour obtenir une approximation de $\\pi$ à $10^{-p}$ près, il suffit donc de choisir $n$ tel que\n",
    "\n",
    "$$\\frac 4 {2n+1}\\frac 1 {5^{2n}}\\le 10^{-p}$$\n",
    "\n",
    "ou encore\n",
    "\n",
    "$$5^{2n}\\ge \\frac{4 \\times 10^p}{2n+1}$$\n",
    "\n",
    "Comme nous avons compris que le facteur $2n+1$ ne change pas la face du monde, il suffit donc de choisir $n$ tel que\n",
    "\n",
    "$$5^{2n}\\ge 4 \\times 10^p$$\n",
    "\n",
    "ou encore\n",
    "\n",
    "$$n\\ge \\frac {p+ \\log 4}{2\\log 5}$$ "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Voulons-nous 1000 décimales de $\\pi$ ?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(1000 + math.log(4, 10)) / (2 * math.log(5, 10))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Voici le Machin."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def machin(p):\n",
    "    x15 = N(Rational(1, 5), p + 10)\n",
    "    x1239 = N(Rational(1, 239), p + 10)\n",
    "    return 4 * (4 * approx_atan(x15, p) - approx_atan(x1239, p))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "machin(1000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Vérifions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pi.evalf(1000) - machin(1000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On peut lancer `machin(2000)` ou `machin(3000)` mais n'espérons pas obtenir un million de chiffres de $\\pi$ par cette méthode !"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Remarque__ : John Machin n'avait pas d'ordinateur. Il a donc calculé 100 décimales de $\\pi$ __à la main__."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Le record actuel du nombre de décimales calculées de $\\pi$ est 22 trillions ... Comment atteindre une telle précision en un temps raisonnable ? Nous allons examiner deux algorithmes permettant de calculer extrêmement rapidement des valeurs approchées de $\\pi$.\n",
    "\n",
    "L'intérêt d'avoir deux algorithmes totalement distincts est de pouvoir __COMPARER__ les résultats renvoyés, et donc de __CONFIRMER__ la valeur obtenue. Nous allons commencer par l'algorithme de Salamin et Brent."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. L'algorithme de Salamin-Brent"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On va construire dans ce qui suit deux suites récurrentes. Une certaine quantité définie à partir de ces suites converge très rapidement vers $\\pi$. Nous allons voir précisément ce qu'il faut entendre par \"rapidement\"."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.1 Suites récurrentes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On définit les suites $(a_n)$ et $(b_n)$ par récurrence : "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Tout d'abord, $a_0=1$, $b_0=\\frac 1 {\\sqrt 2}$. \n",
    "\n",
    "- Pour tout $n\\ge 0$, $a_{n+1}=\\frac {a_n+b_n} 2$, $b_{n+1}=\\sqrt{a_nb_n}$.\n",
    "\n",
    "Les suites $(a_n)$ et $(b_n)$ sont donc formées en calculant des moyennes arithmétiques et des moyennes géométriques. Elles convergent _très rapidement_ vers ce que l'on appelle la moyenne __arithmético-géométrique__ de $1$ et $\\frac 1 {\\sqrt 2}$, que nous noterons ci-dessous $AGM(1, \\frac 1 {\\sqrt 2})$.\n",
    "\n",
    "Définissons ensuite $c_{n}=a_n^2-b_n^2$ puis\n",
    "\n",
    "$$\\pi_n=\\frac{2a_n^2}{1-2\\sum_{k=0}^n2^kc_k^2}$$\n",
    "\n",
    "__Proposition__ :\n",
    "\n",
    "$$0\\le \\pi - \\pi_n\\le\\frac{\\pi^2 2^{n + 4}e^{-\\pi2^{n+1}}}{AGM(1,\\frac{\\sqrt 2}{2})^2}$$\n",
    "\n",
    "__Démonstration__ : J'admets cette proposition. La démonstration en est abordable mais un peu longue."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.2 L'erreur commise"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Nous tenons clairement là un algorithme intéressant pour calculer des approximations de $\\pi$. Le majorant de l'erreur tend très rapidement vers 0, ce qui nous permet d'espérer de très bonnes approximations en très peu d'itérations. Combien d'itérations, au juste ? Regardons cela."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Voici tout d'abord une fonction `agm` qui calcule la moyenne arithmético-géométrique de deux réels $a$ et $b$ à $10^{-10}$ près. Inutile d'être précis, nous voulons juste majorer $\\pi-\\pi_n$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def agm(a, b):\n",
    "    u = a\n",
    "    v = b\n",
    "    while abs(v - u) > 1e-10:\n",
    "        u, v = (u + v) / 2, math.sqrt(u * v)\n",
    "    return u"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "agm(1, 1 / math.sqrt(2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Combien de termes sont nécessaires pour approcher $\\pi$ à $10^{-p}$ près ?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def nb_termes(p):\n",
    "    n = 0\n",
    "    z = math.log(math.pi ** 2 / (agm(1, math.sqrt(2) / 2) - 1e-10), 10)\n",
    "    while math.pi * 2 ** (n + 1) * math.log(math.e, 10) - (n + 4) * math.log(2, 10) < p + math.log(z, 10):\n",
    "        n = n + 1\n",
    "    return n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "100000 chiffres de $\\pi$ ?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nb_termes(10 ** 5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "16 itérations suffiront. Et 1 million de chiffres ?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nb_termes(10 ** 6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Il suffit de faire __TROIS__ itérations de plus. Et 22 trillions de décimales ?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nb_termes(22 * 10 ** 12)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Eh oui, 22 trillions de décimales en 43 itérations ..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.3 La boucle magique"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La fonction `salamin` est évidente. Elle calcule une approximation de $\\pi$ à $10^{-p}$ près. La fonction prend également un paramètre opitonnel $n$. Si $n$  vaut `None`, la fonction appelle `nb_termes` pour décider du nombre d'itérations à effectuer. Sinon, elle effectue $n$ itérations.\n",
    "\n",
    "Le paramètre optionnel `dbg` affiche en particulier le numéro de l'itération en cours. Cela permet de s'assurer que les calculs avancent, et fait patienter celui qui est assis devant son ordinateur.\n",
    "\n",
    "__Remarquez__ le $p+10$ dans les approximations numériques. On prévoit 10 chiffres significatifs supplémentaires dans les calculs en espérant que cela suffira pour compenser les erreurs d'arrondi."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def salamin(p, n=None, dbg=False):\n",
    "    if n == None:\n",
    "        n = nb_termes(p) + 1\n",
    "    if dbg: print('Nombre itérations :', n)\n",
    "    a = N(1, p + 10)\n",
    "    b = 1 / N(sqrt(2), p + 10)\n",
    "    s = N(1, p + 10) / 2\n",
    "    for k in range(1, n + 1):\n",
    "        if dbg: print(k)\n",
    "        a, b = (a + b) / 2, sqrt(a * b)\n",
    "        c = a * a - b * b\n",
    "        s = s - 2 ** k * c\n",
    "        #gc.collect()\n",
    "    return N(2 * a ** 2 / s, p + 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Voici un petit graphique montrant le nombre de décimales obtenues en fonction du nombre d'itérations. Les calculs sont faits avec 10000 chiffres après la virgule. En gros le nombre de décimales exactes __double__ à chaque itération : l'algorithme de Salamin-Brent a une vitesse de convergence __quadratique__."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams['figure.figsize'] = (12, 6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "s = []\n",
    "p = 10000\n",
    "Npi = N(pi, p + 10)\n",
    "ks = list(range(1, 15))\n",
    "for k in ks:\n",
    "    s.append(salamin(p, k) - Npi)\n",
    "s = [-log(abs(x), 10) for x in s]\n",
    "plt.plot(ks, s, 'k')\n",
    "plt.plot(ks, s, 'or')\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Nous y sommes : calculons $\\pi$ avec $10^5$ chiffres après la virgule. Si vous avez installé `gmpy` vous pouvez mettre sans risque $10^6$ ."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p = 10 ** 5\n",
    "x = salamin(p, None, True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Écrivons la valeur obtenue dans un fichier !"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f = open('pi_salamin_' + str(p) + '.txt', 'w')\n",
    "f.write(str(x))\n",
    "f.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Bonne lecture :-)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sur ma machine, j'obtiens $10^6$ décimales en environ trois secondes. Pour $10^7$ décimales, une quarantaine de secondes suffisent. Et pour $10^8$ (cent millions de décimales), une quinzaine de minutes s'imposent. Le fichier obtenu fait alors 100 Méga-Octets et mon éditeur de textes ne peut pas l'ouvrir. En plus ma machine utilise pour le calcul toute la mémoire disponible (4 Giga-Octets). Inutile donc de tenter le milliard de décimales :-)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.4 On a fini, alors ?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Euh, pas tout à fait. Avez-vous lu le fichier `salamin_100000.txt` ? Avez-vous bien vérifié chaque décimale de $\\pi$ ? N'y a-t-il pas d'erreur ? Nous pourrions évidemment demander à `sympy` une évaluation et comparer avec notre résultat. Mais que faire lorsqu'on tente de battre un record est qu'on est le premier à calculer $N$ décimales de $\\pi$ ? La solution qui s'impose est de refaire le calcul avec un autre algorithme puis de comparer les deux résultats. C'est ce que nous allons faire."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. L'algorithme de Borwein"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "L'algorithme de Borwein et Borwein (Jonathan et Peter) est un algorithme __quartique__, c'est à dire que le nombre de décimales calculées __quadruple__ à chaque itération. Il fonctionne comme suit :\n",
    "\n",
    "On pose $a_0=6-4\\sqrt 2$ et $y_0=\\sqrt 2-1$.\n",
    "\n",
    "Puis, pour tout entier $k\\ge 0$, on pose\n",
    "\n",
    "$$y_{k+1}=\\frac{1-(1-y_k^4)^{1/4}}{1+(1-y_k^4)^{1/4}}$$\n",
    "\n",
    "et\n",
    "\n",
    "$$a_{k+1}=a_k(1+y_{k+1})^4-2^{2k+3}y_{k+1}(1+y_{k+1}+y_{k+1}^2)$$\n",
    "\n",
    "Alors $\\frac 1 {a_k}$ converge quartiquement vers $\\pi$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La fonction `borwein` ci-dessous calcule $p$ décimales de $\\pi$. N'ayant pas réussi à trouver dans la littérature un majorant explicite de l'erreur, je lui fais effectuer une sortie de boucle lorsque deux valeurs successives $a_k$ et $a_{k+1}$ sont égales à $10^{-p}$ près. Remarquez, encore une fois, les 10 décimales de \"sécurité\" pour éviter les erreurs d'arrondi."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def borwein(p, dbg=False):\n",
    "    a = N(6 - 4 * sqrt(2), p + 10)\n",
    "    y = N(sqrt(2) - 1, p + 10)\n",
    "    k = 0\n",
    "    while True:\n",
    "        if dbg: print(k)\n",
    "        z = sqrt(sqrt(1 - y ** 4))\n",
    "        y = (1 - z) / (1 + z)\n",
    "        a1 = a * (1 + y) ** 4 - 2 ** (2 * k + 3) * y * (1 + y + y * y)\n",
    "        if N(a1, p) == N(a, p): break\n",
    "        a = a1\n",
    "        k = k + 1\n",
    "        #gc.collect()\n",
    "    return N(1 / a, p + 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "p = 10 ** 5\n",
    "x = borwein(p, True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f = open('pi_borwein_' + str(p) + '.txt', 'w')\n",
    "f.write(str(x))\n",
    "f.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Bon, et maintenant ? Nous disposons de deux fichiers, `pi_salamin_100000.txt` et `pi_borwein_100000.txt`. Ces deux fichiers devraient être identiques.\n",
    "\n",
    "Si votre ordinateur tourne sous Linux ou MacOS, exécutez la cellule ci-dessous. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!diff -s pi_borwein_100000.txt pi_salamin_100000.txt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Cela nous remplit d'une joie sans bornes :-)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Si vous êtes sous Windows, essayez la commande `FC` (je ne l'ai pas testée)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5. Peut-on faire mieux ?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Il existe des algorithmes, toujours basés sur des récurrences, dont la convergence est cubique, quartique (algorithme de Borwein), ..., nonique. Nonique, cela veut dire qu'à chaque itération le nombre de décimales exactes est multiplié par 9. Je ne les présenterai pas ici, mais une petite recherche sur Internet vous en donnera une idée. \n",
    "\n",
    "Apparemment, la formule tendance pour obtenir les derniers records est la formule de Chudnovsky :\n",
    "\n",
    "$$\\frac 1 \\pi=12\\sum_{k=0}^\\infty \\frac{(-1)^k(6k)!(13591409+545140134k)}{(3k)!(k!)^3640320^{3k+3/2}}$$\n",
    "\n",
    "Chaque terme de la formule permet d'obtenir 14 chiffres corrects supplémentaires. De plus, utilisée de façon extrêmement rusée, cette formule permet d'utiliser un résultat obtenu à une certaine précision pour obtenir un résultat à une meilleure précision."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Depuis quelques années, les records dans le calcul des décimales de $\\pi$ sont obtenus par des particuliers sur des PC survitaminés. Le record actuel est, à ma connaissance, de 22.459.157.718.361 décimales (22 trillions). Le calcul a demandé 105 jours. Pourquoi ce nombre précis de décimales ? Regardez ci-dessous :-)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "math.pi ** math.e"
   ]
  },
  {
   "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
}