{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Opérations sur les polynômes\n", "\n", "Marc Lorenzi\n", "\n", "1er mars 2019" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "from sympy import S\n", "import math\n", "import random" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plt.rcParams['figure.figsize'] = (10, 5)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Voici deux compteurs, ils nous serviront à compter les additions et les multiplications effectuées par les fonctions que nous allons écrire. Si vous voyez apparaître dans le notebook une ligne du type `cnta += 1`, c'est qu'une somme (ou une différence) de coefficients de polynômes a été faite à la ligne précédente. Les fonctions incrémentant ces compteurs sont clairement signalées par une première ligne du type `global cnta`. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "cnta = 0\n", "cntm = 0" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 1. Introduction" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 1.1 Représentation des polynômes en Python" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Dans ce notebook, nous nous intéressons aux diverses \"opérations\" (au sens large) que l'on peut effectuer dans l'algèbre $\\mathbb R[X]$ des polynômes à coefficients réels. Un polynôme $A\\in\\mathbb R[X]$ s'écrit de façon unique\n", "\n", "$$A=\\sum_{k=0}^\\infty a_kX^k$$\n", "\n", "où les $a_k$ sont des réels tous nuls sauf un nombre fini. Les $a_k$ sont les __coefficients__ du polynôme $A$.\n", "\n", "Nous représenterons un tel polynôme par la liste de ses coefficients. Soyons plus précis :\n", "\n", "- Si $A=0$ est le polynôme nul, nous représenterons $A$ par la liste vide $[]$.\n", "- Si $A\\ne 0$, soit $n$ le plus grand entier naturel tel que $a_n\\ne 0$ (il n'a pas un nom, cet entier ?). Nous représenterons alors $A$ par la liste $[a_0,a_1,\\ldots,a_n]$.\n", "\n", "Il pourrait arriver qu'à la suite de certaines opérations nous obtenions un \"polynôme\" représenté par une liste finissant par des zéros. Il nous faudra alors __normaliser__ cette liste pour respecter nos conventions de représentation des polynômes. La fonction `normaliser` effectue ce travail. Elle prend en paramètre une liste et renvoie cette liste purgée des zéros finaux.\n", "\n", "Il pourra également nous arriver de représenter des polynômes sous forme non normalisée, c'est à dire par des listes finissant par des zéros. Nous le signalerons lorsque cela arrive." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def normaliser(A):\n", " n = len(A) - 1\n", " while n >= 0 and A[n] == 0: n = n - 1\n", " return A[:n + 1]" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "normaliser([1, 2, 0, 1, 0, 0])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 1.2 Évaluer un polynôme en un point" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Soit $A$ un polynôme représenté par la liste $[a_0, a_1, \\ldots, a_{n-1}]$. Soit $x\\in \\mathbb R$. Comment calculer $A(x)$ ? Ceci peut être fait en $O(n)$ opérations. L'idée est simple. Il suffit de remarquer que\n", "\n", "$$A(x)=a_{n-1}x^{n-1}+a_{n-2}x^{n-2}+\\ldots+a_0= (a_{n-1}x^{n-2}+a_{n-2}x^{n-3}+\\ldots+a_1)x + a_0$$\n", "\n", "En réitérant l'opération on obtient __l'algorithme de Horner__ :\n", "\n", "$$A(x)=(\\ldots((a_{n-1}x+a_{n-2})x+a_{n-3})x+\\ldots)x+a_0$$" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def evaluer(A, x):\n", " s = 0\n", " n = len(A)\n", " for k in range(len(A) - 1, -1, -1):\n", " s = s * x + A[k]\n", " return s" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "evaluer([1, 2, 1, 1], 2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "La boucle `for` effectue $n$ itérations. À chaque itération on effectue une addition et une multiplication, soit, au total, $2n$ opérations. Inutile d'en dire plus." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 1.3 Le degré d'un polynôme" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Soit $A=\\sum_{k=0}^\\infty a_kX^k\\in\\mathbb R[X]$.\n", "\n", "- Si $A\\ne 0$, le degré de $A$ est le plus grand entier $k$ tel que $a_k\\ne 0$.\n", "- Le degré du polynôme nul est $-\\infty$.\n", "\n", "La fonction `degre` renvoie le degré du polynôme $A$. Enfin, presque. Si $A$ est le polynôme nul la fonction renvoie $-1$. Il faudra nous en souvenir dans la suite. Je rappelle que $A$ est implicitement supposé normalisé. Si ce n'est pas le cas, la fonction `degre` renvoie un résultat inexact." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def degre(A):\n", " return len(A) - 1" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "degre([1, 3, 2, 1])" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "degre([])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 1.4 Représentations lisibles d'un polynôme" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "La représentation d'un polynôme sous forme de liste est parfaitement lisible. Mais bon, pour faire plaisir aux fanatiques de l'indéterminée $X$, voici une fonction `vers_chaine`. Elle prend un polynôme $A$ en paramètre et renvoie une chaîne de caractères qui est la représentation \"classique\" du polynôme $A$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def vers_chaine(A):\n", " n = degre(A)\n", " if n == -1: return \"0\"\n", " s = ''\n", " for k in range(n + 1):\n", " fmt = '%s*X**%s' \n", " if A[k] > 0: s += ('+' + fmt) % (A[k], k)\n", " elif A[k] < 0: s += fmt % (A[k], k)\n", " return s" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "vers_chaine([1, -3, 0, 1])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Prévoyons aussi une représentation sous forme de graphique, pour les artistes. Cela ne nous coûte pas grand chose." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def plot_poly(A):\n", " plt.plot(A, '-k')\n", " plt.plot(A, 'or')\n", " plt.grid()" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true }, "outputs": [], "source": [ "plot_poly([1, 3, 7, -1, 1, 2])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 1.5 Coefficients" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Nous allons dans la suite passer notre temps à accéder aux coefficients d'un polynôme. Facile, me direz-vous, le coefficient de degré $k$ du polynôme $A$, c'est $A[k]$ ? Eh bien pas exactement. Si $k\\le d^o A$, c'est vrai. Mais sinon ? Histoire d'éviter des problèmes, écrivons illico une fonction `coef`. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def coef(A, k):\n", " if k <= degre(A): return A[k]\n", " else: return 0" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Et tant que nous y sommes ..." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def coef_dominant(A):\n", " n = degre(A)\n", " if n < 0: raise Exception('Polynôme nul')\n", " else: return A[-1]" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "coef_dominant([1, 3, 2, 7])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 2. Somme, produit par un scalaire" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 2.1 Somme et différence" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Soient $A=\\sum_{k=0}^\\infty a_k X^k$ et $B=\\sum_{k=0}^\\infty b_k X^k$ deux polynômes. La somme de $A$ et $B$ est le polynôme\n", "\n", "$$A+B=\\sum_{k=0}^\\infty (a_k + b_k) X^k$$\n", "\n", "On définit de même $A-B$. Rappelons que\n", "\n", "$$d^o(A+B)\\le\\max(d^oA,d^oB)$$\n", "\n", "et de même pour $A-B$. Les fonctions `somme` et `difference` ne présentent pas de difficulté. Il faut juste ne pas oublier de normaliser le résultat avant de le renvoyer.\n", "\n", "Ces deux fonctions prennent un paramètre optionnel `norm`. Si `norm` vaut `True`, ce qui est vrai par défaut, le polynôme renvoyé est normalisé." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def somme(A, B, norm=True):\n", " global cnta\n", " n = max(degre(A), degre(B))\n", " R = (n + 1) * [0]\n", " for k in range(n + 1):\n", " R[k] = coef(A, k) + coef(B, k)\n", " cnta += 1\n", " if norm: return normaliser(R)\n", " else: return R" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "somme([1, 2, 1], [2, 3, -1])" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "somme([], [1, 2, 3])" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def difference(A, B, norm=True):\n", " global cnta\n", " n = max(degre(A), degre(B))\n", " R = (n + 1) * [0]\n", " for k in range(n + 1):\n", " R[k] = coef(A, k) - coef(B, k)\n", " cnta += 1\n", " if norm: return normaliser(R)\n", " else: return R" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "difference([1, 2, 1], [2, 3, -1])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 2.2 Produit par un scalaire" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Soit $A=\\sum_{k=0}^\\infty a_k X^k$un polynôme. Soit $t\\in\\mathbb R$. La produit de $t$ et $A$ est le polynôme\n", "\n", "$$tA=\\sum_{k=0}^\\infty (ta_k) X^k$$\n", "\n", "Aucune difficulté pour écrire une fonction `prod_scal`. Gardons juste à l'esprit que toutes nos fonctions doivent renvoyer des polynômes normalisés. Il faut donc traiter à part le cas où $t=0$. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def prod_scal(t, A):\n", " global cntm\n", " if t == 0: return []\n", " else:\n", " cntm += len(A)\n", " return [t * c for c in A]" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "prod_scal(2, [1, 3, 2])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 3. Produit, puissances" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 3.1 Produit de deux polynômes" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Soient $A=\\sum_{k=0}^\\infty a_k X^k$ et $B=\\sum_{k=0}^\\infty b_k X^k$ deux polynômes. Le produit de $A$ et $B$ est le polynôme\n", "\n", "$$AB=\\sum_{k=0}^\\infty c_k X^k$$\n", "\n", "où $c_k=\\sum_{i=0}^k a_ib_{k-i}$.\n", "\n", "Rappelons que $d^o(AB)=d^oA+d^oB$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Écrivons tout d'abord une fonction `somme_deg` qui ajoute des degrés. Euh, n'est-ce pas tout simplement \"$+$\" ? Non, à cause du polynôme nul. À cause des conventions que nous avons faites, nous devons faire en sorte que $-1$ soit absorbant pour l'addition (si j'ose dire)." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def somme_deg(m, n):\n", " if m == -1 or n == -1: return -1\n", " else: return m + n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def produit(A, B):\n", " global cnta, cntm\n", " p = somme_deg(degre(A), degre(B))\n", " C = (p + 1) * [0]\n", " for k in range(p + 1):\n", " for i in range(k + 1):\n", " C[k] += coef(A, i) * coef(B, k - i)\n", " cnta += 1\n", " cntm += 1\n", " return C" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "cnta, cntm = 0, 0\n", "print(produit([1, 2, 1], [1, 3, 3, 1]))\n", "print(cnta, cntm)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Quelle est la complexité de la fonction `produit` en nombre d'opérations sur des coefficients de polynômes ? Soient $m=d^o A$ et $n = d^o B$. Supposons $A$ et $B$ non nuls. Le degré de $C=AB$ est $p=m+n$, nous avons donc $p+1$ coefficients à calculer. Pour $k=0,1,\\ldots,p$, on a\n", "\n", "$$c_k=\\sum_{i=0}^k a_ib_{k-i}$$\n", "\n", "Le calcul de $c_k$ nécessite donc $k+1$ additions et $k+1$ multiplications. Vous me direz qu'il n'y a que $k$ additions, mais regardez le code de `produit` ... \n", "\n", "Le nombre total d'additions ou de multiplications à effectuer est donc\n", "\n", "$$\\sum_{k=0}^{p} (k+1)=\\frac 1 2 (p+1)(p+2)$$\n", "\n", "Le nombre total d'opérations est donc $(p+1)(p+2)$.\n", "\n", "Sur l'exemple ci-dessus, $m=2$, $n=3$ et $p=5$. D'où $\\frac 1 2(p+1)(p+2)=\\frac 1 26\\times 7=21$, ce qui est bien la valeur des compteurs en fin de calcul." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Notre algorithme de multiplication a donc une complexité en $O((m+n)^2)$. Pour des polynômes de degrés raisonables, cela ne pose pas de problème. En revanche, cela n'est pas acceptable lorsque nos polynômes possèdent de grands degrés. Nous verrons à la fin de ce notebook un algorithme dû à Anatolii Alexevich Karatsuba meilleur que notre algorithme \"naïf\".\n", "\n", "__Remarque__ : \"naïf\", en algorithmique, cela veut dire la même chose que \"j'utilise les définitions de mon cours de maths\" :-). \n", "\n", "__Exercice__ : Les mathématiciens sont-ils naïfs ?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 3.3 Puissances" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Maintenant que nous savons calculer un produit, calculer une puissance ne pose pas de problème. La fonction `puissance` prend en paramètres un polynôme $A$ et un entier $n$. Elle renvoie $A^n$. Cette fonction utilise l'algorithme d'exponentiation rapide. Si vous ne savez pas ce que sait, allez consulter le notebook à ce sujet." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def puissance(A, n):\n", " if n == 0: return [1]\n", " else:\n", " B = puissance(produit(A, A), n // 2)\n", " if n % 2 == 0: return B\n", " else: return produit(A, B)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "puissance([1, 1], 30)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "De façon plus graphique, voici $(1+X)^{30}$. Si vous reconnaissez des coefficients binomiaux ne soyez pas surpris." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plot_poly(puissance([1, 1], 30))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 4. Composée" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$\\mathbb R[X]$ est également muni d'une opération de composition. Soient $A=\\sum_{k=0}^\\infty a_k X^k$ et $B$ deux polynômes. La composée de $A$ et $B$ est\n", "\n", "$$A\\circ B=\\sum_{k=0}^\\infty a_k B^k$$\n", "\n", "Si $B$ n'est pas constant, on a $d^o(A\\circ B)=d^oA\\times d^o B$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def composee(A, B):\n", " C = []\n", " for k in range(degre(A) + 1):\n", " C = somme(C, prod_scal(coef(A, k), puissance(B, k)))\n", " return C" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "composee([1, 2, 3, 1], [1, 1, 2])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 5. Division" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 5.1 La division euclidienne" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "__Proposition__ : Soient $A, B\\in\\mathbb R[X]$, $B\\ne 0$. Il existe un unique couple $(Q,R)$ de polynômes tel que $A=BQ+R$ et $d^o R= m:\n", " c = coef_dominant(R) / coef_dominant(B)\n", " cntm += 1\n", " D = decaler([c], n - m)\n", " R = difference(R, decaler(prod_scal(c, B), n - m))\n", " Q = somme(Q, D)\n", " n = degre(R)\n", " return (Q, R)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "La fonction `decaler` multiplie $A$ par $X^n$, c'est à dire rajoute $n$ zéros en tête de la liste $A$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def decaler(A, n):\n", " return (n * [0]) + A" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "decaler([1, 2, 3], 5)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "A = S([1, 4, 6, 4, 1])\n", "B = [1, 3, 2]\n", "Q, R = division(A, B)\n", "print(Q, R)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Petite vérification ..." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "somme(produit(B, Q), R)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "__Remarque__ : À quoi sert le $S$ dans l'exemple ci-dessus ? À effectuer des calculs exacts. Sans celui-ci, les calculs seraient faits en flottants." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 5.2 Complexité" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Combien d'opérations sur des coefficients de polynômes doit-on faire lors d'une division euclidienne ? Soient $A$, $B$ de degrés respectifs $n$ et $m$, où $n\\ge m$. Le pire cas se produit (preuve ?) lorsque, à chaque itération, le degré du reste chute seulement de 1. On effectue alors $n-m+1$ itérations que nous allons numéroter $k=0,1,\\ldots,n-m$.\n", "\n", "- Itération 0 : le calcul de $c$ coûte une division. $D$ ne coûte rien, simple décalage. Pour calculer $R$, on calcule d'abord $cB$ : $m+1$ multiplications. Puis on fait la différence $R-cDB$ : $n+1$ soustractions. Le calcul de $Q$ demande $n-m+1$ additions, la plupart, avouons-le, triviales. Total : $1+(m+1)+(n+1)+(n-m+1)=2n+4$ opérations.\n", "\n", "- Itération 1 : $R$ est maintenant de degré $n-1$. le calcul de $c$ coûte une division. $D$ est toujours gratuit. Pour calculer $R$, on calcule d'abord $cB$ : $m+1$ multiplications. Puis on fait la différence $R-cDB$ : $n$ soustractions. Le calcul de $Q$ demande toujours $n-m+1$ additions. Total : $1+(m+1)+n+(n-m+1)=2n+3$ opérations.\n", "\n", "- Itération 2 : même raisonnement, $2n+2$ opérations.\n", "\n", "- ...\n", "\n", "- Itération $k$ : $2n+4-k$ opérations.\n", "\n", "Au total, on effectue donc en pire cas $\\sum_{k=0}^{n-m}(2n+4-k)=\\frac 1 2(n-m+1)(3n+m+8)$ opérations.\n", "\n", "Petite vérification ?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "cnta, cntm = 0, 0\n", "A = 101 * [1]\n", "B = [1, 2, 3, 4, 1]\n", "Q, R = division(A, B)\n", "print(cnta + cntm)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "n = 100\n", "m = 4\n", "print((n - m + 1) * (3 * n + m + 8) / 2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Nous sommes comblés.\n", "\n", "__Remarque__ : Notons $C_\\times$ et $C_\\div$ les complexités de la multiplcation et de la division euclidienne. Nous avons\n", "\n", "$$C_\\times = (m+n+1)(m+n+2)$$\n", "\n", "et, en pire cas,\n", "\n", "$$C_\\div=\\frac 1 2(n-m+1)(3n+m+8)$$\n", "\n", "Histoire de comparer des choses comparables, fixons $m$ et faisons tendre $n$ vers l'infini. Il vient\n", "\n", "$$C_\\times\\sim n^2 \\quad{\\rm\\ et}\\quad C_\\div\\sim \\frac 3 2 n^2$$\n", "\n", "__Faire une division est à peine 50% plus coûteux que faire une multiplication.__" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 5.3 Exercices" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Nous avons maintenant défini les opérations présentes dans un tout bon cours de maths sur les polynômes. L'existence d'une division euclidienne dans l'anneau $\\mathbb R[X]$ fait de celui-ci un __anneau euclidien__. Dans un tel anneau on dispose des notions de pgcd et de ppcm. Le théorème de Bézout y est valide.\n", "\n", "1. Écrire une fonction `pgcd` prenant en paramètres deux polynômes $A$ et $B$ et renvoyant leur pgcd. Cette fonction implémentera bien entendu le célèbre algorithme d'Euclide.\n", "2. Écrire une fonction `ppcm`. C'est immédiat à partir de `pgcd`.\n", "3. Écrire une fonction `bezout` prenant en paramètres deux polynômes $A$ et $B$ et renvoyant un triplet $(U,V,\\Delta)$ de polynômes tel que $UA+VB=\\Delta$ et $\\Delta = A\\land B$. Il suffit de reprendre de façon un peu plus poussée la fonction de l'exercice 1." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 6. L'algorithme de Karatsuba" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Dans cette section, nous ne supposons pas forcément les polynômes représentés sous forme normalisée. Nous appellerons __longueur__ d'un polynôme la longueur d'une liste qui le représente. Cette quantité est mal définie, nous devrions plutôt dire la longueur __d'une__ représentation du polynôme par une liste. Mais faisons court ..." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 6.1 L'idée" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Soit $p\\in\\mathbb N^*$. Soient $A, B$ deux polynômes de longueur $n=2p$. Écrivons $A=A'+A''X^p$ et $B=B'+B''X^p$ où les polynômes $A',A'',B',B''$ sont de longueur $p$. On a\n", "\n", "$$AB=(A'+A''X^p)(B'+B''X^p)=A'B'+(A'B''+A''B')X^p+A''B''X^{2p}$$\n", "\n", "Pour multiplier $A$ et $B$ il suffit donc de faire __quatre__ multiplications de polynômes de longueur $p$ (je vous laisse compter). On ne compte pas les multiplications par des puissances de $X$, celles-ci peuvent être effectuées par de simples décalages. Peut-on faire mieux ? Oui !\n", "\n", "Posons $C=(A'+A'')(B'+B'')$. On a $C=A'B'+A''B''+A'B''+A''B'$. Ainsi, si l'on pose $C'=A'B'$ et $C''=A''B''$, on a\n", "\n", "$$AB=C'+(C-C'-C'')X^p+C''X^{2p}$$\n", "\n", "On constate donc que __trois__ multiplications suffisent.\n", "\n", "Ceci suggère donc un algorithme de type \"diviser pour régner\" pour multiplier deux polynômes. On coupe les polynômes en deux, on effectue trois produits de polynômes de longueurs \"deux fois plus petites\", on recolle les morceaux. Il n'y a plus qu'à écrire du code. Mais avant tout, voyons en quoi nous y gagnons.\n", "\n", "__Remarque__ : Si $n=2p+1$ (où $p\\ge 0$) est un nombre impair, on peut décider, par exemple, de prendre les polynômes $A$ et $B$ \"prime\" de taille $p$ et les polynômes $A$ et $B$ \"seconde\" de taille $p+1$. Cela nous donne bien entendu les mêmes formules." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 6.2 Complexité" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Notons $M_n$ le nombre de multiplications de coefficients à effectuer pour calculer avec l'algorithme de Karatsuba le produit de deux polynômes de longueur $n$. \n", "\n", "- Nous n'avons pas parlé d'initialisation pour l'instant, mais Pour $n=0$ il n'y a rien à faire. Ainsi, $M_0 = 0$.\n", "- Pour $n=1$, on a $M_1=1$. Multiplier deux polynômes constants ce n'est pas bien délicat.\n", "- Soit $p\\ge 1$. On a clairement $M_{2p}=3M_p$.\n", "- Enfin, soit $p\\ge 0$. On a $M_{2p+1} = 2M_{p+1} + M_p$. Je vous laisse regarder pourquoi.\n", "\n", "Je vais ici me contenter de regarder les polynômes dont la longueur est une puissance de 2. Que vaut $u_n=M_{2^n}$ pour $n\\in\\mathbb N$ ?\n", "\n", "- Pour $n=0$, on a $2^n=1$ : $u_0=M_1=1$.\n", "- Pour $n\\ge 0$, $u_{n+1}=M_{2^{n+1}}=3M_{2^n}=3u_n$.\n", "\n", "Ainsi, la suite $(u_n)_{n\\ge 0}$ est géométrique de raison 3, de premier terme 1. Pour tout entier $n$, $u_n=3^n$. \n", "\n", "Si $n$ est une puissance de 2, nous avons donc \n", "\n", "$$M_n=3^{\\lg n}=n^{\\lg 3}$$\n", "\n", "où $\\lg$ désigne le logarithme en base 2." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Pour l'algorithme naïf, nous avions besoin (voir paragraphe 3.2) de $n(2n-1)\\sim 2n^2$ opérations. Y a-t-on gagné ? oui, car $\\lg 3 < \\lg 4 = 2$, et donc $n^{\\lg 3}=o(n^2)$. L'algorithme de Karatsuba est donc asymptotiquement meilleur que l'algorithme naïf.\n", "\n", "Euh, me direz-vous, et les additions dans tout ça ? Il faudrait peut-être les compter aussi ? Vous avez raison, nous le ferons, mais attendons d'avoir écrit la fonction `karatsuba`." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 6.3 Complexité : esquisse du cas général" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Si $n$ n'est pas une puissance de 2, $M_n$ n'est certainement pas égal à $n^{\\lg 3}$ vu que ce nombre n'est en général pas un entier, alors que $M_n$ en est un. On peut cependant démontrer que l'ordre de grandeur est le bon : $M_n=O(n^{\\lg 3})$ lorsque $n$ tend vers l'infini. Je ne le ferai pas ici, je me contenterai d'une petite illustration. Nous avons $M_0=0$, $M_1=1$, et la récurrence\n", "\n", "$$M_{2p}=3M_p, M_{2p+1}=2M_{p+1}+M_p$$\n", "\n", "Voici donc la fonction `nbmult_kara` qui prend en paramètre un entier $n$ et renvoie $M_n$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def nbmult_kara(n):\n", " if n == 0: return 0\n", " elif n == 1: return 1\n", " else:\n", " p = n // 2\n", " if n % 2 == 0: return 3 * nbmult_kara(p)\n", " else: return 2 * nbmult_kara(p + 1) + nbmult_kara(p)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print([nbmult_kara(n) for n in range(20)])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Ci dessous, \n", "\n", "- en noir le tracé de $M_n-n^{\\lg 3}$.\n", "- en rouge le tracé de $\\frac 2 5 n^{\\lg 3}$. Le facteur $\\frac 2 5$, choisi par tâtonnements, n'est sans doute pas optimal.\n", "- en bleu l'axe $Ox$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "l = math.log(3, 2)\n", "xs = range(0, 1025)\n", "ys = [nbmult_kara(n) - n ** l for n in xs]\n", "zs = [0.4 * n ** l for n in xs]\n", "ws = [0 for n in xs]\n", "plt.plot(xs, ys, 'k')\n", "plt.plot(xs, zs, 'r')\n", "plt.plot(xs, ws, 'b')\n", "plt.grid()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "La courbe noire mériterait sans doute un notebook à elle seule ... on peut en tout cas espérer (?) la\n", "\n", "__Conjecture non prouvée__ : Pour tout $n\\in\\mathbb N$, \n", "\n", "$$n^{\\lg 3}\\le M_n\\le \\frac 7 5 n^{\\lg 3}$$\n", "\n", "__Exercice__ : Pour $n=2^p+2^{p-1}$, $p\\ge 1$, montrer que\n", "\n", "$$M_n=Cn^{\\lg 3}$$\n", "\n", "où\n", "\n", "$$C=\\frac 7 3 \\left(\\frac 2 3\\right)^{\\lg 3}\\simeq 1.227$$\n", "\n", "Exemple motivant pour l'exercice :" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "7 / 3 * (2 / 3 * (1024 + 512)) ** math.log(3, 2)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "nbmult_kara(1024 + 512)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Bon, assez tergiversé, codons l'algorithme de Karatsuba." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 6.4 L'algorithme" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "pour ne pas rendre le code inutilement compliqué, nous allons supposer que la fonction `karatsuba` que nous allons écrire prend en paramètres deux polynômes $A$ et $B$ de même longueur. Nous verrons plus loin comment procéder lorsque les longueurs sont distinctes. $A$ et $B$ sont donc représentés par des listes de même taille." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "__Remarque__ : Un tout petit détail avant de nous lancer. Pour que nos polynômes soient bien représentés par des listes de même taille, il ne faut pas normaliser les sommes et les différences, ce qui explique dans le code ci-dessous la présence du paramètre `False`. Dans le cas $n=1$ (polynômes de degré 0), j'avoue avoir un peu abusé : par souci de cohérence la fonction renvoie une liste de taille 2, qui est en fait un polynôme de degré 0, non normalisé." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def karatsuba(A, B):\n", " global cntm\n", " n = len(A)\n", " if n == 0: return []\n", " elif n == 1: \n", " cntm += 1\n", " return [coef(A, 0) * coef(B, 0), 0]\n", " else:\n", " p = n // 2\n", " A1, A2 = (A[:p], A[p:])\n", " B1, B2 = (B[:p], B[p:])\n", " C1 = karatsuba(A1, B1)\n", " C2 = karatsuba(A2, B2)\n", " C = karatsuba(somme(A1, A2, False), somme(B1, B2, False))\n", " D = difference(difference(C, C1, False), C2, False)\n", " U = decaler(D, p)\n", " V = decaler(C2, 2 * p)\n", " return somme(somme(C1, U, False), V, False)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Maintenant, testons. Prenons $A=(1+X)^{15}$. Calculons $A^2$ avec `karatsuba` et `produit`, et comparons." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "A = puissance([1, 1], 15)\n", "cnta, cntm = 0, 0\n", "B = karatsuba(A, A)\n", "print(B)\n", "print(cnta, cntm)\n", "cnta, cntm = 0, 0\n", "C = produit(A, A)\n", "print(cnta, cntm)\n", "print(difference(B, C))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Tout a l'air de fonctionner. De plus :\n", "\n", "- `karatsuba` a effectué 845 additions et 81 multiplications.\n", "- `produit` a effectué 496 additions et 496 multiplications." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 6.5 Combien d'opérations, précisément ?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Nous avons déjà compté le nombre de multiplications de coefficients effectuées par l'algorithme de Karatsuba. Mais combien d'additions de coefficients, au juste ? Notons $A_n$ le nombre d'additions nécessaires pour le produit de deux polynômes de longueur $n$.\n", "\n", "Prenons pour $A$ et $B$ deux polynômes de longueur paire $n=2p$, où $p\\ge 1$.\n", "\n", "Quelles sont les longueurs des polynômes qui interviennent dans la fonction `karatsuba` ? Faisons un petit briefing.\n", "\n", "- $A_1, A_2, B_1, B_2$ : longueur $p$.\n", "- $C_1, C_2$ : longueur $2p$.\n", "- $C$ : longueur $2p$. Rappelons que nous ne normalisons pas les sommes.\n", "- $D$ : longueur $2p$.\n", "- $U$ : longueur $3p$.\n", "- $V$ : longueur $4p$.\n", "\n", "- Calcul de $C$ : 2 appels à `somme`. Ces appels se font sur des polynômes de longueur $p$ : il nous faut faire pour cela $2p$ additions de coefficients.\n", "\n", "- Calcul de $D$ : 2 appels à `difference`. Comme les polynômes $C$, $C_1$ et $C_2$ sont de longueur $2p$, il nous faut faire $2\\times 2p=4p$ additions de coefficients.\n", "\n", "- Nous avons également deux appels à `somme` à la dernière ligne de la fonction `karatsuba`. $U$ est de longueur $3p$ et $V$ est de longueur $4p$. Cela fait donc encore $7n$ additions de coefficients.\n", "\n", "- Enfin, il ne faut pas oublier les additions effectuées lors des 3 appels récursifs : cela nous en fait encore $3A_p$ additions.\n", "\n", "Finalement,\n", "\n", "$$A_{2p}=3A_p+13p$$\n", "\n", "Et, bien entendu, $A_0=A_1=0$.\n", "\n", "Le 13, avouons-le, nous fait un peu peur. En plus, 13 ça porte malheur. Mais pas d'inquiétude, on n'est pas superstitieux :-).\n", "\n", "__Exercice__ : Déterminer une relation reliant, pour $p\\ge 0$, les quantités $A_{2p+1}$, $A_p$ et $A_{p+1}$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Que vaut précisément $A_n$ ? Nous allons regarder uniquement le cas où $n$ est une puissance de 2.\n", "\n", "Posons $v_n=A_{2^n}$. On a $v_0=A_1=0$ et, pour tout $n\\in\\mathbb N$,\n", "\n", "$$v_{n+1}=3v_n+13\\times 2^n$$\n", "\n", "__Proposition__ : Pour tout entier naturel $n$, $v_n=13(3^n-2^n)$.\n", "\n", "__Démonstration__ : Récurrence sur $n$. Pour $n=0$ c'est clair. Soit $n\\in\\mathbb N$. Supposons la propriété vraie pour $n$. On a alors\n", "\n", "$$v_{n+1}=3v_n+13\\times 2^n= 39(3^n-2^n)+13\\times 2^n=13(3^{n+1}-2^{n+1})$$\n", "\n", "Lorsque $n$ est une puissance de 2, nous avons donc\n", "\n", "$$A_n = 13(n^{\\lg 3}-n)$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Rappelons nous le nombre de multiplications : $M_n=n^{\\lg 3}$. Le nombre total d'opérations effectuées par l'algorithme de Karatsuba lorsque $n$ est une puissance de 2 est donc\n", "\n", "$$K_n=14n^{\\lg 3}-13n$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Pour $n = 4$, par exemple :" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "13 * (3 ** 4 - 2 ** 4)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "On retrouve bien les 845 additions vues un peu plus haut." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 6.6 Karatsuba vs naïf. Qui gagne ?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Alors, quel est le meilleur algorithme ? Karatsuba ou l'algorithme naïf ? Pour de grandes valeurs de $n$ c'est clairement l'algorithme de Karatsuba. Mais grandes comment ? Grandes beaucoup ou grandes pas beaucoup ? Faisons un petit calcul.\n", "\n", "Voici deux fonctions renvoyant le nombre d'opérations effectuées par nos algorithmes pour des polynômes de longueur $n$ puissance de 2." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def ops_karatsuba(n):\n", " return 14 * n ** math.log(3, 2) - 13 * n\n", "\n", "def ops_naif(n):\n", " return 2 * n * (2 * n - 1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Faisons la différence, évaluons, et regardons ce que l'on obtient. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def phi(n):\n", " return ops_naif(n) - ops_karatsuba(n)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "s = [phi(k) for k in range(21)]\n", "print(s)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "On voit que $\\varphi(k)$ devient positif pour $k\\ge 13$ (le 13 est une pure coïncidence). L'algorithme de Karatsuba effectue donc moins d'opérations que l'algorithme naïf pour des polynômes de degré supérieur ou égal à 12. Petit dessin pour ceux qui ne veulent pas compter jusqu'à 13 ?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plt.plot(s, '-or')\n", "plt.plot(s, '-k')\n", "plt.grid()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Dans les faits, notre fonction `karatsuba` n'est pas très efficace. Pour ne citer que quelques uns de ses défauts :\n", "\n", "- appels récursifs entraînant des recopies de listes.\n", "- les appels à `decaler` ne sont pas gratuits.\n", "- Des additions et des soustractions sur des polynômes contenant beaucoup de zéros. Clairement, le nombre d'additions non triviales est bien moins important que notre estimation.\n", "\n", "Telle quelle, notre fonction ne devient plus efficace que la fonction naïve que pour des polynômes de degré environ 100 (testez !). Allez, un dernier test avec un gros polynôme ..." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 6.7 Un dernier test" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Prenons $A=1+X+X^2+\\ldots+X^{4095}$ et calculons $A^2$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "A = 4096 * [1]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "D'abord, $A^2$ avec `karatsuba`. Environ 14 secondes sur ma machine." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "cnta, cntm = 0, 0\n", "B = karatsuba(A, A)\n", "print(cnta, cntm)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Puis $A^2$ avec `produit`. Environ 1 minute sur ma machine. Notre fonction `karatsuba` n'est donc pas totalement sans intérêt :-)." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "cnta, cntm = 0, 0\n", "B = produit(A, A)\n", "print(cnta, cntm)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "__Exercice__ : Calculez explicitement $A^2$. Réponse : 1, 2, 3, 4, et ensuite ?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 6.8 Et si les polynômes n'ont pas la même longueur ?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Que faire ? Il suffit par exemple de \"dénormaliser\" le polynôme de plus petite longueur pour que sa longueur soit égale à celle du second. Cela ne ne changera pas l'ordre de grandeur de la complexité, $O(n^{\\lg 3})$, où $n$ est la plus grande des longueurs des deux polynômes." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Voici notre fonction finale `produit_kara`. Elle renvoie, histoire de faire propre, un polynôme normalisé." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def produit_kara(A, B):\n", " m, n = len(A), len(B)\n", " if m == 0 or n == 0: return []\n", " else:\n", " A = A + (n - m) * [0]\n", " return normaliser(karatsuba(A, B))" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print(produit_kara([1, 2, 3], [3, 4, 5, 6, 7]))\n", "print(produit([1, 2, 3], [3, 4, 5, 6, 7]))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Nous __savons__ que `produit_kara` fonctionne. Mais comme un certain Matthieu, nous voudrions voir pour croire. Histoire de clôre le dossier, calculons 1000 produits de polynômes \"aléatoires\" avec `produit_kara` et `produit`, et comparons les résultats." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def random_poly():\n", " n = random.randint(0, 20)\n", " s = n * [0]\n", " for k in range(n):\n", " s[k] = random.randint(-100, 100)\n", " return s" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print(random_poly())" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def test(n):\n", " b = True\n", " for k in range(n):\n", " A = random_poly()\n", " B = random_poly()\n", " C = produit_kara(A, B)\n", " D = produit(A, B)\n", " E = difference(C, D)\n", " if E != []:\n", " b = False\n", " break\n", " print(b)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "test(1000)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Nous avons de bonnes raisons, théoriques et pratiques, de penser que `produit_kara` fonctionne bien :-)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 6.9 Peut-on faire mieux ?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Oui, en utilisant la transformée de Fourier rapide (FFT, voir le notebook à ce sujet). Pour multiplier deux polynômes de longueur $n$ :\n", "\n", "- On calcule leurs transformées de Fourier en ajustant auparavant les polynômes à la longueur $2n$ : $O(n\\lg n)$ opérations.\n", "- On multiplie les transformées, qui sont des listes de longueur $2n$, __terme à terme__ : $2n$ opératiions.\n", "- On calcule la transformée de Fourier inverse du résultat : $O(n\\lg n)$ opérations.\n", "\n", "On obtient alors le produit des deux polynômes. Cet algorithme permet de multiplier deux polynômes de longueur $n$ en $O(n\\lg n)$ opérations. Il devient efficace par rapport à l'algorithme de Karatsuba lorsque les polynômes sont de degré élevé.\n", "\n", "L'algorithme basé sur la FFT est, à ma connaissance, le meilleur algorithme connu pour multiplier deux polynômes." ] } ], "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 }