Algorithmes de calcul formel et numérique

B. Parisse
Institut Fourier
UMR 5582 du CNRS
Université de Grenoble

Giac/Xcas est un logiciel libre de calcul formel dont une caractéristique est de nécessiter peu de ressources sans sacrifier les performances (en particulier sur les calculs polynomiaux). Ce document décrit une partie des algorithmes de calcul formel et numérique qui y sont impleémentés, l’objectif à long terme est de couvrir l’essentiel des algorithmes implémentés. Ce n’est pas le manuel d’utilisation de Xcas, ni un manuel de programmation ou d’exercices illustrés avec Xcas (voir le menu Aide, Manuels : Référence calcul formel, Programmation, Exercices, Amusements...). Ce texte regroupe donc des résultats mathématiques qui ont été ou sont utilisés dans Giac (ou sont susceptibles de l’être), ils sont en général accompagnés de preuves et souvent d’illustrations avec Xcas.
Pour plus d’informations sur Giac/Xcas, cf. :
www-fourier.ujf-grenoble.fr/~parisse/giac_fr.html

N.B.: La version HTML de ce document comporte des champs de saisie interactifs, ceux-ci apparaissent comme des commandes “mortes” dans la version PDF (elles sont exécutées une fois pour toutes par la version non interactive de giac). La version HTML est optimisée pour le navigateur Firefox. Elle est générée avec hevea.inria.fr de Luc Maranget, ou le fork de Yannick Chevallier pour le support mathjax, ainsi qu’une version modifiée de itex2MML de Jacques Distler pour la conversion en MathML. Si vous avez une machine très puissante, vous pouvez exécuter toutes les commandes interactives en cliquant sur le bouton Exécuter. En-dessous de ce bouton se trouve la console de l’interpréteur du logiciel de calcul formel.

Table des matières

Chapitre 1  Plan et index

L’index commence page suivante dans la version PDF.

Quelques conseils de lecture :

Index

Chapitre 2  Trousse de survie Xcas

Cette section peut être vue comme un tutoriel très abrégé pour rapidement prendre en main Xcas par des exemples au niveau fin de licence master de mathématique et préparation aux concours de recrutement d’enseignants. Le lecteur pourra consulter le tutoriel calcul formel (menu Xcas, Aide, Débuter en calcul formel, tutoriel) pour plus de détails ou/et à un niveau mathématique moins élevé.

2.1  Utilisation comme super-calculatrice

2.2  Calcul exact

2.2.1  Arithmétique

2.2.2  Algèbre linéaire exacte

2.3  Calcul scientifique

2.3.1  Analyse numérique

2.3.2  Algèbre linéaire numérique

Chapitre 3  Calculer sur ordinateur

3.1  Représentation des entiers

Proposition 1   Division euclidienne de deux entiers : si aa et bb sont deux entiers, a≥0,b>0a \geq 0, b>0, il existe un unique couple (q,r)(q,r) tel que a=bq+r,r∈[0,b[ a = bq +r , \quad r \in [0, b[

Preuve : On prend pour qq le plus grand entier tel que a−bq≥0a-bq \geq 0.
Exemple :

La division euclidienne permet d’écrire un nombre entier, en utilisant une base bb et des caractères pour représenter les entiers entre 0 et b−1b-1. Nous écrivons les nombres entiers en base b=10b=10 avec comme caractères les chiffres de 0 à 9. Les ordinateurs utilisent des circuits binaires pour stocker les informations, il est donc naturel d’y travailler en base 2 en utilisant comme caractères 0 et 1 ou en base 16 en utilisant comme caractères les chiffres de 0 à 9 et les lettres de A à F. En général, pour trouver l’écriture d’un nombre en base bb (par exemple b=2b=2), on effectue des divisions euclidienne successives par bb du nombre puis de ses quotients successifs jusqu’à ce que le quotient soit 0 et on accolle les restes obtenus (premier reste à droite, dernier reste à gauche). Inversement, pour retrouver un entier dd à partir de son écriture d n...d 0d_n...d_0, on traduit les divisions euclidiennes successives en d = (...((d nb+d n−1)b+d n−2)...+d 1)b+d 0 = d nb n+d n−1b n−1+...+d 0 \begin{matrix} d &=&( ... ((d_n b +d_{n-1})b + d_{n-2})...+d_1)b+d_0\\ &=& d_n b^n + d_{n-1} b^{n-1} + ... + d_0 \end{matrix} Par exemple, vingt-cinq s’écrit en base 16 0x19 car 25 divisé par 16 donne quotient 1, reste 9


En base 2, on trouverait 0b11001 car 25=2 4+2 3+125=2^4+2^3+1.


On peut effectuer les opérations arithmétiques de base (+,-,*, division) directement en base 2 (ou 16). Par exemple la table de l’addition est 0+0=0, 0+1=1+0=1 et 1+1=0 je retiens 1, donc : 

  01001111
+ 01101011
----------
  10111010

Exercice : comment passe-t-on simplement de la représentation d’un nombre en base 2 à un nombre en base 16 et réciproquement ?

Les microprocesseurs peuvent effectuer directement les opérations arithmétiques de base sur les entiers “machine” (déclinés en plusieurs variantes selon la taille et la possibilité d’avoir un signe). Noter que la division de deux entiers aa et bb n’a pas la même signification que la division de deux réels, comme elle ne tomberait pas forcément juste, on calcule le quotient et le reste de la division euclidienne.

Ces entiers machines permettent de représenter de manière exacte des petits entiers relatifs par exemple un entier machine signé sur 4 octets est compris entre [−2 31,2 31−1][-2^{31},2^{31}-1].

Ces entiers machines permettent de faire très rapidement du calcul exact sur les entiers, mais à condition qu’il n’y ait pas de dépassement de capacité, par exemple pour des entiers 32 bits, 2 30+2 30+2 30+2 302^{30}+2^{30}+2^{30}+2^{30} renverra 0. Ils sont utilisables avec tous les langages de programmation traditionnels.

Les logiciels de calcul formel et certains logiciels de programmation permettent de travailler avec des entiers de taille beaucoup plus grande, ainsi qu’avec des rationnels, permettant de faire du calcul exact, mais on paie cette exactitude par un temps de calcul plus long, de plus pas mal de méthodes numériques ne gagnent rien à faire des calculs intermédiaires exacts. Néanmoins, l’utilisation d’un logiciel de calcul formel permettra dans certains cas d’illustrer certains phénomènes dus au calcul approché.

3.2  Les réels

On se ramène d’abord au cas des réels positifs, en machine on garde traditionnellement un bit pour stocker le signe du réel à représenter.

3.2.1  Virgule fixe et flottante.

La première idée qui vient naturellement serait d’utiliser un entier et de déplacer la virgule d’un nombre fixe de position, ce qui revient à mulitplier par une puissance (négative) de la base. Par exemple en base 10 avec un décalage de 4, 1234.5678 serait représenté par 12345678 et 1.2345678 par 12345 (on passe de l’entier au réel par multiplication par 10 −410^{-4}. L’inconvénient d’une telle représentation est qu’on ne peut pas représenter des réels grands ou petits, comme par exemple le nombre d’Avogadro, la constante de Planck, etc.

D’où l’idée de ne pas fixer la position de la virgule, on parle alors de représentation à virgule flottante ou de nombre flottant : on représente un nombre par deux entier, l’un appelé mantisse reprend les chiffres significatifs du réel sans virgule, l’autre l’exposant, donne la position de la virgule. Attention, le séparateur est un point et non une virgule dans la grande majorité des logiciels scientifiques. On sépare traditionnellement la mantisse de l’exposant par la lettre e. Par exemple 1234.5678 peut être représenté par 12345678e-8 (mantisse 12345678, exposant -8) mais aussi par 1234567800e-10.

Naturellement, sur un ordinateur, il y a des limites pour les entiers représentant la mantisse mm et l’exposant ee. Si on écrit les nombres en base bb, la mantisse mm s’écrira avec un nombre nn fixé de chiffres (ou de bits en base 2), donc m∈[0,b n[m \in [0,b^n[. Soit un réel xx représenté par x=mb e,m∈[0,b n[ x=mb^e, \quad m \in [0,b^n[ Si m∈[0,b n−1[m\in [0,b^{n-1}[, alors on peut aussi écrire x=m′b e−1x=m' b^{e-1} avec m′=mb∈[0,b n[m'=mb \in [0,b^n[, quelle écriture faut-il choisir? Intuitivement, on sent qu’il vaut mieux prendre m′m' le plus grand possible, car cela augmente le nombre de chiffres significatifs (alors que des 0 au début de mm ne sont pas significatifs). Ceci est confirmé par le calcul de l’erreur d’arrondi pour représenter un réel. En effet, si xx est un réel non nul, il ne s’écrit pas forcément sous la forme mb emb^e, on doit l’arrondir, par exemple au plus proche réel de la forme mb emb^e. La distance de xx à ce réel est inférieure ou égale à la moitié de la distance entre deux flottants consécutifs, mb emb^e et (m+1)b e(m+1)b^e, donc l’erreur d’arrondi est inférieure ou égale à b e/2b^e/2. Si on divise par x≥mb ex \geq mb^e, on obtient une erreur relative d’arrondi majorée par 1/(2m)1/(2m). On a donc intérêt à prendre mm le plus grand possible pour minimiser cette erreur. Quitte à mulitplier par bb, on peut toujours se ramener (sauf exceptions, cf. ci-dessous), à m∈[b n−1,b n[m \in [b^{n-1},b^n[, on a alors une erreur d’arrondi relative majorée par 12b n−1 \frac{1}{2b^{n-1}}

On appelle flottant normalisé un flottant tel que m∈[b n−1,b n[m \in [b^{n-1},b^n[. Pour écrire un réel sous forme de flottant normalisé, on écrit le réel en base bb, et on déplace la virgule pour avoir exactement nn chiffres non nuls avant la virgule et on arrondit (par exemple au plus proche). L’exposant est égal au décalage effectué. Notez qu’en base 2, un flottant normalisé commence forcément par 1, ce qui permet d’économiser un bit dans le stockage.

Ainsi, l’erreur d’arrondi commise lorsqu’on représente un réel (connu exactement) par un double normalisé est une erreur relative inférieure à de 2 −532^{-53} (b=2b=2 et n=52+1n=52+1 pour les doubles).

Exemples :

Il existe une exception à la possibilité de normaliser les flottants, lorsqu’on atteint la limite inférieure de l’exposant ee. Soit en effet e me_m le plus petit exposant des flottants normalisés et considérons les flottants x=b e m(1+1/b)x=b^{e_m}(1+1/b) et y=b e my=b^{e_m}. Ces flottants sont distincts, mais leur différence n’est plus représentable par un flottant normalisé. Comme on ne souhaite pas représenter x−yx-y par 0, (puisque le test x==yx==y renvoie faux), on introduit les flottants dénormalisés , il s’agit de flottants dont l’exposant est l’exposant minimal représentable sur machine et dont la mantisse appartient à [0,b n−1[[0,b^{n-1}[. Par exemple 0 est représenté par un flottant dénormalisé de mantisse 0 (en fait 0 a deux reprśentation, une de signe positif et une de signe négatif).

Enfin, on utilise traditionnellement une valeur de l’exposant pour représenter les nombres plus grands que le plus grand réel reprśentable sur machine (traditionnellement appelé plus ou moins infini) et les erreurs (par exemple 0./0. ou racine carrée d’un nombre réel négatif, traditionnellement appelé NaN, Not a Number).

Exercice : quels sont les nombres réels représentables exactement en base 10 mais pas en base 2 ? Si on écrit 1/101/10 en base 2 avec 53 bits de précision, puis que l’on arrondit avec 64 bits de précision, ou si on écrit 1/101/10 en base 2 avec 64 bits de précision, obtient-on la même chose ?

Les ordinateurs reprśentent généralement les flottants en base 2 (cf. la section suivante pour plus de précisions), mais cette représentation n’est pas utilisée habituellement par les humains, qui préfèrent compter en base 10. Les ordinateurs effectuent donc la conversion dans les routines d’entrée-sortie. Le format standard utilisé pour saisir ou afficher un nombre flottant dans un logiciel scientifique est composé d’un nombre à virgule flottante utilisant le point comme séparateur décimal (et non la virgule) suivi si nécessaire de la lettre e puis de l’exposant, par exemple 1.23e-5 ou 0.0000123. Dans les logiciels de calcul formel, pour distinguer un entiers représentés par un entier d’un entier représenté par un flottant on écrit l’entier suivi de .0 par exemple 23.0.

Remarque :
Les microprocesseurs ayant un mode BCD peuvent avoir un format de représentation des flottants en base 10, les nombres décimaux comme par exemple 0.3 peuvent être représentés exactement. Certains logiciels, notamment maple, utilisent par défaut des flottants logiciels en base 10 sur des microprocesseurs sans mode BCD, ce qui entraine une baisse de rapidité importante pour les calculs numériques (on peut partiellement améliorer les performances en utilisant evalhf en maple).

3.2.2  Les flottants au format double

Cette section développe les notions de la section précédente pour les flottants machine selon la norme IEEE-754, utilisables dans les langage de programmation usuels, elle peut être omise en première lecture. La représentation d’un double en mémoire se compose de 3 parties : le bit de signe s=±1s=\pm 1 sur 1 bit, la mantisse M∈[0,2 52[M \in [0,2^{52}[ sur 52 bits, et l’exposant e∈[0,2 11[e \in [0, 2^{11}[ sur 11 bits. Pour les nombres “normaux”, l’exposant est en fait compris entre 1 et 2 11−22^{11}-2, le nombre représenté est le rationnel (1+M2 52)2 e+1−2 10 (1+\frac{M}{2^{52}}) 2^{e+1-2^{10}} Pour écrire un nombre sous cette forme, il faut d’abord chercher par quel multiple de 2 il faut le diviser pour obtenir un réel rr dans [1,2[[1,2[, ce qui permet de déterminer l’exposant ee. Ensuite on écrit la représentation en base 2 de r−1∈[0,1[r-1 \in [0,1[. Exemples :

On observe que la représentation en base 2 de 6.4 a du être arrondie (car elle est infinie en base 2) bien qu’elle soit exacte (finie) en base 10. Seuls les entiers et les rationnels dont le dénominateur est une puissance de 2 peuvent être représentés exactement. Ceci entraine des résultats qui peuvent surprendre comme par exemple le fait que 0.5 - 5*0.1 n’est pas nul.

Des représentations spéciales (avec e=0e=0 ou e=2 11−1e=2^{11}-1) ont été introduites pour représenter ±∞\pm \infty (pour les flottants plus grands en valeur absolue que le plus grand flottant représentable), et pour représenter les nombres non nuls plus petits que le plus petit flottant représentable de la manière exposée ci-dessus (on parle de flottants dénormalisés), ainsi que le nombre NaN (Not a Number) lorsqu’une opération a un résultat indéfini (par exemple 0/0).

Remarque : Sur les processeurs compatibles avec les i386, le coprocesseur arithmétique i387 gère en interne des flottants avec 80 bits dont 64 bits de mantisse. Sur les architectures 64 bits (x86 ou AMD), le jeu d’instruction SSE permet de travailler avec des flottants de 128 bits. Le compilateur gcc permet d’utiliser ces flottants longs avec le type long double ou les types __float80 et __float128 en utilisant un drapeau de compilation du type -msse

3.2.3  Opérations sur les flottants

Les opérations arithmétiques de base sur les flottants se font de la manière suivante :

3.2.4  Erreurs

La représentation des nombres réels par des doubles présente des avantages, les opérations arithmétiques sont faites au plus vite par le microprocesseur. Les coprocesseurs arithmétiques (intégrés sur les microprocesseurs de PC) proposent même le calcul des fonctions usuelles (trigonométriques, racine carrée, log et exp) sur le type double et utilisent des formats de représentation interne ayant plus de 64 bits pour les doubles, ce qui permet de limiter les erreurs d’arrondi. Par contre, des erreurs vont être introduites, on parle de calcul approché par opposition au calcul exact sur les rationnels. En effet, la représentation doit d’abord arrondir tout réel qui n’est pas un rationnel dont le dénominateur est une puissance de 2. Ensuite chaque opération va entrainer une propagation de ces erreurs et va y ajouter une erreur d’arrondi sur le résultat. Enfin, l’utilisation du type double peut provoquer un dépassement de capacité (par exemple 100!*100!).

Pour diminuer ces erreurs et les risques de dépassement de capacité, il existe des types flottants multiple précision, qui permettent de travailler avec un nombre fixé à l’avance de décimales et une plage d’exposants plus grande. Les calculs sont plus longs mais les erreurs plus faibles. Attention, il s’agit toujours de calcul approché! De plus, pour des quantités dont la valeur est déterminée de manière expérimentale, la source principale de propagation d’erreurs est la précision des quantités initiales, il ne sert souvent à rien d’utiliser des types flottants multiprécision car les erreurs dus à la représentation (double) sont négligeables devant les erreurs de mesure. Dans ce cas, il est pertinent lorsqu’on évalue f(x)f(x) avec xx mal connu de calculer aussi f′(x)f'(x), en effet : f(x(1+h))=f(x)+xhf′(x)+O(h 2)f(x(1+h))= f(x)+xh f'(x) + O(h^2) l’erreur relative sur f(x)f(x) est donc au premier ordre multipliée par |xf′(x)f(x)||\frac{xf'(x)}{f(x)}| Par exemple, l’erreur relative sur e xe^x est au premier ordre l’erreur relative sur xx multipliée par |x||x|.


3.2.5  Erreur absolue, relative, arrondi propagation des erreurs.

On a vu précédemment que pour représenter un réel, on devait l’arrondir, ce qui introduit une erreur même si le réel est connu exactement (par exemple 1/10). Voyons comment se propagent les erreurs dans les opérations arithmétiques de base : on distingue l’addition, la multiplication et l’inversion. La soustraction se ramène à l’addition car le calcul de l’opposé n’introduit aucune erreur nouvelle. Pour l’addition, si |x−x 0|≤ε 0|x -x_0| \leq \varepsilon_0 et si |y−y 0|≤ε 1|y-y_0| \leq \varepsilon_1 alors par l’inégalité triangulaire (|a+b|≤|a|+|b||a+b|\leq |a|+|b|), on a : |(x+y)−(x 0+y 0)|≤|x−x 0|+|y−y 0|≤ε 0+ε 1 |(x+y)-(x_0+y_0)| \leq |x-x_0| + | y-y_0 | \leq \varepsilon_0 + \varepsilon_1 on dit que les erreurs absolues s’additionnent.

Définition 2   L’erreur absolue est définie comme un majorant de la valeur absolue de la différence entre le nombre réel et son représentant double : |x−x 0|≤ε |x-x_0| \leq \varepsilon

Mais comme il faut représenter x 0+y 0x_0+y_0 en machine, on doit ajouter une erreur d’arrondi, qui est proportionnelle à la valeur absolue de x 0+y 0x_0+y_0 d’où la notion d’erreur relative :

Définition 3   L’erreur relative est égale à l’erreur absolue divisée par la valeur absolue du nombre |x−x 0|≤ε|x 0| |x-x_0| \leq \varepsilon |x_0|

Remarquons au passage que les erreurs de mesure expérimentales sont pratiquement toujours des erreurs relatives.

Donc lorsqu’on effectue une addition (ou une soustraction) de deux réels sur machine, on doit additionner les deux erreurs absolues sur les opérandes et ajouter une erreur d’arrondi (relative de 2 −532^{-53}, à titre d’exercice, on pourra vérifier que cette erreur d’arrondi est majorée par l’erreur absolue de la somme x+yx+y dès l’instant où xx et yy ont eux-même une erreur d’arrondi).

Lorsqu’on effectue une multiplication de deux nombres x,yx,y dont les représentants x 0,y 0x_0,y_0 sont non nuls, on a |xy−x 0y 0x 0y 0|=|xx 0yy 0−1|=|(xx 0−1)(yy 0−1)+(xx 0−1)+(yy 0−1)| \left| \frac{xy-x_0 y_0}{x_0 y_0} \right| = \left| \frac{x}{x_0} \frac{y}{y_0} -1 \right| = \left| (\frac{x}{x_0}-1)(\frac{y}{y_0} -1)+(\frac{x}{x_0}-1)+(\frac{y}{y_0} -1) \right| l’erreur relative est donc la somme des erreurs relatives et du produit des erreurs relatives (on peut souvent négliger le produit devant la somme). Il faut aussi y ajouter une erreur relative d’arrondi de 2 −532^{-53} sur x 0y 0x_0 y_0.

On observe que la multiplication est une opération posant moins de problèmes que l’addition, car on manipule toujours des erreurs relatives, par exemple si l’erreur relative sur deux doubles xx et yy non nuls est de 2 −532^{-53}, alors l’erreur relative sur xyxy sera de 2 −53+2 −53+2 −106+2 −53≈3×2 −53 2^{-53} + 2^{-53} + 2^{-106} + 2^{-53} \approx 3 \times 2^{-53} Lorsque l’erreur relative sur les données est grande devant 2 −532^{-53}, l’erreur relative d’arrondi final est négligeable, on peut alors dire que les erreurs relatives s’additionnent pour un produit (c’est aussi vrai pour un quotient: exercice!). Par contre, si on additionne deux nombres dont le représentant de la somme est proche de 0, la somme des erreurs absolues peut devenir non négligeable par rapport à la somme des représentants, entrainant une erreur relative très grande. Par exemple si xx est représenté par x 0=1+2 −52x_0=1+2^{-52} avec une erreur d’arrondi de 2 −532^{-53} et yy par y 0=−1y_0=-1 avec la même erreur d’arrondi, l’addition de xx et yy renvoie 2 −522^{-52} avec une erreur absolue de 2*2 −532 * 2^{-53} (ici il n’y a pas d’arrondi lorsqu’on fait la somme). C’est une erreur relative de 11 (qui domine largement l’erreur d’arrondi) ce qui signifie que dans la mantisse, seul le premier bit sur les 52 a un sens, la perte de précision est très grande.

Une autre conséquence importante est que l’addition de réels sur machine n’est pas une opération associative, par exemple (2.0 −53+2.0 −53)+1.0→1+2 −52 (2.0^{-53}+2.0^{-53})+1.0 \rightarrow 1+2^{-52} alors que 2.0 −53+(2.0 −53+1.0)→1 2.0^{-53}+(2.0^{-53}+1.0) \rightarrow 1 Dans Xcas, il n’y a que 48 bits de mantisse :

Si on a plusieurs termes à additionner, il faut commencer par additionner entre eux les termes les plus petits, pour que les petits termes ne soient pas absorbés un à un dans les erreurs d’arrondi (les petits ruisseaux font les grands fleuves).

Exercice : pour calculer la valeur numérique d’une dérivée de fonction, il vaut mieux calculer (f(x+h)−f(x−h))/(2h)(f(x+h)-f(x-h))/(2h) que (f(x+h)−f(x))/h(f(x+h)-f(x))/h car le terme d’erreur est en O(h 2)O(h^2) et non en O(h)O(h). Attention toutefois à ne pas prendre hh trop petit, sinon x+h=xx+h=x en flottants et même si x+h≠xx+h \neq x, l’erreur absolue sur f(x+h)−f(x−h)f(x+h)-f(x-h) est (au moins) d’ordre ε|f(x)|\varepsilon |f(x)|, donc l’erreur relative est d’ordre ε/h|f(x)|\varepsilon/h |f(x)|. Par exemple pour h=1e-8 le reste est en O(h 2)O(h^2) donc de l’ordre des erreurs d’arrondi mais l’erreur relative sur f(x+h)−f(x−h)f(x+h)-f(x-h) est d’ordre ε/h\epsilon/h largement supérieure (en flottants double-précision). On choisira plutôt hh tel que ε/h\epsilon/h soit proche de h 2h^2, donc de l’ordre de 1e-5, qui fournira une valeur approchée avec une erreur relative de l’ordre de 1e-10. Exemple : calcul de la dérivée numérique de exp(sin(x))\exp(\sin(x)) en x=1x=1

Remarquons néanmoins que les erreurs calculées ici sont des majorations des erreurs réelles (ou si on préfère l’erreur obtenue dans le pire des cas), statistiquement les erreurs sur les résultats sont moindres, par exemple si on effectue nn calculs susceptibles de provoquer des erreurs indépendantes suivant une même loi d’espérance nulle, la moyenne des erreurs divisée par l’écart-type de la loi tend vers une loi normale centrée réduite. De manière plus déterministe, on a l’inégalité de Bienaymé-Tchebyshev P(|X|>α)≤nσ 2α 2 P(|X|>\alpha) \leq \frac{n\sigma^2}{\alpha^2} où XX est la variable aléatoire somme des nn erreurs, α\alpha l’erreur et nσ 2n\sigma^2 la variance de la somme nn erreurs supposées indépendantes, cette probabilité tend vers 0 pour nn grand si α\alpha est d’ordre nn, et ne tend pas vers 0 si α\alpha est de l’ordre de n\sqrt{n}. Exemple : somme de n=400n=400 nombres répartis sur [−1,1][-1,1] selon la loi uniforme (représentant des erreurs), on divise par n=\sqrt{n}=20, on effectue plusieurs tirages (par exemple 500) on trace l’histogramme et on compare avec la loi normale de moyenne nulle (l’espérance de la somme) et d’écart-type celui de la loi uniforme.

Attention, si on effectue la somme de nn réels ∑ jx j\sum_j x_j, les erreurs d’arrondis ne satisfont pas à ces hypothèses. En effet, l’erreur d’arrondi à chaque opération est une erreur relative, l’erreur absolue correspondante est ε|x 1+x 2|\epsilon |x_1+x_2| puis ε|x 1+x 2+x 3|\epsilon |x_1+x_2+x_3| puis ... ε|x 1+x 2+...+x n|\epsilon |x_1+x_2+...+x_n|, que l’on peut majorer par ε((n−1)|x 1|+(n−2)|x 2|+...+|x n||)\epsilon ((n-1)|x_1|+(n-2)|x_2|+...+|x_n||) La majoration de l’erreur d’arrondi dépend donc de l’ordre des termes, on a intérêt à sommer en commençant par les termes les plus petits en valeur absolue. Mais on peut faire mieux, il est possible de corriger les erreurs d’arrondi dans une somme avec le programme suivant pour une liste (on peut bien sur adapter à la somme d’une expression dépendant d’une variable entière sans stocker de liste) : 

Somme(l):={
  local x,s,c;
  s:=0.0;
  c:=0.0;
  pour x in l faire
    c += (x-((s+x)-s));
    s += x;
  fpour;
  print(c);
  return s+c;
}:;

onload
En effet, cc devrait valoir 0 sans erreurs d’arrondis, avec les erreurs d’arrondis, on a le premier calcul s+xs+x qui donnera une erreur opposée à celui du calcul de ss à la ligne suivante, le 2ième calcul effectué (s+x)−s(s+x)-s donne une erreur absolue en ε|x|\epsilon |x| au pire (car c’est une erreur relative par rapport à (s+x)−s(s+x)-s), et la 3ième erreur d’arrondi est négligeable (puisque la somme vaut 0). On a donc une erreur absolue sur s+cs+c qui est au premier ordre au pire en O(ε∑|x i|)O(\epsilon \sum|x_i|), bien meilleure que la majoration ε((n−1)|x 1|+(n−2)|x 2|+...+|x n||)\epsilon ((n-1)|x_1|+(n-2)|x_2|+...+|x_n||) calculée précédemment.

Par exemple
à comparer avec
(le calcul de SS est fait en exact, celui de sum(1. /j,j,1,n) est approché sans correction).

En conclusion, il est souvent très difficile de calculer une majoration rigoureuse de l’erreur pour des calculs (sauf les plus simples), et cette majoration est en général bien trop pessimiste. Lorsqu’on doute de la précision d’un calcul, un test peu couteux consiste à refaire ce calcul en utilisant des flottants en précision plus grande et tester si le résultat varie en fonction du nombre de chiffres significatifs utilisés, ou faire varier légèrement les données et observer la sensibilité du résultat. Si on veut travailler en toute rigueur sans pour autant calculer les erreurs à priori, il faut utiliser un logiciel utilisant des intervalles pour représenter les réels (section suivante)

3.3  L’arithmétique d’intervalle.

Certains systèmes de calcul formel peuvent manipuler directement des intervalles réels, par exemple par l’intermédiaire de la bibliothèque C MPFI. Les opérations arithmétiques sur des intervalles renvoient alors le meilleur intervalle possible contenant toutes les valeurs possibles lorsque les opérandes parcourent leurs intervalles respectifs. Exemple en Xcas (version 1.1.1 et ultérieures) : [-1..2]*[-1..2] renvoie [-2..4]. Attention ici on parcourt toutes les valeurs possibles de xy,x∈[−1,2],y∈[−1,2]xy, \ x \in [-1,2], y \in [-1,2]. Ce qui est différent du carré d’un intervalle ou plus généralement de l’évaluation d’un polynôme en un intervalle, horner(x^2,[-1..2]) renvoie ainsi [0..4].

Les fonctions disponibles sont souvent moins riches qu’en arithmétique flottante, le calcul d’une fonction non monotone sur un intervalle peut s’avérer délicat, alors que si la fonction est monotone, il suffit de calculer l’image des deux bornes de l’intervalle. Pour les polynômes, Xcas décompose les coefficients en deux parties P=P +−P −P=P_+-P_- en fonction du signe, puis utilise la monotonie de P +P_+ et P −P_- sur ℝ +\mathbb{R}^+ et ℝ −\mathbb{R}^- respectivement.

L’arithmétique d’intervalle dans ℂ\mathbb{C} est beaucoup plus difficile à mettre en oeuvre puisqu’il n’y a plus d’ordre ni de monotonie, on doit alors s’en remettre à des estimations sur les parties réelles et imaginaires qui ne tiendront pas compte du phénomène ci-dessus sur la différence entre xy,x∈[−1,2],y∈[−1,2]xy, \ x \in [-1,2], y \in [-1,2] et x 2,x∈[−1,2]x^2, \ x \in [-1,2].

3.4  Calcul exact et approché, types, évaluation.

Dans les langages de programmation traditionnel (C, Pascal,...), il existe déjà des types permettant une représentation exacte des données (type entier) ou une représentation approchée (type flottant). Mais ces types de donnée de base occupent une taille fixe en mémoire, le type entier est donc limité à un intervalle d’entiers (par exemple [0,2 32−1][0,2^{32}-1] pour un entier non signé sur une machine utilisant un processeur 32 bits) alors que le type flottant peut représenter des nombres réels, mais est limité à une précision en nombre de digits de la mantisse et de l’exposant (par exemple 12 chiffres significatifs et un exposant compris entre -499 et 499).

En calcul formel, on souhaite pouvoir calculer rigoureusement d’une part, et avec des paramètres dont la valeur n’est pas connue d’autre part ; il faut donc s’affranchir de ces limites :

Enfin, il faut pouvoir évaluer un objet (en particulier symbolique) : par exemple évaluer sin(x)\sin(x) lorsqu’on assigne une valeur à xx. Dans cet exemple, on voit qu’il faut d’abord remplacer xx par sa valeur avant de lui appliquer la fonction sinus. C’est le mécanisme général de l’évaluation, mais il y a quelques exceptions où on souhaite empêcher l’évaluation d’un ou plusieurs arguments d’une fonction avant l’évaluation de la fonction. Par exemple si on veut calculer la valeur numérique d’une intégrale par des méthodes de quadrature, on ne souhaitera pas rechercher une primitive de la fonction à intégrer. Dans le jargon, on parle alors de “quoter” un argument (l’origine du terme vient probablement de la notation ' du langage Lisp). Certaines fonctions doivent toujours quoter leurs arguments (par exemple la fonction qui permet de purger le contenu d’un paramètre), on parle parfois d’autoquotation.

3.5  Forme normale et reconnaissance du 0.

Une fois défini ces types de base représentant les nombres d’un système de calcul formel, il faut pouvoir comparer ces nombres, en particulier décider si deux représentations distinctes correspondent au même nombre ou, ce qui revient au même, par soustraction décider quand un nombre est nul. Par exemple 4/24/2 et 2 représentent le même nombre. Lorsqu’on dispose d’un algorithme permettant de représenter un nombre d’une manière unique, on parle de forme normale. C’est par exemple le cas pour les nombres rationnels, la forme normale usuelle est la fraction irréductible de dénominateur positif. C’est aussi le cas pour les fractions rationnelles de polynômes à coefficients entiers représentées par une fraction irréductible, avec au dénominateur un coefficient de plus haut degré positif. Malheureusement, il n’est pas toujours possible de trouver une forme normale pour diverses raisons théoriques ou pratiques :

En résumé, au mieux on a une forme normale, au pire on risque de ne pas reconnaître un zéro, entre les deux on peut ne pas avoir de forme normale mais être capable de reconnaître à coup sûr une expression nulle (par contre, si le système de calcul formel détermine qu’une expression est nulle, alors elle l’est).

Il n’existe pas d’algorithme solution pour le problème de la reconnaissance du zéro pour une classe d’expressions "assez générale". Heureusement, dans la plupart des cas pratiques on sait résoudre ce problème, en se ramenant le plus souvent au cas des polynômes et fractions rationnelles. Par exemple, pour simplifier une expression trigonométrique, on remplace les fonctions trigonométriques sin(x),cos(x),tan(x)\sin(x), \cos(x), \tan(x) par leur expression en fonction de t=tan(x/2)t=\tan(x/2), on est ainsi ramené à une fraction rationnelle en tt que l’on écrit sous forme normale.

Les polynômes ont un rôle central dans tout système de calcul formel puisque sauf dans les cas les plus simples (fractions d’entiers par exemple), la simplification d’expressions fait appel à un moment ou à un autre à des calculs de PGCD de polynômes. Le PGCD de polynômes est un algorithme très sollicité auquel nous consacrerons une section. En effet, l’application brutale de l’algorithme d’Euclide pose des problèmes d’efficacité ce qui a obligé à inventer des méthodes plus efficaces. Anticipons rapidement sur un exemple qui montre l’un des problèmes majeurs des algorithmes de calcul formel, l’explosion en taille (ici des coefficients des restes successifs). Voici donc les restes successifs lorsqu’on applique l’algorithme d’Euclide pour calculer le PGCD de P(x)=(x+1) 7−(x−1) 6P(x)=(x+1)^{7}-(x-1)^{6} avec sa dérivée (les deux polynômes sont premiers entre eux) : 7×(x+1) 6−6×(x−1) 5 16249×x 5+−39049×x 4+106049×x 3+−78049×x 2+47449×x+−7849 157780729×x 4+−5076402187×x 3+290864729×x 2+−101528729×x+28028729 149×(14003282645×x 3+−7328882645×x 2+11333523703×x+−73288818515) 12187×(21618163768324669921×x 2+−5554368469444669921×x+3019170248644669921) 1907235×(469345063045455129411872×x+−47641670106615129411872) 5497465490623352995840209648836272383412129 \begin{matrix} 7\* (x+1)^{6}-6\* (x-1)^{5} & &\\ \frac{162}{49} \* x^{5}+\frac{-390}{49} \* x^{4}+\frac{1060}{49} \* x^{3}+\frac{-780}{49} \* x^{2}+\frac{474}{49} \* x+\frac{-78}{49}& &\\ \frac{157780}{729} \* x^{4}+\frac{-507640}{2187} \* x^{3}+\frac{290864}{729} \* x^{2}+\frac{-101528}{729} \* x+\frac{28028}{729}& &\\ \frac{1}{49} \* (\frac{1400328}{2645} \* x^{3}+\frac{-732888}{2645} \* x^{2}+\frac{1133352}{3703} \* x+\frac{-732888}{18515})& &\\ \frac{1}{2187} \* (\frac{2161816376832}{4669921} \* x^{2}+\frac{-555436846944}{4669921} \* x+\frac{301917024864}{4669921})& &\\ \frac{1}{907235} \* (\frac{469345063045455}{129411872} \* x+\frac{-47641670106615}{129411872})& &\\ \frac{5497465490623352995840}{209648836272383412129} \end{matrix} Le lecteur voulant tester d’autres exemples pourra utiliser le programme Xcas suivant : 

pgcdderiv(a):={
  local b,r,res;
  b:=diff(a,x);
  res:=NULL;
  for (;b!=0;){
    res:=res,b;
    r:=rem(a,b);
    a:=b;
    b:=r;
  }
  return(res);
}



3.6  Valeur générique des variables et hypothèses

Lorsqu’on utilise un symbole sans lui affecter de valeurs en mathématiques on s’attend à une discussion en fonction du paramètre représenté par ce symbole. Ce qui nécessiterait de créer un arborescence de calculs (on retrouve ici les problèmes d’explosion évoqués dans la section précédente). La plupart des systèmes de calcul formel contournent la difficulté en supposant que le paramètre possède une valeur générique (par exemple la solution de (t 2−1)x=t−1(t^2-1)x=t-1 sera x=1/(t+1)x=1/(t+1)) ou choisissent une branche pour les fonctions possédant un point de branchement (par exemple pour résoudre x 2=tx^2=t en fonction de tt). Certains systèmes demandent de manière interactive à l’utilisateur si la variable est par exemple positive ou différente de 1 mais cela s’oppose à un traitement automatique. On peut aussi anticiper ce type de décision en faisant des hypothèses sur une paramètre, la plupart des systèmes de calcul formel actuel proposent cette possibilité.

3.7  Structures de données

On a vu plus haut qu’on souhaitait manipuler des entiers de taille non fixe, des réels de précision fixe ou non, des fractions, des nombres complexes, des extensions algébriques, des paramètres, des expressions symboliques. La plupart des systèmes proposent un type générique qui recouvre ces divers types de scalaire. On peut par exemple utiliser un type structuré comportant un champ type et la donnée ou un pointeur sur la donnée (avec dans ce cas un pointeur sur un compteur de références de la donnée pour pouvoir la détruire dès qu’elle n’est plus référencée1). En programmation orientée objet, on utiliserait plutôt un type abstrait dont dérivent ces différents scalaires et le polymorphisme.

Il faut aussi un type pour les vecteurs, les matrices et les listes. Il faut prendre garde à la méthode utilisée par le système lorsqu’on modifie un élément d’un vecteur, matrice ou liste : soit on effectue une copie de tout l’objet en modifiant l’élément, soit on modifie l’élément de l’objet original. La première méthode (par valeur) est plus aisée à comprendre pour un débutant mais la seconde méthode (par référence) est bien plus efficace.

On peut se poser la question de savoir s’il faut inclure ces types dans le type générique ; en général la réponse est affirmative, une des raisons étant que les interpréteurs qui permettront de lire des données dans un fichier texte sont en général basé sur le couple de logiciels lex(flex)/yacc(bison) qui ne peut compiler qu’à destination d’un seul type. Ceci permet également d’unifier en un seul type symbolique les fonctions ayant un ou plusieurs arguments en voyant plusieurs arguments comme un vecteur d’arguments. Les fonctions sont le plus souvent elle-même incluses dans le type générique permettant ainsi à l’utilisateur de saisir des commandes ou programmes fonctionnels (on peut utiliser une fonction comme argument d’une commande).

Pour des raisons d’efficacité, les systèmes de calcul formel utilisent souvent des représentations particulières pour les polynômes dont on a dit qu’ils jouaient un rôle central. Pour les polynômes à une variable, on peut utiliser la liste des coefficients du polynôme, on parle alors de représentation dense. On peut aussi décider de ne stocker que les coefficients non nuls, on parle alors de représentation creuse (on stocke alors un couple formé par le coefficient et le degré du monôme correspondant). Pour les polynômes à plusieurs variables, on peut les considérer comme des polynômes à une variable à coefficients polynomiaux, on parle alors de représentation récursive. On peut aussi décider de ne pas briser la symétrie entre les variables (pas de variable principale), on parle alors de représentation distribuée, le plus souvent les représentation distribuées sont creuses car les représentations denses nécessitent très vite beaucoup de coefficients. Les méthodes de représentation creuses sont parfois aussi utilisées pour les matrices ayant beaucoup de coefficients nuls.

Voyons maintenant plus précisément sur quelques exemples de logiciels de calcul formel répandus quelles structures de données sont utilisées. Plusieurs éléments entrent en compte dans les choix faits :

Voyons quelques exemples, d’abord Giac, puis des systèmes pour ordinateur où les ressources (par exemple mémoire) sont moins limitées ce qui permet d’utiliser des langages de programmation de plus haut niveau. On termine par les calculatrices formelles HP et TI des années 20002. Ce sont des systèmes plutôt destinés à l’enseignement, soumis à de fortes contraintes en termes de taille mémoire, et destinés à traiter des petits problèmes.

3.7.1  Maple, Mathematica, ...

Ces systèmes ont un noyau fermé, au sens où l’utilisateur n’a pas accès du tout, ou en tout cas pas facilement, aux structures de données de base. Je ne dispose donc pas d’information sur les structures de données utilisées par le noyau.

L’interaction système-utilisateur se fait quasiment toujours en utilisant le langage de programmation propre au système, langage interprété par le noyau du système (ce qui ralentit l’exécution). Ces langages utilisateurs sont essentiellement non typés : on travaille avec des variables du type générique sans pouvoir accéder aux types sous-jacents. On ne bénéficie en général pas des vérifications faites lors de la compilation avec un langage typé, de plus ces systèmes ne sont pas toujours fourni avec de bon outils de mise au point. Enfin ces langages ne sont pas standardisés d’un système à l’autre et il est en général impossible d’utiliser ces systèmes comme des librairies depuis un langage de programmation traditionnel. Leur intérêt principal réside donc dans une utilisation interactive en profitant de la librairie de fonctions accessibles.

3.7.2  Giac/Xcas

Il s’agit du système de calcul formel que j’implémente actuellement sous forme d’une bibliothèque C++ (ce qui permettra aux programmes tiers d’utiliser beaucoup plus facilement du calcul formel qu’avec les systèmes précédents). L’objectif est d’avoir un système facile à programmer directement en C++, proche du langage utilisateur, lui-même compatible avec Maple ou MuPAD, tout cela sans trop perdre en performances comparativement aux librairies spécialisées écrites en C/C++. Ce qui explique un choix de type générique (gen) non orienté objet, avec un champ type et soit une donnée immédiate (pour les nombres flottants par exemple), soit un pointeur vers un objet du type correspondant au champ type pour les données de taille non fixe (on pourrait donc se contenter du langage C, mais le langage C++ permet de redéfinir les opérateurs sur des types utilisateurs ce qui améliore considérablement la lisibilité du code source). Les données dynamiques ne sont pas dupliquées, Giac utilise un pointeur sur un compteur de référence pour détruire ces données lorsqu’elles ne sont plus référencées.

Les entiers en précision arbitraire sont hérités de la bibliothque GMP (écrite en C) du projet GNU. Les flottants en précision arbitraire utiliseront aussi GMP (plus précisément MPFR). Il y a un type fraction, structure C composé d’un champ numérateur et d’un champ dénominateur, et un type nombre complexe.

Les listes, vecteurs, matrices utilisent le type paramétré vector<> de la librairie standard C++ (Standard Template Library). Les objets symboliques sont des structures composés d’un champ sommet qui est une fonction prenant un argument de type gen et renvoyant un résultat de type gen, et d’un champ feuille qui est de type gen. Lorsqu’une fonction possède plusieurs arguments, ils sont rassemblés en une liste formant le champ feuille de l’objet symbolique. Les programmes sont aussi des objets symboliques, dont le champ sommet est la fonction évaluation d’un programme. Les listes sont aussi utilisées pour représenter vecteurs, matrices et polynômes en une variable en représentation dense, on peut y accéder par valeur (:=) ou par référence (=<). Ces polynômes servent eux-mêmes á représenter des éléments d’une extension algébrique de ℚ\mathbb{Q} (vus comme un couple de polynômes P,QP,Q, où QQ est un polynome minimal irréductible à coefficients entiers, autrement dit P,QP,Q vaut P(α)P(\alpha) où Q(α)=0Q(\alpha)=0), ou des éléments d’un corps fini (comme ci-dessus, mais ici QQ est à coefficients dans ℤ/pℤ\mathbb{Z}/p\mathbb{Z} avec pp premier, cf. la commande GF). Giac posséde aussi un type pour les polynômes en représentation creuse distribuée en plusieurs indéterminées (cf. les commandes symb2poly et poly2symb).

L’évaluation d’un objet symbolique se fait en regardant d’abord si la fonction au sommet doit évaluer ou non ses arguments (autoquote), on évalue les arguments si nécessaire puis on applique la fonction.

Une hypthèse sur un paramètre est une valeur spéciale affectée au paramètre, valeur ignorée par la routine d’évaluation.

3.7.3  Calculatrices formelles HP48/49

Les langages utilisés pour programmer ces calculateurs sont l’assembleur et le RPL (Reverse Polish Lisp) adapté à l’écriture de code en mémoire morte très compact.

Le type générique est implémenté avec un champ type appelé prologue (qui est en fait un pointeur sur la fonction chargée d’évaluer ce type d’objet) suivi de la donnée elle-même (et non d’un pointeur sur la donnée, on économise ainsi la place mémoire du compteur de référence).

Le type entier en précision arbitraire est codé par le nombre de digits (sur 5 quartets3) suivi du signe sur un quartet et de la représentation BCD (en base 10) de la valeur absolue de l’entier. Le choix de la représentation BCD a été fait pour optimiser les temps de conversion en chaîne de caractères pour l’affichage. La mémoire vive disponible est de 256K, c’est elle qui limite la taille des entiers et non le champ longueur de l’entier. Il n’y a pas de type spécifique pour les rationnels (on utilise un objet symbolique normal).

Les fonctions internes des HP49/50/40 utilisent le type programme pour représenter les entiers de Gauß (complexes dont la partie réelle et imaginaire est entière). Les nombres algébriques ne sont pas implémentés, sauf les racines carrées (représentée de manière interne par le type programme). Il y a un type spécifique prévu pour les flottants en précision arbitraire, mais l’implémentation des opérations sur ces types n’a pas été intégrée en ROM à ce jour.

Les types listes, programmes et objet symbolique sont composés du prologue (champ type) suivi par la succession d’objets situés en mémoire vive ou de pointeurs sur des objets situés en mémoire en lecture seule (ROM) et se terminent par un pointeur sur une adresse fixe (appelée SEMI). Ces types sont eux-mêmes des objets et peuvent donc être utilisés de manière récursive. La longueur des types listes, programmes, symboliques n’est stockée nulle part, c’est le délimiteur final qui permet de la connaître, ce qui est parfois source d’inefficacité. On utilise de manière interne les listes pour représenter les polynômes denses (avec représentation récursive pour les polynômes à plusieurs variables).

Les calculatrices HP4xG utilisent une pile4, c’est-à-dire une liste de taille non fixée d’objets. On place les objets sur la pile, l’exécution d’une fonction prend ces arguments sur la pile et renvoie un ou plusieurs résultats sur la pile (ce qui est une souplesse du RPN comparé aux langages où on ne peut renvoyer qu’une valeur de retour). Il faut donc donner les arguments avant d’appeler la fonction correspondante. Par exemple pour calculer a+ba+b on tapera a b +. C’est la syntaxe dite polonaise inversée (RPN). Un avantage de cette syntaxe est que le codage d’un objet symbolique par cette syntaxe est évidente, il suffit de stocker la liste précédente {a b +}. Les objets symboliques sont donc représenté par une suite d’objets écrit en syntaxe polonaise inversée. L’évaluation d’un objet symbolique se fait dans l’ordre polonaise inversé : les arguments sont évalués puis les fonctions leur sont appliqués. Pour des raisons d’efficacité, on représente souvent les objets composites (listes, symboliques) par leurs composants placés sur la pile (appelé meta-objets).

Une rigidité de la syntaxe polonaise est que les fonctions ont toujours un nombre fixe d’arguments5, par exemple l’addition a toujours 2 arguments, ainsi a+b+ca+b+c est obtenu par (a+b)+c(a+b)+c ou par a+(b+c)a+(b+c) c’est-à-dire respectivement a b + c + ou a b c + + ce qui brise parfois artificiellement la symétrie de certaines opérations. En polonaise inversée, le système doit de plus jongler avec l’autoquote puisque les arguments sont évalués avant l’opérateur qui éventuellement demanderait à ne pas évaluer ses arguments. À noter l’existence d’une commande QUOTE permettant à l’utilisateur de quoter une sous-expression.

Les hypothèses sur des variables réelles sont regroupées dans une liste stockée dans la variable globale REALASSUME, on peut supposer qu’une variable est dans un intervalle. Il n’y a pas à ce jour de possibilité de supposer qu’une variable est entière (ni à fortiori qu’une variable à une valeur modulo un entier fixé), bien qu’il ait été décidé de réserver la variable globale INTEGERASSUME à cet effet. Il n’y a pas de possibilité de faire des hypothèses ayant une portée locale.

3.7.4  Calculatrices formelles TI92/89/Voyage 200

Le langage utilisé pour programmer ces calculatrices est le langage C (on peut aussi écrire du code en assembleur pour ces calculatrices). On retrouve ici les différents types de données regroupé en un type générique qui est un tableau d’octets (aussi appelé quantum). Le champ type est appelé tag dans la documentation TI. Contrairement à ce qui précède, ce champ type est placé en mémoire à la fin de l’objet, ce qui est possible car la longueur d’un objet est toujours indiquée au début de l’objet. Ceci est fait afin de faciliter l’évaluation (cf. infra).

Les entiers en précision arbitraire sont codés par un tag parmi deux (pour différencier le signe), un octet pour la longueur, puis la valeur absolue de l’entier (en base 256). Ils sont donc limités par le champ longueur à 255 octets, le plus grand entier représentable est 6 (256 255−1)(256^{255}-1). Il existe un tag spécifique pour les rationnels, pour les constantes réelles et entières qui apparaissent par exemple en résolvant une équation. Il existe des tags utilisés de manière interne, par exemple pour les nombres complexes. Il n’y a pas de tag prévu pour les flottants en précision arbitraire. ni pour les nombres algébriques (racines carrées par exemple).

Les listes sont codées par la succession de leurs éléments. En principe elles ne peuvent pas contenir des listes (sauf pour représenter une matrice). Quelques fonctions utilisent les listes pour représenter des polynômes denses à une variable, mais probablement pas pour représenter de manière récursive des polynômes à plusieurs variables (puisque le type liste n’est en principe pas récursif).

Comme les HP, les TI utilisent une pile (non visible par l’utilisateur) appelée expression stack afin de traduire un expression mathématique sous forme d’un texte en un objet symbolique codé exactement comme ci-dessus en syntaxe polonaise. Toutefois, la présence du champ longueur permet d’évaluer un objet symbolique sans perdre en efficacité en partant de l’opérateur final et en redescendant ensuite sur ces arguments, c’est la stratégie adoptée. C’est pour cela que le tag d’identification se trouve à la fin de l’objet. L’utilisation de cette méthode facilite grandement l’autoquotation (on peut toutefois regretter que le système n’ait pas prévu d’instruction permettant à l’utilisateur d’empêcher l’évaluation d’une sous-expression).

On ne peut pas faire d’hypothèse globale sur un paramètre par contre on peut faire des hypothèses de type appartenance à un intervalle ayant une portée locale.

3.8  Algorithmes et complexité.

On va présenter dans la suite quelques algorithmes que l’on peut considérer comme classiques dans le domaine du calcul formel. Avant d’implémenter ce type d’algorithmes, on a besoin des algorithmes de base en arithmétique.

La plupart des problèmes posés en calcul formel nécessitent des calculs dont la taille croit de manière exponentielle voire doublement exponentielle en fonction de la taille des données et ce même si le résultat est lui aussi de taille petite. Un exemple est la réduction des systèmes de plusieurs équations polynomiales (bases de Groebner).

3.8.1  Algorithmes modulaires ou p-adiques

Dans certains cas, l’application de théories mathématiques parfois sophistiquées permet de réduire la complexité (par exemple, M. Van Hoeij a découvert récemment qu’un algorithme très utilisé en théorie des nombres, l’algorithme LLL, permettait d’améliorer la complexité d’une des étapes de la factorisation des polynomes à coefficients entiers sur les entiers). Heureusement, dans de nombreux cas, on peut réduire la complexité (donc le temps de calcul) par des adaptations au problème d’une même idée à condition de faire des hypothèses sur les données (autrement dit en abandonnant la volonté d’implémenter un algorithme très générique, ou tout au moins en spécialisant des algorithmes génériques). Par exemple lorsqu’on travaille avec des entiers (ou des polynômes à coefficients entiers, ou des matrices à coefficients entiers...) on utilise souvent des algorithmes modulaires et pp-adiques. Comme le calcul exact nécessite presque toujours de calculer avec des entiers, ces méthodes ont un rôle central en calcul formel, nous les présentons donc maintenant brièvement. Dans les prochaines sections, nous utiliserons ce type de méthode, par exemple pour le calcul de PGCD ou la factorisation de polynômes à coefficients entiers.

Les méthodes modulaires consistent à réduire un problème dans ℤ\mathbb{Z} à son équivalent dans ℤ/nℤ\mathbb{Z}/n\mathbb{Z} pour une ou plusieurs valeurs de nn, nombre premier. Le calcul dans ℤ/nℤ\mathbb{Z}/n\mathbb{Z} a l’avantage de se faire avec des entiers dont la taille est bornée. Ensuite à l’aide d’estimations à priori sur la taille des solutions éventuelles du problème initial, on reconstruit la solution au problème initial avec le théorème des restes chinois.

Par exemple, on peut calculer un déterminant d’une matrice à coefficients entiers en cherchant ce déterminant dans ℤ/nℤ\mathbb{Z}/n\mathbb{Z} pour plusieurs nombres premiers nn, dont le produit est deux fois plus grand qu’une estimation à priori de la taille du déterminant (donnée par exemple par l’inégalité d’Hadamard, cf. Cohen, p. 50).

Les méthodes pp-adiques commencent de manière identique par un calcul dans ℤ/nℤ\mathbb{Z}/n\mathbb{Z}, on augmente ensuite la précision de la solution en la «liftant»de ℤ/n kℤ\mathbb{Z}/n^k \mathbb{Z} vers ℤ/n k+1ℤ\mathbb{Z}/n^{k+1}\mathbb{Z} ou vers ℤ/n 2kℤ\mathbb{Z}/n^{2k}\mathbb{Z} (lift linéaire ou lift quadratique), on s’arrête lorsque kk est assez grand (à l’aide d’estimations à priori) et on reconstruit alors la solution initiale. L’étape de «lift»est en général un lemme de Hensel dont on verra quelques exemples dans les prochains articles. L’algorithme commun au lemme de Hensel et au théorème des restes chinois est l’identité de Bézout, que l’on retrouve d’ailleurs un peu partout (par exemple pour le calcul de primitives).

Illustrons cette méthode sur un exemple simple, la recherche de racines rationnelles d’un polynôme P(X)=a dX d+⋯+a 0P(X)=a_d X^d + \cdots + a_0 à coefficients entiers ou polynomiaux, avec a da_d et a 0a_0 non nuls. L’algorithme générique (assez connu) consiste à chercher les diviseurs de a 0a_0 et de a da_d et à tester toutes les fractions de ces diviseurs, on montre en effet aisément que si X=p/qX=p/q fraction irréductible est racine de PP alors qq divise a da_d et pp divise a 0a_0. Cet algorithme est très inefficace si a da_d ou a 0a_0 est un grand entier (car on ne sait pas forcément le factoriser) ou s’il a beaucoup de facteurs premiers (la liste des diviseurs à tester est alors très grande).

Lorsque les coefficients de PP sont entiers, la recherche précédente revient à trouver un facteur à coefficients entiers qX−pqX-p de PP, on peut donc réduire le problème modulo un entier premier nn qui ne divise pas a da_d : si un tel facteur existe dans ℤ\mathbb{Z} alors ce facteur (réduit modulo nn) est un facteur de PP dans ℤ/nℤ\mathbb{Z}/n\mathbb{Z} donc PP admet une racine dans ℤ/nℤ\mathbb{Z}/n\mathbb{Z} (puisque qq est inversible modulo nn car on a choisi nn premier ne divisant pas a da_d). On évalue maintenant PP en les nn éléments de ℤ/nℤ\mathbb{Z}/n\mathbb{Z}. S’il n’y a pas de 0, alors PP n’a pas de racine rationnelle. S’il y a des racines, on va les lifter de ℤ/n kℤ\mathbb{Z}/n^k\mathbb{Z} dans ℤ/n 2kℤ\mathbb{Z}/n^{2k}\mathbb{Z}.

On suppose donc que pour k≥1k\geq 1, il existe un entier p kp_k tel que P(p k)=0(modn k) P(p_k)=0 \pmod{n^k} Il s’agit de trouver un entier xx tel que p k+1=p k+n k×xp_{k+1}=p_k+n^k \* x vérifie P(p k+1)=0(modn 2k) P(p_{k+1})=0 \pmod{n^{2k}} On applique la formule de Taylor à l’ordre 1 pour PP en p kp_k, le reste est nul modulo n 2kn^{2k}, donc : P(p k)+n k×xP′(p k)=0(modn 2k) P(p_k)+ n^k \* x P'(p_k)=0 \pmod{n^{2k}} soit finalement : x=−P(p k)n k×(P′(p k)(modn k)) −1 x=-\frac{P(p_k)}{n^k} \* ( P'(p_k) \pmod{n^k}) ^{-1} On reconnaît au passage la méthode de Newton, pour qu’elle fonctionne il suffit que P′(p k)≠0(modn)P'(p_k) \neq 0 \pmod n ce qui permet de l’inverser modulo n kn^k (et c’est ici qu’intervient l’identité de Bézout). En pratique quand on factorise un polynôme, on commence par retirer les multiplicités, on peut donc supposer que PP est sans facteur multiple dans ℤ\mathbb{Z}. Ceci n’entraîne pas forcément qu’il le reste dans ℤ/nℤ\mathbb{Z}/n\mathbb{Z} ce qui crée une contrainte supplémentaire sur le choix de nn, à savoir que PP et P′P' restent premier entre eux dans ℤ/nℤ\mathbb{Z}/n\mathbb{Z} (il existe forcément de tels nn, par exemple nn premier plus grand que le plus grand entier intervenant dans le calcul du PGCD de PP et P′P' dans ℤ\mathbb{Z}).

Reste donc à revenir dans ℤ\mathbb{Z} à partir d’une racine p kp_k dans ℤ/(n kℤ)\mathbb{Z}/(n^k \mathbb{Z}) (où on peut choisir kk). On va maintenant utiliser la représentation modulaire symétrique : on prend comme représentant modulaire d’un entier zz dans ℤ/n kℤ\mathbb{Z}/n^k\mathbb{Z} l’unique entier congru à zz modulo nn qui est strictement compris entre −n k/2-n^k/2 et n k/2n^k/2 (si nn est pair, la deuxième inégalité est choisie large).

Si qX−pqX-p est un facteur de PP, alors a dX−a dqpa_dX-\frac{a_d}{q}p est encore un facteur de PP (le quotient de PP par a dX−a dqpa_dX-\frac{a_d}{q}p est à coefficients rationnels mais le facteur est à coefficients entiers). Si on a choisi kk tel que n k>2|a da 0|n^k&gt;2|a_d a_0|, l’écriture en représentation modulaire symétrique de a dX−a dqpa_dX-\frac{a_d}{q}p est inchangée, en effet on a des estimations à priori sur les entiers pp et qq : |q|≤|a d||q|\leq |a_d| et |p|≤|a 0||p| \leq |a_0| puisque qq divise a da_d et pp divise a 0a_0. Comme a dX−a dqpa_dX-\frac{a_d}{q}p est égal à a d(X−p k)a_d(X-p_k) dans ℤ/(n kℤ)\mathbb{Z}/(n^k \mathbb{Z}), il nous suffit d’écrire en représentation modulaire symétrique a d(X−p k)=a dX−p′a_d(X-p_k)=a_d X-p'. Pour conclure, on sait que a dX−p′a_d X-p' est un multiple entier de qX−pqX-p. On divise donc le facteur a dX−p′a_d X-p' par le pgcd de a da_d et p′p' et on teste la divisibilité de PP par ce facteur réduit.

Exemple
Considérons le polynôme 2X 3−X 2−X−32 X^3-X^2-X-3 qui est sans facteur carré. On ne peut pas choisir n=2n=2 car on réduirait le degré, pour n=3n=3, on a P′=X−1P'=X-1 qui est facteur de PP, pour n=5n=5, P′=6X 2−2X−1P'=6X^2-2X-1, on vérifie que PP et P′P' sont premiers entre eux (par exemple avec GCDMOD sur une HP49 où on aura fixé la variable MODULO à 5).

On teste ensuite les entiers de -2 à 2 sur PP. Seul -1 est racine modulo 5 (P(−1)=−5P(-1)=-5), on va maintenant lifter p 1=−1p_1=-1.

L’estimation à priori est 2|a d||a 0|=122|a_d||a_0|=12 donc k=2k=2 (5 2=25>125^2=25&gt;12), une itération suffira. On a P′(−1)=7P'(-1)=7, l’inverse de P′(−1)(mod5)P'(-1) \pmod 5 est -2 donc: x=−P(−1)5(−2)=−(−1)×(−2)=−2 x= -\frac{P(-1)}{5} (-2) = -(-1) \* (-2)=-2 et p 2=−1+5×(−2)=−11p_2=-1+5\times(-2)=-11 est racine de PP dans ℤ/25ℤ\mathbb{Z}/25\mathbb{Z}. On calcule ensuite a d(X−p k)=2(X+11)=2X+22=2X−3a_d(X-p_k)=2(X+11)=2X+22=2X-3 en représentation symétrique, le PGCD de 2 et -3 est 1 donc on teste le facteur 2X−32X-3, ici il divise PP donc PP admet un unique facteur entier de degré 1 qui est 2X−32X-3.

3.8.2  Algorithmes déterministes. Algorithmes probabilistes: Las Vegas et Monte-Carlo

L’algorithme p-adique présenté ci-dessus est un algorithme déterministe, il renvoie toujours un résultat certifié et le temps de calcul nécessaire à son exécution ne dépend pas du hasard (sauf si on choisit le nombre premier pp au hasard...). Ce type d’algorithmes est parfois trop long par rapport à d’autres type d’algorithmes utilisant le hasard :

Dans Xcas, certains algorithmes sont de type Monte-Carlo par défaut, notamment le calcul de déterminant de grandes matrices à coefficients entiers ou de bases de Gröbner, et un warning s’affiche alors. La variable proba_epsilon permet de régler le niveau de probabilité d’erreur acceptée, on peut la mettre à 0 pour forcer l’utilisation d’algorithmes déterministes ou de type Las Vegas avec certification du résultat. Si l’on fait des calculs à but expérimental pour établir une conjecture, il n’est pas nécessaire de certifier un calcul et il ne sert à rien de mettre proba_epsilon à 0. Par contre, pour établir une preuve (au sens mathématique du terme) qui nécessite un calcul fait sur machine, on prendra soin de mettre proba_epsilon à 0. On remarquera au passage que ce type de preuve ne peut se faire qu’avec un logiciel open-source, puisqu’il faut aussi pouvoir montrer que l’algorithme utilisé est correctement implémenté.

3.9  Entiers et polynômes.

L’ensemble des entiers relatifs ℤ\mathbb{Z} est muni des opérations + et * qui en font un anneau commutatif intègre.

Si aa et bb sont deux entiers naturels, b≠0b \neq 0, on définit le quotient euclidien de aa par bb comme le plus grand entier naturel qq tel que a−bq≥0a-bq \geq 0 (l’ensemble des entiers naturels tels que a−bq≥0a-bq \geq 0 est borné par aa donc admet un élément maximal), On définit rr le reste euclidien de aa par bb est alors a−bqa-bq. On vérifie que r∈[0,b−1]r \in [0,b-1] (sinon a−b(q+1)a-b(q+1) serait positif ou nul).

3.9.1  Petits et grands entiers

Les microprocesseurs sont capables d’effectuer des opérations arithmétiques de base (toujours +,-, pratiquement toujours * et la division euclidienne) sur des petits entiers compris dans l’intervalle [−2 t−1,2 t−1−1][-2^{t-1},2^{t-1}-1] (entier signé) ou [0,2 t−1][0,2^{t}-1] (entier non signé), pour t=8t=8, t=16t=16, t=32t=32 et t=64t=64 (parfois aussi pour t=128t=128). En cas de dépassement, les calculs sont faits modulo 2 t2^t. Les logiciels de calcul doivent pouvoir travailler avec des entiers de taille plus grande, il faut les stocker en mémoire et implémenter les opérations arithmétiques de base en se ramenant à des petits entiers. Pour faire cela, on va écrire des entiers dans des bases adéquates.

Ecriture en base bb

On se fixe une fois pour toutes un entier b>1b&gt;1. Soit N∈ℕN \in \mathbb{N}, on peut alors écrire NN de manière unique sous la forme N=∑ i=0 nN ib i,b i∈[0,b−1],b n≠0(1) N=\sum_{i=0}^n N_i b^i , \quad b_i \in [0,b-1], b_n \neq 0 \qquad (1) Les N iN_i sont appelés bits si b=2b=2 ou digits si b=10b=10 (parfois aussi pour b>2b&gt;2, par exemple b=16b=16). Pour montrer l’unicité de l’écriture, on suppose qu’il y a deux écritures distinctes de ce type, et on regarde le coefficient d’indice maximal qui est différent (appelons le jj), on a alors (N j−N˜ j)b j=∑ i=0 j−1(N i−N˜ i)b i(N_j-\tilde{N}_j)b^j = \sum_{i=0}^{j-1} (N_i-\tilde{N}_i)b^i quitte à changer de signe, on peut supposer que le membre de gauche est strictement positif, on a alors (N j−N˜ j)b j≥b j(N_j-\tilde{N}_j)b^j \geq b^j Mais le membre de droite se majore par ∑ i=0 j−1(b−1)b j=b j−1\sum_{i=0}^{j-1} (b-1) b^j=b^j-1 absurde.

L’existence se prouve par l’algorithme suivant qui utilise la division euclidienne. 

#
def ecriture_base(N,b):
    L=[]
    while N>0:
        L.append(N % b)
        N=N//b
    return L

onload

i.e. tant que N>0N&gt;0, on calcule le reste de NN par bb, qui donne le coefficient de poids le plus faible de l’écriture de NN en base bb, on ajoute lécriture en base bb du quotient de NN par bb. L’algorithme s’arrête au bout de partie entière de log b(N)+1\log_b(N)+1 itérations. Réciproquement, on vérifie que l’écriture obtenue convient en développant b 0+b(b 1+b(...+b(b n−1+b(b n))...))(2) b_0+b(b_{1}+b(...+b(b_{n-1}+b(b_n))...)) \qquad (2) On observe au passage que l’écriture de NN sous la forme ci-dessus nécessite nn additions et nn multiplications, et est donc plus efficace que le calcul sous la forme développée 1. C’est la méthode dite de Hörner.

Exemple : en base b=2b=2. Pour écrire N=12N=12 en base 2, on calcule le reste de 12 par 2 donc N 0=0N_0=0, le quotient de 12 par 2 vaut 6, on divise par 2, reste 0 donc N 1=0N_1=0, quotient 3, on divise par 2, reste N 2=1N_2=1, quotient 1, on divise par 2, reste N 3=1N_3=1 quotient 0 on s’arrête, donc 12=0b1100. Réciproquement on a bien 0+2×(0+2×(1+2×(1)))=120 +2\times(0+2\times(1+2\times(1)))=12

Exercice : la commande b:=convert(N,base,b) de Xcas effectue la conversion d’un entier NN en la liste des coefficients de son écriture en base bb, la réciproque étant convert(L,base,b) ou horner(L,b). Implémenter des fonctions équivalentes dans votre langage de programmation préféré.

Exercice : comment passe-t-on simplement de la représentation d’un nombre en base 2 à un nombre en base 16 et réciproquement ?

Lien avec les polynômes de A[X]A[X] avec A=ℤA=\mathbb{Z}.

A tout entier naturel NN, on peut associer un polynôme à coefficients entiers qui est son écriture en base bb, les coefficients du polynôme sont dans l’intervalle [0,b−1][0,b-1]. Réciproquement, si les coefficients d’un polynôme sont des entiers compris entre [0,b−1][0,b-1] alors ils correspondent à l’écriture d’un entier en base bb.

On va voir que les opérations arithmétiques de base sur les grands entiers reviennent à effectuer des opérations de base sur les polynômes, avec une petite difficulté supplémentaire, il faut tenir compte de retenues. Les algorithmes naifs pour additionner et multiplier deux polynômes correspondent précisément aux algorithmes enseignés à l’école primaire pour effectuer l’addition ou la multiplication de deux entiers.

3.9.2  Opérations arithmétiques de base

Addition, soustraction

Si N=∑ i=0 nN ib iN=\sum_{i=0}^n N_i b^i et M=∑ i=0 mM ib iM=\sum_{i=0}^m M_i b^i, alors N+M=∑ i=0 max(n,m)(N i+M i)b i,N+M=\sum_{i=0}^{\mbox{max}(n,m)} (N_i+M_i) b^i, il faut donc additionner les digits de l’écriture en base bb comme pour additionner deux polynômes. On a N i+M i∈[0,2b−1]N_i+M_i \in [0,2b-1] , si N i+M i<bN_i+M_i&lt;b, il n’y a rien à faire, sinon on remplace par N i+M i−bN_i+M_i-b et on ajoute 1 (retenue) à N i+1+M i+1N_{i+1}+M_{i+1} qui appartient à [0,2b−2+1][0,2b-2+1], etc. Le nombre d’opérations à effectuer est de l’ordre du maximum du nombre de bits/chiffres/digits de NN et MM.

Multiplication

Si N=∑ i=0 nN ib iN=\sum_{i=0}^n N_i b^i et M=∑ j=0 mM jb iM=\sum_{j=0}^m M_j b^i, alors on se ramène à la multiplication de deux polynômes : N+M=∑ i=0 n∑ j=0 mN iM jb i+jN+M=\sum_{i=0}^n \sum_{j=0}^m N_i M_j b^{i+j} Si b>2b&gt;2, par exemple b=10b=10 ou b=16b=16, N iM jN_i M_j sera très souvent supérieur à bb, il y aura souvent des retenues! On peut regrouper les termes b i+jb^{i+j} ayant le même exposant, en utilisant des décalages pour tenir compte de b ib^i (c’est l’algorithme de multiplication posée en base 10) ou additionner au fur et à mesure N iM jN_i M_j au coefficient actuel de b i+jb^{i+j} (exercice de la feuille de TD). En base 2, si on pose la multiplication comme en base 10, il est conseillé d’effectuer les additions une par une, sinon la prise en compte des retenues est délicate.

Le nombre de multiplications et d’additions est de l’ordre de nmnm (on peut montrer que c’est aussi le cas en tenant compte des retenues, car si une retenue se propage de rr positions dans la boucle intérieure en jj, elle ne pourra pas se propager de plus de 2 positions pour les r−2r-2 itérations suivantes de la boucle intérieure).

Comme le résultat a pour taille n+m+1n+m+1 ou n+m+2n+m+2 bits/chiffres/digits, on peut espérer l’existence d’algorithmes plus efficaces. Nous allons en présenter un dans le cadre plus simple de la multiplication des polynômes (où on n’a pas à gérer les retenues).

Diviser pour régner: multiplication de Karatsuba.

Soient P,QP, Q deux polynômes de degrés strictement inférieur à 2n2n. On suppose que le cout d’une opération arithmétique dans l’ensemble des coefficients vaut 1 et on néglige les autres opérations (on suppose par exemple que l’ensemble des coefficients est fini). On écrit P=A+x nB,Q=C+x nD P=A+x^n B, \quad Q=C+x^n D avec A,B,C,DA,B,C,D de degrés strictement inférieur à nn, on a alors : PQ=AC+x n(AD+BC)+x 2nBDP Q = AC + x^n(AD+BC)+x^{2n} BD Il y a 4 produits de polynômes de degrés <n&lt;n, mais au prix d’additions intermédiaires, on peut se ramener à 3 produits, en effet (A+B)(C+D)−AC−BD=AD+BC(A+B)(C+D)-AC-BD = AD+BC donc pour calculer le facteur de x nx^n il suffit de soustraire à (A+B)(C+D)(A+B)(C+D) les produits ACAC et BDBD que l’on doit calculer par ailleurs. Au premier ordre, le temps nécessaire pour multiplier les 2 polynômes de degré <2n&lt;2n est multiplié par 3, au lieu de 4.

Plus précisément, soit M(n)M(n) le temps nécessaire pour calculer le produit de 2 polynômes par cette méthode, on a alors M(2n)=3M(n)+8nM(2n) = 3M(n)+ 8n où 8n8n représente le nombre d’additions ou de soustractions pour former A+BA+B, C+DC+D, soustraire ACAC et BDBD, et tenir compte du fait que les termes de degré ≥n\geq n de ACAC se combinent avec ceux de degré <2n&lt;2n de AD+BCAD+BC et les termes de degré <3n&lt; 3n de x 2nBDx^{2n}BD avec ceux de degré ≥2n\geq 2n de AD+BCAD+BC. On en déduit u n=M(2 n),u n+1=3u n+8×2 nu_n=M(2^n), \quad u_{n+1}=3u_n+8 \times 2^n cette récurrence se résoud facilement par la commande

ce qui donne M(2 n)=u n=−8⋅2 n+9⋅3 nM(2^n)=u_n=-8\cdot 2^{n}+9\cdot 3^{n}.

Asymptotiquement, M(2 n)≈9⋅3 nM(2^n) \approx 9\cdot 3^{n} ce qui est bien meilleur que la multiplication naive en 2⋅4 n2 \cdot 4^n, mais pour de petites valeurs de nn, la multiplication naive est plus rapide comme l’illustre le graphe ci-dessous (en rouge la multiplication naive, en bleu Karatsuba complet) : 

plot([log(2*4^n),log(-8*2^n+9*3^n)],
n=1..10,display=[red,blue])

onload
On utilise Karatsuba (récursivement) uniquement pour des valeurs de nn suffisamment grandes (théoriquement lorsque 8n8n, le surcout dû aux additions est plus petit que la multiplication économisée, soit 8n<2n 28n&lt;2n^2 soit n>4n&gt;4, en pratique plutôt pour nn de l’ordre de quelques dizaines selon les implémentations, car nous n’avons tenu compte que des opérations arithmétiques).

Exercice : vérifier que la multiplication des entiers longs en langage Python se comporte pour nn grand comme la multiplication de Karatsuba.

Les logiciels de calcul intensif utilisent la plupart du temps des algorithmes plus efficaces lorsque nn est encore plus grand, qui se rapprochent de nln(n)n \ln(n). C’est en particulier le cas de tous les logiciels de calcul formel qui utilisent la libraire GMP écrite en langage C (Xcas, Sagemath, Maple, Mathematica, ...). Ces algorithmes utilisent la transformation de Fourier rapide (FFT).

3.9.3  Opérations de base sur les petits entiers.

Les petits entiers sont représentés en base b=2b=2, les opérations de base sont directement faites par le microprocesseur en se ramenant aux calculs en base 2:

3.9.4  Opérations de base sur les grands entiers

Les grands entiers naturels sont représentés dans une base qui est une puissance de 2, b=2 t′b=2^{t'}, les coefficients de l’écriture étant de petits entiers, ils sont stockés en mémoire à des adresses consécutives. On utilise le plus souvent t′=tt'=t la taille des régistres du microprocesseur, parfois une valeur légèrement inférieure pour simplifier la gestion des retenues.

Les implémentations de la multiplication font intervenir des algorithmes plus efficaces (comme Karatsuba), ainsi que des découpages par blocs visant à minimiser les opérations de lecture/écriture en mémoire au profit de l’utilisation de régistres et des caches de bas niveaux. En effet, avec les processeurs modernes, le cout du calcul d’un produit de 2 grands entiers par la méthode naive est dominé par les accès mémoire et non par les opérations arithmétiques elles-même.

3.10  Euclide

Dans cette section, on se concentre sur tout ce qui tourne autour de la division euclidienne.

3.10.1  Division euclidienne.

On a vu la définition de la division euclidienne pour les entiers naturels a=bq+r,r∈[0,b−1]a=bq+r, r\in[0,b-1]. Pour les entiers relatifs, plusieurs conventions de signe sont possible pour qq et rr, par exemple on peut conserver la convention de signe de reste positif. On peut aussi adopter la convention du reste symétrique, c’est-á-dire que r∈]−|b|2,|b|2]r \in ]-\frac{|b|}{2},\frac{|b|}{2}]. En général, on peut supposer que b>0b&gt;0, ce qui donne le programme suivant : 

def rsym(a,b):
    r = a % b
    if r<0:
        r += b
    if r<=b//2: 
        return r
    return r-b

onload

Optimisation : sur les microprocesseurs modernes, on évite si possible l’utilisation de tests pour optimiser le temps de calcul. En C, si m>0m&gt;0, et rr de signe quelconque sont des entiers 31 bits, unsigned(r)>>31 vaut 0 (si r≥0r\geq 0) ou -1 (sinon) et on peut écrire : 

int rsym(int r,int m){
  r = r % m;
  r += (unsigned(r)>>31)*m; // make positive
  return r-(unsigned((m>>1)-r)>>31)*m;
}

Pour les polynômes en une variable, si AA et BB sont deux polynômes à coefficients dans un corps KK, il existe un unique polynôme QQ et un unique polynôme RR tels que : A=BQ+R,deg(R)<deg(B)A=BQ+R, \quad \mbox{deg}(R)&lt;\ \mbox{deg}(B) L’unicité vient du fait que R−R′=B(Q−Q′)R-R'=B(Q-Q') est incompatible avec les degrés si Q≠Q′Q \neq Q'. L’existence se prouve par récurrence sur la différence de degré entre AA et BB, si elle est négative Q=0,R=AQ=0, R=A, si elle est nulle QQ est le quotient des coefficients dominants A aA_a de AA de degré aa et B bB_b de BB de degré bb, sinon, on pose Q=x a−bA a/B b+...Q=x^{a-b} A_a/B_b+..., ce qui annule le coefficient dominant de AA et on est ramené à une différence de degré diminuée de 1 (au moins).

Remarque 1 : si AA et BB sont à coefficients entiers et si BB est unitaire, alors QQ et RR sont à coefficients entiers. Sinon, on peut définir la pseudo-division de AA par BB en multipliant AA par le coefficient dominant de BB à la puissance la différence des degrés + 1.

Remarque 2 : pour les polynômes à plusieurs variables, il n’existe malheureusement pas de division euclidienne.

On dit que bb divise aa si le reste de la division de aa par bb est nul. Si aa et bb sont des entiers naturels non nuls et si a=bqa=bq alors a≥ba \geq b. Si AA et BB sont des polynômes alors le degré de BB est inférieur ou égal au degré de AA.

3.10.2  Anneau euclidien

L’ensemble des entiers relatifs est un anneau commutatif, mais il a des propriétés supplémentaires :

Cette propriété est également vraie pour les polynômes à coefficients dans un corps.

Observation : les seuls entiers ayant un inverse pour la multiplication sont 1 et -1. Pour les polynômes à coefficients dans un corps KK, il s’agit des polynômes constants non nuls (ceci est une conséquence du degré d’un produit qui est égal à la somme des degrés des facteurs). Ces éléments, appelés inversibles de l’anneau, jouent un role important, en effet si bb divise aa, i.e. a=bqa=bq alors pour tout élément inversible ii, bibi divise aa puisque a=(bi)i −1qa=(bi) i^{-1}q. Dans le cas des entiers, cela signifie qu’on peut négliger le signe pour étudier la divisibilité, autrement dit on peut se limiter aux entiers naturels. Dans le cas des polynômes, on peut multiplier par n’importe quel coefficient non nul, autrement dit on peut se limiter aux polynômes unitaires dont le coefficient dominant est 1.

3.10.3  Le PGCD

Soient aa et bb deux entiers. On se restreint au cas des entiers naturels en enlevant les signes.

Définition 4   On définit le PGCD (plus grand commun diviseur) de deux entiers naturels aa et bb non simultanément nuls comme le plus grand entier naturel dd qui divise simultanément aa et bb.

Si l’un des deux éléments est nul, disons bb, alors la réponse est aa. Sinon, comme 1 divise aa et de bb, l’ensemble des diviseurs communs à aa et bb est non vide, et son plus grand élément est inférieur ou égal au minimum de aa et de bb, le PGCD existe donc bien.

Pour le calculer, on observe que le PGCD de aa et de bb est le même que celui de bb et de rr, le reste de la division euclidienne de aa par bb. En effet si dd est un diviseur commun à aa et bb, alors a=da′a=da' et b=db′b=db' donc r=a−bq=d(a′−qb′)r=a-bq=d(a'-qb'). Réciproquement si dd est un diviseur commun à bb et rr, alors b=db′b=db' et r=dr′r=dr' donc a=bq+r=d(b′q+r′)a=bq+r=d(b'q+r'). Cet algorithme se traduit en un programme récursif : 

def pgcdr(a,b):
     if b=0: return a
     return pgcdr(b,a % b)

onload

A chaque appel récursif, la valeur du reste diminue au moins de 1 puisque r<br&lt;b. Au bout d’un nombre fini d’appels, b=0b=0 et on a le PGCD. Si on utilise le reste symétrique au lieu du reste euclidien positif, le reste perd un bit dans son écriture en base 2 à chaque itération, le nombre d’itérations est donc au plus égal au minimum du nombre de bits de aa et bb.

Exercice : traduire l’algorithme en un programme non récursif.

Le même algorithme fonctionne pour les polynômes, cette fois on cherche un polynôme de degré le plus grand possible qui divise AA et BB. L’algorithme ci-dessus s’adapte en remplaçant la division euclidienne des entiers par celle des polynômes, et le nombre d’appels récursifs est fini car le degré du reste qui diminue de 1 au moins à chaque itération. On obtient alors un PGCD, les autres PGCD sont des multiples non nuls de ce PGCD.

Il existe une version adaptée à l’écriture en base 2 des entiers sur machine, appelée PGCD binaire, qui évite de devoir faire des divisions euclidiennes (les divisions par 2 sont des décalages)

Chaque itération nécessite au plus le nombre nn de bits de l’écriture en base 2 de aa et bb, et il y a au plus 2n2n itérations (car on diminue de 1 le nombre de bits d’au moins un des 2 opérandes). Cet algorithme a donc un coût d’au plus O(n 2)O(n^2).

3.10.4  L’identité de Bézout

On va construire deux entiers relatifs uu et vv tels que au+bv=d=pgcd(a,b)au+bv=d = \mbox{pgcd}(a,b) On considère la matrice (L 1 1 0 a L 2 0 1 b)\left(\begin{array}{rccc} L_1 & 1 & 0 & a \\ L_2 & 0 & 1 & b \end{array}\right) qu’il faut interpréter ligne par ligne comme a×a\times1er coefficient+b×+b\times 2ème coefficient=3ième coefficient. L’algorithme de construction de uu et vv est une version arithmétique de l’algorithme du pivot de Gauss, où on souhaite créer un zéro dans la dernière colonne de la matrice. Mais seules les manipulations de lignes du type L n+2=L n−qL n+1L_{n+2}=L_n-qL_{n+1} avec qq entier sont admises. Le mieux que l’on puisse faire est donc de créer L 3=L 1−qL 2L_3=L_1-qL_2 avec qq quotient euclidien de aa par bb. Il faudra donc plusieurs manipulations de lignes avant d’avoir un 0, et la dernière colonne sera la liste des restes successsifs de l’algorithme d’Euclide pour calculer le PGCD de aa et bb. La ligne précédent la ligne avec un 0 en dernière colonne se termine par le PGCD (dernier reste non nul) et les deux premières colonnes donnent les coefficients de Bézout. Exemple avec a=125a=125 et b=45b=45 (L 1 1 0 125 L 2 0 1 45 L 3=L 1−2L 2 1 −2 35 L 4=L 2−L 3 −1 3 10 L 5=L 3−4L 4 4 −11 5 L 6=L 4−2L 5 −9 25 0 )\left(\begin{array}{rccc} L_1 & 1 & 0 & 125 \\ L_2 & 0 & 1 & 45 \\ L_3=L_1-2L_2 & 1 & -2 & 35\\ L_4=L_2-L_3 & -1 & 3 & 10 \\ L_5=L_3-4L_4 & 4 & -11 & 5 \\ L_6=L_4-2L_5 & -9 & 25 & 0 \\ \end{array}\right)

Algorithme :

Exemple : exécutez la commande ci-dessous qui calcule l’identité de Bézout pour 125 et 45 en affichant les étapes de l’algorithme dans la partie basse de la page HTML

Explication : A chaque itération, l’élément l[2] d’indice 2 de la liste l 1l_1 (ou l 2l_2) prend comme valeur l’un des restes successifs de l’algorithme d’Euclide, et on a un invariant de boucle a*l[0]+b*l[1]=l[2] qui donne l’identité de Bézout lorsque l[2] est le dernier reste non nul. Cet invariant se prouve facilement par récurrence.

Regardons maintenant la taille de uu et de vv. On peut supposer que a≥ba \geq b quitte à échanger aa et bb et on ignore les cas triviaux où aa et/ou bb valent 1.

Proposition 5   Si a≥b≥2a \geq b \geq 2, alors |u|≤b|u| \leq b et |v|≤a|v| \leq a.

Preuve (abrégée) : On considére les trois suites u nu_n (coefficients de la colonne d’indice 0), v nv_n (colonne d’indice 1) et r nr_n (colonne d’indice 2), ces suites vérifient la même relation de récurrence : u n+2=u n−q nu n+1,v n+2=v n−q nv n+1,r n+2=r n−q nr n+1u_{n+2}=u_n-q_nu_{n+1}, \quad v_{n+2}=v_n-q_nv_{n+1}, \quad r_{n+2}=r_n-q_nr_{n+1} où q nq_n est le quotient de la division de r nr_n par r n+1r_{n+1}. On montre par récurrence v nr n+1−v n+1r n=(−1) n+1a,(−1) n+1v ncroissante (stricte pourn≥2)v_n r_{n+1}-v_{n+1} r_n=(-1)^{n+1} a, \quad (-1)^{n+1}v_n \ \mbox{croissante (stricte pour}\ n\geq 2 ) au cran après Bézout, on a |v N+1|pgcd(a,b)=a|v_{N+1}| \mbox{pgcd}(a,b)=a donc la suite |v n||v_n| est majorée par aa et |v|≤a|v| \leq a On en déduit |u|≤b|u| \leq b puisque u=(1−bv)/au=(1-bv)/a.

Cet algorithme fonctionne à l’identique avec les polynômes à coefficients dans un corps, en utilisant la division euclidienne des polynômes.

L’identité de Bézout a de nombreuses conséquences en arithmétique.

3.10.5  Nombres premiers entre eux

Deux entiers sont dits premiers entre eux si leur PGCD vaut 1. Par exemple 22 et 15 sont premiers entre eux.

Théorème 6   (Lemme de Gauss) :
Si aa et bb sont premiers entre eux et si aa divise bcbc alors aa divise cc.

Preuve : par hypothèse, il existe un entier qq tel que bc=qabc=qa et par Bézout, deux entiers u,vu,v tels que au+bv=1au+bv=1, on élimine bb entre ces deux équations ce qui donne qav=bcv=bvc=(1−au)cqav=bcv=bvc=(1-au)c er c=a(qv+uc)c=a(qv+uc) est bien multiple de aa.

Proposition 7   Si aa et bb sont premiers entre eux et divisent tous deux cc alors abab divise cc.

En effet, il existe des entiers q,q′,u,vq,q',u,v tels que c=qa,c=q′b,au+bv=1c=qa, \quad c=q'b, \quad au+bv=1 Donc c=c(au+bv)=acu+bvc=aq′bu+bvqa=ab(q′u+qv)c=c(au+bv)=acu+bvc=aq'bu+bvqa=ab(q'u+qv)

Ces résultats restent vrais avec des polynômes à coefficients dans un corps (où premiers entre eux signifie que le PGCD est constant).

Application : décomposition d’une fraction nab\frac{n}{ab} avec aa et bb premiers entre eux (entiers ou polynômes).
Il existe alors uu et vv tels que au+bv=1au+bv=1, donc UU et VV tels que n=aU+bVn=aU+bV et nab=aU+bVab=aUab+bVab=Ub+Va\frac{n}{ab} = \frac{aU+bV}{ab}=\frac{aU}{ab}+\frac{bV}{ab}=\frac{U}{b}+\frac{V}{a} Exemple 1 : 23×5\frac{2}{3\times 5}, on a 1=2×3−1×51=2\times 3-1 \times 5 donc 2=4×3−2×52=4\times 3-2 \times 5 et 215=4×3−2×53×5=45−23\frac{2}{15} = \frac{4\times 3-2 \times 5}{3 \times 5}=\frac{4}{5}-\frac{2}{3}

Exemple 2 : 2x 2−1=2(x−1)(x+1)\frac{2}{x^2-1}=\frac{2}{(x-1)(x+1)} Les polynômes x−1x-1 et x+1x+1 sont premiers entre eux, on a l’identité de Bézout x−1−(x+1)=−2x-1 - (x+1) = -2 donc 2(x−1)(x+1)=(x+1)−(x−1)(x−1)(x+1)=1x−1−1x+1\frac{2}{(x-1)(x+1)} = \frac{(x+1)-(x-1)}{(x-1)(x+1)} = \frac{1}{x-1} - \frac{1}{x+1} Ce qui permet de calculer une primitive de 2x 2−1\frac{2}{x^2-1}. ∫2x 2−1=ln(|x−1|)−ln(|x+1|)\int \frac{2}{x^2-1}= \ln(|x-1|)-\ln(|x+1|)

3.10.6  Les restes chinois

L’identité de Bézout a une application très importante en calcul formel, elle permet de reconstruire un entier dont on connait le reste de la division par des entiers premiers entre eux :

Théorème 8   (restes chinois)
Soit mm et nn deux entiers naturels premiers entre eux, alors pour tout aa et bb entiers relatifs, il existe un entier cc tel que c−ac-a est divisible par mm et c−bc-b est divisible par nn. L’entier cc n’est pas unique, deux solutions cc et c′c' distinctes diffèrent par un multiple de nmnm.

Remarque : lorsqu’on reconstruit un entier naturel plus petit que nmnm, on choisit l’unique valeur de cc telle que c∈[0,nm[c \in [0,nm[ (c’est le reste de la division par nmnm de n’importe quelle solution), lorsqu’on reconstruit un entier relatif de valeur absolue plus petite que nm/2nm/2, on choisit cc tel que |c|≤nm/2|c|\leq nm/2 (on prend le reste symétrique).

Preuve : Si on a deux solutions cc et c′c' alors c−c′c-c' est divisible par nn et mm qui sont premiers entre eux, donc c−c′c-c' est divisible par nmnm. Cherchons une solution, cela revient à chercher UU et VV tels que c−a=mU,c−b=nVc-a=mU, \quad c-b= nV On a donc a+mU=b+nV⇔−mU+nV=b−aa+mU=b+nV \quad \Leftrightarrow \quad -mU+nV=b-a Or mm et nn sont premiers entre eux, donc par Bézout, il existe uu et vv tels que mu+nv=1mu+nv=1 il suffit de multiplier par b−ab-a pour trouver un couple de valeurs U,VU,V qui convient.

On sait que le couple U,VU,V n’est pas unique, si U′,V′U',V' est une autre solution, alors b−a=−mU+nV=−mU′+nV′b-a=-mU+nV=-mU'+nV' donc n(V−V′)=m(U−U′)n(V-V')=m(U-U') donc nn divise U−U′U-U' et mm divise V−V′V-V' avec le même quotient. Pour déterminer efficacement une valeur de cc qui convient, on calculera vv tel que |v|≤m|v|\leq m par l’algorithme de Bézout, on multipliera par b−ab-a et on prendra pour VV le reste de V(b−a)V(b-a) par mm.

Exemple : trouver cc tel que le reste de cc par 13 est 11 et le reste de cc par 7 est 3. On veut c=11+13U=3+7Vc=11+13U=3+7V donc 13U−7V=−813U-7V=-8 On a l’identité de Bézout : (−1)×13+7×2=1(-1)\times 13 + 7\times2 = 1 on multiplie par -8 donc on peut prendre V=16V=16 (U=8U=8), on préfère prendre le reste par 13 donc V=3V=3, c=24c=24 (et U=1U=1).

Remarque : en calcul formel, on prendra m=p 1...p k−1m=p_1...p_{k-1} un produit de nombres premiers distincts, et n=p kn=p_{k} un nombre premier distinct des précédents, on calcule la valeur de cc modulo chaque nombre premier p jp_j et dès qu’on a assez de nombres premiers p kp_k (i.e. si |c|<p 1...p k/2|c|&lt;p_1...p_k/2) on en déduit cc. On parle d’algorithme modulaire. Par exemple, on peut calculer le déterminant d’une matrice à coefficients entiers de cette manière.

Cette méthode s’applique aussi aux polynômes, mais elle est moins utile, car l’interpolation polynomiale de Lagrange est plus efficace lorsque les nombres premiers p ip_i sont remplacés par des x−α ix-\alpha_i

3.10.7  Nombres premiers, factorisation

Définition 9   Un entier naturel est premier s’il est divisible seulement par deux entiers distincts: 1 et lui-même.

Exemple : 17 est premier, mais pas 18. 1 n’est pas premier.

Un nombre non premier admet un diviseur différent de 1 et lui-même donc une factorisation n=pqn=pq avec p≠1p \neq 1 et q≠1q \neq 1. On peut supposer que p≤qp \leq q, pour savoir si un nombre est premier on peut donc tester que tous les entiers plus petits ou égaux à n\sqrt{n} ne divisent pas nn.

Proposition 10   Si pp et qq sont deux nombres premiers distincts, alors ils sont premiers entre eux.

En effet un diviseur commun d≠1d \neq 1 de pp et qq est égal à pp et à qq par définition de pp et qq premiers, absurde.

On en déduit par récurrence :

Proposition 11   Si q 1,..,q kq_1,..,q_k sont des nombres premiers tous distincts de pp premier, alors pp et q 1...q kq_1...q_k sont premiers entre eux.

On a vu que c’était le cas pour k=1k=1, pour k>1k&gt;1, soit d≠1d \neq 1 un diviseur commun de pp et q 1...q kq_1...q_k, dd est donc égal à pp, donc pp divise q 1...q kq_1...q_k. Comme q k≠pq_k \neq p, pp divise q 1..q k−1q_1..q_{k-1} ce qui contredit l’hypothèse de récurrence.

Théorème 12   Un entier naturel n>1n&gt;1 s’écrit de manière unique (à permutation près) comme produit de nombres premiers.

Preuve de l’existence par récurrence. Pour n=2n=2, on a 2=22=2. Pour n>2n&gt;2, si nn est premier, alors n=nn=n, sinon il admet un plus petit diviseur strictement supérieur à 1, on a n=pqn=pq avec p<np&lt;n premier et q<nq&lt;n à qui on applique l’hypothèse de récurrence.

Exercice : traduire cette preuve en un algorithme.

Unicité: supposons p 1...p k=p′ 1...p′ lp_1...p_k=p'_1...p'_l. Donc p jp_j divise p′ 1...p′ lp'_1...p'_l. Comme p jp_j est premier, il doit être égal à un des p′ lp'_l. On divise par ce facteur commun et on recommence.

Des résultats analogues s’appliquent pour les polynômes à coefficients dans un corps: on parle de facteur irréductible (analogue de nombre premier), divisible seulement par des polynomes constants et des polynomes multiples non nuls d’eux-mêmes. Par exemple tous les polynômes de degré 1 sont irréductibles. Un polynôme s’écrit de manière unique à permutation près et à inversibles près comme produit de facteurs irréductibles.

3.10.8  Remarque : la divisibilité est une relation d’ordre partiel

On observe que :

Une relation réflexive, antisymétrique et transitive est appelée relation d’ordre. Ici l’ordre est partiel, on peut avoir des éléments qui ne sont pas comparables, par exemple 2 ne divise pas 3 et 3 ne divise pas 2. Lorsque deux éléments sont comparables, on parle de relation d’ordre total (par exemple ≤\leq sur les entiers ou les réels ou l’ordre lexicographique sur les chaines de caractères ou les listes).

On verra plus loin un autre type de relation, les relations d’équivalence, où on remplace l’antisymétrie par la symétrie (aRba\ R\ b entraine bRab\ R\ a)

3.10.9  Prolongement : idéaux

Un idéal est un sous-ensemble d’un anneau qui est stable pour l’addition (la somme de deux éléments de l’idéal est encore dans l’idéal), et stable par multiplication par tout élément de l’anneau. Par exemple si on veut résoudre un système de plusieurs équations en une inconnue, l’ensemble des polynômes en une variable qui s’annulent en un point est un idéal.

Dans les anneaux euclidiens, tous les idéaux sont principaux, i.e. engendrés par un élément, on prend un plus petit entier en valeur absolue pour les entiers ou un plus petit polynôme en degré pour les polynômes, la division euclidienne par cet élément doit renvoyer 0 comme reste sinon on aurait un élément dans l’idéal plus petit. Dans l’exemple ci-dessus, on peut se ramener à une équation polynomiale en une variable (c’est un PGCD des polynômes des équations).

Dans les anneaux non euclidiens, comme par exemple les polynômes à plusieurs variables, les idéaux ne sont pas forcément engendrés par un seul élément, ce qui rend la résolution de systèmes polynomiaux à plusieurs inconnues beaucoup plus délicate qu’en une inconnue. Il est en général difficile de vérifier qu’un élément de l’anneau appartient ou non à l’idéal (alors que c’est juste une division dans le cas à une variable). Heureusement, il existe des familles particulières qui engendrent l’idéal pour lesquelles la vérification se fait par une généralisation de la division euclidienne. La recherche efficace de telles familles, appelées bases de Groebner, est un domaine de recherche actif!

3.11  Quelques algorithmes d’arithmétique de base.

3.11.1  Calcul de la racine carrée entière

Étant donné un entier NN, il s’agit de déterminer le plus grand entier nn tel que n 2≤Nn^2\leq N, nn est la racine carrée de NN. On choisit une base bb par exemple b=10b=10 pour un humain ou une puissance de 2 pour une machine, et on écrit NN en base bb, en découpant les chiffres par blocs de 2 en commençant par la droite, par exemple 2 00 00 00. On initialise la racine carrée nn à 0 et son carré cc à 0, on va calculer la racine carrée entière bloc par bloc en commençant par la gauche. Pour calculer le bloc suivant, on multiplie nn par bb et cc par b 2b^2 (c’est un simple décalage de l’écriture en ajoutant un ou deux zéros). Puis on ajoute les nombres impairs successifs 2n+12n+1, (2n+1)+2(2n+1)+2, ... à cc tant que l’on est inférieur à NN tronqué au bloc. Le nombre d’impairs successifs ajouté est ajouté à nn. En pratique, il suffit de conserver N−cN-c tronqué et de lui retrancher les impairs successifs.

Ainsi, pour 2 00 00 00, au 1er bloc 2, on initialise n=c=0n=c=0, on ajoute 2n+1=12n+1=1 à cc qui vaut alors 1 et on s’arrête car 1+3 est supérieur à 2. On passe au 2ième bloc, N−cN-c tronqué vaut 100, nn vaut 10, 2n+12n+1 vaut 21, on retranche donc à 100 successivement 21, 23, 25, 27 et on s’arrête car le reste est 4. Donc nn devient 14, et N−c=4N-c=4. On passe au troisième bloc, N−c=400N-c=400 et n=140n=140 donc 2n+1=2812n+1=281, on retranche de 400 les impairs successifs à partir de 281, ce qui n’est possible qu’une seule fois, cela donne N−c=119N-c=119 et n=141n=141. On passe au dernier bloc, N−c=11900N-c=11900 et n=1410n=1410 donc 2n+1=28212n+1=2821, on soustrait 2821, 2823, 2825, 2827 de 11900, il reste 604 et n=1414n=1414.

Exercice : calculer la quatrième décimale de 2\sqrt{2} de cette manière.

La complexité de cet algorithme est en O(log b(N) 2)O(\log_b(N)^2). En effet, pour calculer un chiffre il faut faire un nombre de soustraction au plus égal à bb, ces soustractions ayant au plus le nombre de chiffres de NN en base bb. (On peut accélérer le calcul à la manière de Karatsuba en choisissant une base bb puissance de 2 (ou 10) de l’ordre de N\sqrt{N} et en divisant pour régner).

isqrt(x):={
  local l,j,k,s,n,N,res;
  l:=revlist(convert(x,base,100));
  res:=seq(0,size(l));
  s:=0;
  N:=0;
  pour k de 0 jusque size(l)-1 faire
    N := (N-s)*100+l[k];
    n:=2*horner(res[0..k],10)+1;
    s:=n; // ajout de la somme des impairs consecutifs
    pour j de 0 jusque 10 faire
      si s>N alors break; fsi;
      n+=2;
      s+=n;
    fpour;
    s -= n;
    res[k]:=j;
  fpour;
  retourne horner(res,10);
}:;



3.11.2  Bezout sur les entiers et les fractions continues

Il existe une variante de l’identité de Bézout présentée ci-dessus pour les entiers. Soient a≥b>0a\geq b&gt;0 deux entiers, on pose (L n)au n−bv n=(−1) nr n(L_n) \quad a u_n - b v_n = (-1)^n r_n où r 0=a,r 1=br_0=a, r_1=b et r n+2r_{n+2} est le reste de la division euclidienne de r nr_n par r n+1r_{n+1} (q n+2q_{n+2} le quotient), u 0=1,u 1=0,v 0=0,v 1=1u_0=1, u_1=0, v_0=0,v_1=1. Comme précedemment, chaque ligne s’obtient par combinaison linéaire des deux précédentes, mais cette fois avec une addition L n+2=L n+q n+2L n+1L_{n+2}=L_n+q_{n+2} L_{n+1} ce qui se traduit par : u n+2=u n+q n+2u n+1,v n+2=v n+q n+2v n+1u_{n+2}=u_n+q_{n+2} u_{n+1}, \quad v_{n+2}=v_n+q_{n+2} v_{n+1} Les suites u nu_n et v nv_n sont alors strictement croissantes (à partir du rang 1 pour u nu_n). Au rang kk du dernier reste non nul on a : au k−bv k=(−1) kr k,r k=d=gcd(a,b)a u_k - b v_k = (-1)^k r_k, \quad r_k=d=\mbox{gcd}(a,b) et au rang suivant : au k+1−bv k+1=0au_{k+1} -b v_{k+1}=0 On montre par récurrence que v nr n+1+v n+1r n=av_n r_{n+1} + v_{n+1} r_n=a et une relation analogue pour u nu_n, on en déduit alors que v k+1=a/dv_{k+1}=a/d et u k+1=b/du_{k+1}=b/d (ce sont les cofacteurs du PPCM de aa et bb), en particulier les coefficients de Bézout vérifient u k<bu_k&lt;b et v k<av_k&lt;a.

On va aussi voir que u n+2/v n+2u_{n+2}/v_{n+2} est la nn-ième réduite du développement en fractions continues de a/ba/b (donc les coefficients de Bézout se lisent sur l’avant-dernière réduite). On introduit la notation [a 0,a 1,..,a n]=a 0+1a 1+1a 2+...a k[a_0,a_1,..,a_n] =a_0+\frac{1}{a_1+\frac{1}{a_2+\frac{...}{a_k}}} pour a 0≥0,a 1>0,...,a n>0a_0 \geq 0, a_1&gt;0, ..., a_n&gt;0. On a alors : ab=[q 2,q 3,..,q k]\frac{a}{b}=[q_2,q_3,..,q_k] En effet : ab=r 0r 1=q 2+r 2r 1=q 2+1r 1r 2=...\frac{a}{b}= \frac{r_0}{r_1}=q_2 +\frac{r_2}{r_1} = q_2 + \frac{1}{\frac{r_1}{r_2}} = ... D’autre part, on montre par récurrence sur n≥1n\geq 1 que si x>0x&gt;0 [q 2,...,q n,x]=v nx+v n−1u nx+u n−1[q_2,..., q_n,x]=\frac{v_{n}x+v_{n-1}}{u_{n}x+u_{n-1}} en effet au rang n=1n=1 [x]=x=v 1x+v 0u 1x+u 0[x]=x=\frac{v_1 x + v_0}{u_1 x+u_0 } et pour l’induction : [q 2,...,q n,x] = [q 2,...,q n−1,q n+1x] = v n−1(q n+1/x)+v n−2u n−1(q n+1/x)+u n−2 = x(v n−1q n+v n−2)+v n−1x(u n−1q n+u n−2)+u n−1 = v nx+v n−1u nx+u n−1 \begin{matrix} [q_2,..., q_n,x] &= & [q_2,..., q_{n-1},q_n+\frac{1}{x}] \\ &=& \frac{v_{n-1}(q_n+1/x)+v_{n-2}}{u_{n-1}(q_n+1/x)+u_{n-2}} \\ &=& \frac{x(v_{n-1}q_n+v_{n-2})+v_{n-1}}{x(u_{n-1}q_n+u_{n-2})+u_{n-1}} \\ & = & \frac{v_{n}x+v_{n-1}}{u_{n}x+u_{n-1}} \end{matrix} Donc au rang n−1n-1 et pour x=q nx=q_{n}, on obtient [q 2,...,q n]=v n+1u n+1[q_2,..., q_n]=\frac{v_{n+1}}{u_{n+1}}

Les fractions continues servent bien entendu aussi et d’abord à approcher les réels par des rationnels. L’algorithme de calcul des termes du développement est le suivant : Soit x≥0x\geq0. On initialise y=xy=x et la liste des a pa_p à vide. Puis on fait une boucle : on ajoute la partie entière de yy à la liste, on calcule la partie fractionnaire de yy, si elle est nulle on s’arrête (dans ce cas x∈ℚx\in \mathbb{Q}), sinon on stocke dans yy l’inverse de cette partie fractionnaire et on recommence. On note classiquement : h −2=0, h −1=1, h p=a ph p−1+h p−2 k −2=1, k −1=0, k p=a pk p−1+k p−2(3) \begin{matrix} h_{-2}=0, & h_{-1}=1, & h_p=a_p h_{p-1}+h_{p-2}\\ k_{-2}=1, & k_{-1}=0, & k_p=a_p k_{p-1}+k_{p-2} \end{matrix} \qquad (3) On a h 0=a 0,h 1=a 1a 0+1,k 0=1,k 1=a 1h_0=a_0, h_1=a_1 a_0+1, k_0=1, k_1=a_1. Les suites h ph_p et k pk_p sont donc positives et strictement croissantes pour p≥1p \geq 1, puisque pour p≥1p \geq 1, a p≥1a_p\geq 1, elles tendent vers l’infini au moins aussi vite que des suites de Fibonacci (à vitesse au moins géométrique donc). On a aussi aisément par récurrence : h pk p−1−h p−1k p=(−1) p+1(4) h_p k_{p-1} - h_{p-1}k_p=(-1)^{p+1} \qquad (4) On montre aussi comme ci-dessus : [a 0,...,a p−1,y]=yh p−1+h p−2yk p−1+k p−2[a_0,...,a_{p-1},y]=\frac{yh_{p-1}+h_{p-2}}{yk_{p-1}+k_{p-2}} On définit x px_p par x=[a 0,...,a p−1,x p]x=[a_0,...,a_{p-1},x_p], en faisant y=x py=x_p on a alors x=x ph p−1+h p−2x pk p−1+k p−2x=\frac{x_ph_{p-1}+h_{p-2}}{x_p k_{p-1}+k_{p-2}} ce qui donne x px_p en fonction de xx et a p=floor(−xk p−2−h p−2xk p−1−h p−1)a_p=\mbox{floor}\left( - \frac{xk_{p-2}-h_{p-2}}{xk_{p-1}-h_{p-1}} \right) En faisant y=a py=a_p on obtient [a 0,...,a p]=h pk p[a_0,...,a_p]=\frac{h_p}{k_p}. On montre ensuite que les suites (h p/k p)(h_p/k_p) pour les indices pairs et impairs sont deux suites adjacentes qui convergent vers xx, et on a h pk p−h p−1k p−1=(−1) p−1k pk p−1(5) \frac{h_p}{k_p} - \frac{h_{p-1}}{k_{p-1}} = \frac{(-1)^{p-1}}{k_p k_{p-1}} \qquad (5) En effet, la dernière égalité est une conséquence immédiate de (4), la croissance ou décroissance des suites d’indice pair ou impair s’en déduit en ajoutant (5) au cran suivant. La convergence vient de la limite infinie de k pk_p en l’infini. On a donc x=a 0+∑ p=0 ∞(−1) p−1k pk p+1,1k p(k p+k p+1)≤|x−h pk p|≤1k pk p+1x=a_0+\sum_{p=0}^\infty \frac{(-1)^{p-1}}{k_p k_{p+1}}, \quad \frac{1}{k_p(k_p+k_{p+1})} \leq |x-\frac{h_p}{k_p}| \leq \frac{1}{k_p k_{p+1}} La convergence est d’autant plus rapide que les k pk_p tendent rapidement vers l’infini, donc si les a pa_p sont plus grands que 1. La convergence la plus lente correspond au cas où tous les a p=1a_p=1 cas du nombre d’or, ou à partir d’un certain rang (nombre de Q[5]Q[\sqrt{5}]).

3.11.3  La puissance rapide itérative

Pour calculer a k(modn)a^k \pmod n, on décompose kk en base 2 k=∑ j=0 Jk j2 j,a k=∏ j=0 Ja k j2 j=∏ j/k j≠0a 2 jk=\sum_{j=0}^J k_j 2^j, \quad a^k = \prod_{j=0}^{J} a^{k_j 2^j} = \prod_{j/k_j \neq 0} a^{2^j} On initialise une variable B à 1, B vaudra a k(modn)a^k \pmod n en fin de calcul, on initialise une variable k à kk. On calcule dans une boucle les carrés successifs de a(modn)a \pmod n que l’on stocke dans une variable A (A vaudra donc successivement a(modn),a 2(modn),a 4(modn),...a \pmod n, a^2 \pmod n, a^{4} \pmod n, ...) et simultanément on teste si k jk_j vaut 1 en prenant le reste de la division par 2 de k (dans ce cas on multuplie B par A modulo nn), on divise ensuite k par 2 au sens du quotient euclidien. 

rapide(a,k,n):={
  local A,B;
  A:=a; B:=1;
  tantque k!=0 faire
    si irem(k,2)==1 alors B:=irem(A*B,n); fsi;
    k:=iquo(k,2);
    A:=irem(A*A,n);
  ftantque;
  return B;
}




Pour tester en mode pas à pas, ajoutez debug( au début de la commande qui précède
debug(rapide(123,4567890,123456))

3.11.4  Application de la puissance rapide: trouver un gńérateur de F p *=ℤ/p *F_p^*=\mathbb{Z}/p^*

Un entier gg est un générateur de F p *F_p^* si et seulement si son ordre est exactement p−1p-1. Si on suppose connue la factorisation de p−1p-1, cela revient à tester que g p−1p i≠1(modp)g^{\frac{p-1}{p_i}} \neq 1 \pmod p pour tous les diviseurs premiers p ip_i de p−1p-1. On essaie des entiers gg compris entre 1 et pp jusqu’à ce que ces tests soient vérifiés. La probabilité de réussite est de ϕ(p−1)/(p−1)\phi(p-1)/(p-1) donc en moyenne le nombre de tests à effectuer est (p−1)/ϕ(p−1)(p-1)/\phi(p-1). On va montrer que cela se comporte au pire en log(log(p))\log(\log(p)).

Le cas le pire à taille de pp fixée est celui où on prend tous les facteurs premiers depuis 2 et où la multiplicité des p ip_i dans la factorisation de p−1p-1 est 1, dans ce cas ϕ(p−1)p−1=∏ i(1−1p i),p−1=∏ ip i\frac{\phi(p-1)}{p-1} = \prod_i (1-\frac{1}{p_i}), \quad p-1=\prod_i p_i On va observer ce qui se passe lorsque on prend de plus en plus de facteurs premiers. Fixons nous un entier kk assez grand, on prend d’abord tous les nombres premiers entre 2 et kk, puis on regarde comment évolue log(∏ i(1−1p i))=∑ ilog(1−1p i)≥−∑ i1p i\log(\prod_i (1-\frac{1}{p_i})) = \sum_i \log(1-\frac{1}{p_i}) \geq -\sum_i \frac{1}{p_i} lorsqu’on ajoute les nombres premiers entre kk et 2k2k, puis entre 2k2k et 4k4k, ... puis entre 2 Nk2^Nk et 2 N+1k2^{N+1}k. Avec le résultat sur la densité des nombres premiers, il va y avoir k/log(k)k/\log(k) nombres premiers dans la première tranche, puis 2k/log(2k)2k/\log(2k), ... puis 2 Nk/log(2 Nk)2^Nk/\log(2^Nk), d’où la minoration log(∏ i(1−1p i))≥−1log(k)−1log(2k)−...−1log(2 Nk)\log(\prod_i (1-\frac{1}{p_i})) \geq -\frac{1}{\log(k)} - \frac{1}{\log(2k)} - ... - \frac{1}{\log(2^Nk)} Quand NN est grand, les log sont équivalents à Nlog(2)N\log(2) et on se retrouve avec la série harmonique, donc un équivalent de la minoration de log(ϕ(p−1)/(p−1))\log(\phi(p-1)/(p-1)) en −log(N)/log(2) -\log(N)/\log(2). Comparons NN avec la taille de pp. On a log(p−1)=∑ ilog(p i)\log(p-1)=\sum_i \log(p_i), si on fait le même découpage en tranche, le log de p−1p-1 est encadré à une constante multiplicative près par ∑ i≤N2 i=2 N+1−1\sum_{i\leq N} 2^i=2^{N+1}-1. Donc NN est en log(log(p))\log(\log(p)) et au final on a une minoration de ϕ(p−1)/(p−1)\phi(p-1)/(p-1) en O(log(log(log(p))))O(\log(log(\log(p)))) donc un double log en la taille de pp.

3.11.5  Application de la puissance rapide: cryptographie RSA

Interprétation du codage et du décodage.

Nous revenons sur la cryptographie RSA, en utilisant les notions vue dans les sections précédentes. On rappelle qu’on se donne un entier n=pqn=pq produit de deux nombres premiers pp et qq. Les entiers pp et qq sont secrets alors que nn est public. On choisit ensuite deux entiers dd et ee tels que de=1+k(p−1)(q−1)de=1 + k(p-1)(q-1), ou encore dd et ee sont inverses modulo l’indicatrice d’Euler de nn : de=1(modφ(n))de=1 \pmod {\varphi(n)} Le codage et le décodage se font avec l’algorithme de la puissance modulaire rapide : a→b=a e(modn)→c=b d(modn)a \rightarrow b=a^e \pmod n \rightarrow c=b^d \pmod n le cout est de O(ln(e)ln(n) 2)O(\ln(e)\ln(n)^2) et O(ln(d)ln(n) 2)O(\ln(d)\ln(n)^2) opérations.

On peut maintenant facilement montrer que ces deux fonctions sont réciproques l’une de l’autre lorsque aa est premier avec nn. En effet, on a a φ(n)=1(modn)a^{\varphi(n)}=1 \pmod n donc : c=a de(modn)=a 1+kφ(n)(modn)=a(a φ(n)) k(modn)=a(modn)c=a^{de} \pmod n = a^{1+k \varphi(n)} \pmod n =a (a^{\varphi(n)})^k \pmod n = a \pmod n

Comment générer des clefs

On choisit pp et qq en utilisant le test de Miller-Rabin. Par exemple

On choisit un couple de clefs privée-publique en utilisant l’identité de Bézout (ou inverse modulaire). Par exemple

Ici, on a pris E de sorte que l’exponentiation modulaire rapide à la puissance E nécessite peu d’opérations arithmétiques (17), comparé au calcul de la puissance d, ceci permet de faire l’opération “publique” plus rapidement (encodage ou vérification d’une signature) si le microprocesseur a peu de ressources (par exemple puce d’une carte bancaire).

Sur quoi repose la sécurité de RSA.

Logarithme discret

Il s’agit de résoudre en xx entier léquation g x=a(modn)g^x=a \pmod n Si aa est dans le sous-groupe multiplicatif engendré par gg, il y a des solutions qui sont égales modulo l’ordre de gg dans ℤ/n *\mathbb{Z}/n^*, et on peut noter x=log g(a)x=\log_g(a) la plus petite solution positive par analogie avec les réeks.

Pour trouver xx, on peut bien sur tester tous les entiers en partant de 1 (en multipliant par gg à chaque itération) jusqu’à trouver aa mais le cout est en O(nlog(n) 2)O(n \log(n)^2).

Si φ(n)\varphi(n) est connu (par exemple si n=pn=p est premier) et si on connait la factorisation de φ(n)\varphi(n) p−1=∏ ip i m ip-1=\prod_i p_i^{m_i} on peut faire mieux, on se ramène d’abord à résoudre plusieurs équations g i x i=a i(modn)g_i^{x_i}=a_i \pmod n, en prenant l’équation de départ à la puissance e i=(p−1)/p i m i=∏ j≠ip j m je_i=(p-1)/p_i^{m_i}=\prod_{j\neq i} p_j^{m_j}, (avec g i=g e ig_i=g^{e_i} et a i=a e ia_i=a^{e_i}), puis on retrouve xx à partir des x ix_i par les restes chinois appliqués à x=x i(modp i m i)x=x_i \pmod{p_i^{m_i}}. Ces équations sont moins couteuses à résoudre, car si gg est un générateur de ℤ/n *\mathbb{Z}/n^*, l’ordre de g ig_i est p i m ip_i^{m_i}. On peut d’ailleurs appliquer une méthode pp-adique pour trouver x ix_i et se ramener à des équations avec un g ig_i d’ordre p ip_i, en commençant par prendre l’équation à la puissance p i m i−1p_i^{m_i-1} ce qui donne x i(modp i)x_i \pmod {p_i} puis on cherche le 2ème terme du développement en puissance croissantes de x ix_i en prenant l’équation à la puissance p i m i−2p_i^{m_i-2}, etc.

On peut donc se ramener à résoudre une équation de logarithme discret avec gg d’ordre p ip_i plus petit (premier ou suffisament petit pour pouvoir faire une recherche exhaustive) et aa dans le sous-groupe engendré par gg.

Lorsque la recherche exhaustive en O(n)O(n) (à des termes logarithmiques près) devient couteuse (par exemple nn de l’ordre de quelques milliers), il est possible d’utiliser une méthode déterministe (pas de bébé, pas de géant) ou probabiliste de type Pollard-rho pour diminuer le cout en O(n)O(\sqrt{n}) (à des termes logarithmiques près).

La méthode pas de bébé, pas de géant décompose x=qm+rx=qm+r où mm est la partie entière d’un majorant de la racine carrée de l’ordre du sous-groupe engendré par gg. On calcule toutes les puissances rr-ième de gg de 0 à mm et on place les paires (g r,r)(g^r,r) dans une liste triée selon la valeur de g rg^r. On compare ensuite les a(g −m) qa(g^{-m})^q pour q de 0 à mm avec les éléments de la table, si on a égalité on en déduit x=qm+rx=qm+r. La méthode de Pollard-rho permet de gagner sur le stockage de la table des pas de bébé.

On peut aussi rechercher des kk tels que g kg^k se factorise sur une base de nombre premiers b ib_i fixée à l’avance g k=∏ ib i m i,kg^k=\prod_i b_i^{m_{i,k}} Ce qui permet de calculer les log en base gg des b ib_i en résolvant un système linéaire dans ℤ/(p−1)\mathbb{Z}/(p-1) (par restes chinois et calcul p ip_i-adique).

On cherche alors une puissance ll telle que ag lag^l se factorise sur la base, ce qui permet de trouver le log de aa.

3.12  Algorithmes rapides sur les polynômes en une variable.

On a présenté à la section 3.9.2 un algorithme de multiplication rapide de polynômes. On verra à la section 25.2 un algorithme plus rapide. En utilisant un algorithme de multiplication rapide de polynômes, on peut accélérer la division euclidienne et l’algorithme d’Euclide.

3.12.1  Inverse modulo x lx^l

Soit AA un polynôme de coefficient constant 1. On cherche BB de degré strictement plus petit que ll tel que AB=1(modx l)AB=1 \pmod{x^l}. Pour l=1l=1, il suffit de poser B=1B=1. Supposons le problème résolu pour ll AB=1(modx l)AB=1 \pmod{x^l} cherchons C=B+x lDC=B+x^lD tel que A(B+x lD)=1(modx 2l)A(B+x^lD)= 1 \pmod{x^{2l}} Ceci équivaut à AD=1−ABx l(modx l)AD=\frac{1-AB}{x^l} \pmod{x^l} donc D=B1−ABx l(modx l)D=B \frac{1-AB}{x^l} \pmod{x^l} on calcule B(1−AB)B(1-AB), on tronque pour avoir un polynôme de degré 2l−12l-1, on divise par x lx^l (on sait que 1−AB1-AB est divisible par x lx^l). On peut aussi calculer CC directement C=B+x lD=2B−AB 2(modx 2l)C=B+x^lD= 2B-AB^2 \pmod{x^{2l}} Pour passer de la précision ll à la précision 2l2l, le coût est celui d’une multiplication de 2 polynômes en précision ll et d’une multiplication en précision 2l2l, donc M(l)+M(2l)M(l)+M(2l). Si ll est une puissance de 2, l=2 nl=2^n, le cout total est ∑ k=1 nM(2 k)+M(2 k−1)\sum_{k=1}^n M(2^k)+M(2^{k-1}) Si la multiplication est Karatsuba, on obtient la même complexité pour l’inverse. Si la multiplication est la FFT en O(k2 k)O(k2^k), on obtient aussi la même complexité pour l’inverse.

3.12.2  Division euclidienne.

Soient AA et BB deux polynômes de degré respectifs nn et mm avec n≥mn \geq m, on pose la division euclidienne : A=BQ+RA=BQ+R Le degré de QQ est n−m+1n-m+1, le cout de l’algorithme naïf est donc proportionnel à m(n−m+1)m(n-m+1). Lorsque n−m+1n-m+1 et mm sont grands, on peut espérer trouver un algorithme plus efficace pour calculer QQ et RR.

On peut ramener le calcul de QQ à celui d’un inverse modulo une puissance de xx en faisant le changement de variables x=1/xx=1/x et de multiplier par x nx^n. Si aa est le polynôme ayant les coefficients de AA écrits dans l’autre sens, on a a(x)=x nA(1x)a(x) = x^n A(\frac{1}{x}) Avec les mêmes conventions de notations pour bb, qq et rr (en considérant RR comme de degré exactement égal à m−1m-1), on a : a=bq+x n−m+1ra=bq+x^{n-m+1} r Donc a=bq(modx n−m+1) a=bq \pmod{x^{n-m+1}}, et q=ainv(b(modx n−m+1))q=a \mbox{inv}(b \pmod{x^{n-m+1}}) ce qui détermine qq (car son degré est n−m+1n-m+1).

On calcule ensuite QQ et R=A−BQR=A-BQ.

3.12.3  Euclide

L’observation de base pour accélérer Euclide est que pour faire une étape d’Euclide, on divise génériquement un polynôme AA de degré nn par un polynôme BB de degré n−1n-1, avec un quotient de degré 1, dont les deux coefficients sont entièrement déterminés par les deux coefficients dominants et sous-dominants de AA et BB.

On peut généraliser cette observation si on pense en terme de l’algorithme d’Euclide étendu (Bézout), le vecteur colonne des deux restes R iR_i et R i+1R_{i+1} successifs à l’étape ii s’obtient en faisant le produit d’une matrice 2,2 par le vecteur colonne AA et BB (R i R i+1)=M i(A B)\begin{pmatrix} R_i \\ R_{i+1} \end{pmatrix} = M_i \begin{pmatrix} A \\ B \end{pmatrix} Les coefficients de M iM_i ne dépendent que des quotients initiaux, ils sont donc de degré au plus ii (environ). Si on change les ll coefficients de plus bas degré de AA ou BB, les coefficients de degré plus grand que (environ) i+li+l de R iR_i et R i+1R_{i+1} ne changent pas. On va prendre tout d’abord i=n/4i=n/4 et l=n/2l=n/2 (environ) pour AA et BB et plus loin dans l’algorithme i=n/4=li=n/4=l. Les termes de de degré >3n/4&gt;3n/4 de R iR_i et R i+1R_{i+1} ne vont pas dépendre des termes de degré plus petit que n/2n/2 de AA et BB, que l’on peut annuler, ce qui revient à mettre X n/2X^{n/2} en facteur ou encore à effectuer le calcul des n/4n/4 premières étapes avec uniquement les n/2n/2 coefficients de plus haut degré de AA et BB. Comme on a fait n/4n/4 étapes environ, on a diminué d’environ n/4n/4 le degré de R iR_i et R i+1R_{i+1}, qui encadrent 3n/43n/4. On refait alors environ n/4n/4 étapes en gardant les termes de degré plus grand que n/4n/4 de R kR_k et R k+1R_{k+1} (k=l=n/4k=l=n/4 ci-dessus), on a alors diminué de n/2n/2 le degré de la paire de restes successifs.

L’algorithme HGCD (half-GCD) calcule la matrice d’Euclide étendu qui permet de passer dans la suite des restes de A,BA,B avec deg(A)=n(A)=n>deg(B)(B) à R k,R k+1R_k,R_{k+1} avec deg(R k)≥n/2>(R_k) \geq n/2 &gt; deg(R k+1(R_{k+1}. Si nn est impair, on pose m=(n+1)/2m=(n+1)/2, si nn est pair, m=n/2m=n/2.

On a donc ramené hgcd en taille nn au calcul de deux hgcd en taille n/2n/2 et des produits de taille au plus nn.

Pour une preuve rigoureuse, on renvoie à Yap.

3.13  Pour en savoir plus.

Sur des aspects plus théoriques :

Sur des aspects plus pratiques, quelques références en ligne, la plupart sont accessibles gratuitement :

3.14  Exercices sur types, calcul exact et approché, algorithmes de bases

Vous pouvez tester directement dans votre navigateur Pour télécharger et installer Xcas sur votre ordinateur, suivre les instructions données sur
http://www-fourier.ujf-grenoble.fr/~parisse/giac_fr.html
Pour lancer xcas sous linux, cherchez Xcas dans le menu Education ou ouvrir un fenêtre terminal et taper la commande
xcas &
Lors de la première exécution, vous devrez choisir entre différents types de syntaxe (compatible C, maple ou TI89). Vous pouvez changer ce choix à tout moment en utilisant le menu Configuration->mode (syntaxe). On vous propose aussi d’ouvrir le tutoriel, qui est également accessible depuis le menu Aide, Débuter en calcul formel.

L’aide en ligne est accessible en tapant ?nom_de_commande. Dans Xcas, vous pouvez aussi taper le début d’un nom de commande puis la touche de tabulation (à gauche du A sur un clavier francais), sélectionner la commande dans la boite de dialogues puis cliquer sur Details pour avoir une aide plus complète dans votre navigateur. Pour plus de détails sur l’interface de Xcas, consultez le manuel (Aide->Interface). Si vous n’avez jamais utilisé de logiciel de calcul formel, vous pouvez commencer par lire le tutoriel (menu Aide->Debuter en calcul formel->tutoriel) et faire certains des exercices proposés (des corrigés sous forme de sessions Xcas sont dans Aide->Debuter en calcul formel->solutions)

Il peut être interessant de tester ces exercices en parallèle avec Xcas et des calculatrices formelles....

  1. À quelle vitesse votre logiciel multiplie-t-il des grands entiers (en fonction du nombre de chiffres)? On pourra tester le temps de calcul du produit de a(a+1)a(a+1) où a=10000!a=10 000!, a=15000!a=15000!, etc. . Même question pour des polynômes en une variable (à générer par exemple avec symb2poly(randpoly(n)) ou avec poly1[op(ranm(.))]).
  2. Comparer le temps de calcul de a n(modm)a^n \pmod m par la fonction powmod et la méthode prendre le reste modulo mm après avoir calculé a na^n.

    Programmez la méthode rapide et la méthode lente. Refaites la comparaison. Pour la méthode rapide, programmer aussi la version itérative utilisant la décomposition en base 2 de l’exposant : on stocke dans une variable locale bb les puissances successives a 2 0(modm),a 2 1(modm),...,a 2 k(modm),...a^{2^0} \pmod m,a^{2^1} \pmod m, ..., a^{2^k} \pmod m, ..., on forme a n(modn)a^n \pmod n en prenant le produit modulo mm de ces puissances successives lorsque le bit correspondant est à 1 (ce qui se détecte par le reste de divisions euclidiennes sucessives par 2, le calcul de bb et du bit correspondant se font dans une même boucle).
  3. Déterminer un entier cc tel que c=1(mod3)c=1 \pmod 3, c=3(mod5)c=3 \pmod 5, c=5(mod7)c=5 \pmod 7 et c=2(mod11)c=2 \pmod{11}.
  4. Calculez dans ℤ/11ℤ\mathbb{Z}/11\mathbb{Z} ∏ a=0 10(x−a) \prod_{a=0}^{10} (x-a)
  5. Algorithmes fondementaux : écrire des programmes implémentant
    1. le pgcd de 2 entiers
    2. l’algorithme de Bézout
    3. l’inverse modulaire en ne calculant que ce qui est nécessaire dans l’algorithme de Bézout
    4. les restes chinois
  6. Construire un corps fini de cardinal 128 (GF), puis factoriser le polynôme x 2−yx^2-y où yy est un élément quelconque du corps fini. Comparer avec la valeur de y\sqrt{y}.
  7. Utiliser la commande type ou whattype ou équivalent pour déterminer la représentation utilisée par le logiciel pour représenter une fraction, un nombre complexe, un flottant en précision machine, un flottant avec 100 décimales, la variable xx, l’expression sin(x)+2\sin(x)+2, la fonction x->sin(x), une liste, une séquence, un vecteur, une matrice. Essayez d’accéder aux parties de l’objet pour les objets composites (en utilisant op par exemple).
  8. Comparer le type de l’objet t si on effectue la commande t[2]:=0; après avoir purgé t ou après avoir affecté t:=[1,2,3] ?
  9. Comparer l’effet de l’affectation dans une liste et dans un vecteur ou une matrice sur votre logiciel (en Xcas, on peut utiliser =< au lieu de := pour stocker par référence).
  10. Le programme ci-dessous calcule la base utilisée pour représenter les flottants.
    Testez-le
    et expliquez.
  11. Déterminer le plus grand réel positif xx de la forme 2 −n2^{-n} (nn entier) tel que (1.0+x)−1.0(1.0+x)-1.0 renvoie 0 sur PC avec la précision par défaut puis avec Digits:=30.
  12. Calculer la valeur de a:=exp(π163)a:=\exp(\pi \sqrt{163}) avec 30 chiffres significatifs, puis sa partie fractionnaire. Proposez une commande permettant de décider si aa est un entier.
  13. Déterminer la valeur et le signe de la fraction rationnelle F(x,y)=13354y 6+x 2(11x 2y 2−y 6−121y 4−2)+112y 8+x2y F(x,y)= \frac{1335}{4} y^6 + x^2 (11x^2 y^2-y^6 -121y^4-2) + \frac{11}{2} y^8 + \frac{x}{2y} en x=77617x=77617 et y=33096y=33096 en faisant deux calculs, l’un en mode approché et l’autre en mode exact. Que pensez-vous de ces résultats? Combien de chiffres significatifs faut-il pour obtenir un résultat raisonnable en mode approché?

  14. Que se passe-t-il si on essaie d’appliquer l’algorithme de la puissance rapide pour calculer (x+y+z+1) k(x+y+z+1)^{k} par exemple pour k=64k=64 ? Calculer le nombre de termes dans le développement de (x+y+z+1) n(x+y+z+1)^n et expliquez.
  15. Programmation de la méthode de Horner
    Il s’agit d’évaluer efficacement un polynôme P(X)=a nX n+...+a 0 P(X) = a_n X^n + ... + a_0 en un point. On pose b 0=P(α)b_0=P(\alpha ) et on écrit : P(X)−b 0=(X−α)Q(X) P(X)-b_0=(X-\alpha )Q(X) où : Q(X)=b nX n−1+...+b 2X+b 1 Q(X) = b_n X^{n-1} + ... +b_2 X + b_1 On calcule alors par ordre décroissant b nb_n, b n−1b_{n-1}, ..., b 0b_0.
    1. Donner b nb_n en fonction de a na_n puis pour i≤n−1i\leq n-1, b ib_i en fonction de a ia_i et b i+1b_{i+1}. Indiquez le détail des calculs pour P(X)=X 3−2X+5P(X)=X^3-2X+5 et une valeur de α\alpha non nulle.
    2. Écrire un fonction horn effectuant ce calcul: on donnera en arguments le polynôme sous forme de la liste de ces coefficients (dans l’exemple [1,0,-2,5]) et la valeur de α\alpha et le programme renverra P(α)P(\alpha ). (On pourra aussi renvoyer les coefficients de QQ).
    3. En utilisant cette fonction, écrire une fonction qui calcule le développement de Taylor complet d’un polynôme en un point.
1
Certains systèmes de calcul formel (calculatrices par exemple) utilisent d’ailleurs des méthodes spécifiques pour gérer le problème de la fragmentation de la mémoire, appelés “garbage collector”. Ce type de méthode est intégré dans des langages comme Lisp ou Java, en C/C++ on trouve des libraries pour cela, par exemple GC de Boehm, incluse dans la distribution de GCC.
2
Les HP Prime utilisent Giac comme noyau de calcul formel, les TI Nspire CAS utilisent sans doute une version actualisée du système utilisé sur les TI 89, 92, Voayge 200.
3
un quartet=un demi octet
4
Plus précisément deux piles, la pile de donnée et la pile gérant le flux d’exécution. Cette dernière n’est pas visible par l’utilisateur
5
Sauf si on utilise comme dernier argument le nombre d’arguments de la fonction ou si on utilise (cf. infra) un tag de début de liste d’arguments
6
Toutefois une adaptation du logiciel utilisant comme quantum de base par exemple 32 bits porterait cette limite à 65536 65535−165536^{65535}-1

Chapitre 4  Les générateurs de nombres pseudo-aléatoires.

4.1  Selon la loi uniforme

Les générateurs d’entiers dans une plage donnée selon la loi uniforme servent en général de base pour générer des nombres aléatoires entiers ou non selon des lois classiques. Ils doivent à la fois être rapides, avoir une période égale à la plage donnée et avoir de bonnes propriétés statistiques.

Xcas utilise un “tiny” Mersenne Twister (de période environ 2 1272^{127}), certaines implantations de Giac utilisent un générateur congruentiel.

4.1.1  Les générateurs congruentiels à 1 cran.

Etant donnés trois entiers a,ca, c et mm on considère la suite u n+1=au n+c(modm)u_{n+1}=au_n+c \pmod m où on choisit (par exemple) comme représentant de u nu_n le reste de la division euclidienne par mm. La valeur de u 0u_0 est appelée seed en anglais, elle est initialisée usuellement soit à 0 (ce qui permet de reproduire des bugs dans un programme dépendant du hasard), soit avec l’horloge système ou tout autre entrée de l’ordinateur (par exemple périphériques).

On supposera que a≠1a\neq 1, le cas a=1a=1 n’est pas très intéressant. On a alors : u n=a nu 0+a n−1a−1c(modm)u_n=a^n u_0 + \frac{a^n-1}{a-1} c \pmod m On cherche à réaliser une période la plus grande possible idéalement mm, mais m−1m-1 peut fort bien convenir, et c’est possible si mm est premier en choisissant aa générateur du groupe cyclique, car on a alors a≠1(modm)a\neq 1 \pmod m et : u n=a n(u 0+ca−1)−ca−1(modm)u_n=a^n (u_0 + \frac{c}{a-1}) - \frac{c}{a-1} \pmod m donc la suite est stationnaire ou prend toutes les valeurs sauf −ca−1- \frac{c}{a-1} .

Exemple : choisir pour mm une puissance de 2 permet d’effectuer la division euclidienne très rapidement, mais cela a un inconvénient assez important : les bits de poids faible de u nu_n ont une périodicité très (trop) petite. Il est alors intéressant de prendre m=2 k±1m=2^k \pm 1, parce que la division euclidienne par mm peut se coder efficacement en base 2, on divise par 2 k2^k (décalage de kk bits) et on ajuste x=(2 k±1)q+r=2 kq+(r±q)x=(2^k \pm 1)q+r=2^k q + (r \pm q). Ainsi pour k=4k=4 et m=2 4+1=17m=2^4+1=17, mm est premier. On peut construire une suite de période 16 en choisissant aa générateur de (ℤ/17ℤ) *(\mathbb{Z}/17\mathbb{Z})^*, par exemple a=3a=3 et c=2c=2 donne la suite 0,2,8,9,12,4,14,10,15,13,7,6,3,11,1,50,2,8,9,12,4,14,10,15,13,7,6,3,11,1,5.

On a le :

Théorème 13   La suite (u n)(u_n) définie ci-dessus est de périodicité maximale mm si et seulement si :
  1. cc et mm sont premiers entre eux
  2. a−1a-1 est divisible par tous les facteurs premiers de mm
  3. a−1a-1 est multiple de 4 si mm l’est.

On observe d’abord que vouloir la périodicité maximale revient à pouvoir supposer que u 0=0u_0=0. Il est donc nécessaire d’avoir cc et mm premiers entre eux, sinon tous les u nu_n sont multiples du pgcd de cc et mm. Ensuite, on pose m=∏p i r im=\prod p_i^{r_i} la décomposition en facteurs premiers de mm et on raisonne modulo chaque premier (par le lemme chinois, la périodicité est le PPCM des périodicités modulo chaque p i r ip_i^{r_i}). Si a≠1(modp) ia\neq 1 \pmod p_i alors a−1a-1 est inversible modulo p ip_i donc modulo p i r ip_i^{r_i} on a u n=a n(u 0+ca−1)+−ca−1u_n=a^n (u_0 + \frac{c}{a-1}) + \frac{-c}{a-1} et la valeur −c/(a−1)-c/(a-1) ne peut pas être atteinte (ou alors la suite est stationnaire). Donc a−1a-1 doit être divisible par tous les facteurs premiers de mm pour avoir la périodicité maximale. Réciproquement, il faut trouver le premier ordre nn tel que (a n−1)/(a−1)=0(modp r)(a^n-1)/(a-1)=0 \pmod{p^r}. On pose a=b+1a=b+1, on a a n−1a−1=(b+1) n−1b=∑ k=1 n(nk)b k−1=n+n(n−1)2b+...\frac{a^n-1}{a-1}=\frac{(b+1)^n-1}{b} = \sum_{k=1}^n \left(^n_k\right) b^{k-1} = n +\frac{n(n-1)}{2}b +... On sait que b=a−1b=a-1 est un multiple de pp, disons b=qpb=qp, on en déduit que pour n=p rn=p^r, on a bien (a n−1)/(a−1)=0(modp r)(a^n-1)/(a-1)=0 \pmod{p^r}, alors que pour n=p r−1n=p^{r-1} et p≠2p\neq 2, (a n−1)/(a−1)=n(modp r)≠0(a^n-1)/(a-1)=n \pmod{p^r} \neq 0. Le même calcul pour p=2p=2 (prise en compte de la division par 2 de n(n−1)n(n-1)) donne la condition b=a−1b=a-1 est multiple de 4 si mm l’est.

On trouvera dans Knuth une discussion détaillée du choix de a,b,ma, b, m.

Exemple : m=2 31−1m=2^{31}-1 est premier, on peut donc construire un générateur congruentiel de période m−1m-1 en choisissant aa générateur de ℤ/mℤ *\mathbb{Z}/m\mathbb{Z}^*. Pour en trouver un, on peut tester aa pris au hasard et voir si a m−1j≠1(modm)a^{\frac{m-1}{j}} \neq 1 \pmod m pour tous les diviseurs premiers de m−1m-1. Par exemple





initialise l’état du générateur. Un appel à r() renvoie un entier entre 1 et m−1m-1, pour avoir un g’enérateur pseudo-aléatoire selon la loi uniforme sur ]0,1[]0,1[, on tape
.

Ainsi

permet de vérifier visuellement si les réels générés sont bien répartis, ou bien

qui détecte des biais invisibles avec le test précédent, par exemple pour
.

4.1.2  Récurrence à kk éléments

Au lieu d’une récurrence u k+1=au k+cu_{k+1}=au_k+c on conserve en mémoire k+1k+1 valeurs successives de la suite et on calcule u n+k+1=a 0u n+...+a ku n+k(modp)u_{n+k+1} = a_0 u_n+...+a_{k}u_{n+k} \pmod p Si on note U nU_n le vecteur (u n,...,u n+k)(u_n,...,u_{n+k}) et AA la matrice companion du polynôme a 0+a 1x+...+a kx ka_0+a_1x+...+a_kx^k, on a U n+1=AU nU_{n+1}=AU_n. Rechercher un générateur de période maximale revient à chercher AA d’ordre le plus grand possible, donc les valeurs propres de AA, i.e. les racines de PP, doivent être racines de l’unité d’ordre le plus grand possible donc p k−1p^k-1. Ce que l’on peut faire en construire un polynôme PP irréductible primitif (cf. la section 16 sur la construction de représentation des corps finis).

4.1.3  Mersenne twister.

Ce sont des générateurs plus performants, avec un état interne en général plus grand, dont l’état initial est généré par un générateur congruentiel. Ils utilisent une relation de récurrence qui ressemble aux générateurs congruentiels, mais au lieu de travailler sur de grands entiers, on découpe l’entier en mots de taille gérée par le CPU, et on fait des opérations de type matriciels avec des opérations bit à bit (ou exclusif par exemple) au lieu d’opérations arithmétiques.

4.2  Selon plusieurs lois classiques

La méthode générale consiste à calculer la distribution cumulée de la loi et à prendre la fonction réciproque d’un réel généré aléatoirement entre 0 et 1 selon la loi uniforme. Lorsqu’on a un nombre discret de valeurs possibles pas trop grand et que l’on veut générer plusieurs nombres selon la même loi, on peut précalculer la distribution cumulée en chaque valeur, et faire une dichotomie pour trouver la valeur de la fonction réciproque du nombre aléatoire généré. Les calculs peuvent être rendus difficiles par des dépassement de capacité des flottants si on utilise des méthodes naives pour estimer les fonction de répartition. On trouvera dans Abramowitz-Stegun diverses formules pour initialiser les méthodes de Newton pour inverser les fonction de répartition courante.

Il existe aussi quelques cas particuliers où on peut obtenir plus facilement un réel selon la loi donnée :

Chapitre 5  Le PGCD de polynômes.

Lorsqu’on travaille avec polynômes en une variable à coefficients sur un corps fini, l’algorithme d’Euclide est un algorithme efficace, pourvu que les degrés ne soient pas trop grands (dans ce cas il faut utiliser des algorithmes de type diviser pour régner, ici par exemple l’algorithme halfgcd qui consiste en gros à calculer d’abord l’identité de Bézout sous forme matricielle pour les polynômes tronqués en enlevant les degrés les plus bas). Mais l’algorithme d’Euclide est inefficace pour calculer le pgcd de deux polynômes à coefficients entiers ou à plusieurs variables. On va présenter ici les algorithmes utilisés habituellement par les systèmes de calcul formel: sous-résultant (PRS), modulaire (GCDMOD), pp-adique (EEZGD) et heuristique (GCDHEU). Le premier est une adaptation de l’algorithme d’Euclide et s’adapte à des coefficients assez génériques. Les trois autres ont en commun d’évaluer une ou plusieurs variables du polynôme (dans ce dernier cas il est nécessaire de bien distinguer le cas de polynômes à plusieurs variables) et de reconstruire le pgcd par des techniques distinctes, la plupart du temps ces algorithmes fonctionnent seulement si les coefficients sont entiers.

Soit donc P P et QQ deux polynômes à coefficients dans un corps. Le pgcd de PP et QQ n’est défini qu’à une constante près. Mais lorsque les coefficients de PP et QQ sont dans un anneau euclidien comme par exemple ℤ\mathbb{Z} ou ℤ[i]\mathbb{Z}[ i ], on appellera pgcd de PP et QQ un polynôme DD tel que P/DP / D et Q/DQ / D soient encore à coefficients dans l’anneau, et que DD soit optimal, c’est-à-dire que si un multiple μD\mu D de DD vérifie P/μDP / \mu D et Q/μDQ / \mu D sont à coefficients dans l’anneau, alors μ\mu est inversible.

La première étape d’un algorithme de calcul de pgcd consiste donc à diviser par son contenu (pgcd des coefficients entiers) chaque polynôme.

Exemple: P=4X 2−4P = 4 X^2 - 4 et Q=6X 2+12X+6Q = 6 X^2 + 12 X + 6. Le polynôme X+1X + 1 est un pgcd de PP et QQ puisqu’il est de degré maximal divisant PP et QQ mais le pgcd de PP et QQ est 2(X+1)2 ( X + 1 ). Remarquons qu’avec notre définition −2(X+1)- 2 ( X + 1 ) convient aussi. Par convention on appelera pgcd dans ℤ[X]\mathbb{Z}[X] le polynôme ayant un coefficient dominant positif.

Définition: On appelle contenu c(P)c ( P ) d’un polynôme PP le pgcd des coefficients de PP. On définit alors la partie primitive de PP: pp(P)=P/c(P)( P ) = P / c ( P ). Si c(P)=1c(P)=1, on dit que PP est primitif.

Proposition : Si AA et BB sont primitifs alors ABAB est primitif.
Sinon, on prend un facteur premier pp du contenu de ABAB, AB=0(modp)AB=0 \pmod p donc A=0A=0 ou B=0B=0 modulo pp, absurde.

Proposition : le contenu de ABAB est le produit des contenus de AA et de BB.
En effet le produit des contenus de AA et BB divise le contenu de ABAB, et AA/contenu de AA est primitif, BB/contenu de BB est primitif donc le produit l’est,

Proposition : Si AA et BB sont primitifs et si BB divise AA dans ℚ[X]\mathbb{Q}[X] alors A/B∈ℤ[X]A/B \in \mathbb{Z}[X].

Preuve : Soit Q=A/B∈ℚ[X]Q=A/B \in \mathbb{Q}[X]. Soit q∈ℕq \in \mathbb{N} le PPCM des dénominateurs des coefficients de QQ et notons P=qQ∈ℤ[X]P=qQ \in \mathbb{Z}[X]. On a A=BQA=BQ donc qA=BPqA=BP donc le contenu de qAqA est le produit du contenu de BB par celui de PP, donc le contenu de P=qQP=qQ est qq, donc Q∈ℤ[X]Q \in \mathbb{Z}[X].

Donc le PGCD de AA et BB, polynômes primitifs de ℤ[X]\mathbb{Z}[X] est obtenu en prenant un PGCD de AA et BB dans ℚ[X]\mathbb{Q}[X], en multipliant par le PPCM des dénominateurs et en rendant le polynôme obtenu primitif (on change le signe du résultat si nécessaire pour avoir un coefficient dominant positif).

On en déduit que : D=pgcd(P,Q)=pgcd(c(P),c(Q))pgcd(pp(P),pp(Q)) D = \mbox{pgcd} ( P, Q ) = \mbox{pgcd} ( c ( P ), c ( Q )) \mbox{pgcd} ( \mbox{pp} ( P ), \mbox{pp} ( Q ))

5.1  Le sous-résultant.

La première idée qui vient à l’esprit pour améliorer l’efficacité de l’algorithme d’Euclide consiste à éviter les fractions qui sont créées par les divisions euclidiennes. On utilise à cet effet la pseudo-division: au lieu de prendre le reste RR de la division euclidienne du polynôme PP par QQ, on prend le reste de la division de Pq δ+1P q^{\delta + 1} par QQ, où qq désigne le coefficient dominant de QQ et δ\delta la différence entre le degré de PP et de QQ.

Exercice: En utilisant votre système de calcul formel préféré, calculez les restes intermédiaires générés dans l’algorithme d’Euclide lorsqu’on utilise la pseudo-division par exemple pour les polynômes P(x)=(x+1) 7−(x−1) 6P ( x ) = ( x + 1 )^7 - ( x - 1 )^6 et sa dérivée.

Une solution avec giac/xcas: 

pgcd(a,b,prs):={ 
 local P,p,Q,q,R,g,h,d,res;
 res:=NULL;
 // convertit a et b en polynomes listes 
 // et extrait la partie primitive   
 P:=symb2poly(a);
 p:=lgcd(P); // pgcd des elements de la liste
 P:=P/p; 
 Q:=symb2poly(b);
 q:=lgcd(Q);
 Q:=Q/q; 
 if (size(P)<size(Q)){ // echange P et Q
  R:=P; P:=Q; Q:=R; 
 } 
 // calcul du contenu du pgcd
 p:=gcd(p,q);
 g:=1;
 h:=1;
 while (size(Q)!=1){
  q:=Q[0]; // coefficient dominant
  d:=size(P)-size(Q);
  R:=rem(q^(d+1)*P,Q);
  if (size(R)==0) return(p*poly12symb(Q/lgcd(Q),x));
  P:=Q;
  Q:=R;
  if (prs==1) Q:=Q/content(Q);
  if (prs==2) Q:=R/(g*h^d);
  res:=res,Q;
  if (prs==2) g:=q; h:=q^d/h^(d-1);
 } 
 return(p,res);
}:;



On s’aperçoit que les coefficients croissent de manière exponentielle (comparer avec ce qui se passe en mettant 1 comme dernier argument). La deuxième idée qui vient naturellement est alors à chaque étape de rendre le reste primitif, donc de diviser RR par le pgcd de ces coefficients. Cela donne un algorithme plus efficace, mais encore assez peu efficace car à chaque étape on doit calculer le pgcd de tous les coefficients, on peut imaginer le temps que cela prendra en dimension 1 et à fortiori en dimension supérieure. L’idéal serait de connaitre à l’avance une quantité suffisamment grande qui divise tous les coefficients du reste.

C’est ici qu’intervient l’algorithme du sous-résultant : après chaque pseudo-division euclidienne, on exhibe un coefficient "magique" qui divise les coefficients du reste (pour tester mettre le dernier argument de pgcd à 2). Ce coefficient n’est pas le pgcd mais il est suffisamment grand pour qu’on évite la croissance exponentielle des coefficients.

Algorithme du sous-résultant

Arguments: 2 polynômes PP et QQ primitifs. Valeur de retour: le pgcd de PP et QQ.

Pour calculer le coefficient "magique" on utilise 2 variables auxiliaires gg et hh initialisées a 1.

Boucle à effectuer tant que QQ est non nul:

Si on sort normalement de la boucle, QQ est nul, on renvoie donc la partie primitive de PP qui est le pgcd cherché.

Pour tester l’algorithme avec xcas, il suffit de décommenter les deux lignes Q:=R/(g*h^d); et g:=q; h:=q^d/h (d-1); ci-dessus.

La preuve de l’algorithme est un peu longue et par ailleurs bien expliquée dans le 2ème tome de Knuth (The Art of Computer Programming, Semi-numerical Algorithms), on y renvoie donc le lecteur intéressé. L’idée générale (et l’origine du nom de l’algorithme) est de considérer la matrice de Sylvester des polynômes de départ PP et QQ (celle dont le déterminant est appelé résultant de PP et QQ) et de traduire les pseudo-divisions qui permettent de calculer les restes successifs du sous-résultant en opération de ligne sur ces matrices. On démontre alors que les coefficients de RR divisés par gh δg h^{\delta} peuvent être interprétés comme des déterminants de sous-matrices de la matrice de Sylvester après réduction et c’est cela qui permet de conclure qu’ils sont entiers.

Par exemple, supposons que P=R 0P=R_0, Q=R 1Q=R_1, R 2R_2... diminuent de 1 en degré à chaque division (c’est le cas générique dans le déroulement de l’algorithme d’Euclide). Dans ce cas, δ=1\delta=1, il s’agit par exemple de montrer que le reste R 3R_3 de Q=R 1Q=R_1 par R 2R_2 est divisible par le carré du coefficient dominant de Q=R 1Q=R_1. Voyons comment on obtient les coefficients de R 3R_3 à partir de la matrice de Sylvester de PP et QQ. Prenons la sous-matrice constituée des 2 premières lignes de PP et des 3 premières lignes de QQ et réduisons-la sous forme échelonnée sans introduire de dénominateur. (p n p n−1 p n−2 p n−3 ... 0 p n p n−1 p n−2 ... q n−1 q n−2 q n−3 q n−4 ... 0 q n−1 q n−2 q n−3 ... 0 0 q n−1 q n−2 ...) \left( \begin{array}{ccccc} p_n & p_{n-1} & p_{n-2} & p_{n-3} & ... \\ 0 & p_n & p_{n-1} & p_{n-2} & ... \\ q_{n-1} & q_{n-2} & q_{n-3} & q_{n-4} & ... \\ 0 & q_{n-1} & q_{n-2} & q_{n-3} & ... \\ 0 & 0 & q_{n-1} & q_{n-2} & ... \end{array} \right) On effectue L 1←q n−1L 1−p nL 3L_1 \leftarrow q_{n-1} L_1 - p_n L_3 et L 2←q n−1L 2−p nL 4L_2 \leftarrow q_{n-1} L_2 - p_n L_4, ce qui correspond à l’élimination du terme en xx du quotient de PP par QQ (0 q n−1p n−1−p nq n−2 ... ... ... 0 0 q n−1p n−1−p nq n−2 ... ... q n−1 q n−2 q n−3 q n−4 ... 0 q n−1 q n−2 q n−3 ... 0 0 q n−1 q n−2 ...) \left( \begin{array}{ccccc} 0 & q_{n-1} p_{n-1} - p_n q_{n-2} & ... & ... & ... \\ 0 & 0 & q_{n-1} p_{n-1} - p_n q_{n-2} & ... & ... \\ q_{n-1} & q_{n-2} & q_{n-3} & q_{n-4} & ... \\ 0 & q_{n-1} & q_{n-2} & q_{n-3} & ... \\ 0 & 0 & q_{n-1} & q_{n-2} & ... \end{array} \right) on effectue ensuite L 1 ← q n−1L 1−(q n−1p n−1−p nq n−2)L 4 L 2 ← q n−1L 2−(q n−1p n−1−p nq n−2)L 5 \begin{matrix} L_1 & \leftarrow &q_{n-1} L_1 - (q_{n-1} p_{n-1} - p_n q_{n-2}) L_4 \\ L_2 & \leftarrow & q_{n-1} L_2 - (q_{n-1} p_{n-1} - p_n q_{n-2}) L_5 \end{matrix} ce qui correspond à l’élimination du terme constant du quotient de PP par QQ, on obtient (0 0 r 2,n−2 ... ... 0 0 0 r 2,n−2 ... q n−1 q n−2 q n−3 q n−4 ... 0 q n−1 q n−2 q n−3 ... 0 0 q n−1 q n−2 ...) \left( \begin{array}{ccccc} 0 & 0 & r_{2,n-2} & ... & ... \\ 0 & 0 & 0 & r_{2,n-2} & ... \\ q_{n-1} & q_{n-2} & q_{n-3} & q_{n-4} & ... \\ 0 & q_{n-1} & q_{n-2} & q_{n-3} & ... \\ 0 & 0 & q_{n-1} & q_{n-2} & ... \end{array} \right) si on enlève les lignes 3 et 4, et les colonnes 1 et 2, on obtient (après échanges de lignes) une sous-matrice de la matrice de Sylvester de QQ et R 2R_2 (q n−1 q n−2 ... r 2,n−2 ... ... 0 r 2,n−2 ...) \left( \begin{array}{ccc} q_{n-1} & q_{n-2} & ... \\ r_{2,n-2} & ... & ... \\ 0 & r_{2,n-2} & ... \end{array} \right) On recommence les opérations de réduction de cette sous-matrice correspondant à la division euclidienne de QQ par R 2R_2, on obtient (0 0 r 3,n−3 r 2,n−2 ... ... 0 r 2,n−2 ...) \left( \begin{array}{ccc} 0 & 0 & r_{3,n-3} \\ r_{2,n-2} & ... & ... \\ 0 & r_{2,n-2} & ... \end{array} \right) puis après suppression des colonnes 1 et 2 et des lignes 2 et 3 la ligne des coefficients de R 3R_3.

Supposons qu’on se limite dès le début de la réduction à ne garder que les colonnes 1 à 4 et une 5-ième colonne parmi les suivantes, on obtient à la fin de la réduction une matrice 1,1 qui contient un des coefficients de R 3R_3 (selon le choix de la 5-ième colonne). Donc ce coefficient est égal au déterminant de la matrice 1,1 qui est égal, au signe près, au déterminant de la matrice 3,3 dont il est issu par notre réduction (en effet, dans la 2ième partie de la réduction, on a multiplié deux fois L 1L_1 par r 2,n−2r_{2,n-2}, mais on doit ensuite diviser le déterminant par r 2,n−2 2r_{2,n-2}^2 pour éliminer les colonnes 1 et 2). Quant au déterminant de la matrice 3,3, il se déduit du déterminant de la matrice 5,5 par multiplication par q n−1 4q_{n-1}^4 (2 lignes ont été multipliées 2 fois par q n−1q_{n-1}) et division par q n−1 2q_{n-1}^2 (élimination des colonnes 1 et 2). Au final, tout coefficient de R 3R_3 est égal au produit d’un déterminant 5,5 extrait de la matrice de Sylvester de PP et QQ par q n−1 2q_{n-1}^2, qui est justement le coefficient “magique” par lequel on divise le reste de R 1=QR_1=Q par R 2R_2 lors de l’algorithme du sous-résultant.

5.2  Le pgcd en une variable.

5.2.1  Le pgcd heuristique.

On suppose ici que les coefficients sont entiers ou entiers de Gauss. On peut donc se ramener au cas où les polynômes sont primitifs.

L’idée consiste à évaluer PP et QQ en un entier zz et à extraire des informations du pgcd gg des entiers P(z)P ( z ) et Q(z)Q ( z ). Il faut donc un moyen de remonter de l’entier gg à un polynôme GG tel que G(z)=gG ( z ) = g. La méthode consiste à écrire en base zz l’entier gg, avec une particularité dans les divisions euclidiennes successives on utilise le reste symétrique (compris entre −z/2- z / 2 et z/2z / 2). Cette écriture donne les coefficients d’un polynôme GG unique. On extrait ensuite la partie primitive de ce polynôme GG. Lorsque zz est assez grand par rapport aux coefficients des polynômes PP et QQ, si pp(G)\mbox{pp} ( G ) divise PP et QQ, on va montrer que le pgcd de PP et de QQ est D=pp(G)D = \mbox{pp} ( G ).

On remarque tout d’abord que d:=D(z)d : = D ( z ) divise gg. En effet DD divise PP et QQ donc pour tout entier (ou entier de Gauss) zz, D(z)D ( z ) divise P(z)P ( z ) et Q(z)Q ( z ). Il existe donc une constante aa telle que g=ad g = a d On a aussi pp(G)\mbox{pp} ( G ) divise DD. Il existe donc un polynôme CC tel que : D=pp(G)C D = \mbox{pp} ( G ) C Nous devons prouver que CC est un polynôme constant. On suppose dans la suite que ce n’est pas le cas. Evaluons l’égalité précédente au point zz, on obtient d=gc(G)C(z) d = \frac{g}{c ( G )} C ( z ) Finalement 1=ac(G)C(z) 1 = \frac{a}{c ( G )} C ( z ) La procédure de construction de GG nous donne une majoration de ces coefficients par |z|/2| z | / 2, donc de c(G)c ( G ) par |z|/2| z | / 2, donc C(z)C ( z ) divise un entier de module plus petit que |z|/2| z | / 2, donc |C(z)|≤|z|2 | C ( z ) | \leq \frac{| z |}{2} On considère maintenant les racines complexes z 1,….,z nz_1, \ldots ., z_n du polynôme CC (il en existe au moins une puisqu’on a supposé CC non constant). On a: C(X)=c n(X−z 1)….(X−z n) C ( X ) = c_n ( X - z_1 ) \ldots . ( X - z_n ) Donc, comme c nc_n est un entier (ou entier de Gauss) non nul, sa norme est supérieure ou égale à 1 et : |C(z)|≥∏ j=1 n(|z|−|z j|) | C ( z ) | \geq \prod^n_{j = 1} ( | z | - | z_j | ) Il nous reste à majorer les racines de CC pour minorer |C(z)|| C ( z ) |. Comme CC divise DD il divise PP et QQ donc les racines de CC sont des racines communes à PP et QQ. On va appliquer le:

Lemme 14   Soit x une racine complexe d’un polynôme P=a nX n+….+a 0P = a_n X^n + \ldots . + a_0.

Alors  |x|<|P||a n|+1,|P|=max 0≤i≤n(|a i|) | x | &lt; \frac{| P |}{| a_n |} + 1, | P | = \max_{0 \leq i \leq n} ( | a_i | )

Application du lemme à C(X)C(X) : on a 1/|c n|≤11/|c_n|\leq 1 donc si on a choisi zz tel que |z|≥2min(|P|,|Q|)+2| z | \geq 2 \min ( | P |, | Q | ) + 2, alors pour tout jj, |z j|<|z|/2| z_j | &lt; | z | / 2 donc |C(z)|>(|z|2) n | C ( z ) | &gt; \left( \frac{| z |}{2} \right)^n qui contredit notre majoration de |C(z)|| C ( z ) |.

Théorème 15   Soit PP et Q deux polynômes à coefficients entiers. On choisit un entier z tel que |z|≥2min(|P|,|Q|)+2| z | \geq 2 \min ( | P |, | Q | ) + 2, si la partie primitive du polynôme GG reconstruit à partir du pgcd de P(z)etP ( z ) \mbox{et}Q(z) par écriture en base zz (avec comme reste euclidien le reste symétrique) divise PP et QQ alors c’est le pgcd de PP et QQ.

Pour finir la démonstration du théorème, il nous faut encore montrer le lemme. On a −a nx n=a n−1x n−1+….+a 0 - a_n x^n = a_{n - 1} x^{n - 1} + \ldots . + a_0 Donc |a n||x| n≤|P|(1+….+|x| n−1)=|P||x| n−1|x|−1 | a_n | | x |^n \leq | P | ( 1 + \ldots . + | x |^{n - 1} ) = | P | \frac{| x |^n - 1}{| x | - 1} Ici on peut supposer que |x|≥1| x | \geq 1, sinon le lemme est démontré, donc |x|−1| x | - 1 est positif et |a n|(|x|−1)≤|P||x| n−1|x| n⇒|x|−1<|P||a n| | a_n | ( | x | - 1 ) \leq | P | \frac{| x |^n - 1}{| x |^n} \Rightarrow | x | - 1 &lt; \frac{| P |}{| a_n |} Remarques

Exemple 16   Si P 0=6(X 2−1)P_0 = 6 ( X^2 - 1 ) et Q 0=4(X 3−1)Q_0 = 4 ( X^3 - 1 ).

Le contenu de P 0P_0 est 6, celui de Q 0Q_0 est 4.
On a donc pgcd des contenus = 2, P=X 2−1,Q=X 3−1P = X^2 - 1, Q = X^3 - 1. La valeur initiale de zz est donc 2*1+2=42 \ast 1 + 2 = 4. On trouve P(4)=15,Q(4)=63P ( 4 ) = 15, Q ( 4 ) = 63. Le pgcd entier de 15 et 63 est 3 que nous écrivons symétriquement en base 4 sous la forme 3=1*4−13 = 1 \ast 4 - 1, donc G=X−1G = X - 1, sa partie primitive est X−1X - 1. On teste si X−1X - 1 divise PP et QQ, c’est le cas, donc c’est le pgcd de PP et QQ et le pgcd de P 0P_0 et Q 0Q_0 est 2(X−1)2 ( X - 1 ).

Algorithme gcdheu
En arguments deux polynômes P 0P_0 et Q 0Q_0 à coefficients entiers ou entiers de Gauss. Retourne le pgcd de P 0P_0 et Q 0Q_0 ou faux en cas d’échec.

  1. Calculer le contenu de P 0P_0 et Q 0Q_0. Vérifier que les coefficients sont entiers de Gauss sinon retourner faux.
  2. Extraire la partie primitive PP de P 0P_0 et QQ de Q 0Q_0, calculer le pgcd cc des contenus de P 0P_0 et Q 0Q_0
  3. Déterminer z=2min(|P|,|Q|)+2z = 2 \min ( | P |, | Q | ) + 2.
  4. Début de boucle: initialisation du nombre d’essais à 1, test d’arrêt sur un nombre maximal d’essais, avec changement de zz entre deux itérations (par exemple z←2zz \leftarrow 2 z).
  5. Calculer le pgcd gg de P(z)P ( z ) et Q(z)Q ( z ) puis son écriture symétrique en base zz dont on extrait la partie primitive GG.
  6. Si GnedivisepasG \mbox{ne} \mbox{divise} \mbox{pas}PP passer à l’itération suivante. De même pour QQ.
  7. Retourner cGc G
  8. Fin de la boucle
  9. Retourner faux.

On remarque au passage qu’on a calculé le quotient de PP par GG et le quotient de QQ par GG lorsque la procédure réussit. On peut donc passer à la procédure gcdheu deux paramètres supplémentaires par référence, les deux polynômes que l’on affectera en cas de succès, ce qui optimise la simplification d’une fraction de 2 polynômes.

5.2.2  Le pgcd modulaire

On part du fait que si DD est le pgcd de PP et QQ dans ℤ\mathbb{Z} (ou ℤ[i])\mathbb{Z} [ i ] ) alors après réduction modulo un nombre premier nn qui ne divise pas les coefficients dominants de PP et QQ, DD divise le pgcd GG de PP et QQ dans ℤ/nℤ\mathbb{Z} / n \mathbb{Z} (par convention, le pgcd dans ℤ/nℤ\mathbb{Z} / n \mathbb{Z} est normalisé pour que son coefficient dominant vaille 1). Comme on calcule GG dans ℤ/nℤ\mathbb{Z} / n \mathbb{Z}, les coefficients des restes intermédiaires de l’algorithme d’Euclide sont bornés, on évite ainsi la croissance exponentielle des coefficients. Il faudra ensuite reconstruire DD à partir de GG.

On remarque d’abord que si on trouve G=1,G = 1, alors PP et QQ sont premiers entre eux. En général, on peut seulement dire que le degré de GG est supérieur ou égal au degré de DD. En fait, le degré de GG est égal au degré de DD lorsque les restes de l’algorithme d’Euclide (calculé en effectuant des pseudo-divisions, cf. l’exercice 1) ont leur coefficient dominant non divisible par nn. Donc plus nn est grand, plus la probabilité est grande de trouver GG du bon degré.

Dans la suite, nous allons déterminer une borne bb à priori majorant les coefficients de DD. On utilisera ensuite la même méthode que dans l’algorithme modulaire de recherche de racines évidentes: on multiplie GG dans ℤ/nℤ\mathbb{Z} / n \mathbb{Z} par le pgcd dans ℤ\mathbb{Z} des coefficients dominants pp et qq de PP et QQ. Soit D˜=pgcd(p,q)G\tilde{D} = \mbox{pgcd} ( p, q ) G le résultat écrit en représentation symétrique. Si n≥bpgcd(p,q)n \geq b \mbox{pgcd} ( p, q ) et si GG est du bon degré, on montre de la même manière que D=D˜D = \tilde{D}. Comme on ne connait pas le degré de DD, on est obligé de tester si D˜\tilde{D} divise PP et QQ. Si c’est le cas, alors D˜\tilde{D} divise DD donc D˜=D\tilde{D} = D puisque degre(D˜)=degre(G)≥degre(D)\mbox{degre} ( \tilde{D} ) = \mbox{degre} ( G ) \geq \mbox{degre} ( D ). Sinon, nn est un nombre premier malchanceux pour ce calcul de pgcd (degre(G)≥degre(D)\mbox{degre} ( G ) \geq \mbox{degre} ( D )), il faut essayer un autre premier.

Remarque: On serait tenté de dire que les coefficients de DD sont bornés par le plus grand coefficient de PP. C’est malheureusement faux, par exemple (X+1) 2( X + 1 )^2 dont le plus grand coefficient est 2 divise (X+1) 2(X−1)( X + 1 )^2 ( X - 1 ) dont le plus grand coefficient (en valeur absolue) est 1.

Soit P=∑p iX iP = \sum p_i X^i un polynôme à coefficients entiers. On utilise la norme euclidienne |P| 2=∑|p i| 2(6) | P |^2 = \sum | p_i |^2 \qquad (6) On établit d’abord une majoration du produit des racines de norme supérieure à 1 de PP à l’aide de |P| | P |^{}. Ensuite si DD est un diviseur de PP, le coefficient dominant dd de DD divise le coefficient dominant pp de PP et les racines de DD sont aussi des racines de PP. On pourra donc déterminer une majoration des polynômes symétriques des racines de DD et donc des coefficients de DD.

Lemme 17   Soit A=∑ j=0 aa jX jA = \sum_{j = 0}^a a_j X^j un polynôme et α∈ℂ\alpha \in \mathbb{C}. Alors |(X−α)A|=|(α¯X−1)A| | ( X - \alpha ) A | = | ( \overline{\alpha} X - 1 ) A |

Pour prouver le lemme 17, on développe les produits de polynômes. On pose a −1=a a+1=0a_{-1} = a_{a + 1} = 0 et on note ℜ\Re la partie réelle. |(X−α)A| 2=∑ j=0 a+1|a j−1−αa j| 2=∑ j=0 a+1|a j−1| 2+|α| 2|a j| 2−2ℜ(a j−1αa j¯) | ( X - \alpha ) A |^2 = \sum_{j = 0}^{a + 1} | a_{j - 1} - \alpha a_j |^2 = \sum_{j = 0}^{a + 1} | a_{j - 1} |^2 + | \alpha |^2 | a_j |^2 - 2 \Re ( a_{j - 1} \overline{\alpha a_j} ) |(α¯X−1)A| 2=∑ j=0 a+1|α¯a j−1−a j| 2=∑ j=0 a+1|α| 2|a j−1| 2+|a j| 2−2ℜ(α¯a j−1a j¯) | ( \overline{\alpha} X - 1 ) A |^2 = \sum_{j = 0}^{a + 1} | \overline{\alpha} a_{j - 1} - a_j |^2 = \sum_{j = 0}^{a + 1} | \alpha |^2 | a_{j - 1} |^2 + | a_j |^2 - 2 \Re ( \overline{\alpha} a_{j - 1} \overline{a_j} ) Les deux donnent bien le même résultat.

Soit P(X)=p∏(X−α j)P ( X ) = p \prod ( X - \alpha_j ) la factorisation de PP sur ℂ\mathbb{C}. On introduit le polynôme P˜=p∏ j/|α j|≥1(X−α j)∏ j/|α j|<1(α j¯X−1) \tilde{P} = p \prod_{j / | \alpha_j | \geq 1} ( X - \alpha_j ) \prod_{j / | \alpha_j | &lt; 1} ( \overline{\alpha_j} X - 1 ) qui d’après le lemme a la même norme que PP. La norme de PP majore donc le coefficient constant de P˜\tilde{P} d’où: ∏ j/|α j|≥1|α j|≤|P||p|(7) \prod_{j / | \alpha_j | \geq 1} | \alpha_j | \leq \frac{| P |}{| p |} \qquad (7) On remarque que (7) reste vraie si on considère les racines δ j\delta_j de norme plus grande que 1 d’un diviseur DD de PP puisque le produit porte alors sur un sous-ensemble. On écrit maintenant l’expression des coefficients d jd_j de DD à l’aide des racines δ j\delta_j de DD: |d m−j|=|d||∑ choixdejracinesparmilesmracinesdeD∏ δ k∈racineschoisiesδ k| | d_{m - j} | = | d | \left| \sum_{\mbox{choix} \mbox{de} j \mbox{racines} \mbox{parmi} \mbox{les} m \mbox{racines} \mbox{de} D} \quad \prod_{\delta_k \in \mbox{racines} \mbox{choisies}} \delta_k \right| Pour majorer |d m−j|| d_{m - j} |, on commence par majorer |δ k|| \delta_k | par β k=max(1,|δ k|)\beta_k = \max ( 1, | \delta_k | ). On est donc ramené à majorer σ j,m(β)=∑ choixdejparmimvaleursβ k∏ β k∈choixβ k \sigma_{j, m} ( \beta ) = \sum_{\mbox{choix} \mbox{de} j \mbox{parmi} m \mbox{valeurs} \beta_k} \quad \prod_{\beta_k \in \mbox{choix}} \beta_k avec pour hypothèse une majoration de M=∏ k=1 mβ kM = \prod_{k = 1}^m \beta_k donnée par la relation (7). Pour cela, on cherche le maximum de σ j,m(β)\sigma_{j, m} ( \beta ) sous les contraintes MM fixé et β k≥1\beta_k \geq 1.

On va montrer que le maximum ne peut être atteint que si l’un des β k=M\beta_k = M (et tous les autres β k=1)\beta_k = 1 ). Sinon, quitte à réordonner supposons que les β k\beta_k sont classés par ordre croissant. On a donc β m−1≠1\beta_{m - 1} \neq 1, on pose β k˜=β k\widetilde{\beta_k} = \beta_k pour k≤m−2k \leq m - 2, β˜ m−1=1\tilde{\beta}_{m - 1} = 1 et β˜ m=β m−1β m\tilde{\beta}_m = \beta_{m - 1} \beta_m. Comparons σ j,m(β)\sigma_{j, m} ( \beta ) et σ j,nm(β˜)\sigma_{j, \mbox{nm}} ( \tilde{\beta} ). Si le choix de jj parmi mm comporte k=m−1k = m - 1 et k=mk = m, le produit est inchangé. Sinon on a la somme de deux produits, l’un contenant k=m−1k = m - 1 et l’autre k=mk = m. On compare donc B(β m−1+β m)B ( \beta_{m - 1} + \beta_m ) et B(1+β m−1β m)B ( 1 + \beta_{m - 1} \beta_m ) avec B=∏ β k∈resteduchoixβ kB = \prod_{\beta_k \in \mbox{reste} \mbox{du} \mbox{choix}} \beta_k. Comme 1+β m−1β m≥β m−1+β m 1 + \beta_{m - 1} \beta_m \geq \beta_{m - 1} + \beta_m puisque la différence est le produit (1−β m)(1−β m−1)(1-\beta_m)(1-\beta_{m-1}) de deux nombres positifs, on arrive à la contradiction souhaitée.

Ensuite on décompose les choix de σ m,j\sigma_{m, j} en ceux contenant MM et des 1 et ceux ne contenant que des 1, d’où la majoration σ j,m(β)≤(m−1 j−1)M+(m−1 j) \sigma_{j, m} ( \beta ) \leq \left(\begin{array}{c} m - 1\\ j - 1 \end{array}\right) M + \left(\begin{array}{c} m - 1\\ j \end{array}\right) et finalement |d m−j|≤|d|((m−1 j−1)|P||p|+(m−1 j))(8) | d_{m - j} | \leq | d | \left( \left(\begin{array}{c} m - 1\\ j - 1 \end{array}\right) \frac{| P |}{| p |} + \left(\begin{array}{c} m - 1\\ j \end{array}\right) \right) \qquad (8) On peut en déduire une majoration indépendante de jj sur les coefficients de DD, en majorant |d|| d | par |p|| p | (puisque dd divise pp) et les coefficients binomiaux par 2 m−12^{m - 1} (obtenue en développant (1+1) m−1( 1 + 1 )^{m - 1}). D’où le

Théorème 18   (Landau-Mignotte) Soit PP un polynôme à coefficients entiers (ou entiers de Gauss) et DD un diviseur de PP de degré mm. Si |P|| P | désigne la norme euclidienne du vecteur des coefficients de PP et pp le coefficient dominant de PP alors les coefficients d jd_j de DD satisfont l’inégalité |d j|≤2 m−1(|P|+|p|)(9) | d_j | \leq 2^{m - 1} ( | P | + | p | ) \qquad (9)

Avec cette estimation, on en déduit que si nn est un premier plus grand que min(2 degre(P)−1(|P|+|p|),2 degre(Q)−1(|Q|+|q|)),(10) \min \left( 2^{\mbox{degre} ( P ) - 1} ( | P | + | p | ), 2^{\mbox{degre} ( Q ) - 1} ( | Q | + | q | ) \right), \qquad (10) alors le pgcd trouvé dans ℤ/nℤ\mathbb{Z} / n \mathbb{Z} va se reconstruire en un pgcd dans ℤ\mathbb{Z} si son degré est le bon.

Malheureusement la borne précédente est souvent très grande par rapport aux coefficients du pgcd et calculer dans ℤ/nℤ\mathbb{Z} / n \mathbb{Z} s’avèrera encore inefficace (surtout si le pgcd est 1). Cela reste vrai même si on optimise un peu la majoration (10) en repartant de (8).

L’idée est donc de travailler modulo plusieurs nombres premiers plus petits et reconstruire le pgcd des 2 polynômes à coefficients entiers à partir des pgcd des polynômes dans ℤ/nℤ\mathbb{Z} / n \mathbb{Z} et du théorème des restes chinois. En pratique on prend des nombres premiers inférieurs à la racine carrée du plus grand entier hardware de la machine (donc plus petits que 2 162^{16} sur une machine 32 bits) ce qui permet d’utiliser l’arithmétique hardware du processeur sans risque de débordement.

Algorithme du PGCD modulaire en 1 variable:

En argument: 2 polynômes primitifs PP et QQ à coefficients entiers. Le résultat renvoyé sera le polynôme pgcd.

Variable auxiliaire: un entier NN initialisé à 1 qui représente le produit des nombres premiers utilisés jusqu’ici et un polynôme HH initialisé à 0 qui représente le pgcd dans ℤ/Nℤ\mathbb{Z} / N \mathbb{Z}.

Boucle infinie :

  1. Chercher un nouveau nombre premier nn qui ne divise pas les coefficients dominants pp et qq de PP et QQ
  2. Calculer le pgcd GG de PP et QQ dans ℤ/nℤ\mathbb{Z} / n \mathbb{Z}. Si GG=1, renvoyer 1.
  3. Si H=0H = 0 ou si le degré de GG est plus petit que le degré de HH, recopier GG dans HH et nn dans NN, passer à la 6ème étape
  4. Si le degré de GG est plus grand que celui de HH passer à l’itération suivante
  5. Si le degré de GG est égal au degré de HH, en utilisant le théorème des restes chinois, calculer un polynôme H˜\tilde{H} tel que H˜=H\tilde{H} = H modulo NN et H˜=G\tilde{H} = G modulo nn. Recopier H˜\tilde{H} dans HH et nNn N dans NN.
  6. Ecrire pgcd(p,q)H\mbox{pgcd} ( p, q ) H en représentation symétrique. Soit H˜\tilde{H} le résultat rendu primitif. Tester si H˜\tilde{H} divise PP et QQ. Si c’est le cas, renvoyer H˜\tilde{H}, sinon passer à l’itération suivante.

Finalement on n’a pas utilisé bb, la borne de Landau-Mignotte. On peut penser que l’étape 6 ne devrait être effectuée que lorsque NN est plus grand que pgcd(p,q)b\mbox{pgcd} ( p, q ) b. En pratique, on effectue le test de l’étape 6 plus tôt parce que les coefficients du pgcd sont rarement aussi grand que bb. Mais pour éviter de faire le test trop tôt, on introduit une variable auxiliaire H′H' qui contient la valeur de HH de l’itération précédente et on ne fait le test que si H′=HH' = H (ou bien sûr si on a dépassé la borne).

Remarque:

L’algorithme ci-dessus fonctionne également pour des polynômes à plusieurs variables.

Exemple 1:

Calcul du pgcd de (X+1) 3(X−1) 4( X + 1 )^3 ( X - 1 )^4 et (X 4−1) ( X^4 - 1 )^{}. Prenons pour commencer n=2n = 2. On trouve comme pgcd X 4+1X^4 + 1 (en effet −1=1- 1 = 1 donc on cherchait le pgcd de (X+1) 7( X + 1 )^7 et de X 4+1=(X+1) 4X^4 + 1 = ( X + 1 )^4). On teste si X 4+1X^4 + 1 divise PP et QQ, ce n’est pas le cas donc on passe au nombre premier suivant. Pour n=3n = 3, on trouve X 2−1X^2 - 1. Donc n=2n = 2 n’était pas un bon nombre premier pour ce calcul de pgcd puisqu’on a trouvé un pgcd de degré plus petit. On teste si X 2−1X^2 - 1 divise PP et QQ, c’est le cas ici donc on peut arrêter, le pgcd cherché est X 2−1X^2-1.

Exemple 2 :

Calcul du pgcd de (X+1) 3(X−1) 4( X + 1 )^3 ( X - 1 )^4 et (X 4−1) 3( X^4 - 1 )^3. Pour n=2n = 2, on trouve un polynôme de degré 7. Pour n=3n = 3, on trouve X 6−1X^6 - 1 donc n=2n = 2 était une mauvaise réduction. Comme X 6−1X^6 - 1 ne divise pas PP et QQ, on passe à n=5n = 5. On trouve X 6+2X 4−2X 2−1X^6 + 2 X^4 - 2 X^2 - 1. On applique le théorème des restes chinois qui va nous donner un polynôme dans ℤ/15ℤ\mathbb{Z} / 15 \mathbb{Z}. On cherche donc un entier congru à 2 modulo 5 et à 0 modulo 3, -3 est la solution (écrite en représentation symétrique), donc le polynôme modulo 15 est X 6−3X 4+3X 2−1=(X 2−1) 3X^6 - 3 X^4 + 3 X^2 - 1 = ( X^2 - 1 )^3. Ce polynôme divise PP et QQ, c’est donc le pgcd de PP et de QQ.

5.3  Le pgcd à plusieurs variables.

5.3.1  Le pgcd heuristique.

On suppose comme dans le cas à une variable que les polynômes sont primitifs, donc qu’on a simplifié les polynômes par le pgcd entier de leurs coefficients entiers.

Le principe est identique à celui du PGCD à 1 variable, on évalue les deux polynômes PP et QQ de kk variables X 1,….,X kX_1, \ldots ., X_k en un X k=zX_k = z et on calcule le pgcd gg des 2 polynômes P(z)P ( z ) et Q(z)Q ( z ) de k−1k - 1 variables. On remonte ensuite à un polynôme GG par écriture symétrique en base zz de gg et on teste si pp(G)\mbox{pp} ( G ) divise PP et QQ. Il s’agit à nouveau de montrer que si zz est assez grand, alors pp(G)\mbox{pp} ( G ) est le pgcd cherché. On sait que d=D(z)d = D ( z ) divise gg. Il existe donc un polynôme aa de k−1k - 1 variables tel que g=adg = a d. On sait aussi que pp(G)\mbox{pp} ( G ) divise DD, donc il existe un polynôme CC de kk variables tel que D=C*pp(G).D = C \ast \mbox{pp} ( G ) . On évalue en zz et on obtient d=C(z)g/c(G)d = C ( z ) g / c ( G ), où c(G)c ( G ) est un entier, donc c(G)=a*C(z) c ( G ) = a \ast C ( z ) Comme c(G)c ( G ) est un entier, aa et C(z)C ( z ) sont des polynômes constants. Comme précédemment, on a aussi |C(z)|≤|z|/2| C ( z ) | \leq | z | / 2 puisque |c(G)|≤|z|/2| c ( G ) | \leq | z | / 2.

En pratique, cet algorithme nécessite le calcul récursif de pgcd sans garantie de réussite. On l’évite donc s’il y a beaucoup de variables (la limite est par exemple de 5 pour MuPAD).

5.3.2  Le pgcd modulaire multivariables.

Ici, on travaille modulo X n−α X_n - \alpha_{}, où X 1,….,X nX_1, \ldots ., X_n désignent les variables des polynômes. On considère donc deux polynômes PP et QQ comme polynômes de la variables X nX_n avec des coefficients dans ℤ[X 1,….,X n−1]\mathbb{Z} [ X_1, \ldots ., X_{n - 1} ]. On évalue en X n=αX_n = \alpha, on obtient deux polynômes en n−1n - 1 variables dont on calcule le pgcd (récursivement).

Il s’agit de reconstruire le pgcd par interpolation. Tout d’abord, on a une borne évidente sur le degré du pgcd par rapport à la variable X nX_n, c’est le minimum δ\delta des degrés par rapport à X nX_n des polynômes PP et QQ. A première vue, il suffit donc d’évaluer les polynômes en δ+1\delta + 1 points α\alpha.

Il faut toutefois prendre garde aux mauvaises évaluations et à la normalisation des pgcd avant d’interpoler. En effet, si D(X 1,….,X n)D ( X_1, \ldots ., X_n ) désigne le pgcd de PP et QQ et G(X 1,….,X n−1)G ( X_1, \ldots ., X_{n - 1} ) le pgcd de P(X 1,….,X n−1,α)P ( X_1, \ldots ., X_{n - 1}, \alpha ) et de Q(X 1,….,X n−1,α)Q ( X_1, \ldots ., X_{n - 1}, \alpha ), on peut seulement dire D(X 1,….,X n−1,α)D ( X_1, \ldots ., X_{n - 1}, \alpha ) divise GG. Plusieurs cas sont donc possibles lorsqu’on évalue en un nouveau point α\alpha:

On voit que les mauvaises évaluations se détectent simplement par les degrés. Pour la normalisation, on utilise une petite astuce: au lieu de reconstruire lepgcdD\mbox{le} \mbox{pgcd} D, on va reconstruire un multiple du pgcd DD (ce multiple appartiendra à ℤ[X n])\mathbb{Z} [ X_n ] ). On voit maintenant PP et QQ comme des polynômes en n−1n - 1 variables X 1,….,X n−1X_1, \ldots ., X_{n - 1} à coefficients dans ℤ[X n]\mathbb{Z} [ X_n ]. Alors lcoeff(D)(D), le coefficient dominant de DD (relativement à l’ordre lexicographique sur les variables X 1,...,X n−1X_1,...,X_{n-1}), est un polynôme en X nX_n qui divise le coefficient dominant de PP et de QQ donc divise le coefficient dominant du pgcd des coefficients dominants de PP et de QQ. On va donc reconstruire le polynôme : D′=DΔ(X n)lcoeff(D)(X n),Δ(X n)=pgcd(lcoeff(P)(X n),lcoeff(Q)(X n)) D' = D \frac{\Delta ( X_n )}{\mbox{lcoeff} ( D ) ( X_n )}, \Delta ( X_n ) = \mbox{pgcd} ( \mbox{lcoeff} ( P ) ( X_n ), \mbox{lcoeff} ( Q ) ( X_n )) c’est-à-dire DD multiplié par un polynôme qui ne dépend que de X nX_n.

Revenons à GG en un point α\alpha de bonne évaluation. C’est un multiple entier de D(X 1,….,X n−1,α)D ( X_1, \ldots ., X_{n - 1}, \alpha ): G=βD(X 1,….,X n−1,α) G = \beta D ( X_1, \ldots ., X_{n - 1}, \alpha ) Donc, comme polynômes de X 1,...,X n−1X_1,...,X_{n-1} à coefficients dans ℤ[X n]\mathbb{Z}[X_n] ou dans ℤ\mathbb{Z}, lcoeff(G)=βlcoeff(D) |X n=α\mbox{lcoeff} ( G ) = \beta \mbox{lcoeff} ( D )_{| X_n = \alpha}. Comme lcoeff(D)\mbox{lcoeff} ( D ) divise Δ(X n)\Delta ( X_n ), il en est de même en X n=αX_n = \alpha donc lcoeff(G)(G) divise βΔ(α)\beta \Delta(\alpha). On en déduit que Δ(α)G \Delta ( \alpha) G qui est divisible par Δ(α)β \Delta (\alpha) \beta est divisible par lcoeff(G)\mbox{lcoeff} ( G ). On va donc considérer le polynôme Δ(α)G/lcoeff(G) \Delta (\alpha) G / \mbox{lcoeff} ( G ) : ses coefficients sont entiers et son coefficient dominant est Δ(α)=lcoeff(D′(X 1,….,X n−1,α))\Delta ( \alpha) = \mbox{lcoeff}(D'( X_1, \ldots ., X_{n - 1}, \alpha )) donc Δ(α)G/lcoeff(G)=D′(X 1,….,X n−1,α) \Delta (\alpha) G / \mbox{lcoeff} ( G )= D'( X_1, \ldots ., X_{n - 1}, \alpha )

Algorithme du pgcd modulaire à plusieurs variables (interpolation dense):

Arguments: 2 polynômes primitifs PP et QQ de nn variables X 1,….,X nX_1, \ldots ., X_n à coefficients entiers. Renvoie le pgcd de PP et QQ.

  1. Si n=1n = 1, renvoyer le pgcd de PP et QQ en une variable.
  2. Test rapide de pgcd trivial par rapport à X nX_n. On cherche des n−1n - 1-uplets α\alpha tels que P(α,X n)P ( \alpha, X_n ) et Q(α,X n)Q ( \alpha, X_n ) soient de même degré que PP et QQ par rapport à la variable X nX_n. On calcule le pgcd GG de ces 2 polynômes en une variable. Si le pgcd est constant, alors on retourne le pgcd des coefficients de PP et QQ.
  3. On divise PP et QQ par leur contenu respectifs vu comme polynômes en X 1,….,X n−1X_1, \ldots ., X_{n - 1} à coefficients dans ℤ[X n]\mathbb{Z} [ X_n ], on note C(X n)C ( X_n ) le pgcd des contenus. On calcule aussi le pgcd Δ(X n)\Delta ( X_n ) des coefficients dominants de PP et de QQ.
  4. On initialise D′D' le pgcd reconstruit à 0, I(X n)I ( X_n ) le polynôme d’interpolation à 1, δ=(δ 1,...,δ n−1)\delta=(\delta_1,...,\delta_{n-1}) la liste des degrés partiels du pgcd par rapport à X 1,….,X n−1X_1, \ldots ., X_{n - 1} au minimum des degrés partiels de PP et QQ par rapport à X 1,….,X n−1X_1, \ldots ., X_{n - 1}, ee le nombre d’évaluation à 0 et EE l’ensemble des points d’interpolation à la liste vide.
  5. Boucle infinie:
    • Faire α\alpha=entier aléatoire n’appartenant pas à EE jusqu’à ce que degre(P(X 1,….,X n−1,α))=degre X n(P(X 1,….,X n) degre(Q(X 1,….,X n−1,α))=degre X n(Q(X 1,….,X n)) \begin{matrix} \mbox{degre}(P ( X_1, \ldots ., X_{n - 1}, \alpha ))=\mbox{degre}_{X_n} ( P ( X_1, \ldots ., X_n ) & & \\ \mbox{degre} ( Q ( X_1, \ldots ., X_{n - 1}, \alpha )) = \mbox{degre}_{X_n} ( Q ( X_1, \ldots ., X_n )) & & \end{matrix}
    • Calculer le pgcd G(X 1,….,X n−1)G ( X_1, \ldots ., X_{n - 1} ) en n−1n - 1 variables de P(X 1,….,X n−1,α)P ( X_1, \ldots ., X_{n - 1}, \alpha ) et Q(X 1,….,X n−1,α)Q ( X_1, \ldots ., X_{n - 1}, \alpha ).
    • Si degre (G) i<δ i\mbox{degre}_{} ( G )_i &lt; \delta_i pour un indice au moins. Si degre(G)≤δ\mbox{degre} ( G ) \leq \delta, on pose δ=degre(G)\delta = \mbox{degre} ( G ), D′=GΔ(α)lcoeff(G)D' = G \frac{\Delta ( \alpha )}{\mbox{lcoeff} ( G )}, I=X n−αI = X_n - \alpha, e=1e = 1 et E=[α]E = [ \alpha ], sinon on pose δ=min(δ,degre(G)),D′=0,I=1,e=0,E=[]\delta = \min ( \delta, \mbox{degre} ( G )), D' = 0, I = 1, e = 0, E = [ ]. On passe à l’itération suivante.
    • Si degre(G)>δ\mbox{degre} ( G ) &gt; \delta, on passe à l’itération suivante.
    • Si degre(G)=δ\mbox{degre} ( G ) = \delta, on interpole:
      • G:=GΔ(α)lcoeff(G)G := G \frac{\Delta ( \alpha )}{\mbox{lcoeff} ( G )}
      • D′:=D′+I(X n)∏ α j∈E(α−α j)(G−D′(X 1,….,X n−1,α))D' := D' + \frac{I ( X_n )}{\prod_{\alpha_j \in E} ( \alpha - \alpha_j )} ( G - D' ( X_1, \ldots ., X_{n - 1}, \alpha ))
      • I:=I*(X n−α)I := I \ast ( X_n - \alpha )
      • e:=e+1e := e + 1 et ajouter α\alpha à EE
      • Si ee est strictement plus grand que le minimum des degrés partiels de PP et QQ par rapport à X nX_n, on pose D˜\tilde{D} la partie primitive de D′D' (vu comme polynôme à coefficients dans ℤ[X n]\mathbb{Z} [ X_n ]), on teste si PP et QQ sont divisibles par D˜\tilde{D}, si c’est le cas, on renvoie D=C(X n)D˜D = C ( X_n ) \tilde{D}

On observe que dans cet algorithme, on fait le test de divisibilite de D˜\tilde{D} par PP et QQ. En effet, même après avoir évalué en suffisamment de points, rien n’indique que tous ces points sont des points de bonne évaluation. En pratique cela reste extrêmement improbable. En pratique, on teste la divisibilité plus tôt, dès que D′D' n’est pas modifié par l’ajout d’un nouveau point à la liste des α j\alpha_j.

Il existe une variation de cet algorithme, appelé SPMOD (sparse modular), qui suppose que seuls les coefficients non nuls du pgcd en n−1n - 1 variables sont encore non nuls en nn variables (ce qui a de fortes chances d’être le cas). L’étape d’interpolation est alors remplacée par la résolution d’un sous-système d’un système de Vandermonde. Cette variation est intéressante si le nombre de coefficients non nuls en n−1n - 1 variables est petit devant le degré. Si elle échoue, on revient à l’interpolation dense.

Notons enfin qu’on peut appliquer cette méthode lorsque les coefficients de PP et QQ sont dans ℤ/nℤ\mathbb{Z} / n \mathbb{Z} mais il faut alors vérifier qu’on dispose de suffisamment de points d’interpolation. Ce qui en combinant avec l’algorithme modulaire à une variable donne un algorithme doublement modulaire pour calculer le pgcd de 2 polynômes à coefficients entiers. C’est cette méthode qu’utilise par exemple MuPAD (en essayant d’abord SPMOD puis l’interpolation dense).

Exemple:

Dans cet exemple, on donne FF et GG sous forme factorisée, le but étant de faire comprendre l’algorithme. En utilisation normale, on n’exécuterait cet algorithme que si FF et GG étaient développés.

P=((x+1)y+x 2+1)(y 2+xy+1),Q=((x+1)y+x 2+1)(y 2−xy−1)P = (( x + 1 ) y + x^2 + 1 ) ( y^2 + x y + 1 ), Q = (( x + 1 ) y + x^2 + 1 ) ( y^2 - x y - 1 ).

Prenons xx comme variable X 1X_1 et yy comme variable X 2X_2. Les coefficients dominants de PP et QQ sont respectivement yy et −y- y donc Δ=y\Delta = y.

En y=0y = 0, P(x,0)=x 2+1P ( x, 0 ) = x^2 + 1 n’est pas du bon degré.

En y=1y = 1, P(x,1)=(x+x 2+2)(x+2)P ( x, 1 ) = ( x + x^2 + 2 ) ( x + 2 ) et Q(x,1)=(x+x 2+2)(−x)Q ( x, 1 ) = ( x + x^2 + 2 ) ( - x ) sont du bon degré. Leur pgcd est G=x 2+x+2G = x^2 + x + 2, Δ(1)=1\Delta ( 1 ) = 1, donc D′=x 2+x+1D' = x^2 + x + 1. On teste la divisibilité de PP par D′D', le teste échoue.

En y=2y = 2, P(x,2)=(x 2+2x+3)(2x+5)P ( x, 2 ) = ( x^2 + 2 x + 3 ) ( 2 x + 5 ) et Q(x,2)=(x 2+2x+3)(−2x+3)Q ( x, 2 ) = ( x^2 + 2 x + 3 ) ( - 2 x + 3 ) donc G=x 2+2x+3G = x^2 + 2 x + 3, Δ(2)=2\Delta ( 2 ) = 2. On interpole: D′=x 2+x+2+y−12−1(2(x 2+2x+3)−(x 2+x+2))=y(x 2+3x+4)−(2x+2) D' = x^2 + x + 2 + \frac{y - 1}{2 - 1} ( 2 ( x^2 + 2 x + 3 ) - ( x^2 + x + 2 )) = y ( x^2 + 3 x + 4 ) - ( 2 x + 2 ) On teste la divisibilité de PP par D′D', le test échoue.

En y=3y = 3, P(x,3)=(x 2+3x+4)(3x+10)P ( x, 3 ) = ( x^2 + 3 x + 4 ) ( 3 x + 10 ) et Q(x,3)=(x 2+3x+4)(−3x+8)Q ( x, 3 ) = ( x^2 + 3 x + 4 ) ( - 3 x + 8 ) donc G=x 2+3x+4G = x^2 + 3 x + 4, Δ(3)=3\Delta ( 3 ) = 3. On interpole: D′ = y(x 2+3x+4)−(2x+2)+ (y−2)(y−1)(3−2)(3−1)(3(x 2+3x+4)−(3(x 2+3x+4)−(2x+2))) \begin{matrix} D' &= &y ( x^2 + 3 x + 4 ) - ( 2 x + 2 ) + \\ & & \frac{( y - 2 ) ( y - 1 )}{( 3 - 2 ) ( 3 - 1 )} \left( 3 ( x^2 + 3 x + 4 ) - ( 3 ( x^2 + 3 x + 4 ) - ( 2 x + 2 )) \right) \end{matrix} donc D′=y(x 2+3x+4)−(2x+2)+(y−2)(y−1)2(−2x−2)=x 2y+xy 2+y 2+y D' = y ( x^2 + 3 x + 4 ) - ( 2 x + 2 ) + \frac{( y - 2 ) ( y - 1 )}{2} ( - 2 x - 2 ) = x^2 y + x y^2 + y^2 + y On divise D′D' par son contenu et on trouve x 2+xy+y+1x^2 + x y + y + 1 qui est bien le pgcd de PP et QQ.

5.3.3  EZGCD.

Il s’agit d’une méthode pp-adique. On évalue toutes les variables sauf une, on calcule le pgcd en une variable et on remonte au pgcd variable par variable (EEZGCD) ou toutes les variables simultanément (EZGCD) par un lemme de Hensel. Il semble qu’il est plus efficace de remonter les variables séparément.

Soit donc FF et GG deux polynômes primitifs dépendant des variables X 1,…,X nX_1, \ldots, X_n de pgcd DD, on fixe une des variables qu’on appelera X 1X_1 dans la suite. Soient lcoeff(F)\mbox{lcoeff} ( F ) et lcoeff(G)\mbox{lcoeff} ( G ) les coefficients dominants de FF et GG par rapport à X 1X_1. On évalue FF et GG en un n−1n - 1 uplet bb tel que le degré de FF et GG par rapport à X 1X_1 soit conservé après evaluation en bb. On suppose que D b(X 1)=pgcd(F(b),G(b))D_b ( X_1 ) = \mbox{pgcd} ( F ( b ), G ( b )) a le même degré que D(b)D ( b ). On a donc l’égalité: (F*lcoeff(F))(b)=(D blcoeff(F(b))lcoeff(D b))*(F(b)D blcoeff(F)(b)lcoeff(F(b)D b)) ( F \ast \mbox{lcoeff} ( F )) ( b ) = \left( D_b \frac{\mbox{lcoeff} ( F ( b ))}{\mbox{lcoeff} ( D_b )} \right) \ast \left( \frac{F ( b )}{D_b} \frac{\mbox{lcoeff} ( F ) ( b )}{\mbox{lcoeff} ( \frac{F ( b )}{D_b} )} \right) et de même en remplaçant FF par GG.

Pour pouvoir lifter cette égalité (c’est-à-dire généraliser à plusieurs variables), il faut que D bD_b et F(b)D b\frac{F ( b )}{D_b} soient premiers entre eux. Sinon, on peut essayer de lifter l’égalité analogue avec GG. En général, on montre qu’il existe un entier jj tel que D bD_b et F(b)+jG(b)D b\frac{F ( b ) + j G ( b )}{D_b} soient premiers entre eux. En effet, sinon au moins un des facteurs irréductibles de D bD_b va diviser F(b)+jG(b)D b\frac{F ( b ) + j G ( b )}{D_b} pour deux valeurs distinctes de jj et va donc diviser à la fois F(b)D b\frac{F ( b )}{D_b} et G(b)D b\frac{G ( b )}{D_b} en contradiction avec la définition de D b=pgcd(F(b),G(b))D_b = \mbox{pgcd} ( F ( b ), G ( b )). On lifte alors l’égalité obtenue en remplaçant FF par (F+kG)( F + k G ) ci-dessus. Dans la suite, on suppose qu’on peut prendre j=0j = 0 pour alléger les notations.

On va aussi supposer que b=0b = 0. Sinon, on fait un changement d’origine sur les polynômes FF et GG pour que b=0b = 0 convienne, on calcule le pgcd et on lui applique la translation d’origine opposée.

On adopte ensuite la notation suivante: si kk est un entier, on dit qu’un polynôme PP est un O(k)O ( k ) si la valuation de PP vu comme polynôme en X 2,….,X nX_2, \ldots ., X_n à coefficients dans ℤ[X 1]\mathbb{Z} [ X_1 ] est supérieure ou égale à k k^{}, ou de manière équivalente si P(X 1,hX 2,….,hX n)=O h→0(h k) P ( X_1, h X_2, \ldots ., h X_n ) = O_{h \rightarrow 0} ( h^k ) L’égalité à lifter se réécrit donc: Flcoeff(F)=P 0Q 0+O(1) F \mbox{lcoeff} ( F ) = P_0 Q_0 + O ( 1 ) où P 0=P_0 =D blcoeff(F(b))lcoeff(D b)D_b \frac{\mbox{lcoeff} ( F ( b ))}{\mbox{lcoeff} ( D_b )} et Q 0=F(b)D blcoeff(F)(b)lcoeff(F(b)D b)Q_0 = \frac{F ( b )}{D_b} \frac{\mbox{lcoeff} ( F ) ( b )}{\mbox{lcoeff} ( \frac{F ( b )}{D_b} )} sont premiers entre eux et de degré 0 par rapport aux variables X 2,….,X nX_2, \ldots ., X_n. Cherchons P 1=O(1)P_1 = O ( 1 ) et Q 1=O(1)Q_1 = O ( 1 ) de degré 1 par rapport aux variables X 2,….,X nX_2, \ldots ., X_n tels que Flcoeff(F)=(P 0+P 1)(Q 0+Q 1)+O(2) F \mbox{lcoeff} ( F ) = ( P_0 + P_1 ) ( Q_0 + Q_1 ) + O ( 2 ) Il faut donc résoudre Flcoeff(F)−P 0Q 0=P 0Q 1+Q 0P 1+O(2) F \mbox{lcoeff} ( F ) - P_0 Q_0 = P_0 Q_1 + Q_0 P_1 + O ( 2 ) On peut alors appliquer l’identité de Bézout qui permet de déterminer des polynômes P 1P_1 et Q 1Q_1 satisfaisant l’égalité ci-dessus (avec comme reste O(2)O ( 2 ) nul) puisque P 0P_0 et Q 0Q_0 sont premiers entre eux. De plus, on choisit P 1P_1 et Q 1Q_1 tels que degre X 1P 1≤degre X 1(F)−degre (Q 0)=degre(P 0)\mbox{degre}_{X_1} P_1 \leq \mbox{degre}_{X_1} ( F ) - \mbox{degre}_{} ( Q_0 ) = \mbox{degre} ( P_0 ) et degre X 1(Q 1)≤degre(Q 0)\mbox{degre}_{X_1} ( Q_1 ) \leq \mbox{degre} ( Q_0 ) et lcoeff X 1(P 0+P 1)+O(2)=lcoeff X 1(Q 0+Q 1)+O(2)=lcoeff X 1(F)\mbox{lcoeff}_{X_1} ( P_0 + P_1 ) + O ( 2 ) = \mbox{lcoeff}_{X_1} ( Q_0 + Q_1 ) + O ( 2 ) = \mbox{lcoeff}_{X_1} ( F ). On tronque ensuite P 1P_1 et Q 1Q_1 en ne conservant que les termes de degré 1 par rapport à X 2,….,X nX_2, \ldots ., X_n.

On trouve de la même manière par récurrence P kP_k et Q kQ_k homogènes de degré kk par rapport à X 2,….,X kX_2, \ldots ., X_k, de degré par rapport à X 1X_1 respectivement inférieur aux degrés de Q 0Q_0 et de P 0P_0 et tels que Flcoeff(F)=(P 0+….+P k)(Q 0+….+Q k)+O(k+1)(11) F \mbox{lcoeff} ( F ) = ( P_0 + \ldots . + P_k ) ( Q_0 + \ldots . + Q_k ) + O ( k + 1 ) \qquad (11) et lcoeff(F)=lcoeff X 1(P 0+….+P k)+O(k+1)=lcoeff X 1(Q 0+….+Q k)+O(k+1)\mbox{lcoeff} ( F ) = \mbox{lcoeff}_{X_1} ( P_0 + \ldots . + P_k ) + O ( k + 1 ) = \mbox{lcoeff}_{X_1} ( Q_0 + \ldots . + Q_k ) + O ( k + 1 ).

Si on est bien en un point de bonne évaluation et si kk est plus grand que le degré total (par rapport aux variables X 2,….,X nX_2, \ldots ., X_n) du polynôme Flcoeff(F)F \mbox{lcoeff} ( F ) on va vérifier que P 0+….+P k=Dlcoeff(F)lcoeff(D)P_0 + \ldots . + P_k = D \frac{\mbox{lcoeff} ( F )}{\mbox{lcoeff} ( D )}. En effet, si on a deux suites de polynômes PP et P′P' et QQ et Q′Q' satisfaisant (11) avec les même termes de degré zéro P 0P_0 et Q 0Q_0, alors en prenant la différence, on obtient: (P 0+P 1…+P k)(Q 0+Q 1…+Q k)=(P 0+P 1′…+P k′)(Q 0+Q 1′…+Q k′)+O(k+1) ( P_0 + P_1 \ldots + P_k ) ( Q_0 + Q_1 \ldots + Q_k ) = ( P_0 + P_1' \ldots + P_k' ) ( Q_0 + Q_1' \ldots + Q_k' ) + O ( k + 1 ) On égale alors les termes homogènes de degré jj, pour j=1j = 1, on obtient P 0(Q 1−Q 1′)=Q 0(P 1−P 1′)P_0 ( Q_1 - Q_1' ) = Q_0 ( P_1 - P_1' ), donc Q 0Q_0 divise Q 1−Q 1′Q_1 - Q_1' qui est de degré strictement inférieur au degré de Q 0Q_0 par rapport à X 1X_1 (car on a l’inégalité large et les termes de plus haut degré sont égaux), donc Q 1=Q 1′Q_1 = Q_1' et P 1=P 1′P_1 = P_1'. On montre de la même manière que Q j=Q j′Q_j = Q_j' et P j=P j′P_j = P_j'. L’écriture est donc unique, c’est donc l’écriture en polynôme homogène de degré croissant de Dlcoeff(F)lcoeff(D)D \frac{\mbox{lcoeff} ( F )}{\mbox{lcoeff} ( D )} que l’on reconstruit.

Cet algorithme permet donc de reconstruire DD, il suffit de tester à chaque étape si P 0+….+P kP_0 + \ldots . + P_k divise Flcoeff(F)F \mbox{lcoeff} ( F ). On appelle cette méthode de remontée lemme de Hensel linéaire. Il existe une variante dite lemme de Hensel quadratique qui consiste à passer de O(k)O ( k ) à O(2k)O ( 2 k ). Elle nécessite toutefois un calcul supplémentaire, celui de l’identité de Bézout à O(2k)O ( 2 k ) près pour les polynômes P 0+….+P k−1P_0 + \ldots . + P_{k - 1} et Q 0+….+Q k−1Q_0 + \ldots . + Q_{k - 1}. Ce calcul se fait également par lifting.

Algorithme EZGCD (Hensel linéaire)

Arguments: 2 polynômes FF et GG à coefficients entiers et primitifs. Renvoie le pgcd de FF et GG ou false.

  1. Evaluer FF et GG en (X 2,….,X n)=(0,….,0)( X_2, \ldots ., X_n ) = ( 0, \ldots ., 0 ), vérifier que les coefficients dominants de FF et de GG ne s’annulent pas. Calculer le pgcd D bD_b de F(0)F ( 0 ) et de G(0)G ( 0 ). Prendre un autre point d’évaluation au hasard qui n’annule pas les coefficients dominants de FF et de GG et vérifier que le pgcd a le même degré que D bD_b. Sinon, renvoyer false (on peut aussi faire une translation d’origine de FF et de GG en un autre point mais cela diminue l’efficacité de l’algorithme).
  2. On note lcF\mbox{lcF} et lcG\mbox{lcG} les coefficients dominants de FF et de GG par rapport à X 1X_1.
  3. Si degre(F)≤degre(G)\mbox{degre} ( F ) \leq \mbox{degre} ( G ) et degre(D b)=degre(G)\mbox{degre} ( D_b ) = \mbox{degre} ( G ) et FF divise GG renvoyer FF
  4. Si degre(G)<degre(F)\mbox{degre} ( G ) &lt; \mbox{degre} ( F ) et degre(D b)=degre(F)\mbox{degre} ( D_b ) = \mbox{degre} ( F ) et GG divise FF renvoyer GG
  5. Si degre(F)=degre(D b)\mbox{degre} ( F ) = \mbox{degre} ( D_b ) ou si degre(G)=degre(D b)\mbox{degre} ( G ) = \mbox{degre} ( D_b ) renvoyer false
  6. Boucle infinie sur jj entier initialisé à 0, incrémenté de 1 à chaque itération: si pgcd(D b,F(0)+jG(0)D b)=C\mbox{pgcd} ( D_b, \frac{F ( 0 ) + j G ( 0 )}{D_b} ) = C constant, alors arrêter la boucle
  7. Lifter l’égalité (F+jG)(lcF+jlcG)(0)=(D b(lcF+jlcG)(0)lcoeff(D b))*….( F + j G ) ( \mbox{lcF} + j \mbox{lcG} ) ( 0 ) = \left( D_b \frac{( \mbox{lcF} + j \mbox{lcG} ) ( 0 )}{\mbox{lcoeff} ( D_{b )}} \right) \ast \ldots . par remontée de Hensel linéaire ou quadratique. Si le résultat est false, renvoyer false. Sinon renvoyer le premier polynôme du résultat divisé par son contenu vu comme polynôme en X 1X_1 à coefficients dans ℤ[X 2,….,X n]\mathbb{Z} [ X_2, \ldots ., X_n ].

Remontée de Hensel linéaire:

Arguments: FF un polynôme, lcF\mbox{lcF}=lcoeff(F)(F) son coefficient dominant, P 0P_0 un facteur de F(0)F ( 0 ) ayant comme coefficient dominant lcF(0)\mbox{lcF} ( 0 ) et dont le cofacteur Q 0Q_0 est premier avec P 0P_0.

Renvoie deux polynômes PP et QQ tels que FlcF=PQF \mbox{lcF} = P Q et P(0)=P 0P ( 0 ) = P_0 et lcoeff(P)=lcoeff(Q)=lcF\mbox{lcoeff} ( P ) = \mbox{lcoeff} ( Q ) = \mbox{lcF}.

  1. Soit G=FlcFG = F \mbox{lcF}, , Q 0=G(0)/P 0Q_0 = G ( 0 ) / P_0, P=P 0P = P_0, Q=Q 0Q = Q_0.
  2. Déterminer les deux polynômes UU et VV de l’identité de Bézout (tels que P 0U+Q 0V=dP_0 U + Q_0 V = d où dd est un entier).
  3. Boucle infinie avec un compteur kk initialisé à 1, incrémenté de 1 à chaque itération
    • Si k>degre X 2,….,X n(G)k &gt; \mbox{degre}_{X_2, \ldots ., X_n} ( G ), renvoyer false.
    • Si PP divise GG, renvoyer PP et G/PG / P.
    • Soit H=G−PQ=O(k)H = G - P Q = O ( k ). Soit u=UHdu = U \frac{H}{d} et v=VHdv = V \frac{H}{d}, on a P 0u+Q 0v=HP_0 u + Q_0 v = H
    • Remplacer vv par le reste de la division euclidienne de vv par P 0P_0 et uu par le reste de la division euclidienne de uu par Q 0Q_0. La somme des deux quotients est égale au quotient euclidien de HH par P 0Q 0P_0 Q_0, c’est-à-dire au coefficient dominant de HH divisé par le produit des coefficients dominants de P 0P_0 et Q 0Q_0 (qui sont égaux) donc on a l’égalité: P 0u+Q 0v=H−lcoeff(H)lcoeff(P 0) 2P 0Q 0 P_0 u + Q_0 v = H - \frac{\mbox{lcoeff} ( H )}{\mbox{lcoeff} ( P_0 )^2} P_0 Q_0
    • Soit α=(lcoeff(F)−lcoeff(P))/lcoeff(P 0)\alpha = ( \mbox{lcoeff} ( F ) - \mbox{lcoeff} ( P )) / \mbox{lcoeff} ( P_0 ) et β=(lcoeff(F)−lcoeff(Q))/lcoeff(P 0)\beta = ( \mbox{lcoeff} ( F ) - \mbox{lcoeff} ( Q )) / \mbox{lcoeff} ( P_0 ). On ajoute αP 0\alpha P_0 à vv, ainsi lcoeff(P +v)=lcoeff(F)+O(k+1)\mbox{lcoeff} ( P_{} + v ) = \mbox{lcoeff} ( F ) + O ( k + 1 ) et βQ 0\beta Q_0 à uu, ainsi lcoeff(Q +u)=lcoeff(F)+O(k+1)\mbox{lcoeff} ( Q_{} + u ) = \mbox{lcoeff} ( F ) + O ( k + 1 )

      Remarque: on montre alors que α+β=lcoeff(H)lcoeff(P 0Q 0)+O(k+1)\alpha + \beta = \frac{\mbox{lcoeff} ( H )}{\mbox{lcoeff} ( P_0 Q_0 )} + O ( k + 1 ) donc P 0u+Q 0v=H+O(k+1)P_0 u + Q_0 v = H + O ( k + 1 ) en utilisant les propriétés : lcoeff(F)=lcoeff(P)+O(k)=lcoeff(Q)+O(k)=lcoeff(P 0)+O(1) \mbox{lcoeff} ( F ) = \mbox{lcoeff} ( P ) + O ( k ) = \mbox{lcoeff} ( Q ) + O ( k ) = \mbox{lcoeff} ( P_0 ) + O ( 1 )

    • Réduire uu et vv en éliminant les termes de degré strictement supérieur à kk par rapport à X 2,….,X nX_2, \ldots ., X_n. S’il reste un coefficient non entier, renvoyer false
    • Remplacer PP par P+vP + v et QQ par Q+uQ + u, passer à l’itération suivante.

Exemple:

F=((x+1)y+x 2+1)(y 2+xy+1),G=((x+1)y+x 2+1)(y 2−xy−1)F = (( x + 1 ) y + x^2 + 1 ) ( y^2 + x y + 1 ), G = (( x + 1 ) y + x^2 + 1 ) ( y^2 - x y - 1 )

On a F(0,y)=(y+1)(y 2+1)F ( 0, y ) = ( y + 1 ) ( y^2 + 1 ) et G(0,y)=(y+1)(y 2−1)G ( 0, y ) = ( y + 1 ) ( y^2 - 1 ), le pgcd est donc D b=(y+1)D_b = ( y + 1 ). On remarque que D bD_b est premier avec le cofacteur de FF mais pas avec le cofacteur de GG. Si on évalue en un autre point, par exemple x=1x = 1, on trouve un pgcd D 1D_1 de même degré, donc 0 est vraissemblablement un bon point d’évaluation (ici on en est sûr puisque le pgcd de FF et GG se calcule à vue...). On a lcoeff(F)=x+1\mbox{lcoeff} ( F ) = x + 1, on va donc lifter G=((x+1)y+x 2+1)(y 2+xy+1)(x+1)=PQG = (( x + 1 ) y + x^2 + 1 ) ( y^2 + x y + 1 ) ( x + 1 ) = P Q où P 0=(y+1)P_0 = ( y + 1 ) et Q 0=(y 2+1)Q_0 = ( y^2 + 1 ).

On calcule les polynômes de l’identité de Bézout U=(1−y)U = ( 1 - y ) et V=1V = 1 avec d=2d = 2, puis à l’ordre k=1k = 1: H=G−P 0Q 0=(2y 3+2y 2+3y+1)x+O(2) H = G - P_0 Q_0 = ( 2 y^3 + 2 y^2 + 3 y + 1 ) x + O ( 2 ) donc u=reste(UH/d,Q 0)=xyu = \mbox{reste} ( U H / d, Q_0 ) = x y et v=reste(VH/d,P 0)=−xv = \mbox{reste} ( V H / d, P_0 ) = - x.

Donc Q 1=xy+αQ 0Q_1 = x y + \alpha Q_0 avec α=(x+1−1)/lcoeff(P 0)=x\alpha = ( x + 1 - 1 ) / \mbox{lcoeff} ( P_0 ) = x et Q 0+Q 1=(y 2+1)(x+1)+xyQ_0 + Q_1 = ( y^2 + 1 ) ( x + 1 ) + x y. De même, P 1=−x+βP 0P_1 = - x + \beta P_0, avec β=(x+1−1)/lcoeff(P 0)=x\beta = ( x + 1 - 1 ) / \mbox{lcoeff} ( P_0 ) = x donc P 0+P 1=(y+1)(x+1)−xP_0 + P_1 = ( y + 1 ) ( x + 1 ) - x. On remarque que P 0+P 1P_0 + P_1 et Q 0+Q 1Q_0 + Q_1 sont bien à O(2)O ( 2 ) près les facteurs de Flcoeff(F)F \mbox{lcoeff} ( F ): P=(x+1)y+x 2+1=P 0+P 1+O(2),Q=(x+1)(y 2+xy+1)=Q 0+Q 1+O(2) P = ( x + 1 ) y + x^2 + 1 = P_0 + P_1 + O ( 2 ), \ Q = ( x + 1 ) ( y^2 + x y + 1 ) = Q_0 + Q_1 + O ( 2 ) Une deuxième itération est nécessaire. On calcule H=G−(P 0+P 1)(Q 0+Q 1)=(2y 2+y+1)x 2+O(3) H = G - ( P_0 + P_1 ) ( Q_0 + Q_1 ) = ( 2 y^2 + y + 1 ) x^2 + O ( 3 ) puis reste(UH/d,Q 0)=yx 2\mbox{reste} ( U H / d, Q_0 ) = y x^2 et reste(VH/d,P 0)=x 2\mbox{reste} ( V H / d, P_0 ) = x^2. Ici les coefficients α\alpha et β\beta sont nuls car lcoeff(F)\mbox{lcoeff} ( F ) n’a pas de partie homogène de degré 2. On trouve alors P=P 0+P 1+P 2P = P_0 + P_1 + P_2 et Q=Q 0+Q 1+Q 2Q = Q_0 + Q_1 + Q_2. Pour calculer le pgcd, il suffit de calculer la partie primitive de PP vu comme polynôme en yy, ici c’est encore PP car le contenu de PP est 1 (remarque: pour QQ le contenu est x+1x + 1).
On trouve donc PP comme pgcd.

5.4  Quel algorithme choisir?

Il est toujours judicieux de faire une évaluation en quelques n−1n - 1 uplets pour traquer les pgcd triviaux. (E)EZGCD sera efficace si (0,...,0) est un point de bonne évaluation et si le nombre de remontées nécessaires pour le lemme de Hensel est petit donc pour les pgcd de petit degré, GCDMOD est aussi efficace si le degré du pgcd est petit. Le sous-résultant est efficace pour les pgcd de grand degré car il y a alors peu de divisions euclidiennes à effectuer et les coefficients n’ont pas trop le temps de croitre. SPMOD est intéressant pour les polynômes creux de pgcd non trivial creux. GCDHEU est intéressant pour les problèmes relativement petits.

Avec des machines multiprocesseurs, on a probablement intérêt à lancer en parallèle plusieurs algorithmes et à s’arrêter dès que l’un deux recontre le succès.

5.5  Pour en savoir plus.

Parmi les références citées dans le premier article, ce sont les livres de Knuth, H. Cohen, et Davenport-Siret-Tournier qui traitent des algorithmes de pgcd. On peut bien sûr consulter le source de son système de calcul formel lorsqu’il est disponible :

Sur le web on trouve quelques articles en lignes sur le sujet en cherchant les mots clefs GCDHEU, EZGCD, SPMOD sur un moteur de recherche, il y a par exemple une description un peu différente du pgcd heuristique sur:
www.inf.ethz.ch/personal/gonnet/CAII/HeuristicAlgorithms/node1.html
et un article de comparaison de ces algorithmes par Fateman et Liao (dont la référence bibliographique est Evaluation of the heuristic polynomial GCD. in: ISSAC pages 240–247, 1995). Quelques autres références :

Chapitre 6  Le résultant

6.1  Définition

Il s’agit d’un point de vue d’algèbre linéaire sur le PGCD. Considérons deux polynômes AA et BB à coefficients dans un corps, de degrés pp et qq et de pgcd DD et l’identité de Bézout correspondante : AU+BV=D(12) A U + B V =D \qquad (12) avec degré(U)<q(U)&lt;q et degré(V)<p(V)&lt;p. Imaginons qu’on cherche UU et VV en oubliant qu’il s’agit d’une identité de Bézout, en considérant simplement qu’il s’agit d’un problème d’algèbre linéaire de p+qp+q équations (obtenues en développant et en identifiant chaque puissance de XX de 0 à p+q−1p+q-1) à p+qp+q inconnues (les pp coefficients de VV et les qq coefficients de UU) On sait que AA et BB sont premiers entre eux si et seulement si ce problème d’algèbre linéaire a une solution pour D=1D=1. Donc si le déterminant du système est non nul, alors AA et BB sont premiers entre eux. Réciproquement si AA et BB sont premiers entre eux, le système a une solution unique non seulement avec comme second membre 11 mais avec n’importe quel polynôme de degré inférieur p+qp+q, donc le déterminant du système est non nul.

Définition:
On appelle résultant de AA et BB le déterminant de ce système (12). Il s’annule si et seulement si AA et BB ne sont pas premiers entre eux (ont au moins une racine commune). On appelle matrice de Sylvester la transposée de la matrice du système (les inconnues étant par ordre décroissant les coefficients de UU et VV) M(A,B)=(A a A a−1 … … A 0 0 … 0 0 A a … … A 1 A 0 … 0 ⋮ ⋮ 0 0 … A 0 B b B b−1 … B 0 0 0 … 0 ⋮ ⋮ 0 0 … B 0) M(A,B)=\left( \begin{array}{cccccccc} A_a & A_{a-1} & \ldots & \ldots & A_0 & 0 & \ldots & 0 \\ 0 & A_a & \ldots & \ldots & A_1 & A_0 & \ldots & 0 \\ \vdots & & & & & & & \vdots \\ 0 & 0 & \ldots & & & & & A_0 \\ B_b & B_{b-1} & \ldots & B_0 & 0 & 0 & \ldots & 0 \\ \vdots & & & & & & & \vdots \\ 0 & 0 & \ldots & & & & & B_0 \end{array} \right) (cette matrice contient b=b=degré(B)(B) lignes de coefficients du polynôme AA et a=a=degré(A)(A) lignes de coefficients du polynôme BB, les coefficients diagonaux sont les A aA_a et B 0B_0)

Remarques
Le résultant s’exprime polynomialement en fonction des coefficients des polynômes AA et BB donc aussi en fonction des coefficients dominants de AA et BB et des racines α 1,...,α a\alpha_1,..., \alpha_a de AA et β 1,...,β b\beta_1,...,\beta_b BB, or si on fait varier les racines de BB on annulera le résultant si l’une d’elle coincide avec une racine de AA, ceci montre que le résultant est divisible par le produit des différences des racines β j−α i\beta_j-\alpha_i de AA et BB. On montre que le quotient est A a bB b aA_a^b B_b^a en regardant le coefficient dominant du résultant en degré total par rapport aux β j\beta_j : dans le déterminant il faut prendre le produit des termes diagonaux pour avoir le degré maximal en les β j\beta_j. On peut aussi l’écrire sous la forme resultant(A,B)=A a b∏ i=1 aB(α i)\mbox{resultant}(A,B)=A_a^b \prod_{i=1}^a B(\alpha_i)

Soit PP un polynôme de degré nn et coefficient dominant p np_n. Le coefficient dominant de P′P' est np nnp_n, un multiple de p np_n, le résultant de PP et P′P' est donc divisible par p np_n, on appelle le quotient discriminant. En terme des racines r ir_i de PP, on a disc(P)=resultant(P,P′)p n=p n n−2∏ i=1 nP′(r i)=p n 2n−2∏ 1≤i<j≤n(r i−r j) 2\mbox{disc}(P)=\frac{\mbox{resultant}(P,P')}{p_n} =p_n^{n-2} \prod_{i=1}^n P'(r_i) = p_n^{2n-2} \prod_{1\leq i&lt;j\leq n} (r_i-r_j)^2 Ce résultat a un intérêt pour par exemple minorer à priori l’écart entre 2 racines d’un polynôme à coefficients entiers.

6.2  Applications

Revenons au cas où nos polynômes sont à coefficients dans un anneau contenu dans un corps, par e