{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# L'exponentiation rapide\n", "\n", "Marc Lorenzi\n", "\n", "2 avril 2023" ] }, { "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'] = (8, 3)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Comment calculer le millionième nombre de Fibonacci en moins de 3 secondes ? Lisez et vous saurez." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 1. L'algorithme" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 1.1 La vision naïve du calcul des puissances" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Soit $\\mathcal M$ un monoïde de neutre $e$ dont la loi est notée multiplicativement. Pour tout $x\\in\\mathcal M$, on définit par récurrence sur $n$ la $n$ième puissance de $x$, $x^n$ par :\n", "\n", "- $x^0=e$\n", "- Pour tout $n\\in\\mathbb N^*$, $x^{n}=x^{n-1}x$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Pour éviter les complications inutiles, supposons pour le moment que $x$ est un nombre (entier, réel, complexe, peu importe). La définition mathématique des puissances nous donne sans aucun doute un algorithme pour les calculer. La fonction `puissance_naive` prend en paramètre un nombre $x$ et un entier naturel $n$. Elle renvoie $x^n$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def puissance_naive(x, n):\n", " if n == 0: return 1\n", " else: return puissance_naive(x, n - 1) * x" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "puissance_naive(2, 11)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Python n'aime pas un trop grand nombre d'appels récursifs imbriqués : essayez donc ci-dessus de calculer $2^{3000}$. Mais il est très facile de réécrire notre fonction avec une simple boucle. On en profite pour ajouter un compteur `c` qui enregistre le nombre de multiplications effectuées par la fonction." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def puissance_naive(x, n):\n", " p = 1\n", " c = 0\n", " for k in range(n):\n", " p = p * x\n", " c += 1\n", " return (p, c)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "puissance_naive(2, 3000)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Rien de surprenant, il faut faire 3000 multiplications pour calculer $2^{3000}$. Imaginons maintenant que l'on veuille calculer $2^{10^{10}}$. Inutile d'essayer, il faudrait faire dix milliards de multiplications, qui plus est avec des entiers énormes. Alors comment calculer des puissances avec d'énormes exposants ? Il nous faut une idée géniale." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 1.2 L'idée géniale" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Mettons que je veuille calculer $2^{16}$. Je peux remarquer que $2^{16}=((2^2)^2)^2)^2$. Chaque élévation au carré comptant pour une multiplication, cela fera au total 4 multiplications, et pas 16. On tient quelque chose.\n", "\n", "Maintenant, me direz vous, c'était facile parce que 16 est une puissance de 2. Et si je veux calculer $2^{13}$ ? Eh bien j'écris $2^{13}=2^82^42^1=(2^4)^22^42^1$. Il suffira de 5 multiplications. Sur l'exemple, on voit. Mais comment faire dans le cas général ?\n", "\n", "Prenons d'abord les choses \"à l'envers\", c'est à dire récursivement. Je veux écrire une fonction puissance qui sache calculer $x^n$. Si $n=0$ c'est évident. Sinon, $n=2p$ ou $n=2p+1$, où $p$ est un entier. Je calcule $y=(x^2)^p$. Comment ? Eh bien en calculant $x\\times x$ puis en appelant récursivement la fonction puissance avec le paramètre $p 0$, $C_{2p}=C_p+1$\n", "- Pour tout $p\\ge 0$, $C_{2p+1}=C_p+2$\n", "\n", "__Démonstration__ : Les deux premiers points sont laissés au lecteur. Traitons le cas de $C_{n}$, où $n=2p+1$, $p\\ge 0$. Comme $n>0$, le premier test échoue. La fonction calcule donc $x^2$, ce qui fait une multiplication. Puis elle s'appelle récursivement avec $p$ comme second paramètre. L'appel récursif effectue $C_p$ multiplications. Enfin, comme $n$ est impair, le second test échoue et le calcul de $xy$ effectue encore une multiplication. Au total : $C_p + 2$ multiplications, comme prévu. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Écrivons une fonction `complexite` qui calcule la complexité de notre exponentiation rapide. Oui, je sais, ça donne un peu le vertige, surtout si je commence à me demander quelle est la complexité de la fonction `complexite` :-)." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def complexite(n):\n", " if n == 0: return 0\n", " else:\n", " c = complexite(n // 2)\n", " if n % 2 == 0: return c + 1\n", " else: return c + 2" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print([(n, complexite(n)) for n in range(1, 20)])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "__Exercice obligé__ : Soit $A_n$ le nombre d'additions effectuées lors de l'appel à `complexite(n)`. Montrer que $A_n=C_n$. En d'autres termes, la complexité de `complexite` est la même que celle de la fonction qui a pour complexité `complexite` :-)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 1.4 Encadrements de la complexité" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Montrons tout d'abord une majoration de $C_n$, en accord avec ce que nous avons \"vu\" sur le graphe un peu plus haut." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Proposition.** Pour tout $n\\in\\mathbb N$, $C_n\\le 2\\lg(n+1)$, où $\\lg$ désigne le logarithme en base 2.\n", "\n", "**Démonstration.** On fait une récurrence forte sur $n$. Pour $n=0$ c'est évident puisque \n", "\n", "$$C_0=0\\le 2\\lg (0+1)=0$$\n", "\n", "Soit $n>0$. Supposons l'inégalité vraie pour tous les entiers strictement inférieurs à $n$.\n", "\n", "- Cas 1 : $n = 2p$, où $p>0$. On a \n", "\n", "$$C_n=C_p+1\\le 2\\lg(p+1)+1$$\n", "\n", "par l'hypothèse de récurrence. On a gagné si on montre que \n", "\n", "$$2\\lg(p+1)+1\\le\\lg(2p+1)$$\n", "\n", "ou encore \n", "\n", "$$\\lg\\frac{2p+1}{p+1}\\ge\\frac 1 2$$\n", "\n", "Soit $f:x\\mapsto \\frac{2x+1}{x+1}$, définie sur $\\mathbb R_+$. $f$ est dérivable, et pour tout $x\\ge 0$ on a \n", "\n", "$$f'(x)=\\frac 1{(x+1)^2}>0$$\n", "\n", "Notre fonction est donc strictement croissante sur $\\mathbb R_+$. Ainsi, pour tout entier $p\\ge 1$, $f(p)\\ge f(1)=\\frac 3 2$. Donc, \n", "\n", "$$2\\lg\\frac{2p+1}{p+1}\\ge 2\\lg\\frac 3 2\\ge 2\\lg\\sqrt 2= 1$$\n", "\n", "- Cas 2: $n=2p+1$ où $p\\ge 0$. On a \n", "\n", "$$C_n=C_p + 2\\le 2\\lg(p+1)+2=2\\lg(2(p+1))$$\n", "\n", "(eh oui, $\\lg 2 = 1$ !). Bref, \n", "\n", "$$C_n\\le 2\\lg(2p+2)=2\\lg(n+1)$$\n", "\n", "Ce cas était plus facile que l'autre." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Passons à la minoration." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Proposition.** Pour tout entier $n\\ge 1$, $C_n\\ge\\lg n + 2$.\n", "\n", "**Démonstration.** Encore une fois ce sera une récurrence forte. Mais à partir de $n=1$, parce que $\\lg 0$ ça fait mauvais effet. Pour $n=1$, on a $C_1=2=\\lg 1 +2$, donc tout va bien.\n", "\n", "Soit $n>1$. Supposons l'inégalité vraie pour tous les entiers strictement inférieurs à $n$ et strictement positifs.\n", "\n", "- Cas 1 : $n=2p$ où $p>0$. On a \n", "\n", "$$C_n=C_p+1\\ge\\lg p + 2 + 1=\\lg(2p)+2=\\lg n + 2$$\n", "\n", "- Cas 2 : $n = 2p+1$ où $p> 0$. On a \n", "\n", "$$C_n=C_p+2\\ge \\lg p + 2 + 2$$\n", "\n", "On a gagné si on montre que \n", "\n", "$$\\lg p + 2 \\ge \\lg (2p+1)$$\n", "\n", "ou encore \n", "\n", "$$\\lg\\frac{2p+1}{p}\\le 2$$\n", "\n", "On pose, pour $x>0$, \n", "\n", "$$f(x)=\\frac{2x+1}{x}$$\n", "\n", "$f$ est dérivable, de dérivée \n", "\n", "$$f'(x)=-\\frac 1 {x^2}<0$$\n", "\n", "Notre fonction est donc strictement décroissante sur $\\mathbb R_+^*$. Ainsi, pour tout entier $p\\ge 1$, $f(p)\\le f(1)=3$. Donc, \n", "\n", "$$\\lg\\frac{2p+1}{p}\\le \\lg 3\\le \\lg 4 = 2$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Bilan.** Pour tout $n\\ge 1$, \n", "\n", "$$2+\\lg n \\le C_n \\le \\lg(n+1)$$\n", "\n", "Le nombre de multiplications effectuées par la fonction d'exponentiation rapide est **logarithmique** en $n$. D'où l'appellation « rapide »." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 1.5 La valeur exacte de la complexité" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Il se trouve que l'on peut calculer $C_n$ de façon exacte (si l'on peut dire)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Proposition.** Pour tout $n\\in\\mathbb N$, $C_n=\\ell(n)+\\gamma(n)$, où $\\ell(n)$ est le nombre de chiffres de l'écriture binaire de $n$, et $\\gamma(n)$ est le nombre de 1 dans cette même écriture (On convient que $\\ell(0)=0$).\n", "\n", "**Démonstration.** En base 2, multiplier par 2, c'est juste ajouter un 0 ! On remarque donc que\n", "\n", "- $\\ell(0)=\\gamma(0)=0$.\n", "- Pour tout $p\\ge 1$, $\\ell(2p)=\\ell(p)+1$ et $\\gamma(2p)=\\gamma(p)$.\n", "- Pour tout $p\\ge 0$, $\\ell(2p+1)=\\ell(p)+1$ et $\\gamma(2p+1)=\\gamma(p)+1$.\n", "\n", "Posons, pour tout entier $n$, $C'_n=\\ell(n)+\\gamma(n)$. On a $C'_0=C_0=0$ et les suites $(C_n)$ et $(C'_n)$ vérifient les mêmes relations de récurrence. On en déduit (par récurrence) que $C_n=C'_n$.\n", "\n", "**Exercice.** Montrer que pour tout $n\\ge 1$, $\\ell(n)=\\lfloor\\lg n\\rfloor+1$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Remarque** :\n", "\n", "- Supposons que $n=2^k, k\\ge 0$ est une puissance de 2. Alors $\\ell(n)=k+1$ et $\\gamma(n)=1$, donc $C_n=k+2=\\lg n + 2$. Ainsi, notre minorant de $C_n$ est optimal, il est égal à $C_n$ pour des valeurs de $n$ aussi grandes que l'on veut.\n", "\n", "- Supposons maintenant que $n=2^k-1, k\\ge 0$. Alors $\\ell(n)=k$ et $\\gamma(n)=k$, donc $C_n=2k=\\lg (n + 1)$. Ainsi, notre majorant de $C_n$ est optimal, il est égal à $C_n$ pour des valeurs de $n$ aussi grandes que l'on veut.\n", "\n", "Nos fonctions en rouge étaient vraiment très bien choisies." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 1.6 Traçons $\\ell$ et $\\gamma$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "On commence par $\\ell$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def ell(n):\n", " if n == 0: return 0\n", " else: return ell(n // 2) + 1" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "[ell(n) for n in range(17)]" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "s = [ell(n) for n in range(129)]\n", "plt.plot(s, 'ok', ms=2)\n", "plt.grid()\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Passons à $\\gamma(n)$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def gamma(n):\n", " if n == 0: return 0\n", " else:\n", " g = gamma(n // 2)\n", " if n % 2 == 0: return g\n", " else: return g + 1" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "[gamma(n) for n in range(17)]" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "s = [gamma(n) for n in range(129)]\n", "plt.plot(s, 'sk', ms=2)\n", "plt.grid()\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 1.7 Une version itérative" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Je n'ai rien contre les fonctions récursives, mais comme je l'ai déjà dit plus haut, Python est un tantinet réticent à effectuer un très grand nombre d'appels récursifs imbriqués. Par défaut, le nombre maximal d'appels récursifs est 256 (ce nombre est cependant modifiable).\n", "\n", "Profitons du sectarisme de Python pour écrire une fonction non récursive qui effectue une exponentiation rapide. La voici." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def puissance_iter(x, n):\n", " z, y, c, m = 1, x, 0, n\n", " while m > 0:\n", " if m % 2 == 1:\n", " z = z * y\n", " c = c + 1\n", " y = y * y\n", " c = c + 1\n", " m = m // 2\n", " return (z, c)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "puissance_iter(2, 3000)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Bon, ça marche pour $2^{3000}$ et le nombre de multiplications est identique à celui de la fonction récursive. Mais voilà pourquoi les théoriciens adorent les fonctions récursives ! Parce que je ne sais pas ce qu'il en est pour vous mais j'ai beau user mes yeux sur cette fonction, je ne comprends pas ce qu'elle fait !!! \n", "\n", "Nous allons donc MONTRER que `puissance_iter(x, n)` renvoie effectivement $x^n$. Pour cela, nous allons mettre en évidence un *invariant de boucle*. C'est quoi ça ? Eh bien c'est \"quelque chose qui ne varie pas au cours de la boucle\". Plus précisément, un invariant de boucle est un objet mathématique dont la valeur avant une itération est la même que sa valeur après cette itération.\n", "\n", "Cet objet aura donc la même valeur après la dernière itération que celle qu'il avait avant la première ... à condition, évidemment que la boucle termine :-).\n", "\n", "Quel genre d'objet peut-on considérer ? Par exemple une propriété, dont la valeur est un booléen. Ou un nombre, dont la valeur est ... un nombre. Ou tout ce que l'on veut." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Ici, nous allons montrer que l'entier $zy^m$ est un invariant de boucle. Soient $m', z', y'$ les valeurs de $m,z,y$ avant une certaine itération, et $m'',z'',y''$ les valeurs de ces mêmes variables après cette itération. Deux cas sont à considérer :\n", "\n", "- Cas 1 : $m'=2p$ où $p>0$. Le test échoue, on a $z''=z'$, $y''=y'^2$ et $m''=p$. Donc, \n", "\n", "$$z''y''^{m''}=z'(y'^2)^p=z'y'^{2p}=z'y'^{m'}$$\n", "\n", "- Cas 2 : $m'=2p+1$ où $p\\ge0$. Le test réussit, on a $z''=z'y'$, $y''=y'^2$ et $m''=p$. Donc, \n", "\n", "$$z''y''^{m''}=z'y'(y'^2)^p=z'y'^{2p+1}=z'y'^{m'}$$\n", "\n", "Conséquence : la valeur de $zy^m$ AVANT la toute première itération est la même que celle APRÈS la toute dernière itération.\n", "\n", "- Valeur avant première itération : $1x^n=x^n$.\n", "- Valeur après dernière itération : $zy^0=z$.\n", "\n", "Ainsi, la valeur de $z$ après la dernière itération est $x^n$. Et cela tombe bien puisque notre fonction renvoie justement $z$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Remarque.** Quand je parle de la dernière itération, je suis très optimiste : la fonction termine-t-elle ? Mais oui, rassurons-nous. La boucle `while` est exécutée tant que $m>0$, c'est à dire $\\ell(m)>0$. Avant la première itération, $\\ell(m)=\\ell(n)$. À chaque itération, $m$ est divisé par 2, donc $\\ell(m)$ est diminué de 1. Ainsi, la valeur de $\\ell(m)$ après $k$ itérations est $\\ell(n)-k$. Conclusion :\n", "\n", "**Proposition.** Lors de l'appel à `puissance_iter(x, n)`, la boucle `while` effectue $\\ell(n)$ itérations.\n", "\n", "**Exercice.** Montrer que lors de la $k$ième itération, la fonction `puissance_iter` effectue 2 multiplications si le $k$ième chiffre de l'écriture binaire de $n$ est 1, et 1 multiplication si ce chiffre est 0.\n", "\n", "On en déduit facilement :\n", "\n", "**Proposition.** Le nombre de multiplications effectuées par la fonction itérative est le même que celui de la fonction récursive." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 1.8 Généralisons" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Souvenez vous de l'introduction de ce notebook où je parlais de puissances dans un monoïde. Était-ce juste pour faire savant ? Pas du tout. Il est très facile d'écrire une fonction générale d'exponentiation dans un monoïde quelconque. Il suffit de rajouter quelques paramètres.\n", "\n", "Un monoïde est caractérisé par\n", "\n", "- Son opération\n", "- Son élément neutre $e$\n", "\n", "Passons donc ces deux objets en paramètres ! Voici la version itérative." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def puissance(x, n, mul, e):\n", " z, y, m = e, x, n\n", " while m > 0:\n", " if m % 2 == 1: z = mul(z, y)\n", " y = mul(y, y)\n", " m = m // 2\n", " return z" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "puissance(7, 10 ** 1000, lambda x, y: x + y, 0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Eh oui, car un multiple n'est jamais qu'une puissance additive.\n", "\n", "Un autre exemple ? L'élévation à une puissance dans l'anneau des entiers modulo $p$, $\\mathbb Z/p\\mathbb Z$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def puissance_mod(x, n, p):\n", " return puissance(x, n, lambda x, y: (x * y) % p, 1)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "p = 2 ** 1279 - 1\n", "print(p)\n", "puissance_mod(3, (p - 1) // 2, p)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Ce calcul nous montre que le nombre 3 n'est pas un résidu quadratique modulo le nombre de Mersenne $M_{1279}=2^{1279}-1$. Voir le notebook à ce sujet ..." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Un dernier exemple, pour montrer l'extrême généralité de tout cela. L'élévation à la puissance $n$ dans le monoïde $\\mathbb R^{\\mathbb R}$ des fonctions de $\\mathbb R$ vers $\\mathbb R$ muni de la composition des applications. Et le neutre ? L'identité, bien sûr. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def puissance_fonc(f, n):\n", " g = puissance(f, n, lambda f, g:lambda x:f(g(x)), lambda x:x)\n", " return g" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Élevons la fonction $f:x\\mapsto\\sqrt{1+x}$ à la puissance 1000000. C'est à dire $f\\circ f\\ldots\\circ f$, un million de fois." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fmuch = puissance_fonc(lambda x:sqrt(1 + x), 1000000)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fmuch" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Question à 1.000.000 euros** : Remplacez 1000000 ci-dessus par $10^{1000}$, et appuyez sur Entrée. Le calcul termine en zéro seconde. Exponentiation rapide, certes, mais là vous allez devoir expliquer. Parce que pour faire un parallèle avec les nombres, le nombre de chiffres de $2^{10^{1000}}$ est très très supérieur au nombre d'atomes de l'univers. Que calcule _exactement_ Python lorsqu'on appuie sur Entrée ?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Calculons $f^{1000000}(0.3)$. Surtout n'oubliez pas de remettre 1000000 dans la cellule du dessus." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fmuch(0.3)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Tiens, le nombre d'or ? Ne le réveillons pas." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 2. Une illustration : les nombres de Fibonacci" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "À propos du nombre d'or ... Il existe de nombreux algorithmes pour calculer les nombres de Fibonacci. Le sujet de ce notebook étant l'exponentiation rapide, nous allons voir comment des calculs de puissances dans un anneau $\\mathbb A$ judicieusement choisi permettent d'obtenir très efficacement la valeur de ces nombres." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 2.1 Un anneau en or" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Notons $\\Phi$ la racine positive de l'équation $x^2=x+1$. On a $\\Phi=\\frac{1+\\sqrt 5} 2$. Cela implique que, tout comme $\\sqrt 5$, $\\Phi$ est irrationnel. Et surtout, que $\\Phi^2=\\Phi + 1$.\n", "\n", "Soit $\\mathbb A=\\{a+b\\Phi, a, b\\in\\mathbb Z\\}$. On montre facilement que $\\mathbb A$ est un sous-anneau de $\\mathbb R$. Par l'irrationalité de $\\Phi$, tout élement $z$ de $\\mathbb A$ s'écrit __de façon unique__ sous la forme $z=a+b\\Phi$, avec $a,b\\in\\mathbb Z$.\n", "\n", "Soient $z=a+b\\Phi$ et $z'=c+d\\Phi$ deux éléments de notre anneau. On a \n", "\n", "$zz'=(a+b\\Phi)(c+d\\Phi)=ac + (ad+bc)\\Phi + bd\\Phi^2$. \n", "\n", "Mais $\\Phi^2=\\Phi + 1$. En remplaçant, on obtient\n", "\n", "$$(a+b\\Phi)(c+d\\Phi)=(ac+bd) + (ad+bc+bd)\\Phi$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Que faire de tout cela en Python ? Modélisons les éléments de l'anneau $\\mathbb A$ par des couples d'entiers. Plus précisément, le couple $(a,b)$ représente l'élément $a+b\\Phi$. Par exemple, le neutre de l'anneau pour la multiplication, qui est bien entendu le nombre 1, est modélisé par le couple $(1,0)$. Et le nombre d'or $\\Phi=0+1\\Phi$ est modélisé par le couple $(0,1)$.\n", "\n", "La multiplication de deux éléments de notre anneau est immédiate à coder :" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def mul(x, y):\n", " a, b = x\n", " c, d = y\n", " return (a * c + b * d, a * d + b * c + b * d)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "De là, l'élévation de $z\\in\\mathbb A$ à la puissance $n$ est également immédiate, grâce à notre fonction générale d'exponentiation rapide." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def puissphi(z, n):\n", " return puissance(z, n, mul, (1, 0))" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "puissphi((0,1), 10)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 2.2 Quel rapport avec les nombres de Fibonacci ?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Le calcul précédent nous montre que $\\Phi^{10}=34+55\\Phi$. Évidemment ce 34 et ce 55 nous rappellent quelque chose." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "__Définition__ : la suite de Fibonacci $(F_n)_{n\\ge 0}$ est définie par\n", "\n", "- $F_0=0$\n", "- $F_1=1$\n", "- $\\forall n\\in \\mathbb N^*, F_{n+1}=F_{n-1} + F_n$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Les premières valeurs de la suite sont 0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55 ... " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Proposition.** pour tout $n\\in\\mathbb N^*$, on a $\\Phi^n=F_{n-1}+F_n\\Phi$.\n", "\n", "**Démonstration.** Faisons une récurrence simple sur $n$. \n", "\n", "Pour $n=1$ c'est immédiat. \n", "\n", "Soit donc $n\\ge 1$, supposons $\\Phi^n=F_{n-1}+F_n\\Phi$. \n", "\n", "On a alors $\\Phi^{n+1}=\\Phi^n\\Phi=(F_{n-1}+F_n\\Phi)\\Phi=F_{n-1}\\Phi+F_n\\Phi^2$.\n", "\n", "Mais $\\Phi^2=1+\\Phi$. On remplace et on obtient $\\Phi^{n+1}=F_n+(F_{n-1}+F_n)\\Phi=F_n+F_{n+1}\\Phi$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Nous voici donc en possession d'un algorithme pour calculer le $n$ième nombre de Fibonacci avec un nombre _logarithmique_ d'opérations : élever $\\Phi$ à la puissance $n$ et récupérer la partie $\\Phi$-maginaire, si j'ose dire. La fonction `fibonacci` fait le travail." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def fibonacci(n): return puissphi((0, 1), n)[1]" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "for k in range(20):\n", " print('%5d' % k, end='')\n", "print()\n", "for k in range(20):\n", " print('%5d' % fibonacci(k), end='')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 2.3 Le moment de vérité" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Combien de temps pour calculer le millionième nombre de Fibonacci ? Évaluez la cellule ci-dessous. " ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true }, "outputs": [], "source": [ "f = fibonacci(10 ** 6)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Pari gagné. Pour information, ce nombre possède 208988 chiffres." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import sys\n", "sys.set_int_max_str_digits(300000)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print(len(str(f)))" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "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.10.8" } }, "nbformat": 4, "nbformat_minor": 2 }