{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# La millionième décimale de $\\pi$ (II)\n", "\n", "Marc Lorenzi\n", "\n", "16 février 2019" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from sympy import *\n", "init_printing()\n", "x = Symbol('x')\n", "import math" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Nous allons dans ce notebook nous intéresser à la question suivante : \n", "\n", "Est-il possible de calculer la millionième décimale de $\\pi$ \n", "\n", "- __SANS calculer les décimales précédentes__, et\n", "- __En effectuant uniquement des calculs en précision standard__ ?\n", "\n", "En 1997, David Bailey, Peter Borwein et Simon Plouffe ont démontré que ceci était possible, _en base 16_. L'aricle original se trouve facilement sur Internet, et il est tout à fait lisible. Pour les lecteurs motivés, en voici la référence :\n", "\n", "D. Bailey, P. Borwein, S. Plouffe, _On the Rapid Computation of Various Polylogarithmic Constants_, Mathematics of Computation, Vol 66, Number 218, April 1997, Pages 903-913.\n", "\n", "\n", "On peut donc déterminer le $d$ième chiffre de $\\pi$ en base 16, ainsi que quelques chiffres suivants, en un temps quasi-linéaire et en calculant en précision standard. Nous allons tout d'abord établir la formule qui permet ce calcul, puis faire le dit calcul ..." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 1. Une intégrale" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 1.1 Intégrons une fraction rationnelle" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Considérons la fraction rationnelle suivante." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "F = (-8 * x ** 5 - 4 * sqrt(2) * x ** 4 - 8 * x ** 3 + 4 * sqrt(2)) / (1 - x ** 8)\n", "F" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Le dénominateur de $F$ est immédiat à factoriser." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "factor(1 - x ** 8)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "La décomposition de $F$ en éléments simples est" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "apart(F, x)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "On en déduit facilement une primitive de $F$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "G = integrate(F, x)\n", "G" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Calculons $\\int_0^{\\sqrt(2)/2} F(x)dx$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "G.subs(x, sqrt(2) / 2) - G.subs(x, 0)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "init_printing(pretty_print=False)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 1.2 Échange d'une intégrale et d'une somme" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Faisons maintenant un petit calcul ... On a pour tout $x\\in[0,\\frac{\\sqrt 2} 2]$\n", "\n", "$$\\frac{-8x^5-4\\sqrt 2 x^4-8x^3+4\\sqrt 2}{1-x^8}=(-8x^5-4\\sqrt 2 x^4-8x^3+4\\sqrt 2)\\sum_{i=0}^\\infty x^{8i}$$\n", "\n", "La série convergeant uniformément sur $[0,\\frac{\\sqrt 2}2]$, on peut échanger la série et l'intégrale :\n", "\n", "$$\\int_0^{\\sqrt 2 / 2}F(x)dx=\\sum_{i=0}^\\infty\\left( -8\\int_0^{\\sqrt 2 / 2}x^{8i+5}dx-4\\sqrt 2 \\int_0^{\\sqrt 2 / 2}x^{8i+4}dx-8\\int_0^{\\sqrt 2 / 2}x^{8i+3}dx+4\\sqrt 2\\int_0^{\\sqrt 2 / 2}x^{8i}dx\\right)$$\n", "\n", "Laissons le reste au lecteur. Le membre de gauche vaut $\\pi$ et les intégrales du membre de droite sont faciles à calculer. On obtient\n", "\n", "$$\\sum_{i=0}^\\infty \\frac 1 {16^i}\\left(\\frac{4}{8i+1}-\\frac{2}{8i+4}-\\frac{1}{8i+5}-\\frac{1}{8i+6}\\right)=\\pi$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "C'est cette dernière formule qui va nous permettre de calculer des décimales de $\\pi$. Remarquez les puissances de 16. Si nous voulons calculer le $(d+1)$-ième chiffre hexadécimal de $\\pi$, \n", "\n", "- Nous multiplions $\\pi$ par $16^d$.\n", "- Nous ne gardons que la partie fractionnaire du résultat.\n", "- Nous prenons le premier chiffre hexadécimal de ce que nous venons de trouver." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 2. L'algorithme de Bailey, Borwein et Plouffe (BBP)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 2.1 Partie fractionnaire" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Notons, pour tout réel $x\\ge 0$, $x\\bmod 1=x-\\lfloor x\\rfloor$ la partie fractionnaire de $x$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def frac(x):\n", " return x - math.floor(x)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Pour $q=1,4,5,6$, posons \n", "\n", "$$S_q = \\sum_{i=0}^\\infty \\frac 1 {16^i}\\frac{1}{8i+q}$$\n", "\n", "Soit $d\\ge 0$. On a\n", "\n", "$$(16^d S_q)\\bmod 1=\\sum_{i=0}^d \\frac{16^{d-i}}{8i+q} \\bmod 1+\\sum_{i=d+1}^\\infty \\frac{16^{d-i}}{8i+q}\\bmod 1$$\n", "\n", "ou encore\n", "\n", "$$(16^d S_q)\\bmod 1=\\sum_{i=0}^d \\frac{16^{d-i}\\bmod (8i + q)}{8i+q} \\bmod 1+\\sum_{i=d+1}^\\infty \\frac{16^{d-i}}{8i+q}\\bmod 1$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Notons dans ce qui suit\n", "\n", "$$S_{dq}=\\sum_{i=0}^d \\frac{16^{d-i}\\bmod (8i + q)}{8i+q} \\bmod 1$$\n", "\n", "et\n", "\n", "$$S'_{dq}=\\sum_{i=d+1}^\\infty \\frac{16^{d-i}}{8i+q}\\bmod 1$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 2.2 Le calcul de $S_{dq}$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Ce calcul demande l'évaluation de $16^{d-i}$ modulo $8i+q$. Pour cela, on utilise une exponentiation rapide. La fonction `puissance` prend en paramètres trois entiers $x,n,p$. Elle renvoie $x^n\\bmod p$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def puissance(x, n, p):\n", " z = 1\n", " m = n\n", " y = x\n", " while m != 0:\n", " if m % 2 == 1:\n", " 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": [ "puissance(2, 10000000000, 10000)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "La somme $S_{dq}$ est alors facilement calculable." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def somme1(d, q):\n", " s = 0\n", " for i in range(d + 1):\n", " s = frac(s + puissance(16, d - i, 8 * i + q) / (8 * i + q))\n", " return s" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "somme1(10 ** 5, 1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 2.3 Le calcul de $S'_{dq}$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Pour le calcul de $S'_{dq}$ on se contente d'additionner un certain nombre de ses termes, $N$ disons, en virgule flottante. Bref, on approche $S'_{dq}$ par \n", "\n", "$$\\sum_{i=d+1}^{d+N} \\frac{16^{d-i}}{8i+q}$$\n", "\n", "Combien prendre pour $N$ ? Eh bien l'erreur commise est\n", "\n", "$$\\sum_{i=d+N+1}^\\infty \\frac{16^{d-i}}{8i+q}\\le \\sum_{i=d+N+1}^\\infty 16^{d-i}=\\frac 1 {16^{N+1}}\\sum_{k=0}^\\infty\\frac 1 {16^k}=\\frac{1}{15\\times 16^{N}}$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Si l'on veut une erreur inférieure à $10^{-16}$, valeur raisonnable pour des calculs sur des flottants, il suffit donc de choisir $N$ tel que $15\\times 16^{N}\\ge 10^{16}$. Passons au logarithme décimal :\n", "\n", "$$N\\log 16+\\log 15\\ge 16$$\n", "\n", "c'est à dire\n", "\n", "$$N\\ge \\frac{16-\\log 15}{\\log 16}$$" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "(16 - math.log(15, 10)) / math.log(16, 10)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Soyons généreux et prenons 15 termes pour le calcul de la somme en question." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def somme2(d, q):\n", " s = 0\n", " for i in range(d + 1, d + 16):\n", " s = frac(s + 16 ** (d - i) / (8 * i + q))\n", " return s" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "somme2(10 ** 5, 1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 2.4 BBP, enfin" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Tout a été fait, il n'y a plus qu'à additionner. La fonction `somme`, tout d'abord, renvoie $S_{dq}+S'_{dq}$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def somme(d, q):\n", " return somme1(d, q) + somme2(d, q)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "somme(10 ** 5, 1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "La fonction BBP renvoie une valeur approchée de $(16^d\\pi)\\bmod 1$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def BBP(d):\n", " s = frac(4 * somme(d, 1) - 2 * somme(d, 4) - somme(d, 5) - somme(d, 6))\n", " return s" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "BBP(10 ** 5 - 1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 3. Fiabilité des résultats" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 3.1 Comptons les opérations" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "La valeur renvoyée par BBP est-elle correcte ? Après tout nous faisons des calculs en virgule flottante ! Python utilise pour ces calculs le format IEE754, qui garantit une erreur d'au plus 1 bit sur le dernier chiffre pour chaque opération. Combien avons nous effectué d'opérations en virgule flottante ?\n", "\n", "Nous effectuons 4 sommes de $d$ termes. Pour chaque terme on fait une division. Cela nous donne donc $8d-4$ opérations : $4d$ divisions et $4(d-1)$ additions.\n", "\n", "Nous effectuons également 4 sommes de 15 termes : pour chaque terme il y a un calacul de puissance en flottants et une division. Au total $8\\times 15 + 56=176$ opérations en virgule flottante.\n", "\n", "Au total, nous faisons donc $8d+176$ opérations en virgule flottante. Rajoutons quelques additions (il y a 8 sommes à ajouter), une multiplication par 4 et une multiplication par 2. Nous arrivons $8d+185$ opérations (en ai-je oublié ?).\n", "\n", "Imaginons qu'à chaque opération l'erreur de 1 bit s'ajoute aux précédentes. Ceci est très pessimiste ! Nous aurons à la fin de nos calcul un nombre de bits erronés de l'ordre de $\\log_2(8d+185)$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def bits_erreur(d):\n", " return math.log(8 * d + 185, 2)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "bits_erreur(10 ** 6)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Python calcule sur 52 bits. Combien de bits exacts de façon sûre ?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def bits_exacts(d):\n", " return 52 - bits_erreur(d)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "bits_exacts(10 ** 6)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Sachant qu'un chiffre en hexadécimal est constitué de 4 bits, combien de chiffres hexadécimaux exacts ?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def hexa_exacts(d):\n", " return math.floor(bits_exacts(d)) // 4" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "hexa_exacts(10 ** 6)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "En lançant `BBP(10 ** 6 - 1)` nous aurons donc de façon sûre les chiffres 1000000 à 1000006 du nombre $\\pi$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 3.2 Convertissons en hexadécimal" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "La fonction `vers_hexa` prend un réel $x\\in[0,1[$ en paramètre et renvoie une chaîne de caractères qui contient les dix premiers chiffres du développement hexadécimal de $x$. Lisez la fonction, elle ne pose aucune difficulté." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def vers_hexa(x):\n", " s = ''\n", " for k in range(10):\n", " y = 16 * x\n", " c = math.floor(y)\n", " s = s + chiffre_hexa(c)\n", " x = frac(y)\n", " return s" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "La fonction `vers_hexa` utilise `chiffre_hexa` qui prend en paramètre un entier $c$ entre 0 et 15 et renvoie le \"chiffre\" hexadécimal correspondant." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def chiffre_hexa(c):\n", " if c <= 9: return str(c)\n", " elif c == 10: return 'A'\n", " elif c == 11: return 'B'\n", " elif c == 12: return 'C'\n", " elif c == 13: return 'D'\n", " elif c == 14: return 'E'\n", " elif c == 15: return 'F'" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 3.3 Nous y sommes" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### 3.3.1 Cent mille" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Quelle est le cent-millième chiffre hexadécimal de $\\pi$ ?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "vers_hexa(BBP(10 ** 5 - 1))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "C'est un 5. Et les chiffres suivants ? Sont-ils corrects ?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "hexa_exacts(10 ** 5)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "On peut donc garantir les 8 chiffres `535EA16C`." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### 3.3.1 Un million" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Quelle est le millionième chiffre hexadécimal de $\\pi$ ? Soyez patients, environ une minute de calcul sur ma machine. Python n'est pas fort en boucles `for`, et on ne fait que des boucles `for` :-)." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "vers_hexa(BBP(10 ** 6 - 1))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "C'est un 2. Et les chiffres suivants ? Combien sont corrects ?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "hexa_exacts(10 ** 6)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "On peut donc garantir les 7 chiffres `26C65E5`. Au cas où vous vous poseriez la question, il s'avère que les chiffres `2CB` sont également corrects (voir l'article BBP). Mais on ne peut pas les __garantir__ avec cet algorithme." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### 3.3.1 Dix millions" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Quelle est le dix-millionième chiffre hexadécimal de $\\pi$ ? Soyez infiniment patients. Allez prendre un petit café, cela va prendre 10 fois plus de temps que pour un million (donc, dix minutes environ)." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "vers_hexa(BBP(10 ** 7 - 1))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "C'est un 1. Et les chiffres suivants ? Sont-ils corrects ?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "hexa_exacts(10 ** 7)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "On peut donc garantir les 6 chiffres `17AF58`. Il s'avère que les chiffres `63` sont également corrects (voir l'article BBP). Les deux suivants sont `EF` alors que notre fonction renvoie les chiffres erronés `F0`." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 4. Et en base 10 ?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Bien évidemment, tout le monde meurt d'envie de savoir si une telle formule existe pour la base 10. Borwein, Borwein et Galway ont montré que ce genre de formule ne peut pas exister. Je resterai vague là-dessus, et me contenterai de citer leur article :\n", "\n", "\" _Of particular interest, we show that $\\pi$ has no Machin-type BBP arctangent formula when $b\\ne 2$. To the best of our knowledge, when there is no Machin-type BBP formula for a constant then no BBP formula of any form is known for that constant_ \".\n", "\n", "La référence de l'article est :\n", "\n", "J. Borwein, D. Borwein, W. Galway, _Finding and Excluding $b$-ary Machin-Type Individual Digit Formulae_, Canadian Journal of Mathematics, Vol. 56 (5), 2004, Pages 897-925.\n", "\n", "Il convient pour terminer de remarquer que des algorithmes permettant le calcul des décimales individuelles de $\\pi$ existent mais, à l'heure actuelle, ces algorithmes sont moins efficaces que ceux qui consistent à calculer toutes les décimales. " ] } ], "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 }