Factorisation de Horner

Introduction

Évaluer un polynôme est une situation extrêmement courante en informatique, ce qui explique l’importance que l’on accorde à l’algorithme de William George Horner, somme toute très simple. Il ne faut pas perdre de vue qu’in fine, les algorithmes sont exécutés sur des machines dont les opérations élémentaires se réduisent essentiellement à des combinaisons d’additions et de multiplications. Or, ces deux lois de composition sont précisément celles qui structurent l’expression d’un polynôme.

Le procédé de Hörner, en réécrivant le polynôme sous une forme factorisée et hiérarchisée, met en évidence cette structure de manière particulièrement élégante. Il offre une manière naturelle et systématique d’évaluer un polynôme.

Au-delà de son apparente simplicité, cette méthode illustre parfaitement la manière dont une expression symbolique est transformée en processus calculatoire rigoureux et fiable. Elle met en lumière le lien étroit entre l’algèbre formelle et la pratique algorithmique : comprendre la structure d’un objet mathématique conduit à une mise en œuvre rationnelle.

Enfin, il faut souligner que la portée de ce type d’algorithme dépasse très largement le cadre de la simple évaluation d’un polynôme donné. Dans de nombreux domaines — qu’il s’agisse d’approximer des fonctions analytiques (par des développements limités), de décrire des formes géométriques (courbes de Bézier, splines), ou encore de modéliser des phénomènes physiques ou statistiques —, on ramène les objets étudiés à des expressions polynomiales pour bénéficier de la simplicité et de l’efficacité du calcul qui leur est associé. Le schéma de Hörner s’inscrit ainsi au cœur d’une idée fondamentale : dans la pratique du calcul scientifique comme dans celle de la modélisation numérique, une grande partie des problèmes se ramènent, tôt ou tard, à un calcul polynomial.

Présentation du problème

On se donne un anneau intègre \(A\) 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\geqslant N\then a_n=0. \end{equation} 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\) et les valeurs \((a_i)_{i=0}^k\) sont appelées coefficients du polynôme. Si un seul coefficient \(a_d\) n'est pas nul, le polynôme est appelé monôme de degré \(d\). La structure adaptée pour coder un polynôme de degré \(n\) est donc typiquement une liste de longueur \(n+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 deux lois de composition \(+\) et \(\times\) de l'anneau \(A,\) on munit l'ensemble des polynômes à coefficients dans \(A\) de trois lois, une addition et une multiplication internes \(+\) et \(\times,\) (on conserve les mêmes symboles par commodité) 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)_{i=0}^n\) et \((b_j)_{j=0}^m\) et \(a\in A\) : \begin{align*} \text{Addition interne :}&\quad P(X)+Q(X):=\sum_{\ell\in I_{\max\{n,m\}}}({\color{orange}a_\ell+b_\ell})X^\ell.\\ \text{Produit interne :}&\quad P(X)\times Q(X):=\sum_{\ell\in I_{n+m}} \Big(\sum_{\overset{\scriptstyle(i,j)\in I_n\times I_m}{i+j=\ell}} {\color{orange}a_ib_j}\Big)X^\ell.\\ \text{Produit externe :}&\quad 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.\) On l'appelle anneau des polynômes à coefficients dans \(A\) et on le note \(A[X]\). On peut à présent définir la fonction polynomiale \(P:A\to A\) associée à un polynôme \(P\) de coefficients \((a_i)_{i=0}^n\) 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: valeur
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 ne peut être majorée, 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