Factorisation de Horner

Présentation du problème.

On se donne un anneau intègre \(A\) — il n'admet aucun diviseur de \(0\) — muni de ses lois additive \(+\) et multiplicative \(\times\) d'éléments neutres respectifs notés \(0\) et \(1\) et supposés distincts.

Un polynôme à coefficients dans \(A\) est défini comme une suite \(({\color{orange}a_k})_{k\in\N}\) d'éléments de \(A\) nulle à partir d'un certain rang : \begin{equation} \exists N\in\N\ \ \forall n\in\N\quad n\geq N\then a_n=0. \end{equation} Les valeurs \(a_k\) sont appelés les coefficients du polynôme. Si \(\forall k\in\N\ a_k=0,\) le polynôme est appelé polynôme nul. On appelle degré d'un polynôme non-nul, le plus grand indice \(k\in\N\) tel que \(a_k\neq 0\). Si un seul terme \(a_k\) n'est pas nul, le polynôme est appelé monôme de degré \(k\). La structure adaptée pour coder un polynôme de degré \(d\) est donc typiquement une liste de longueur \(d+1\).

On distingue structurellement un polynôme de sa fonction polynomiale car contrairement à \(\R\) ou \(\C,\) pour des anneaux finis \(A,\) il n'y a pas nécessairement bijection entre les polynômes et les fonctions po­ly­no­mia­les. Nous avons déjà rencontré ce problème de représentation avec les expressions polynomiales différentes d'une même fonction booléenne (à coefficients dans un anneau fini, le corps \(\F_2=\{0,1\}\)), justifiant l'introduction des formes normales canoniques.

On exprime généralement un polynôme \(P\) de degré \(n\) à l'aide d'une somme formelle faisant intervenir une indéterminée, à savoir un symbole \(X\) : \begin{equation}\label{exp} P(X) = {\color{orange}a_0} + {\color{orange}a_1}X + {\color{orange}a_2}X² + \cdots + {\color{orange}a_{n-1}}X^{n-1} + {\color{orange}a_n}X^n. \end{equation} À l'aide des opérations de l'anneau \(A,\) on munit l'ensemble des polynômes à coefficients dans \(A\) de trois opérateurs, une addition et une multiplication internes \(+\) et \(\times,\) et une multiplication externe \(\cdot\) sur \(A\). Notons \(I_k:=\ab{0}{k}\). Soit \(P\) et \(Q\) deux polynômes de degrés respectifs \(n\) et \(m\) et de coefficients \((a_i)\) et \((b_j)\) et \(a\in A\) : \begin{align*} P(X)+Q(X)&:=\sum_{\ell\in I_{\max\{n,m\}}}({\color{orange}a_\ell+b_\ell})X^\ell,\\ P(X)\times Q(X)&:=\sum_{\ell\in I_{nm}} \Big(\sum_{\overset{\scriptstyle(i,j)\in I_n\times I_m}{i+j=\ell}} {\color{orange}a_ib_j}\Big)X^\ell\\ a\cdot P(X)&:=\sum_{i\in I_n}(a.{\color{orange}a_i})X^i. \end{align*}

Attention, si la somme \({\color{orange}a_\ell+b_\ell}\) ci-dessus est bien définie mathématiquement pour tout \(\ell\in\N\) puisque par définition les coefficients de rang supérieur au degré d'un polynôme sont nuls, elle ne l'est plus informatiquement que pour \(\ell>\min\{n,m\}\) une fois les polynômes \(P\) et \(Q\) codés par des structures de listes de tailles respectives \(n\) et \(m\).

On vérifie que l'ensemble des polynômes est ainsi équipé d'une structure d'anneau héritée de \(A\) appelé anneau des polynômes à coefficients dans \(A\) et noté \(A[X]\). On peut à présent définir la fonction polynomiale \(P:A\rg A\) associée à un polynôme \(P\) de coefficients \((a_i)\) et de degré \(n\) par : \begin{equation} x\mapsto \sum_{i=0}^n{\color{orange}a_i}x^i. \end{equation}

Nous pouvons à présent considérer la question suivante : Comment calculer efficacement l'image \(P(x)\) d'un élément \(x\in A\) par une fonction polynomiale associée à un polynôme \(P\) de degré \(n\) ?

L'écriture du polynôme suggère de procéder itérativement en sommant chacun des monômes évalué en \(x,\) ce qui implique de réaliser à chaque fois une exponentiation, ce qui peut être évité comme nous le verrons à la section suivante.

L'algorithme de Horner.

La méthode dite de la factorisation de Horner (qui n'est pas réellement une factorisation au sens usuel) consiste à écrire le polynôme \(P(X)\) sous la forme \({\color{steelblue}P_1(X)}.X+a_0\) et à répéter cette opération pour le polynôme \(P_1(X)\) de degré \(n-1\) et ainsi de suite jusqu'au polynôme constant \({\color{orange}P_{n}(X)}={\color{orange}a_n}\) qui achève le pro­ces­sus : \begin{align*} P(X) &= \big({\color{steelblue}a_nX^{n-1} + a_{n-1}X^{n-2} + \cdots + a_2X+a_1}\big)X+a_0\\ &= \big((a_nX^{n-2} + a_{n-1}X^{n-3} + \cdots + a_2)X+a_1\big)X+a_0\\ &\ \,\vdots\\ &= \Big(\big(\ldots({\color{orange}a_n}X + a_{n-1})X + a_{n-2})X + \cdots + a_2\big)X+a_1\Big)X+a_0 \end{align*}

On peut construire cette suite de polynômes \(P_k(X)\) par la relation de récurrence suivante, avec \(P_0(X):=P(X)\) et \begin{align}\label{eq:rec} \forall k\in\ab{0}{n-1}\quad P_{k+1}(X):=\frac{P_{k}(X)-a_{k}}{X}. \end{align}

Avec la définition récurrente \((\ref{eq:rec})\) des polynômes \(P_k(X)\) ci-dessus, vérifiez que \(P_n(X)=a_n.\)
Vérifions le résultat à l'aide d'une récurrence finie. On considère le prédicat \(H(k)\) défini par \[ P_k(X) = \sum_{i=k}^na_iX^{i-k}. \] On sait que \(P_0(X)=P(X)\) donc \(H(0)\) est vrai. Soit \(k\in\ab{0}{n},\) supposons \(H(k)\). On a \begin{align*} P_{k+1}(X) &=X^{-1}(P_k(X)-a_k)\quad\text{d'après \((\ref{eq:rec})\)}\\ &=X^{-1}\left({\color{olive}\left(\sum_{i=k}^na_iX^{i-k}\right)}-a_k\right)\quad\text{Hypothèse de récurrence}\\ &=X^{-1}\left({\color{olive}\left(\sum_{i=k+1}^na_iX^{i-k}\right)+a_k}-a_k\right)\\ &=X^{-1}\left(X\left(\sum_{i=k+1}^na_iX^{i-(k+1)}\right)\right)\\ &=\sum_{i=k+1}^na_iX^{i-(k+1)}. \end{align*} L'hérédité est donc vérifiée. Ainsi \(H(k)\) est vrai pour tout \(k\in\ab{0}{n}\) et en particulier pour \(k=n,\) c'est-à-dire \(P_n(X)=a_n\).

L'algorithme itératif de Horner consiste à remonter le processus défini par la relation de récurrence \((\ref{eq:rec})\) en partant du polynôme constant \(P_n=a_n\). On évalue alors à chaque étape une fonction affine dont le résultat devient le coefficient du monôme \(X\) pour la suivante.

ALGORITHME Horner(P,x):valeur
DONNEES
   P[0:n]: liste de n + 1 valeurs dans un anneau A
   x: réel
VARIABLES
   R: valeur
   i: entier
DEBUT
   R ← 0
   i ← 0
   TQ (i ≤ n) FAIRE
      R ← R * x + P[n - i]
      i ← i + 1
   FTQ
   RENVOYER R
fin
Démontrez que l'algorithme de Horner s'arrête puis démontrez sa justesse. Indication : considérez le prédicat \(P(i)\) suivant : \begin{equation} \label{eq:sumrec} \text{après}\ i\ \text{itérations de la boucle}\ \ R=\sum_{k=0}^{i-1}a_{n-k}x^{i-k-1}. \end{equation} On rappelle qu'une somme sur un ensemble d'indexation vide est nulle par convention.
Arrêt : la variable \(i\) est initialisée à \(0\) avant la boucle et n'est modifiée qu'à la ligne #12 où elle est incrémentée, on a une suite strictement croissante et comme \(\N\) est archimédien, elle atteindra la valeur \(n+1\) faisant échouer l'entrée dans la boucle. Justesse : le prédicat \(P(i)\) est vrai avant la boucle quand \(i=0,\) en effet \(R=0\). Supposons que la proposition \(P(i)\) soit vraie en entrant dans la boucle, i.e. \[ R=\sum_{k=0}^{i-1}a_{n-k}x^{i-k-1}. \] Alors après l'instruction #12, on a \begin{align*} R&=x(\sum_{k=0}^{i-1}a_{n-k}x^{i-k-1})+a_{n-i}\\ &=\sum_{k=0}^{i}a_{n-k}x^{i-k}. \end{align*} et on a donc \(P(i+1)\). Quand on sort de la boucle, on a \(i=n+1\) et en remplaçant dans \((\ref{eq:sumrec}),\) on obtient le résultat désiré.

Complexité.

Pour les mêmes raisons que pour l'exponentiation, c'est le coût en nombre de multiplications qui va nous intéresser, nous faisons donc l'hypothèse que les valeurs sont bornées ce qui nous évite d'intégrer le coût de ces opérations (le tp sur le calcul multiprécision a pour objet de mettre en évidence que cette hypothèse n'est pas toujours pertinente). Étudions brièvement le cas de l'algorithme naïf ci-dessous au préalable:

ALGORITHME EvalNaive(P,x):valeur
donnees
   P: liste de n + 1 valeurs dans un anneau A
   x: réel
variables
   R: valeur
   i: entier
debut
   R ← P[0]
   i ← 1
   TQ (i ≤ n) FAIRE
      R ← R + P[i] * Exp(x,i)
      i ← i + 1
   FTQ
   renvoyer R
fin

Cet algorithme est basé sur une boucle qui calcule chaque terme \(a_kX^k\) du polynôme \(P(X)\) et l'additionne aux résultats qui précèdent. L'évaluation de ce terme demande exactement \(k\) produits dont \(k-1\) pour l'exponentiation Exp (en supposant qu'elle se fait naïvement elle aussi) et \(1\) pour multiplier le coefficient. On a donc au total:

\(\displaystyle \sum_{i=0}^n i=\frac{n(n+1)}{2}\) multiplications.

Le nombre de multiplications effectuées dans l'algorithme de Horner est \(n+1\) puisqu'il n'y a qu'une seule opération de multiplication dans la boucle et que l'on y passe \(n+1\) fois. Finalement, on passe d'une complexité quadratique en nombre de multiplications avec l'algorithme naïf à une complexité linéaire. L'exercice montre que l'utilisation de l'algorithme square & multiply donnerait une complexité en \(n\log_2 n,\) donc intermédiaire entre les deux.

En supposant que l'algorithme naïf utilise non pas l'exponentiation naïve mais l'algorithme Square & Multiply pour calculer \(X^k,\) quelle serait le nombre de multiplications réalisées ?
Dans ce cas il faudrait évaleur les puissances \(X^k\) pour tous les entiers \(k\in\ab{1}{n},\) le coût de l'exponentiation \(X^k\) étant égal à \(\theta(\log(k)\). Il faut donc évaluer \[ \sum_{k=1}^n\Theta(\log(k)). \] Grâce à une minoration et une majoration par une intégrale de la fonction \(\ln(x)\) dont une primitive est \(x\ln(x)-x,\) on vérifie aisément que le nombre de produits est en \(\Theta(n\log(n))\).
Le mathématicien Victor Yakovlevich Pan démontre en 1966 que cet algorithme est optimal, au sens où il est impossible d'évaluer une fonction polynomiale de degré \(n\) en effectuant moins de \(n\) produits.

Travaux pratiques

On se propose de comparer trois algorithmes pour évaluer une fonction polynomiale \(P\) en \(a\). Chacun des trois algorithmes suivants doit renvoyer non seulement la valeur \(P(a)\) mais également le nombre de multiplications effectuées pour l'obtenir. Le langage est le \(C\).

Écrivez une fonction float Eval-Naif(float *P, uint n, float a) qui calcule la valeur de la fonction polynomiale \(P\) en \(a\) de manière directe en utilisant l'algorithme d'exponentiation naïf.
Écrivez une fonction float Eval-SM(float *P, uint n, float a) qui calcule la valeur de la fonction polynomiale \(P\) en \(a\) de manière directe en utilisant la fonction d'exponentiation vue dans l'algorithme Square & Multiply.
Écrivez une fonction float Eval-Horner(float *P, uint n, float a) qui calcule la valeur de la fonction polynomiale \(P\) en \(a\) en implantant la méthode de Horner décrite plus haut.
Écrivez une fonction float *Creer-Poly(uint n) qui renvoie un polynôme de degré \(n,\) donc une liste de \(n+1\) termes non-nuls. Cette liste servira pour l'évaluation empirique de la complexité des trois algorithme, les coefficients peuvent être tirés au hasard ou fixés arbitrairement du moment qu'aucun d'entre eux n'est égal à 0.
Comparez le nombre de multiplications nécessaires aux trois algorithmes pour calculer \(P(a)\) en fonction du degré \(n\) du polynôme \(P(X)\). Pour cela, tracez les fonctions de coût en multiplications en fonction du degré \(n\) du polynôme \(P\) de ces trois algorithmes sur un graphique. Utilisez la commande gnuplot.

Pour cet exercice, les couples \((n,T(n))\) doivent être rangés sur deux colonnes dans un fichier texte (sans parenthèse, ni virgule, l'espacement faisant office de séparateur entre \(n\) et \(T(n)\)). Pour tracer les trois courbes simultanément dans la même fenêtre graphique, sous l'environnement gnuplot, lancez la commande plot NomFichier1, puis replot NomFichier2 et replot NomFichier3.

Les données des trois fonctions de complexité sont à générer dans des fichiers textes à l'aide d'un programme, pas à la main… De la même manière, écrivez le script gnuplot dans un fichier intitulé comparer-eval.gnu dont la dernière commande sera pause -1 afin de conserver l'affichage des trois courbes à l'écran. Le lancement du script se fait dans un terminal : gnuplot comparer-eval.gnu