Factorisation de Hörner

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. L'anneau est commutatif et n'admet aucun diviseur de \(0\).

Un polynôme à coefficients dans \(A\) est défini comme une suite \(({\color{yellow}a_k})_{k\in\N}\) d'éléments de \(A\) nulle à partir d'un certain rang : \[\exists N\in\N\ \ \forall n\in\N\quad n\geq N\then a_n=0. \] 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 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{yellow}a_0} + {\color{yellow}a_1}X + {\color{yellow}a_2}X² + \cdots + {\color{yellow}a_{n-1}}X^{n-1} + {\color{yellow}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\) : 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_{l=0}^{\max\{n,m\}}({\color{yellow}a_l+b_l})X^l,\\ P(X)\times Q(X)&:=\sum_{l=0}^{nm} \left(\ \sum_{\overset{\scriptstyle i\in\ab{0}{n}\ j\in\ab{0}{m}}{(i+j)=l}} {\color{yellow}a_ib_j}\right)X^l,\\ a\cdot P(X)&:=\sum_{i=0}^{n}(a.{\color{yellow}a_i})X^i. \end{align*}

Attention, la somme \({\color{yellow}a_l+b_l}\) ci-dessus est bien définie pour tout \(l\in\N\) puisque les coefficients de rang supérieur au degré d'un polynôme sont nuls, mais ne l'est plus que pour \(l>\min\{n,m\}\) une fois les polynômes \(P\) et \(Q\) codés par des 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{yellow}a_i}x^i. \end{equation}

On distingue structurellement un polynôme de sa fonction polynomiale car pour les anneaux finis \(A\), il n'y a pas nécessairement bijection entre les polynômes et les fonctions polynomiales.

On veut 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 calculer à chaque fois une exponentielle ce qui peut être évité comme nous le verrons à la section suivante.

L'algorithme.

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{#46F}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{yellow}P_{n}(X)}={\color{yellow}a_n}\) qui achève le pro­ces­sus : \begin{align*} P(X) &= \big({\color{#46F}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{#FF0}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, où le premier terme \(P_0(X)\) de la suite est défini par le polynôme \(P(X)\) : \begin{align}\label{rec} P_k(X):= \begin{cases} P(X),&\text{si}\ k=0.\\ \frac{1}{X}(P_{k-1}(X)-a_{k-1}),&\text{si}\ 1\leq k \leq n. \end{cases} \end{align}

Avec la définition récurrente des polynômes \(P_k\) ci-dessus, vérifiez que \(P_n(X)=a_n\).

L'algorithme itératif de Horner consiste à remonter le processus défini par la relation de récurrence \((\ref{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. 

Horner(P,x):valeur
DONNÉES
   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} 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: 

EvalNaive(P,x):valeur
données
   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
   retourner 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 ?

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 Hörner 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 automatiquement dans des fichiers textes à l'aide d'un programme, pas à la main! De la même manière, on écrira 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 sous un terminal avec la commande gnuplot comparer-eval.gnu.