Algorithmes de calcul formel et numérique

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

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. Deux versions HTML sont proposées, une utilisant www.mathjax.org pour un rendu plus fidèle des formules, mais qui nécessite un temps de chargement plus long, l’autre sans. Ces fichiers HTML ont été générés avec hevea.inria.fr de Luc Maranget, et le fork de Yannick Chevallier pour le support mathjax. Vous pouvez exécuter toutes les commandes interactives en cliquant sur le bouton Exécuter, le champ suivant est la console de l’interpréteur du logiciel de calcul formel.

Résumé: 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 poténomi&y). 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

Table des matières

1  Index, plan

L’index commence page suivante dans la version PDF.

Quelques conseils de lecture :

Index

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

3  Calculer sur ordinateur

3.1  Représentation des entiers

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

Preuve : On prend pour q le plus grand entier tel que abq ≥ 0.
Exemple :

La division euclidienne permet d’écrire un nombre entier, en utilisant une base b et des caractères pour représenter les entiers entre 0 et b−1. Nous écrivons les nombres entiers en base b=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 b (par exemple b=2), on effectue des divisions euclidienne successives par b 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 d à partir de son écriture dn...d0, on traduit les divisions euclidiennes successives en

 d=( ... ((dn b +dn−1)b + dn−2)...+d1)b+d0
 =dn bn + dn−1 bn−1 + ... + d0

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=24+23+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 a et b 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 [−231,231−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, 230+230+230+230 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−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 m et l’exposant e. Si on écrit les nombres en base b, la mantisse m s’écrira avec un nombre n fixé de chiffres (ou de bits en base 2), donc m ∈ [0,bn[. Soit un réel x représenté par

x=mbe,    m ∈ [0,bn

Si m∈ [0,bn−1[, alors on peut aussi écrire x=mbe−1 avec m′=mb ∈ [0,bn[, quelle écriture faut-il choisir? Intuitivement, on sent qu’il vaut mieux prendre m′ le plus grand possible, car cela augmente le nombre de chiffres significatifs (alors que des 0 au début de m ne sont pas significatifs). Ceci est confirmé par le calcul de l’erreur d’arrondi pour représenter un réel. En effet, si x est un réel non nul, il ne s’écrit pas forcément sous la forme mbe, on doit l’arrondir, par exemple au plus proche réel de la forme mbe. La distance de x à ce réel est inférieure ou égale à la moitié de la distance entre deux flottants consécutifs, mbe et (m+1)be, donc l’erreur d’arrondi est inférieure ou égale à be/2. Si on divise par xmbe, on obtient une erreur relative d’arrondi majorée par 1/(2m). On a donc intérêt à prendre m le plus grand possible pour minimiser cette erreur. Quitte à mulitplier par b, on peut toujours se ramener (sauf exceptions, cf. ci-dessous), à m ∈ [bn−1,bn[, on a alors une erreur d’arrondi relative majorée par

1
2bn−1

On appelle flottant normalisé un flottant tel que m ∈ [bn−1,bn[. Pour écrire un réel sous forme de flottant normalisé, on écrit le réel en base b, et on déplace la virgule pour avoir exactement n 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−53 (b=2 et n=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 e. Soit en effet em le plus petit exposant des flottants normalisés et considérons les flottants x=bem(1+1/b) et y=bem. Ces flottants sont distincts, mais leur différence n’est plus représentable par un flottant normalisé. Comme on ne souhaite pas représenter xy par 0, (puisque le test x==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,bn−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/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/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=± 1 sur 1 bit, la mantisse M ∈ [0,252[ sur 52 bits, et l’exposant e ∈ [0, 211[ sur 11 bits. Pour les nombres “normaux”, l’exposant est en fait compris entre 1 et 211−2, le nombre représenté est le rationnel

(1+
M
252
) 2e+1−210 

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 r dans [1,2[, ce qui permet de déterminer l’exposant e. Ensuite on écrit la représentation en base 2 de r−1 ∈ [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=0 ou e=211−1) ont été introduites pour représenter ± ∞ (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) avec x mal connu de calculer aussi f′(x), en effet :

f(x(1+h))= f(x)+xh f′(x) + O(h2)

l’erreur relative sur f(x) est donc au premier ordre multipliée par

|
xf′(x)
f(x)
|

Par exemple, l’erreur relative sur ex est au premier ordre l’erreur relative sur x multipliée par |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 |xx0| ≤ ε0 et si |yy0| ≤ ε1 alors par l’inégalité triangulaire (|a+b|≤ |a|+|b|), on a :

|(x+y)−(x0+y0)| ≤ |xx0| + | yy0 | ≤  ε0 + ε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 :
|xx0| ≤ ε 

Mais comme il faut représenter x0+y0 en machine, on doit ajouter une erreur d’arrondi, qui est proportionnelle à la valeur absolue de x0+y0 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
|xx0| ≤ ε |x0

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−53, à titre d’exercice, on pourra vérifier que cette erreur d’arrondi est majorée par l’erreur absolue de la somme x+y dès l’instant où x et y ont eux-même une erreur d’arrondi).

Lorsqu’on effectue une multiplication de deux nombres x,y dont les représentants x0,y0 sont non nuls, on a




xyx0 y0
x0 y0
 





x
x0
 
y
y0
 −1 





(
x
x0
−1)(
y
y0
 −1)+(
x
x0
−1)+(
y
y0
 −1) 


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−53 sur x0 y0.

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 x et y non nuls est de 2−53, alors l’erreur relative sur xy sera de

2−53 + 2−53 + 2−106 + 2−53 ≈ 3 × 2−53 

Lorsque l’erreur relative sur les données est grande devant 2−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 x est représenté par x0=1+2−52 avec une erreur d’arrondi de 2−53 et y par y0=−1 avec la même erreur d’arrondi, l’addition de x et y renvoie 2−52 avec une erreur absolue de 2 * 2−53 (ici il n’y a pas d’arrondi lorsqu’on fait la somme). C’est une erreur relative de 1 (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 

alors que

2.0−53+(2.0−53+1.0) → 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(xh))/(2h) que (f(x+h)−f(x))/h car le terme d’erreur est en O(h2) et non en O(h). Attention toutefois à ne pas prendre h trop petit, sinon x+h=x en flottants et même si x+hx, l’erreur absolue sur f(x+h)−f(xh) est (au moins) d’ordre ε |f(x)|, donc l’erreur relative est d’ordre ε/h |f(x)|. Par exemple pour h=1e-8 le reste est en O(h2) donc de l’ordre des erreurs d’arrondi mais l’erreur relative sur f(x+h)−f(xh) est d’ordre є/h largement supérieure (en flottants double-précision). On choisira plutôt h tel que є/h soit proche de h2, 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)) en x=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 n 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

X est la variable aléatoire somme des n erreurs, α l’erreur et nσ2 la variance de la somme n erreurs supposées indépendantes, cette probabilité tend vers 0 pour n grand si α est d’ordre n, et ne tend pas vers 0 si α est de l’ordre de √n. Exemple : somme de n=400 nombres répartis sur [−1,1] selon la loi uniforme (représentant des erreurs), on divise par √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 n réels ∑j xj, 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 є |x1+x2| puis є |x1+x2+x3| puis ... є |x1+x2+...+xn|, que l’on peut majorer par

є ((n−1)|x1|+(n−2)|x2|+...+|xn||)

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;
}:;


En effet, c devrait valoir 0 sans erreurs d’arrondis, avec les erreurs d’arrondis, on a le premier calcul s+x qui donnera une erreur opposée à celui du calcul de s à la ligne suivante, le 2ième calcul effectué (s+x)−s donne une erreur absolue en є |x| au pire (car c’est une erreur relative par rapport à (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+c qui est au premier ordre au pire en O(є ∑|xi|), bien meilleure que la majoration є ((n−1)|x1|+(n−2)|x2|+...+|xn||) calculée précédemment.

Par exemple
à comparer avec
(le calcul de S 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]. 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 en fonction du signe, puis utilise la monotonie de P+ et P sur + et respectivement.

L’arithmétique d’intervalle dans 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] et x2, x ∈ [−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,232−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) lorsqu’on assigne une valeur à x. Dans cet exemple, on voit qu’il faut d’abord remplacer x 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/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) par leur expression en fonction de t=tan(x/2), on est ainsi ramené à une fraction rationnelle en t 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)6 avec sa dérivée (les deux polynômes sont premiers entre eux) :

7 (x+1)6−6 (x−1)5  
162
49
  x5+
−390
49
  x4+
1060
49
  x3+
−780
49
  x2+
474
49
  x+
−78
49
  
157780
729
  x4+
−507640
2187
  x3+
290864
729
  x2+
−101528
729
  x+
28028
729
  
1
49
  (
1400328
2645
  x3+
−732888
2645
  x2+
1133352
3703
  x+
−732888
18515
)
  
1
2187
  (
2161816376832
4669921
  x2+
−555436846944
4669921
  x+
301917024864
4669921
)
  
1
907235
  (
469345063045455
129411872
  x+
−47641670106615
129411872
)
  
5497465490623352995840
209648836272383412129

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 (t2−1)x=t−1 sera x=1/(t+1)) ou choisissent une branche pour les fonctions possédant un point de branchement (par exemple pour résoudre x2=t en fonction de t). 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 (vus comme un couple de polynômes P,Q, où Q est un polynome minimal irréductible à coefficients entiers, autrement dit P,Q vaut P(α) où Q(α)=0), ou des éléments d’un corps fini (comme ci-dessus, mais ici Q est à coefficients dans /p avec p 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+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+c est obtenu par (a+b)+c ou par 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 (256255−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 p-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 à son équivalent dans Z/n pour une ou plusieurs valeurs de n, nombre premier. Le calcul dans /n 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 pour plusieurs nombres premiers n, 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 p-adiques commencent de manière identique par un calcul dans /n, on augmente ensuite la précision de la solution en la «liftant»de /nk vers /nk+1 ou vers /n2k (lift linéaire ou lift quadratique), on s’arrête lorsque k 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)=ad Xd + ⋯ + a0 à coefficients entiers ou polynomiaux, avec ad et a0 non nuls. L’algorithme générique (assez connu) consiste à chercher les diviseurs de a0 et de ad et à tester toutes les fractions de ces diviseurs, on montre en effet aisément que si X=p/q fraction irréductible est racine de P alors q divise ad et p divise a0. Cet algorithme est très inefficace si ad ou a0 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 P sont entiers, la recherche précédente revient à trouver un facteur à coefficients entiers qXp de P, on peut donc réduire le problème modulo un entier premier n qui ne divise pas ad : si un tel facteur existe dans alors ce facteur (réduit modulo n) est un facteur de P dans /n donc P admet une racine dans /n (puisque q est inversible modulo n car on a choisi n premier ne divisant pas ad). On évalue maintenant P en les n éléments de /n. S’il n’y a pas de 0, alors P n’a pas de racine rationnelle. S’il y a des racines, on va les lifter de /nk dans /n2k.

On suppose donc que pour k≥ 1, il existe un entier pk tel que

P(pk)=0 (mod nk ) 

Il s’agit de trouver un entier x tel que pk+1=pk+nk x vérifie

P(pk+1)=0 (mod n2k ) 

On applique la formule de Taylor à l’ordre 1 pour P en pk, le reste est nul modulo n2k, donc :

P(pk)+ nk  x P′(pk)=0 (mod n2k ) 

soit finalement :

x=−
P(pk)
nk
  ( P′(pk) (mod nk )) −1 

On reconnaît au passage la méthode de Newton, pour qu’elle fonctionne il suffit que P′(pk) ≠ 0 (mod n ) ce qui permet de l’inverser modulo nk (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 P est sans facteur multiple dans . Ceci n’entraîne pas forcément qu’il le reste dans /n ce qui crée une contrainte supplémentaire sur le choix de n, à savoir que P et P′ restent premier entre eux dans /n (il existe forcément de tels n, par exemple n premier plus grand que le plus grand entier intervenant dans le calcul du PGCD de P et P′ dans ).

Reste donc à revenir dans à partir d’une racine pk dans /(nk ) (où on peut choisir k). On va maintenant utiliser la représentation modulaire symétrique : on prend comme représentant modulaire d’un entier z dans /nk l’unique entier congru à z modulo n qui est strictement compris entre −nk/2 et nk/2 (si n est pair, la deuxième inégalité est choisie large).

Si qXp est un facteur de P, alors adXad/qp est encore un facteur de P (le quotient de P par adXad/qp est à coefficients rationnels mais le facteur est à coefficients entiers). Si on a choisi k tel que nk>2|ad a0|, l’écriture en représentation modulaire symétrique de adXad/qp est inchangée, en effet on a des estimations à priori sur les entiers p et q : |q|≤ |ad| et |p| ≤ |a0| puisque q divise ad et p divise a0. Comme adXad/qp est égal à ad(Xpk) dans /(nk ), il nous suffit d’écrire en représentation modulaire symétrique ad(Xpk)=ad Xp′. Pour conclure, on sait que ad Xp′ est un multiple entier de qXp. On divise donc le facteur ad Xp′ par le pgcd de ad et p′ et on teste la divisibilité de P par ce facteur réduit.

Exemple
Considérons le polynôme 2 X3X2X−3 qui est sans facteur carré. On ne peut pas choisir n=2 car on réduirait le degré, pour n=3, on a P′=X−1 qui est facteur de P, pour n=5, P′=6X2−2X−1, on vérifie que P et 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 P. Seul -1 est racine modulo 5 (P(−1)=−5), on va maintenant lifter p1=−1.

L’estimation à priori est 2|ad||a0|=12 donc k=2 (52=25>12), une itération suffira. On a P′(−1)=7, l’inverse de P′(−1) (mod 5 ) est -2 donc:

x= −
P(−1)
5
 (−2) = −(−1)  (−2)=−2 

et p2=−1+5×(−2)=−11 est racine de P dans /25. On calcule ensuite ad(Xpk)=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−3, ici il divise P donc P admet un unique facteur entier de degré 1 qui est 2X−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 p 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  Quelques algorithmes d’arithmétique de base.

3.9.1  Exemple de multiplication rapide : l’algorithme de Karatsuba

Soient P, Q deux polynômes de degrés strictement inférieur à 2n. On suppose que le cout d’une opération arithmétique dans le corps des coefficients vaut 1 et on néglige les autres opérations (on suppose par exemple que le corps des coefficients est un corps fini). On écrit

P=A+xn B,    Q=C+xn D

avec A,B,C,D de degrés strictement inférieur à n, on a alors :

P Q = AC + xn(AD+BC)+x2n BD

Il y a 4 produits de polynômes de degrés <n, mais au prix d’additions intermédiaires, on peut se ramener à 3 produits, en effet

(A+B)(C+D)−ACBD = AD+BC

donc pour calculer le cofacteur de xn il suffit de soustraire à (A+B)(C+D) les produits AC et BD que l’on calcule par ailleurs. Soit 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)+ 8n

où 8n représente le nombre d’additions ou de soustractions pour former A+B, C+D, soustraire AC et BD, et tenir compte des "retenues" (les termes de degré ≥ n de AC se combinent avec ceux de degré <2n de AD+BC et les termes de degré < 3n de x2nBD avec ceux de degré ≥ 2n de AD+BC). On en déduit

un=M(2n),    un+1=3un+8 × 2n 

cette récurrence se résoud facilement par la commande
rsolve(u(n+1)=3*u(n)+8*2^n,u(n),u(0)=1)
qui donne M(2n)=un=−8· 2n+9· 3n.

Asymptotiquement, M(2n) ≈ 9· 3n ce qui est bien meilleur que la multiplication naive en 2 · 4n, mais pour de petites valeurs de n, la multiplication naive est plus rapide, on utilise Karatsuba (récursivement) uniquement pour des valeurs de n suffisamment grandes (théoriquement lorsque 8n, le surcout dû aux additions est plus petit que la multiplication économisée, soit 8n<2n2 soit n>4, en pratique plutôt pour n de l’ordre de quelques dizaines selon les implémentations, car nous n’avons tenu compte que des opérations arithmétiques).

3.9.2  Calcul de la racine carrée entière

Étant donné un entier N, il s’agit de déterminer le plus grand entier n tel que n2N, n est la racine carrée de N. On choisit une base b par exemple b=10 pour un humain ou une puissance de 2 pour une machine, et on écrit N en base b, 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 n à 0 et son carré c à 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 n par b et c par b2 (c’est un simple décalage de l’écriture en ajoutant un ou deux zéros). Puis on ajoute les nombres impairs successifs 2n+1, (2n+1)+2, ... à c tant que l’on est inférieur à N tronqué au bloc. Le nombre d’impairs successifs ajouté est ajouté à n. En pratique, il suffit de conserver Nc tronqué et de lui retrancher les impairs successifs.

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

Exercice : calculer la quatrième décimale de √2 de cette manière.

La complexité de cet algorithme est en O(logb(N)2). En effet, pour calculer un chiffre il faut faire un nombre de soustraction au plus égal à b, ces soustractions ayant au plus le nombre de chiffres de N en base b. (On peut accélérer le calcul à la manière de Karatsuba en choisissant une base b puissance de 2 (ou 10) de l’ordre de √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.9.3  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 ab>0 deux entiers, on pose

(Ln)    a un − b vn = (−1)n rn 

r0=a, r1=b et rn+2 est le reste de la division euclidienne de rn par rn+1 (qn+2 le quotient), u0=1, u1=0, v0=0,v1=1. Comme précedemment, chaque ligne s’obtient par combinaison linéaire des deux précédentes, mais cette fois avec une addition

Ln+2=Ln+qn+2 Ln+1

ce qui se traduit par :

un+2=un+qn+2 un+1,    vn+2=vn+qn+2 vn+1

Les suites un et vn sont alors strictement croissantes (à partir du rang 1 pour un). Au rang k du dernier reste non nul on a :

a uk − b vk = (−1)k rk,    rk=d=gcd(a,b)

et au rang suivant :

auk+1 −b vk+1=0

On montre par récurrence que

vn rn+1 + vn+1 rn=a

et une relation analogue pour un, on en déduit alors que vk+1=a/d et uk+1=b/d (ce sont les cofacteurs du PPCM de a et b), en particulier les coefficients de Bézout vérifient uk<b et vk<a.

On va aussi voir que un+2/vn+2 est la n-ième réduite du développement en fractions continues de a/b (donc les coefficients de Bézout se lisent sur l’avant-dernière réduite). On introduit la notation

[a0,a1,..,an] =a0+
1
a1+
1
a2+
...
ak

pour a0 ≥ 0, a1>0, ..., an>0. On a alors :

a
b
=[q2,q3,..,qk]

En effet :

a
b
r0
r1
=q2 +
r2
r1
 = q2 +
1
r1
r2
 = ...

D’autre part, on montre par récurrence sur n≥ 1 que si x>0

[q2,..., qn,x]=
vnx+vn−1
unx+un−1

en effet au rang n=1

[x]=x=
v1 x + v0
u1 x+u0 

et pour l’induction :

  [q2,..., qn,x]=
[q2,..., qn−1,qn+
1
x
 =
vn−1(qn+1/x)+vn−2
un−1(qn+1/x)+un−2
 
 =
x(vn−1qn+vn−2)+vn−1
x(un−1qn+un−2)+un−1
 
 =
vnx+vn−1
unx+un−1

Donc au rang n−1 et pour x=qn, on obtient

[q2,..., qn]=
vn+1
un+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≥0. On initialise y=x et la liste des ap à vide. Puis on fait une boucle : on ajoute la partie entière de y à la liste, on calcule la partie fractionnaire de y, si elle est nulle on s’arrête (dans ce cas x∈ ), sinon on stocke dans y l’inverse de cette partie fractionnaire et on recommence. On note classiquement :

     
h−2=0, h−1=1, hp=ap hp−1+hp−2    (1)
k−2=1, k−1=0, kp=ap kp−1+kp−2     (2)

On a h0=a0, h1=a1 a0+1, k0=1, k1=a1. Les suites hp et kp sont donc positives et strictement croissantes pour p ≥ 1, puisque pour p ≥ 1, ap≥ 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 :

  hp kp−1 − hp−1kp=(−1)p+1     (3)

On montre aussi comme ci-dessus :

[a0,...,ap−1,y]=
yhp−1+hp−2
ykp−1+kp−2

On définit xp par x=[a0,...,ap−1,xp], en faisant y=xp on a alors x=xphp−1+hp−2/xp kp−1+kp−2 ce qui donne xp en fonction de x et

ap=floor


− 
xkp−2hp−2
xkp−1hp−1
 


En faisant y=ap on obtient [a0,...,ap]=hp/kp. On montre ensuite que les suites (hp/kp) pour les indices pairs et impairs sont deux suites adjacentes qui convergent vers x, et on a

 
hp
kp
 − 
hp−1
kp−1
 = 
(−1)p−1
kp kp−1
    (4)

En effet, la dernière égalité est une conséquence immédiate de (3), la croissance ou décroissance des suites d’indice pair ou impair s’en déduit en ajoutant (4) au cran suivant. La convergence vient de la limite infinie de kp en l’infini. On a donc

x=a0+
p=0
(−1)p−1
kp kp+1
,     
1
kp(kp+kp+1)
 ≤ |x
hp
kp
| ≤ 
1
kp kp+1

La convergence est d’autant plus rapide que les kp tendent rapidement vers l’infini, donc si les ap sont plus grands que 1. La convergence la plus lente correspond au cas où tous les ap=1 cas du nombre d’or, ou à partir d’un certain rang (nombre de Q[√5]).

3.9.4  La puissance rapide itérative

Pour calculer ak (mod n ), on décompose k en base 2

k=
J
j=0
 kj 2j,    ak = 
J
j=0
 akj 2j  = 
 
j/kj ≠ 0
 a2j 

On initialise une variable B à 1, B vaudra ak (mod n ) en fin de calcul, on initialise une variable k à k. On calcule dans une boucle les carrés successifs de a (mod n ) que l’on stocke dans une variable A (A vaudra donc successivement a (mod n ), a2 (mod n ), a4 (mod n ), ...) et simultanément on teste si kj vaut 1 en prenant le reste de la division par 2 de k (dans ce cas on multuplie B par A modulo n), 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;
}




En mode pas à pas :


3.10  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.11  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) où a=10 000!, 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 an (mod m ) par la fonction powmod et la méthode prendre le reste modulo m après avoir calculé an.

    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 b les puissances successives a20 (mod m ),a21 (mod m ), ..., a2k (mod m ), ..., on forme an (mod n ) en prenant le produit modulo m 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 b et du bit correspondant se font dans une même boucle).
  3. Déterminer un entier c tel que c=1 (mod 3 ), c=3 (mod 5 ), c=5 (mod 7 ) et c=2 (mod 11 ).
  4. Calculez dans /11
    10
    a=0
     (xa)

  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 x2yy est un élément quelconque du corps fini. Comparer avec la valeur de √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 x, l’expression 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. Voici un programme qui calcule la base utilisée pour représenter les flottants.
    Base():={
      local A,B;
      A:=1.0; B:=1.0;
      while evalf(evalf(A+1.0)-A)-1.0=0.0 do  A:=2*A; od;
      while evalf(evalf(A+B)-A)-B<>0 do B:=B+1; od;
      return B;
    } :;
    

    Testez-le
    et expliquez.
  11. Déterminer le plus grand réel positif x de la forme 2n (n entier) tel que (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) avec 30 chiffres significatifs, puis sa partie fractionnaire. Proposez une commande permettant de décider si a est un entier.
  13. Déterminer la valeur et le signe de la fraction rationnelle
    F(x,y)= 
    1335
    4
     y6 + x2 (11x2 y2y6 −121y4−2) + 
    11
    2
     y8 + 
    x
    2y
    en x=77617 et y=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 par exemple pour k=64 ? Calculer le nombre de termes dans le développement de (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) = an Xn + ... + a0 
    en un point. On pose b0=P(α ) et on écrit :
    P(X)−b0=(X−α )Q(X
    où :
    Q(X) = bn Xn−1 + ... +b2 X + b1 
    On calcule alors par ordre décroissant bn, bn−1, ..., b0.
    1. Donner bn en fonction de an puis pour in−1, bi en fonction de ai et bi+1. Indiquez le détail des calculs pour P(X)=X3−2X+5 et une valeur de α 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 α et le programme renverra P(α ). (On pourra aussi renvoyer les coefficients de Q).
    3. En utilisant cette fonction, écrire une fonction qui calcule le développement de Taylor complet d’un polynôme en un point.

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 2127), 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, c et m on considère la suite

un+1=aun+c (mod m ) 

où on choisit (par exemple) comme représentant de un le reste de la division euclidienne par m. La valeur de u0 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≠ 1, le cas a=1 n’est pas très intéressant. On a alors :

un=an u0 + 
an−1
a−1
 c (mod m )

On cherche à réaliser une période la plus grande possible idéalement m, mais m−1 peut fort bien convenir, et c’est possible si m est premier en choisissant a générateur du groupe cyclique, car on a alors a≠ 1 (mod m ) et :

un=an (u0 + 
c
a−1
) − 
c
a−1
  (mod m )

donc la suite est stationnaire ou prend toutes les valeurs sauf − c/a−1 .

Exemple : choisir pour m 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 un ont une périodicité très (trop) petite. Il est alors intéressant de prendre m=2k ± 1, parce que la division euclidienne par m peut se coder efficacement en base 2, on divise par 2k (décalage de k bits) et on ajuste x=(2k ± 1)q+r=2k q + (r ± q). Ainsi pour k=4 et m=24+1=17, m est premier. On peut construire une suite de période 16 en choisissant a générateur de (/17)*, par exemple a=3 et c=2 donne la suite 0,2,8,9,12,4,14,10,15,13,7,6,3,11,1,5.

On a le :

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

On observe d’abord que vouloir la périodicité maximale revient à pouvoir supposer que u0=0. Il est donc nécessaire d’avoir c et m premiers entre eux, sinon tous les un sont multiples du pgcd de c et m. Ensuite, on pose m=∏piri la décomposition en facteurs premiers de m et on raisonne modulo chaque premier (par le lemme chinois, la périodicité est le PPCM des périodicités modulo chaque piri). Si a≠ 1 (mod p )i alors a−1 est inversible modulo pi donc modulo piri on a

un=an (u0 + 
c
a−1
) + 
c
a−1
  

et la valeur −c/(a−1) ne peut pas être atteinte (ou alors la suite est stationnaire). Donc a−1 doit être divisible par tous les facteurs premiers de m pour avoir la périodicité maximale. Réciproquement, il faut trouver le premier ordre n tel que (an−1)/(a−1)=0 (mod pr ). On pose a=b+1, on a

an−1
a−1
=
(b+1)n−1
b
 = 
n
k=1

n
 
k

bk−1 = n +
n(n−1)
2
b +... 

On sait que b=a−1 est un multiple de p, disons b=qp, on en déduit que pour n=pr, on a bien (an−1)/(a−1)=0 (mod pr ), alors que pour n=pr−1 et p≠ 2, (an−1)/(a−1)=n (mod pr ) ≠ 0. Le même calcul pour p=2 (prise en compte de la division par 2 de n(n−1)) donne la condition b=a−1 est multiple de 4 si m l’est.

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

Exemple : m=231−1 est premier, on peut donc construire un générateur congruentiel de période m−1 en choisissant a générateur de /m*. Pour en trouver un, on peut tester a pris au hasard et voir si am−1/j ≠ 1 (mod m ) pour tous les diviseurs premiers de m−1. Par exemple





initialise l’état du générateur. Un appel à r() renvoie un entier entre 1 et m−1, pour avoir un g’enérateur pseudo-aléatoire selon la loi uniforme sur ]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 à k éléments

Au lieu d’une récurrence uk+1=auk+c on conserve en mémoire k+1 valeurs successives de la suite et on calcule

un+k+1 = a0 un+...+akun+k (mod p )

Si on note Un le vecteur (un,...,un+k) et A la matrice companion du polynôme a0+a1x+...+akxk, on a Un+1=AUn. Rechercher un générateur de période maximale revient à chercher A d’ordre le plus grand possible, donc les valeurs propres de A, i.e. les racines de P, doivent être racines de l’unité d’ordre le plus grand possible donc pk−1. Ce que l’on peut faire en construire un polynôme P 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 :

5  Le PGCD de polynômes.

Comme on l’a remarqué dans le premier article, l’algorithme d’Euclide est inefficace pour calculer le pgcd de deux polynômes à coefficients entiers. On va présenter ici les algorithmes utilisés habituellement par les systèmes de calcul formel: sous-résultant (PRS), modulaire (GCDMOD), p-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 et Q deux polynômes à coefficients dans un corps. Le pgcd de P et Q n’est défini qu’à une constante près. Mais lorsque les coefficients de P et Q sont dans un anneau euclidien comme par exemple ℤ ou ℤ[ i ], on appellera pgcd de P et Q un polynôme D tel que P / D et Q / D soient encore à coefficients dans l’anneau, et que D soit optimal, c’est-à-dire que si un multiple µ D de D vérifie P / µ D et Q / µ D sont à coefficients dans l’anneau, alors µ 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 = 4 X2 − 4 et Q = 6 X2 + 12 X + 6. Le polynôme X + 1 est un pgcd de P et Q puisqu’il est de degré maximal divisant P et Q mais le pgcd de P et Q est 2 ( X + 1 ). Remarquons qu’avec notre définition − 2 ( X + 1 ) convient aussi. Par convention on appelera pgcd dans [X] le polynôme ayant un coefficient dominant positif.

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

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

Proposition : le contenu de AB est le produit des contenus de A et de B.
En effet le produit des contenus de A et B divise le contenu de AB, et A/contenu de A est primitif, B/contenu de B est primitif donc le produit l’est,

Proposition : Si A et B sont primitifs et si B divise A dans [X] alors A/B ∈ [X].

Preuve : Soit Q=A/B ∈ [X]. Soit q ∈ le PPCM des dénominateurs des coefficients de Q et notons P=qQ ∈ [X]. On a A=BQ donc qA=BP donc le contenu de qA est le produit du contenu de B par celui de P, donc le contenu de P=qQ est q, donc Q ∈ [X].

Donc le PGCD de A et B, polynômes primitifs de [X] est obtenu en prenant un PGCD de A et B dans [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 ( PQ ) = pgcd ( c ( P ), c ( Q )) pgcd ( pp ( P ), 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 R de la division euclidienne du polynôme P par Q, on prend le reste de la division de P qδ + 1 par Q, où q désigne le coefficient dominant de Q et δ la différence entre le degré de P et de Q.

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 )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 R 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 P et Q primitifs. Valeur de retour: le pgcd de P et Q.

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

Boucle à effectuer tant que Q est non nul:

Si on sort normalement de la boucle, Q est nul, on renvoie donc la partie primitive de P 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 P et Q (celle dont le déterminant est appelé résultant de P et Q) 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 R divisés par g hδ 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=R0, Q=R1, R2... 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, il s’agit par exemple de montrer que le reste R3 de Q=R1 par R2 est divisible par le carré du coefficient dominant de Q=R1. Voyons comment on obtient les coefficients de R3 à partir de la matrice de Sylvester de P et Q. Prenons la sous-matrice constituée des 2 premières lignes de P et des 3 premières lignes de Q et réduisons-la sous forme échelonnée sans introduire de dénominateur.







pnpn−1pn−2pn−3... 
0pnpn−1pn−2... 
qn−1qn−2qn−3qn−4... 
0qn−1qn−2qn−3... 
00qn−1qn−2...  
 





On effectue L1qn−1 L1pn L3 et L2qn−1 L2pn L4, ce qui correspond à l’élimination du terme en x du quotient de P par Q







0qn−1 pn−1 − pn qn−2......... 
00qn−1 pn−1 − pn qn−2...... 
qn−1qn−2qn−3qn−4... 
0qn−1qn−2qn−3... 
00qn−1qn−2...  
 





on effectue ensuite

L1qn−1 L1 − (qn−1 pn−1 − pn qn−2)  L4 
L2qn−1 L2 − (qn−1 pn−1 − pn qn−2)  L5

ce qui correspond à l’élimination du terme constant du quotient de P par Q, on obtient







00r2,n−2...... 
000r2,n−2... 
qn−1qn−2qn−3qn−4... 
0qn−1qn−2qn−3... 
00qn−1qn−2...  
 





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 Q et R2




 qn−1qn−2... 
 r2,n−2...... 
 0r2,n−2... 
 


On recommence les opérations de réduction de cette sous-matrice correspondant à la division euclidienne de Q par R2, on obtient




 00r3,n−3 
 r2,n−2...... 
 0r2,n−2... 
 


puis après suppression des colonnes 1 et 2 et des lignes 2 et 3 la ligne des coefficients de R3.

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 R3 (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 L1 par r2,n−2, mais on doit ensuite diviser le déterminant par r2,n−22 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 qn−14 (2 lignes ont été multipliées 2 fois par qn−1) et division par qn−12 (élimination des colonnes 1 et 2). Au final, tout coefficient de R3 est égal au produit d’un déterminant 5,5 extrait de la matrice de Sylvester de P et Q par qn−12, qui est justement le coefficient “magique” par lequel on divise le reste de R1=Q par R2 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 P et Q en un entier z et à extraire des informations du pgcd g des entiers P ( z ) et Q ( z ). Il faut donc un moyen de remonter de l’entier g à un polynôme G tel que G ( z ) = g. La méthode consiste à écrire en base z l’entier g, avec une particularité dans les divisions euclidiennes successives on utilise le reste symétrique (compris entre − z / 2 et z / 2). Cette écriture donne les coefficients d’un polynôme G unique. On extrait ensuite la partie primitive de ce polynôme G. Lorsque z est assez grand par rapport aux coefficients des polynômes P et Q, si pp ( G ) divise P et Q, on va montrer que le pgcd de P et de Q est D = pp ( G ).

On remarque tout d’abord que d : = D ( z ) divise g. En effet D divise P et Q donc pour tout entier (ou entier de Gauss) z, D ( z ) divise P ( z ) et Q ( z ). Il existe donc une constante a telle que

g = a d 

On a aussi pp ( G ) divise D. Il existe donc un polynôme C tel que :

D = pp ( G ) C 

Nous devons prouver que C 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 z, on obtient

d = 
g
c ( G )
 C ( z ) 

Finalement

1 = 
a
c ( G )
 C ( z ) 

La procédure de construction de G nous donne une majoration de ces coefficients par | z | / 2, donc de c ( G ) par | z | / 2, donc C ( z ) divise un entier de module plus petit que | z | / 2, donc

C ( z ) | ⩽ 
z |
2
 

On considère maintenant les racines complexes z1, … ., zn du polynôme C (il en existe au moins une puisqu’on a supposé C non constant). On a:

C ( X ) = cn ( X − z1 ) … . ( X − zn ) 

Donc, comme cn est un entier (ou entier de Gauss) non nul, sa norme est supérieure ou égale à 1 et :

C ( z ) | ⩾ 
n
j = 1
 ( | z | − | zj | ) 

Il nous reste à majorer les racines de C pour minorer | C ( z ) |. Comme C divise D il divise P et Q donc les racines de C sont des racines communes à P et Q. On va appliquer le:

Lemme 5   Soit x une racine complexe d’un polynôme P = an Xn + … . + a0.

Alors 

| x | < | P |/| an | + 1, | P | = max0 ⩽ in ( | ai | ) 

Application du lemme à C(X) : on a 1/|cn|≤ 1 donc si on a choisi z tel que | z | ⩾ 2 min( | P |, | Q | ) + 2, alors pour tout j, | zj | < | z | / 2 donc

C ( z ) | > 


z |
2
 


n



 
 

qui contredit notre majoration de | C ( z ) |.

Théorème 6   Soit P et Q deux polynômes à coefficients entiers. On choisit un entier z tel que | z | ⩾ 2 min( | P |, | Q | ) + 2, si la partie primitive du polynôme G reconstruit à partir du pgcd de P ( z ) etQ(z) par écriture en base z (avec comme reste euclidien le reste symétrique) divise P et Q alors c’est le pgcd de P et Q.

Pour finir la démonstration du théorème, il nous faut encore montrer le lemme. On a

− an xn = an − 1 xn − 1 + … . + a0 

Donc

an | | x |n ⩽ | P | ( 1 + … . + | x |n − 1 ) = | P |
x |n − 1
x | − 1
 

Ici on peut supposer que | x | ⩾ 1, sinon le lemme est démontré, donc | x | − 1 est positif et

an | ( | x | − 1 ) ⩽ | P | 
x |n − 1
x |n
⇒ | x | − 1 < 
P |
an |
 

Remarques

Exemple 7   Si P0 = 6 ( X2 − 1 ) et Q0 = 4 ( X3 − 1 ).

Le contenu de P0 est 6, celui de Q0 est 4.
On a donc pgcd des contenus = 2,
P = X2 − 1, Q = X3 − 1. La valeur initiale de z est donc 2 ∗ 1 + 2 = 4. On trouve P ( 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 − 1, donc G = X − 1, sa partie primitive est X − 1. On teste si X − 1 divise P et Q, c’est le cas, donc c’est le pgcd de P et Q et le pgcd de P0 et Q0 est 2 ( X − 1 ).

Algorithme gcdheu
En arguments deux polynômes P0 et Q0 à coefficients entiers ou entiers de Gauss. Retourne le pgcd de P0 et Q0 ou faux en cas d’échec.

  1. Calculer le contenu de P0 et Q0. Vérifier que les coefficients sont entiers de Gauss sinon retourner faux.
  2. Extraire la partie primitive P de P0 et Q de Q0, calculer le pgcd c des contenus de P0 et Q0
  3. Déterminer z = 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 z entre deux itérations (par exemple z ← 2 z).
  5. Calculer le pgcd g de P ( z ) et Q ( z ) puis son écriture symétrique en base z dont on extrait la partie primitive G.
  6. Si G ne divise pasP passer à l’itération suivante. De même pour Q.
  7. Retourner c G
  8. Fin de la boucle
  9. Retourner faux.

On remarque au passage qu’on a calculé le quotient de P par G et le quotient de Q par G 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 D est le pgcd de P et Q dans ℤ (ou ℤ [ i ] ) alors après réduction modulo un nombre premier n qui ne divise pas les coefficients dominants de P et Q, D divise le pgcd G de P et Q dans ℤ / n ℤ (par convention, le pgcd dans ℤ / n ℤ est normalisé pour que son coefficient dominant vaille 1). Comme on calcule G dans ℤ / n ℤ, 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 D à partir de G.

On remarque d’abord que si on trouve G = 1, alors P et Q sont premiers entre eux. En général, on peut seulement dire que le degré de G est supérieur ou égal au degré de D. En fait, le degré de G est égal au degré de D 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 n. Donc plus n est grand, plus la probabilité est grande de trouver G du bon degré.

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

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

Soit P = ∑pi Xi un polynôme à coefficients entiers. On utilise la norme euclidienne

P |2 = pi |2     (5)

On établit d’abord une majoration du produit des racines de norme supérieure à 1 de P à l’aide de | P |. Ensuite si D est un diviseur de P, le coefficient dominant d de D divise le coefficient dominant p de P et les racines de D sont aussi des racines de P. On pourra donc déterminer une majoration des polynômes symétriques des racines de D et donc des coefficients de D.

Lemme 8   Soit A = ∑j = 0a aj Xj un polynôme et α ∈ ℂ. Alors
| ( X − α ) A | = | ( α X − 1 ) A

Pour prouver le lemme 8, on développe les produits de polynômes. On pose a−1 = aa + 1 = 0 et on note ℜ la partie réelle.

| ( X − α ) A |2 = ∑j = 0a + 1 | aj − 1 − α aj |2 = 
a + 1
j = 0
 | aj − 1 |2 + | α |2 | aj |2 − 2 ℜ ( aj − 1 
α  aj
 ) 
| ( α X − 1 ) A |2 = 
a + 1
j = 0
 | 
α
 aj − 1 − aj |2 = 
a + 1
j = 0
 | α |2 | aj − 1 |2 + | aj |2 − 2 ℜ ( 
α
  aj − 1   
aj
 ) 

Les deux donnent bien le même résultat.

Soit P ( X ) = p ∏( X − αj ) la factorisation de P sur ℂ. On introduit le polynôme

P = p 
 
j / | αj | ⩾ 1
 ( X − αj )
 
j / | αj | < 1
 (  
αj
 X − 1 ) 

qui d’après le lemme a la même norme que P. La norme de P majore donc le coefficient constant de P d’où:

  ∏j / | αj | ⩾ 1 | αj | ⩽ | P |/| p |     (6)

On remarque que (6) reste vraie si on considère les racines δj de norme plus grande que 1 d’un diviseur D de P puisque le produit porte alors sur un sous-ensemble. On écrit maintenant l’expression des coefficients dj de D à l’aide des racines δj de D:

dm − j | = | d | 


 
choix de j racines parmi les m racines de D
     
 
δk ∈ racines choisies
 δk 


Pour majorer | dmj |, on commence par majorer | δk | par βk = max( 1, | δk | ). On est donc ramené à majorer

σjm ( β ) = 
 
choix de j parmi m valeurs βk
    
 
βk ∈ choix
 βk  

avec pour hypothèse une majoration de M = ∏k = 1m βk donnée par la relation (6). Pour cela, on cherche le maximum de σj, m ( β ) sous les contraintes M fixé et βk ⩾ 1.

On va montrer que le maximum ne peut être atteint que si l’un des βk = M (et tous les autres βk = 1 ). Sinon, quitte à réordonner supposons que les βk sont classés par ordre croissant. On a donc βm − 1 ≠ 1, on pose βk = βk pour km − 2, βm − 1 = 1 et βm = βm − 1 βm. Comparons σj, m ( β ) et σj, nm ( β ). Si le choix de j parmi m comporte k = m − 1 et k = m, le produit est inchangé. Sinon on a la somme de deux produits, l’un contenant k = m − 1 et l’autre k = m. On compare donc B ( βm − 1 + βm ) et B ( 1 + βm − 1 βm ) avec B = ∏βk ∈ reste du choix βk. Comme

1 + βm − 1 βm ⩾ βm − 1 + βm 

puisque la différence est le produit (1−βm)(1−βm−1) de deux nombres positifs, on arrive à la contradiction souhaitée.

Ensuite on décompose les choix de σm, j en ceux contenant M et des 1 et ceux ne contenant que des 1, d’où la majoration

σjm ( β ) ⩽ 

     m − 1
     j − 1


M + 

     m − 1
     j


et finalement

dm − j | ⩽ | d | 




    m − 1
    j − 1


P |
p |
 + 

    m − 1
    j





    (7)

On peut en déduire une majoration indépendante de j sur les coefficients de D, en majorant | d | par | p | (puisque d divise p) et les coefficients binomiaux par 2m − 1 (obtenue en développant ( 1 + 1 )m − 1). D’où le

Théorème 9   (Landau-Mignotte) Soit P un polynôme à coefficients entiers (ou entiers de Gauss) et D un diviseur de P de degré m. Si | P | désigne la norme euclidienne du vecteur des coefficients de P et p le coefficient dominant de P alors les coefficients dj de D satisfont l’inégalité
dj | ⩽ 2m − 1 ( | P | + | p | )     (8)

Avec cette estimation, on en déduit que si n est un premier plus grand que

min( 2degre ( P ) − 1 ( | P | + | p | ), 2degre ( Q ) − 1 ( | Q | + | q | ) ),      (9)

alors le pgcd trouvé dans ℤ / n ℤ va se reconstruire en un pgcd dans ℤ 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 ℤ s’avèrera encore inefficace (surtout si le pgcd est 1). Cela reste vrai même si on optimise un peu la majoration (9) en repartant de (7).

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 ℤ 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 216 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 P et Q à coefficients entiers. Le résultat renvoyé sera le polynôme pgcd.

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

Boucle infinie :

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

Finalement on n’a pas utilisé b, la borne de Landau-Mignotte. On peut penser que l’étape 6 ne devrait être effectuée que lorsque N est plus grand que 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 b. Mais pour éviter de faire le test trop tôt, on introduit une variable auxiliaire H′ qui contient la valeur de H de l’itération précédente et on ne fait le test que si H′ = 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 et ( X4 − 1 ). Prenons pour commencer n = 2. On trouve comme pgcd X4 + 1 (en effet − 1 = 1 donc on cherchait le pgcd de ( X + 1 )7 et de X4 + 1 = ( X + 1 )4). On teste si X4 + 1 divise P et Q, ce n’est pas le cas donc on passe au nombre premier suivant. Pour n = 3, on trouve X2 − 1. Donc n = 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 X2 − 1 divise P et Q, c’est le cas ici donc on peut arrêter, le pgcd cherché est X2−1.

Exemple 2 :

Calcul du pgcd de ( X + 1 )3 ( X − 1 )4 et ( X4 − 1 )3. Pour n = 2, on trouve un polynôme de degré 7. Pour n = 3, on trouve X6 − 1 donc n = 2 était une mauvaise réduction. Comme X6 − 1 ne divise pas P et Q, on passe à n = 5. On trouve X6 + 2 X4 − 2 X2 − 1. On applique le théorème des restes chinois qui va nous donner un polynôme dans ℤ / 15 ℤ. 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 X6 − 3 X4 + 3 X2 − 1 = ( X2 − 1 )3. Ce polynôme divise P et Q, c’est donc le pgcd de P et de Q.

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 P et Q de k variables X1, … ., Xk en un Xk = z et on calcule le pgcd g des 2 polynômes P ( z ) et Q ( z ) de k − 1 variables. On remonte ensuite à un polynôme G par écriture symétrique en base z de g et on teste si pp ( G ) divise P et Q. Il s’agit à nouveau de montrer que si z est assez grand, alors pp ( G ) est le pgcd cherché. On sait que d = D ( z ) divise g. Il existe donc un polynôme a de k − 1 variables tel que g = a d. On sait aussi que pp ( G ) divise D, donc il existe un polynôme C de k variables tel que D = C ∗ pp ( G ) . On évalue en z et on obtient d = C ( z ) g / c ( G ), où c ( G ) est un entier, donc

c ( G ) = a ∗ C ( z ) 

Comme c ( G ) est un entier, a et C ( z ) sont des polynômes constants. Comme précédemment, on a aussi | C ( z ) | ⩽ | z | / 2 puisque | c ( G ) | ⩽ | 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 Xn − α, où X1, … ., Xn désignent les variables des polynômes. On considère donc deux polynômes P et Q comme polynômes de la variables Xn avec des coefficients dans ℤ [ X1, … ., Xn − 1 ]. On évalue en Xn = α, on obtient deux polynômes en n − 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 Xn, c’est le minimum δ des degrés par rapport à Xn des polynômes P et Q. A première vue, il suffit donc d’évaluer les polynômes en δ + 1 points α.

Il faut toutefois prendre garde aux mauvaises évaluations et à la normalisation des pgcd avant d’interpoler. En effet, si D ( X1, … ., Xn ) désigne le pgcd de P et Q et G ( X1, … ., Xn − 1 ) le pgcd de P ( X1, … ., Xn − 1, α ) et de Q ( X1, … ., Xn − 1, α ), on peut seulement dire D ( X1, … ., Xn − 1, α ) divise G. Plusieurs cas sont donc possibles lorsqu’on évalue en un nouveau point α:

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 le pgcd D, on va reconstruire un multiple du pgcd D (ce multiple appartiendra à ℤ [ Xn ] ). On voit maintenant P et Q comme des polynômes en n − 1 variables X1, … ., Xn − 1 à coefficients dans ℤ [ Xn ]. Alors lcoeff(D), le coefficient dominant de D (relativement à l’ordre lexicographique sur les variables X1,...,Xn−1), est un polynôme en Xn qui divise le coefficient dominant de P et de Q donc divise le coefficient dominant du pgcd des coefficients dominants de P et de Q. On va donc reconstruire le polynôme :

D′ = D 
Δ ( Xn )
lcoeff ( D ) ( Xn )
, Δ ( Xn ) = pgcd ( lcoeff ( P ) ( Xn ), lcoeff ( Q ) ( Xn )) 

c’est-à-dire D multiplié par un polynôme qui ne dépend que de Xn.

Revenons à G en un point α de bonne évaluation. C’est un multiple entier de D ( X1, … ., Xn − 1, α ):

G = β D ( X1, … ., Xn − 1, α ) 

Donc, comme polynômes de X1,...,Xn−1 à coefficients dans ℤ[Xn] ou dans ℤ, lcoeff ( G ) = β lcoeff ( D )| Xn = α. Comme lcoeff ( D ) divise Δ ( Xn ), il en est de même en Xn = α donc lcoeff(G) divise β Δ(α). On en déduit que Δ ( α) G qui est divisible par Δ (α) β est divisible par lcoeff ( G ). On va donc considérer le polynôme Δ (α) G / lcoeff ( G ) : ses coefficients sont entiers et son coefficient dominant est

Δ ( α) = lcoeff(D′( X1, … ., Xn − 1, α ))

donc

Δ (α) G  / lcoeff ( G )= D′( X1, … ., Xn − 1, α )

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

Arguments: 2 polynômes primitifs P et Q de n variables X1, … ., Xn à coefficients entiers. Renvoie le pgcd de P et Q.

  1. Si n = 1, renvoyer le pgcd de P et Q en une variable.
  2. Test rapide de pgcd trivial par rapport à Xn. On cherche des n − 1-uplets α tels que P ( α, Xn ) et Q ( α, Xn ) soient de même degré que P et Q par rapport à la variable Xn. On calcule le pgcd G de ces 2 polynômes en une variable. Si le pgcd est constant, alors on retourne le pgcd des coefficients de P et Q.
  3. On divise P et Q par leur contenu respectifs vu comme polynômes en X1, … ., Xn − 1 à coefficients dans ℤ [ Xn ], on note C ( Xn ) le pgcd des contenus. On calcule aussi le pgcd Δ ( Xn ) des coefficients dominants de P et de Q.
  4. On initialise D′ le pgcd reconstruit à 0, I ( Xn ) le polynôme d’interpolation à 1, δ=(δ1,...,δn−1) la liste des degrés partiels du pgcd par rapport à X1, … ., Xn − 1 au minimum des degrés partiels de P et Q par rapport à X1, … ., Xn − 1, e le nombre d’évaluation à 0 et E l’ensemble des points d’interpolation à la liste vide.
  5. Boucle infinie:
    • Faire α=entier aléatoire n’appartenant pas à E jusqu’à ce que
            degre(P ( X1, … ., Xn − 1, α ))=degreXn ( P ( X1, … ., Xn )  
            degre ( Q ( X1, … ., Xn − 1, α )) = degreXn ( Q ( X1, … ., Xn ))     
    • Calculer le pgcd G ( X1, … ., Xn − 1 ) en n − 1 variables de P ( X1, … ., Xn − 1, α ) et Q ( X1, … ., Xn − 1, α ).
    • Si degre ( G )i < δi pour un indice au moins. Si degre ( G ) ⩽ δ, on pose δ = degre ( G ), D′ = G Δ ( α )/lcoeff ( G ), I = Xn − α, e = 1 et E = [ α ], sinon on pose δ = min( δ, degre ( G )), D′ = 0, I = 1, e = 0, E = [ ]. On passe à l’itération suivante.
    • Si degre ( G ) > δ, on passe à l’itération suivante.
    • Si degre ( G ) = δ, on interpole:
      • G := G Δ ( α )/lcoeff ( G )
      • D′ := D′ + I ( Xn )/∏αjE ( α − αj ) ( GD′ ( X1, … ., Xn − 1, α ))
      • I := I ∗ ( Xn − α )
      • e := e + 1 et ajouter α à E
      • Si e est strictement plus grand que le minimum des degrés partiels de P et Q par rapport à Xn, on pose D la partie primitive de D′ (vu comme polynôme à coefficients dans ℤ [ Xn ]), on teste si P et Q sont divisibles par D, si c’est le cas, on renvoie D = C ( Xn ) D

On observe que dans cet algorithme, on fait le test de divisibilite de D par P et Q. 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′ n’est pas modifié par l’ajout d’un nouveau point à la liste des αj.

Il existe une variation de cet algorithme, appelé SPMOD (sparse modular), qui suppose que seuls les coefficients non nuls du pgcd en n − 1 variables sont encore non nuls en n 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 − 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 P et Q sont dans ℤ / n ℤ 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 F et G sous forme factorisée, le but étant de faire comprendre l’algorithme. En utilisation normale, on n’exécuterait cet algorithme que si F et G étaient développés.

P = (( x + 1 ) y + x2 + 1 ) ( y2 + x y + 1 ), Q = (( x + 1 ) y + x2 + 1 ) ( y2x y − 1 ).

Prenons x comme variable X1 et y comme variable X2. Les coefficients dominants de P et Q sont respectivement y et − y donc Δ = y.

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

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

En y = 2, P ( x, 2 ) = ( x2 + 2 x + 3 ) ( 2 x + 5 ) et Q ( x, 2 ) = ( x2 + 2 x + 3 ) ( − 2 x + 3 ) donc G = x2 + 2 x + 3, Δ ( 2 ) = 2. On interpole:

D′ = x2 + x + 2 + 
y − 1
2 − 1
 ( 2 ( x2 + 2 x + 3 ) − ( x2 + x + 2 )) = y ( x2 + 3 x + 4 ) − ( 2 x + 2 ) 

On teste la divisibilité de P par D′, le test échoue.

En y = 3, P ( x, 3 ) = ( x2 + 3 x + 4 ) ( 3 x + 10 ) et Q ( x, 3 ) = ( x2 + 3 x + 4 ) ( − 3 x + 8 ) donc G = x2 + 3 x + 4, Δ ( 3 ) = 3. On interpole:

  D=y ( x2 + 3 x + 4 ) − ( 2 x + 2 ) + 
  
y − 2 ) ( y − 1 )
( 3 − 2 ) ( 3 − 1 )
 
3 ( x2 + 3 x + 4 ) − ( 3 ( x2 + 3 x + 4 ) − ( 2 x + 2 )) 

donc

D′ = y ( x2 + 3 x + 4 ) − ( 2 x + 2 ) + 
y − 2 ) ( y − 1 )
2
 ( − 2 x − 2 ) = x2 y + x y2 + y2 + y 

On divise D′ par son contenu et on trouve x2 + x y + y + 1 qui est bien le pgcd de P et Q.

5.3.3  EZGCD.

Il s’agit d’une méthode p-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 F et G deux polynômes primitifs dépendant des variables X1, …, Xn de pgcd D, on fixe une des variables qu’on appelera X1 dans la suite. Soient lcoeff ( F ) et lcoeff ( G ) les coefficients dominants de F et G par rapport à X1. On évalue F et G en un n − 1 uplet b tel que le degré de F et G par rapport à X1 soit conservé après evaluation en b. On suppose que Db ( X1 ) = pgcd ( F ( b ), G ( b )) a le même degré que D ( b ). On a donc l’égalité:

F ∗ lcoeff ( F )) ( b ) = 


Db  
lcoeff ( F ( b ))
lcoeff ( Db )
 


∗ 





F ( b )
Db
 
lcoeff ( F ) ( b )
lcoeff ( 
F ( b )
Db
 )






et de même en remplaçant F par G.

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

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

On adopte ensuite la notation suivante: si k est un entier, on dit qu’un polynôme P est un O ( k ) si la valuation de P vu comme polynôme en X2, … ., Xn à coefficients dans ℤ [ X1 ] est supérieure ou égale à k, ou de manière équivalente si

P ( X1h X2, … ., h Xn ) = Oh → 0 ( hk ) 

L’égalité à lifter se réécrit donc:

F lcoeff ( F ) = P0 Q0 + O ( 1 ) 

P0 =Db lcoeff ( F ( b ))/lcoeff ( Db ) et Q0 = F ( b )/Db lcoeff ( F ) ( b )/lcoeff ( F ( b )/Db ) sont premiers entre eux et de degré 0 par rapport aux variables X2, … ., Xn. Cherchons P1 = O ( 1 ) et Q1 = O ( 1 ) de degré 1 par rapport aux variables X2, … ., Xn tels que

F lcoeff ( F ) = ( P0 + P1 ) ( Q0 + Q1 ) + O ( 2 ) 

Il faut donc résoudre

F lcoeff ( F ) − P0 Q0 = P0 Q1 + Q0 P1 + O ( 2 ) 

On peut alors appliquer l’identité de Bézout qui permet de déterminer des polynômes P1 et Q1 satisfaisant l’égalité ci-dessus (avec comme reste O ( 2 ) nul) puisque P0 et Q0 sont premiers entre eux. De plus, on choisit P1 et Q1 tels que degreX1 P1 ⩽ degreX1 ( F ) − degre ( Q0 ) = degre ( P0 ) et degreX1 ( Q1 ) ⩽ degre ( Q0 ) et lcoeffX1 ( P0 + P1 ) + O ( 2 ) = lcoeffX1 ( Q0 + Q1 ) + O ( 2 ) = lcoeffX1 ( F ). On tronque ensuite P1 et Q1 en ne conservant que les termes de degré 1 par rapport à X2, … ., Xn.

On trouve de la même manière par récurrence Pk et Qk homogènes de degré k par rapport à X2, … ., Xk, de degré par rapport à X1 respectivement inférieur aux degrés de Q0 et de P0 et tels que

F lcoeff ( F ) = ( P0 + … . + Pk ) ( Q0 + … . + Qk ) + O ( k + 1  )      (10)

et lcoeff ( F ) = lcoeffX1 ( P0 + … . + Pk ) + O ( k + 1 ) = lcoeffX1 ( Q0 + … . + Qk ) + O ( k + 1 ).

Si on est bien en un point de bonne évaluation et si k est plus grand que le degré total (par rapport aux variables X2, … ., Xn) du polynôme F lcoeff ( F ) on va vérifier que P0 + … . + Pk = D lcoeff ( F )/lcoeff ( D ). En effet, si on a deux suites de polynômes P et P′ et Q et Q′ satisfaisant (10) avec les même termes de degré zéro P0 et Q0, alors en prenant la différence, on obtient:

P0 + P1 …  + Pk ) ( Q0 + Q1 …  + Qk ) = ( P0 + P1′ …  + Pk′ ) ( Q0 + Q1′ …  + Qk′ ) + O ( k + 1 ) 

On égale alors les termes homogènes de degré j, pour j = 1, on obtient P0 ( Q1Q1′ ) = Q0 ( P1P1′ ), donc Q0 divise Q1Q1′ qui est de degré strictement inférieur au degré de Q0 par rapport à X1 (car on a l’inégalité large et les termes de plus haut degré sont égaux), donc Q1 = Q1′ et P1 = P1′. On montre de la même manière que Qj = Qj′ et Pj = Pj′. L’écriture est donc unique, c’est donc l’écriture en polynôme homogène de degré croissant de D lcoeff ( F )/lcoeff ( D ) que l’on reconstruit.

Cet algorithme permet donc de reconstruire D, il suffit de tester à chaque étape si P0 + … . + Pk divise F 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 ( 2 k ). Elle nécessite toutefois un calcul supplémentaire, celui de l’identité de Bézout à O ( 2 k ) près pour les polynômes P0 + … . + Pk − 1 et Q0 + … . + Qk − 1. Ce calcul se fait également par lifting.

Algorithme EZGCD (Hensel linéaire)

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

  1. Evaluer F et G en ( X2, … ., Xn ) = ( 0, … ., 0 ), vérifier que les coefficients dominants de F et de G ne s’annulent pas. Calculer le pgcd Db de F ( 0 ) et de G ( 0 ). Prendre un autre point d’évaluation au hasard qui n’annule pas les coefficients dominants de F et de G et vérifier que le pgcd a le même degré que Db. Sinon, renvoyer false (on peut aussi faire une translation d’origine de F et de G en un autre point mais cela diminue l’efficacité de l’algorithme).
  2. On note lcF et lcG les coefficients dominants de F et de G par rapport à X1.
  3. Si degre ( F ) ⩽ degre ( G ) et degre ( Db ) = degre ( G ) et F divise G renvoyer F
  4. Si degre ( G ) < degre ( F ) et degre ( Db ) = degre ( F ) et G divise F renvoyer G
  5. Si degre ( F ) = degre ( Db ) ou si degre ( G ) = degre ( Db ) renvoyer false
  6. Boucle infinie sur j entier initialisé à 0, incrémenté de 1 à chaque itération: si pgcd ( Db, F ( 0 ) + j G ( 0 )/Db ) = C constant, alors arrêter la boucle
  7. Lifter l’égalité ( F + j G ) ( lcF + j lcG ) ( 0 ) = ( Db ( lcF + j lcG ) ( 0 )/lcoeff ( Db ) ) ∗ … . 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 X1 à coefficients dans ℤ [ X2, … ., Xn ].

Remontée de Hensel linéaire:

Arguments: F un polynôme, lcF=lcoeff(F) son coefficient dominant, P0 un facteur de F ( 0 ) ayant comme coefficient dominant lcF ( 0 ) et dont le cofacteur Q0 est premier avec P0.

Renvoie deux polynômes P et Q tels que F lcF = P Q et P ( 0 ) = P0 et lcoeff ( P ) = lcoeff ( Q ) = lcF.

  1. Soit G = F lcF, , Q0 = G ( 0 ) / P0, P = P0, Q = Q0.
  2. Déterminer les deux polynômes U et V de l’identité de Bézout (tels que P0 U + Q0 V = dd est un entier).
  3. Boucle infinie avec un compteur k initialisé à 1, incrémenté de 1 à chaque itération
    • Si k > degreX2, … ., Xn ( G ), renvoyer false.
    • Si P divise G, renvoyer P et G / P.
    • Soit H = GP Q = O ( k ). Soit u = U H/d et v = V H/d, on a P0 u + Q0 v = H
    • Remplacer v par le reste de la division euclidienne de v par P0 et u par le reste de la division euclidienne de u par Q0. La somme des deux quotients est égale au quotient euclidien de H par P0 Q0, c’est-à-dire au coefficient dominant de H divisé par le produit des coefficients dominants de P0 et Q0 (qui sont égaux) donc on a l’égalité:
      P0 u + Q0 v = H − 
      lcoeff ( H )
      lcoeff ( P0 )2
       P0 Q0 
    • Soit α = ( lcoeff ( F ) − lcoeff ( P )) / lcoeff ( P0 ) et β = ( lcoeff ( F ) − lcoeff ( Q )) / lcoeff ( P0 ). On ajoute α P0 à v, ainsi lcoeff ( P + v ) = lcoeff ( F ) + O ( k + 1 ) et β Q0 à u, ainsi lcoeff ( Q + u ) = lcoeff ( F ) + O ( k + 1 )

      Remarque: on montre alors que α + β = lcoeff ( H )/lcoeff ( P0 Q0 ) + O ( k + 1 ) donc P0 u + Q0 v = H + O ( k + 1 ) en utilisant les propriétés :

      lcoeff ( F ) = lcoeff ( P ) + O ( k ) = lcoeff ( Q ) + O ( k ) = lcoeff ( P0 ) + O ( 1 ) 
    • Réduire u et v en éliminant les termes de degré strictement supérieur à k par rapport à X2, … ., Xn. S’il reste un coefficient non entier, renvoyer false
    • Remplacer P par P + v et Q par Q + u, passer à l’itération suivante.

Exemple:

F = (( x + 1 ) y + x2 + 1 ) ( y2 + x y + 1 ), G = (( x + 1 ) y + x2 + 1 ) ( y2x y − 1 )

On a F ( 0, y ) = ( y + 1 ) ( y2 + 1 ) et G ( 0, y ) = ( y + 1 ) ( y2 − 1 ), le pgcd est donc Db = ( y + 1 ). On remarque que Db est premier avec le cofacteur de F mais pas avec le cofacteur de G. Si on évalue en un autre point, par exemple x = 1, on trouve un pgcd D1 de même degré, donc 0 est vraissemblablement un bon point d’évaluation (ici on en est sûr puisque le pgcd de F et G se calcule à vue...). On a lcoeff ( F ) = x + 1, on va donc lifter G = (( x + 1 ) y + x2 + 1 ) ( y2 + x y + 1 ) ( x + 1 ) = P QP0 = ( y + 1 ) et Q0 = ( y2 + 1 ).

On calcule les polynômes de l’identité de Bézout U = ( 1 − y ) et V = 1 avec d = 2, puis à l’ordre k = 1:

H = G − P0 Q0 = ( 2 y3 + 2 y2 + 3 y + 1 ) x + O ( 2 ) 

donc u = reste ( U H / d, Q0 ) = x y et v = reste ( V H / d, P0 ) = − x.

Donc Q1 = x y + α Q0 avec α = ( x + 1 − 1 ) / lcoeff ( P0 ) = x et Q0 + Q1 = ( y2 + 1 ) ( x + 1 ) + x y. De même, P1 = − x + β P0, avec β = ( x + 1 − 1 ) / lcoeff ( P0 ) = x donc P0 + P1 = ( y + 1 ) ( x + 1 ) − x. On remarque que P0 + P1 et Q0 + Q1 sont bien à O ( 2 ) près les facteurs de F lcoeff ( F ):

P = ( x + 1 ) y + x2 + 1 = P0 + P1 + O ( 2 ),  Q = ( x + 1 ) ( y2 + x y + 1 ) = Q0 + Q1 + O ( 2 ) 

Une deuxième itération est nécessaire. On calcule

H = G − ( P0 + P1 ) ( Q0 + Q1 ) = ( 2 y2 + y + 1 ) x2 + O ( 3 ) 

puis reste ( U H / d, Q0 ) = y x2 et reste ( V H / d, P0 ) = x2. Ici les coefficients α et β sont nuls car lcoeff ( F ) n’a pas de partie homogène de degré 2. On trouve alors P = P0 + P1 + P2 et Q = Q0 + Q1 + Q2. Pour calculer le pgcd, il suffit de calculer la partie primitive de P vu comme polynôme en y, ici c’est encore P car le contenu de P est 1 (remarque: pour Q le contenu est x + 1).
On trouve donc P comme pgcd.

5.4  Quel algorithme choisir?

Il est toujours judicieux de faire une évaluation en quelques n − 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 :

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 A et B à coefficients dans un corps, de degrés p et q et de pgcd D et l’identité de Bézout correspondante :

  A U + B V =D     (11)

avec degré(U)<q et degré(V)<p. Imaginons qu’on cherche U et V 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+q équations (obtenues en développant et en identifiant chaque puissance de X de 0 à p+q−1) à p+q inconnues (les p coefficients de V et les q coefficients de U) On sait que A et B sont premiers entre eux si et seulement si ce problème d’algèbre linéaire a une solution pour D=1. Donc si le déterminant du système est non nul, alors A et B sont premiers entre eux. Réciproquement si A et B sont premiers entre eux, le système a une solution unique non seulement avec comme second membre 1 mais avec n’importe quel polynôme de degré inférieur p+q, donc le déterminant du système est non nul.

Définition:
On appelle résultant de A et B le déterminant de ce système (11). Il s’annule si et seulement si A et B 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 U et V)

M(A,B)=







AaAa−1A00
0AaA1A0
      ⋮ 
00    A0 
BbBb−1B000
      ⋮ 
00    B0 








(cette matrice contient b=degré(B) lignes de coefficients du polynôme A et a=degré(A) lignes de coefficients du polynôme B, les coefficients diagonaux sont les Aa et B0)

Remarques
Le résultant s’exprime polynomialement en fonction des coefficients des polynômes A et B donc aussi en fonction des coefficients dominants de A et B et des racines α1,..., αa de A et β1,...,βb B, or si on fait varier les racines de B on annulera le résultant si l’une d’elle coincide avec une racine de A, ceci montre que le résultant est divisible par le produit des différences des racines βj−αi de A et B. On montre que le quotient est Aab Bba en regardant le coefficient dominant du résultant en degré total par rapport aux βj : dans le déterminant il faut prendre le produit des termes diagonaux pour avoir le degré maximal en les βj. On peut aussi l’écrire sous la forme

resultant(A,B)=Aab 
a
i=1
 Bi)

Soit P un polynôme de degré n et coefficient dominant pn. Le coefficient dominant de P′ est npn, un multiple de pn, le résultant de P et P′ est donc divisible par pn, on appelle le quotient discriminant. En terme des racines ri de P, on a

disc(P)=
resultant(P,P′)
pn
=pnn−2 
n
i=1
 P′(ri) = pn2n−2 
 
1≤ i<j≤ n
 (rirj)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 exemple ∈ ou un anneau de polynômes [X2,...,Xn] dans son corps de fractions. On remarque alors que l’équation :

AU+BV=C

a une solution dans l’anneau si C est le résultant r de A et B (ou un multiple). En effet, on écrit les solutions comme celles d’un système de Cramer, le dénominateur de chaque inconnue est r, le numérateur est un déterminant ayant les coefficients de C dans une des colonnes, on peut donc y factoriser r et simplifier. On peut le voir directement à partir de la définition du résultant en effectuant sur le déterminant une manipulation de colonnes sur la dernière colonne, on ajoute à cette dernière colonne x fois l’avant-dernière, x2 fois l’avant-avant-dernière etc... La dernière colonne devient








xb−1 A 
... 
A 
xa−1 B 
...
B
 






et en développant le déterminant par rapport à cette dernière colonne, on obtient l’identité de Bézout.

Exemple : le résultant de x+1 et x−1 est 2, donc l’équation

(x+1)U+(x−1)V=2

a une solution U=1 et V=−1 dans [X], par contre

(x+1)U+(x−1)V=1

n’a pas de solution dans [X].

6.2.1  Systèmes polynomiaux

Ceci peut servir à éliminer des inconnues lorsqu’on résoud un système d’équations polynomiales :

P1(X1,...,Xn)=0, ..., Pn(X1,...,Xn)=0

On pose

P11(X1,...,Xn−1)=resultant(P1,Pn,Xn), ...,  Pn−11(X1,...,Xn−1)=resultant(Pn−1,Pn,Xn)

Comme P11, Pn−11, ... sont des combinaisons linéaires des polynômes de départ à coefficients dans l’anneau, si (X1,...,Xn) est solution du système de départ, alors X1,...,Xn−1 est solution du deuxième système. On élimine ainsi les variables les unes après les autres, pour se ramener à une seule équation polynomiale P1n−1(X1)=0, dont on cherche les racines, puis si r1 est une racine de P1n−1, on remonte au système précédent P1n−2(r1,X2)=0, P2n−2(r1,X2)=0, que l’on résoud en cherchant les racines de gcd(P1n−2(r1,X2),P2n−2(r1,X2)), et ainsi de suite jusqu’au système de départ.

Lors des calculs de résultant, il peut arriver que le résultat soit nul si les arguments ne sont pas premiers entre eux, dans ce cas il faut diviser par le PGCD de ces 2 polynômes et traiter le cas du PGCD à part.

Malheureusement, les calculs de résultant deviennent vite impraticables (cf. infra), on ne peut guère traiter par cette méthode que des systèmes 3x3 (ou 4x4 si on est patient). Pour des systèmes plus ambitieux, on utilisera plutôt un calcul de bases de Groebner. Mais le résultant est très bien adapté par exemple à la recherche d’équations cartésiennes d’une courbe ou surface paramétrée par des fractions rationnelles.

6.2.2  Identité de Bézout dans [X].

Lorsque les polynômes A et B sont à coefficients entiers, on peut résoudre l’identité de Bézout en résolvant le système linéaire AU+BV=DD est calculé par un algorithme de calcul de PGCD. C’est un système linéaire de a+b équations en a+b inconnues, que l’on peut résoudre par une méthode p-adique, permettant de calculer efficacement les coefficients rationnels de U et V.

On peut étendre à des polynômes premiers entre eux ayant des coefficients dans une extension algébrique de , en calculant l’identité de Bézout pour Ã=norme(A) et B=norme(B) (norme au sens de norme algébrique, obtenue en multipliant par les conjugués du polynôme, calculé en pratique en prenant le résultant du polynôme avec le polynôme minimal de l’extension algébrique) et en observant que la norme est divisible par le polynôme

àŨ+BṼ=1 ⇒ A 
Ã
A
 Ũ+  B 
B
B
Ṽ=1

On pose alors u le reste de Ã/A Ũ par B et v le reste de B/BṼ par A. Si les polynômes ne sont pas premiers entre eux, on s’y ramène en divisant par leur pgcd.

6.3  Résultant et degrés

Si A et B sont des polynômes en d variables de degré total m et n alors le résultant de A et B par rapport à une des variables, disons la première notée x, est un polynôme en d−1 variables, on va voir qu’on peut majorer son degré total par mn.

Quitte à ajouter une variable d’homogénéisation (appelons-la t), on peut supposer que A et B sont homogènes, par exemple si A=x3+xy+1 on considère At=x3+xyt+t3. Le degré total par rapport aux d−1 variables d’un coefficient Aj de A est alors mj, et pour Bk c’est nk. On développe le déterminant comme somme sur toutes les permutations de a+b éléments, et on regarde le degré total d’un terme par rapport aux d−1 variables, on a donc un produit de ri σ(i) pour i entre 1 et a+b. Pour i entre 1 et b, on est dans les b premières lignes, donc avec des coefficients de A, le degré total de ri σ(i) se déduit de la distance à la diagonale, il vaut ma+i−σ(i) puisque sur la diagonale c’est ma. Pour i entre b+1 et b+a on est dans les a dernières lignes, donc avec des coefficients de B, le degré total est n+σ(i)−i. Le degré total du produit vaut donc

b(ma)+an+
a+b
i=1
 σ(i) −i = b(ma)+an=mn−(ma)(nb)

il est donc au plus mn avec égalité si m=a ou n=b (c’est-à-dire si le degré total est identique au degré partiel en x pour au moins un des deux polynômes).

Lorsqu’on enlève la variable d’homogénéisation (en posant t=1), on peut également perdre un ou plusieurs degrés. Dans le cas de polynômes en 2 variables A(x,y), B(x,y), cela correspond à un point d’intersection à l’infini entre les 2 courbes A(x,y)=B(x,y)=0, en coordonnées homogènes on a t=0 qui est solution, et on remplace t par 0 dans At(x,y,t)=Bt(x,y,t)=0 pour trouver la direction.

Exemple (tiré d’un TP de Frédéric Han) intersection des 2 courbes x*y=4 et y2=(x−3)*(x2−16). On a donc A=xy−4, B=y2−(x−3)(x2−16), m=2, n=3 on définit alors :
A:=x*y-4t^2; B:=y^2*t-(x-3t)*(x^2-16t^2);
On observe que resultant(A,B,x) est bien de degré 6 (car n=b=3), alors que resultant(A,B,y) est de degré 5 (ma, nb). On a donc 5 points d’intersection complexes et un point d’intersection à l’infini correspondant à la racine t=0 du résultant en x de coordonnées homogènes (x,y,t)=(0,1,0). Illustration
solve(subst(resultant(A,B,y),t=1))
implicitplot(subst(A,t=1));implicitplot(subst(B,t=1))

Plus généralement, soit deux courbes algébriques d’équations respectives A(x,y)=0 et B(x,y)=0 de degré totaux m et n et premiers entre eux, alors A et B ont au plus mn points d’intersection (théorème de Bézout). En effet, le résultant en x par exemple est non nul puisque les 2 polynômes sont premiers entre eux, donc est un polynôme en y qui a un nombre fini de racines, puis on cherche les racines en x de gcd(A(.,y),B(.,y)) pour chaque valeur de y racine, il y a donc un nombre fini d’intersections. On peut donc changer de repère et choisir un repère tel que deux points d’intersections distincts aient leurs abscisses distinctes. On refait le même raisonnement, et on utilise la majoration du degré du résultant par rapport à y par mn, on a donc au plus mn valeurs de y, donc au plus mn points d’intersections, puisqu’une valeur de y ne correspond qu’à une valeur de x par choix du repère. Lorsqu’on travaille dans 2, le défaut de nombre de points d’intersection par rapport au majorant mn provient des points à l’infini, à condition de prendre en compte la multiplicité des intersections. Dans 2, on perd aussi les points non réels. Exemple : intersection de (x−2)2+y2=4 et y2=(x−3)*(x2−16).

Le degré du résultant explique pourquoi on ne peut pas résoudre en pratique de grands systèmes polynomiaux avec cet outil d’élimination. Par exemple pour un système de 5 équations en 5 inconnues de degré 5, en éliminant une variable, on passe à 4 équation en 4 inconnues de degré 25, puis à 3 équations en 3 inconnues de degré 252=625, puis 2 équations en 2 inconnues de degré 6252=390625 et enfin un polynôme de degré ... 152587890625. Pour n équations de degré n, on a une majoration par n(2n−1), ainsi pour n=4 on trouve 65536 qui est déjà discutable...

6.4  Lien avec l’algorithme du sous-résultant (calcul de PGCD)

On peut calculer le déterminant avec la suite des restes de divisions euclidiennes de la manière suivante, on part de la pseudo-division de A par B :

Bbab+1 A=BQ+R 

on effectue alors sur chaque ligne contenant les coefficients de A la manipulation de ligne correspondante, c’est-à-dire multiplier la ligne par Bbab+1 et soustraire (q0 fois la ligne de B terminant dans la même colonne+q1 fois la ligne de B terminant une colonne avant+...). Toutes les lignes contenant les coefficients de A ont été remplacées par des lignes contenant les coefficients de R. Ces lignes contiennent k zéros initiaux avec k ≥ 1, ce qui permet de réduire le déterminant à celui de la matrice de Sylvester de R et B (à un coefficient multiplicatif près qui vaut Bbk par rapport au précédent donc Bbkb(ab+1) par rapport au déterminant de départ). On échange ensuite R et B ce qui change éventuellement le signe et on continue en faisant les divisions euclidiennes de l’algorithme du sous-résultant (cf. Knuth où on utilise la matrice de Sylvester pour prouver que l’algorithme du sous-résultant est correct). Rappelons que le sous-résultant définit les suites Ak (A0=A, A1=B), dk le degré de Ak, δk=dkdk+1, gk (g0=1, si k≠ 0, gk coefficient dominant de Ak) hk (h0=1, hk+1=hk1−δk gk+1δk) et

gkδk−1+1 Ak−1 = Ak Qk+1 +  gk−1 hk−1δk−1 Ak+1 
Théorème 10   Le résultant est égal au signe près au coefficient hkk correspond au reste Ak constant (en supposant que le résultant soit non nul).

Preuve
La transcription de l’égalité précédente sur les résultants donne par la méthode ci-dessus :

 gkk−1+1)dkRes(Ak−1,Ak)=gkdk−1dk+1 Res(gk−1 hk−1δk−1 Ak+1,Ak)
 =gkdk−1dk+1 (gk−1 hk−1δk−1)dk Res(Ak+1,Ak)

On en déduit que :

Res(Ak−1,Ak)
gk−1dk hk−1dk−1−1
gkdk−1dk+1−(δk−1+1)dk   hk−1δk−1dk+1−dk−1 Res(Ak+1,Ak

On observe que :

hk−1δk−1dk+1−dk−1 =hk−1k−1−1)(dk−1)
hk−1δk−1−1
dk−1


gkδk−1
hk
 


dk−1



 

donc :

 
Res(Ak−1,Ak)
gk−1dk hk−1dk−1−1
=
gkdk−1dk+1−(δk−1+1)dk  


gkδk−1
hk
 


dk−1



 
  Res(Ak+1,Ak
 =
gkdk−1dk+1dk−δk−1  


1
hk
 


dk−1



 
  Res(Ak+1,Ak
 =
 Res(Ak+1,Ak
 gkdk+1 hkdk−1

Donc en valeur absolue

|
Res(A0,A1)
g0d1 h0d0−1
| = |
Res(Ak−1,Ak)
gk−1dk hk−1dk−1−1
 |

En prenant le rang k tel que Ak est constant, on a dk=0 et le résultant est égal à gkdk−1, on obtient donc :

|Res(A0,A1)|=|
gkdk−1
 hk−1dk−1−1
 |

Comme ici δk−1=dk−1, le terme de droite est |hk|.

Remarque
On peut calculer au fur et à mesure le signe du résultant en tenant compte des degrés de Ak pour inverser l’ordre de Ak−1 et Ak dans le résultant.

Utilisation
La valeur du résultant est très utile pour savoir si 2 polynômes dépendant de paramètres sont premiers entre eux en fonction de la valeur des paramètres. En effet, la fonction gcd d’un logiciel de calcul formel calculera le PGCD par rapport à toutes les variables en incluant les paramètres. En cherchant quand le résultant s’annule en fonction des paramètres on obtient un autre type d’information.

Exemple :
Chercher quand le polynône P=x3+px+q possède une racine multiple en fonction de p et q. On calcule le résultant de P et P′ et on trouve 4p3+27q2, donc P a une racine multiple si et seulement si 4p3+27q2=0.

6.5  Calcul efficace du résultant

On dispose essentiellement de deux stratégies avec des sous-stratégies :

Exemple de comparaison :


7  Localisation des racines

7.1  Majoration

On a vu au lemme 5, que si z est une racine complexe d’un polynôme P de coefficient dominant pn alors

|z| ≤ 1 + 
|P|
|pn|

7.2  Les suites de Sturm.

L’algorithme du sous-résultant appliqué à un polynôme sans racine multiple P et à sa dérivée permet, à condition de changer les signes dans la suite des restes, de connaitre le nombre de racines réelles d’un polynôme à coefficients réels dans un intervalle. Ceci est trè utile pour par exemple simplifier des valeurs absolues de polynômes dans un intervalle.

On définit donc la suite de polynômes A0=P, A1=P′, ..., Ak,0 par :

  Ai = Ai+1 Qi+2 − Ai+2      (12)

avec Ak, le dernier reste non nul, un polynôme constant puisque P n’a pas de racine multiple. On utilise plutot l’algorithme du sous-résultant que l’algorithme d’Euclide, il faut alors s’assurer que les signes de Ai et Ai+2 sont opposés lorsque Ai+1 s’annule quitte à changer le signe de Ai+2 en fonction du signe du coefficient dominant de Ai+1, de la parité de la différence des degrés et du signe du coefficient gh1−δ.

On définit s(a) comme étant le nombre de changements de signes de la suite Ai(a) en ignorant les 0. On a alors le

Théorème 11   Le nombre de racines réelles de A0=P sur l’intervalle ]a,b] est égal à s(a)−s(b).

Preuve
Par continuité de la suite des polynômes, s ne peut varier que si l’un des polynômes s’annule. On considére la suite des signes en un point : elle ne peut contenir deux 0 successifs (sinon toute la suite vaudrait 0 en ce point en appliquant (12), or Ak est constant non nul). Elle ne peut pas non plus contenir +,0,+ ni -,0,- à cause de la convention de signe sur les restes de (12). Donc une racine b de Ai pour 0<i<k, n’influe pas sur la valeur de s au voisinage de b (il y a toujours un changement de signe entre les positions i−1 et i+1). Comme Ak est constant, seules les racines de A0=P sont susceptibles de faire varier s. Comme A1=P′, le sens de variations de A0 au voisinage d’une racine de A0 est déterminé par le signe de A1, donc les possibilités sont -,+ vers +,+ ou +,- vers -,-, ce qui diminue s d’une unité.


7.3  Autres algorithmes

8  Exercices (PGCD, résultant, ...)

8.1  Instructions

Elles sont dans les menus Cmds->Integer et Cmds->Polynomes de Xcas.

8.1.1  Entiers

8.1.2  Polynômes

On peut représenter les polynômes par leur écriture symbolique (par exemple x^2+1), ou par des listes (représentation dense ou creuse, récursive ou distribuée). Xcas propose deux types de représentation, dense à une variable (poly1[ ]), ou distribuée (%%%{ }%%%) et des instructions de conversion (poly2symb et symb2poly) entre représentations. L’intérêt d’une représentation non symbolique est l’efficacité des opérations polynomiales, (et la possibilité de chronométrer des opérations comme le produit de 2 polynômes).

Les instructions qui suivent utilisent la représentation symbolique, certaines acceptent les autres représentations.

Notez aussi que le menu Exemples->poly->pgcd.xws de Xcas contient des exemples de programmes de calcul de pgcd de type Euclide.

8.1.3  Calculs modulo n

Pour travailler dans /n[X] :

8.2  Exercices PGCD

  1. Calculez le pgcd de x202+x101+1 et sa dérivée modulo 3 et modulo 5. Conclusion?
  2. P=51x3−35x2+39x−115 et Q=17x4−23x3+34x2+39x−115. Calculez le pgcd de P et Q modulo 5, 7 et 11. En déduire le pgcd de P et Q par le théorème des restes chinois. Pourquoi ne doit-on pas essayer modulo 17?
  3. Écrire un programme qui détermine le degré probable du pgcd de 2 polynômes en une variable en utilisant le pgcd modulaire (on considère le degré probable déterminé lorsqu’on trouve deux nombres premiers réalisant le minimum des degrés trouvés)
  4. Détaillez l’algorithme du PGCD heuristique pour les polynômes P=(x+1)7−(x−1)6 et sa dérivée. Comparez avec l’algorithme d’Euclide naïf.
  5. Écrire un programme mettant en oeuvre le pgcd heuristique pour des polynômes à une variable.
  6. On veut comprendre comment un logiciel de calcul formel calcule
    x6+2
    (x3+1)2
      dx 
    On se ramène d’abord à une fraction propre (numérateur N de degré inférieur au dénominateur), Soit P=x3+1, calculez le PGCD de P et P′, puis deux polynômes U et V tels que N=UP+VP′ . On décompose alors l’intégrale en deux morceaux :
    N
    P2
    =
    U
    P
      + V 
    P
    P2
      
    Faites une intégration par parties sur le deuxième terme et en déduire la valeur de l’intégrale du départ.
  7. Écrire un programme qui détermine le degré probable du PGCD par rapport à toutes les variables de 2 polynôme à plusieurs variables en utilisant l’évaluation en toutes les variables sauf une.
  8. Calculer le pgcd par une méthode modulaire de (xyx+1)(xy+x2+1) et (xyxy)(xyx+1)
  9. En utilisant uniquement l’instruction de calcul de PGCD déterminer la multiplicité maximale d’un facteur irréductible de x14x13−14x12+12x11+78x10−54x9−224x8+116x7+361x6−129x5−330x4+72x3+160x2−16x−32

8.3  Exercices (résultant)

  1. Pour quelles valeurs de p le polynôme X5+X3pX+1 admet-il une racine multiple?
  2. Résoudre le système en éliminant successivement les variables grâce au résultant :




    a3+b3+c3=
    a2+b2+c2=
    a+b+2c=4
  3. Déterminer l’intersection de xy=4 et y2=(x−3)(x2−16), repésenter graphiquement les courbes. Discuter la multiplicité et le nombre d’intersections.
    Même question pour (x−2)2+y2=4 et y2=(x−3)(x2−16).
  4. Donner le détail des calculs avec Bézout de la décomposition en éléments simples de :
    1
    (x2−1)2(x+2)
    puis calculer le coefficient de xn du développement en séries entières de cette fraction en 0.
  5. Calculer
    1−x2
    1+x4
      dx 
    en utilisant le résultant pour calculer les logarithmes.
  6. Courbe paramétrique dépendant d’un paramètre : on considère la courbe Cm dépendant du réel m :
    x(t)=
    t+m
    t2+1+m2
    ,    y(t)=
    t2
    tm
    Représenter la courbe pour quelques valeurs de m (on pourra utiliser dans un niveau de géométrie, le menu Edit, Ajouter un paramètre pour créer un curseur représentant m, puis plotparam). Donner l’équation cartésienne de Cm. Déterminer les valeurs de m pour lesquelles la courbe admet un point singulier, représenter le graphe dans ce(s) cas et faire l’étude de la courbe.

8.4  Exercice (Bézout modulaire)

Soit A et B deux polynômes à coefficients entiers et premiers entre eux. Soit c* le résultant de A et B, on va calculer les polynômes U et V de l’identité de Bézout

  A U + B V = c ,    deg(U)<deg(B), deg(V)<deg(A)     (13)

par une méthode modulaire.

  1. Montrer, en utilisant les formules de Cramer, que les coefficients de U et V sont des entiers de valeur absolue inférieure ou égale à la borne de Hadamard h de la matrice de Sylvester de A et B (dont le déterminant est c, le résultant de A et B). Calculer h en fonction de la norme euclidienne de A, B et de leurs degrés.
  2. On calcule c* puis on résoud (13) dans /pi Z[X] pour plusieurs nombres premiers pi (choisis si possible inférieurs à √231 pour des raisons d’efficacité), puis on calcule par le théorème des restes chinois (13) dans /∏pi Z[X]. Donner une minoration de ∏i pi faisant intervenir h qui permette de garantir que l’écriture en représentation symétrique de (13) dans /∏pi Z[X] est identique à (13) dans [X].
  3. Application : résoudre de cette manière l’équation de Bézout pour
    A=(X+1)4(X−3),    B=(X−1)4(X+2)
    (vous pouvez utiliser sans justifications l’instruction de calcul de résultant, des coefficients de Bézout dans /piZ[X] et de reste chinois de votre logiciel).
  4. Écrire une fonction mettant en oeuvre cet algorithme.
  5. Que pensez-vous de l’intérêt de cet algorithme par rapport à l’algorithme d’Euclide étendu dans [X]?

8.5  Exercice (Géométrie et résultants).

On cherche une relation algébrique entre les coordonnées de 4 points A,B,C,D qui traduise le fait que ces 4 points sont cocycliques. Cette condition étant invariante par translation, on cherche une relation entre les 6 coordonnées des 3 vecteurs v1=(x1,y1), v2=(x2,y2) et v3=(x3,y3) d’origine A et d’extrémité B, C et D. On peut supposer quitte à translater que le centre du cercle est l’origine, on a donc 5 paramètres : le rayon du cercle R et les 4 angles des points sur le cercle θ0, θ1, θ2 et θ3. La relation cherchée va s’obtenir en éliminant les 5 paramètres des expressions des 6 coordonnées en fonction de ces paramètres.

  1. Exprimer les 6 coordonnées en fonction de R et a=tan(θ0/2), b=tan(θ1/2), c=tan(θ2/2) et d=tan(θ3/2). On obtient ainsi 6 équations, par exemple les deux premières sont de la forme
    x1− F(R,a,b)= 0,    y1− G(R,a,b)= 0 
    F et G sont deux fractions rationnelles.
  2. En réduisant au même dénominateur, calculer 6 polynômes, fonction de x1,y1,x2,y2,x3,y3,R,a,b,c,d, qui doivent s’annuler pour que les points soient cocycliques (Vous pouvez utiliser l’instruction numer pour obtenir le numérateur d’une fraction rationnelle).
  3. Éliminer b des polynômes contenant x1 et y1 et factoriser le polynôme obtenu, faire de même avec c, x2 et y2 et d, x3 et y3, en déduire (en supposant que les points sont tous distincts) 3 polynômes en x1,y1,x2,y2,x3,y3,R,a qui s’annulent.
  4. Éliminer R et a, en déduire la relation cherchée.
  5. Vérifier que cette relation est équivalente à la nullité de la partie imaginaire du birapport des affixes α, β, γ, δ des 4 points :
    ℑ 


    α−β
    α−γ
    δ−γ
    δ−β
     


    = 0

8.6  Décalage entier entre racines.

Soit P un polynôme à coefficients entiers sans racines multiples. On dira que P a la propriété I si deux des racines de P sont décalées d’un entier. En d’autres termes, si r1,...,rn désignent les racines complexes distinctes de P, P possède la propriété I s’il existe au moins un entier parmi les différences rirj pour ij.

  1. Soit
    R(t)=resultantx(P(x),P(x+t)) 
    Montrer que R est à coefficients entiers. Montrer que la propriété I est équivalente à la propriété “R possède une racine entière non nulle”. On va maintenant construire un algorithme déterminant les racines entières du polynôme R.
  2. Après division de R par une puissance de t, on peut supposer que R a un coefficient constant non nul. Après division de R par son contenu, on peut aussi supposer que le contenu de R est 1. En effectuant ensuite une factorisation square-free de R, on peut se ramener au cas où R et R′ sont premiers entre eux. Soit a une racine de R.
    1. Donner une majoration de |a| en fonction du coefficient constant de R.
    2. Soit p un nombre premier ne divisant pas le coefficient dominant de R et tel que R et R′ soient premiers entre eux modulo p. On peut calculer a à partir d’une racine de R modulo p en la “remontant” modulo pk pour k assez grand (algorithme p-adique). Pour quelle valeur de k peut-on reconstruire toutes les racines entières de R ?
    3. Comparer l’algorithme ci-dessus avec les algorithmes suivants : la factorisation de R sur , la recherche numérique des racines complexes de R, la recherche des racines entières de R parmi les diviseurs entiers du coefficient constant de R et leurs opposés.
  3. Une fois les racines entières de R connues, comment peut-on en déduire les facteurs de P dont les racines diffèrent de cet(ces) entier(s)?
  4. Soit
    P(x)=x6+9x5+29x4+41x3+37 x2+59x+31
    Montrer que P a la propriété I. Calculer la ou les racines entières de R et donner la factorisation correspondante de P.
  5. Écrire un programme qui effectue cet algorithme sur un polynôme quelconque. On pourra utiliser la fonction rationalroot de Xcas pour déterminer les racines entières de R.
  6. Application : on cherche à calculer
     
    n
    k=1
     
    −9x2−27x−30
    P(x)
        (14)
    Décomposer cette fraction en éléments simples (donner le détail des calculs en utilisant la factorisation précédente et l’identité de Bezout abcuv en Xcas).
  7. Calculer la somme précédente (14). On pourra remarquer que pour k entier strictement positif, 1/f(x+k)−1/f(x) s’exprime comme une somme de différences 1/f(x+j+1)−1/f(x+j).
  8. Écrire un programme effectuant ce calcul avec une fraction quelconque, lorsque cela est possible.

8.7  Exemple de correction de géométrie et résultant

e1:=x1-R*((1-b^2)/(1+b^2)-(1-a^2)/(1+a^2));
e2:=y1-R*(2b/(1+b^2)-2*a/(1+a^2));
e3:=x2-R*((1-c^2)/(1+c^2)-(1-a^2)/(1+a^2));
e4:=y2-R*(2c/(1+c^2)-2*a/(1+a^2));
e5:=x3-R*((1-d^2)/(1+d^2)-(1-a^2)/(1+a^2));
e6:=y3-R*(2d/(1+d^2)-2*a/(1+a^2));
f1:=factor(resultant(numer(e1),numer(e2),b)/
 (-4)/(a^2+1)^3/R^2);
f2:=factor(resultant(numer(e3),numer(e4),c)/
 (-4)/(a^2+1)^3/R^2);
f3:=factor(resultant(numer(e5),numer(e6),d)/
 (-4)/(a^2+1)^3/R^2);
g1:=factor(resultant(f1,f2,R));
g2:=resultant(f1,f3,R);
r:=factor(resultant(g1/(a^2+1),g2/(a^2+1),a));
eq1:=r[1,3,1];
eq2:=numer(im((-x1-i*y1)/(-x2-i*y2)*
 (x3-x2+i*(y3-y2))/(x3-x1+i*(y3-y1))));
normal(eq1-eq2);

9  Bases de Gröbner.

9.1  Ordre et réduction

L’anneau des polynômes à plusieurs variables n’a pas de division euclidienne. On est donc obligé d’utiliser des outils moins performants. La première chose à faire est de choisir un ordre total sur les monomes, compatible avec la multiplication des monômes (a<b doit entrainer a c<b c) et tel que si un monôme a divise un autre monôme b alors a<b. Exemples d’ordres utilisés fréquemment (ce sont les 3 ordres proposés par les fonctions de Xcas) :

On remarque sur ces 3 exemples qu’il ne peut exister de suite strictement décroissante infinie pour l’ordre >. Lorsque le degré total est le premier critère, c’est évident, puisque le nombre de monomes < à un monome donné est fini. Pour l’ordre lexicographique, on raisonne par l’absurde. On regarde d’abord le premier indice, comme la suite est décroissante, tous les monômes ont un indice inférieur ou égal au premier indice du premier monôme. On peut donc extraire une sous-suite strictement décroissante et infinie de monômes dont le 1er indice est constant. On passe alors au 2ème indice, et ainsi de suite jusqu’au dernier indice qui donne une contradiction. On fait donc dans la suite l’hypothèse qu’il n’existe pas de suite strictement décroissante infinie pour l’ordre >.

On peut alors effectuer une sorte de remplacement de la division euclidienne de A par B, appelée réduction qui consiste à comparer le terme dominant de B au sens de l’ordre (noté LT(B)) aux monomes de A par ordre décroissant, si l’un des monomes de A a toutes ses puissances plus grandes que LT(B), alors on élimine ce terme, disons Ak, en retranchant à A le polynôme Ak/LT(B) B. Ceci ne modifie pas le début de A jusqu’au monôme Ak. Les termes retranchés peuvent eux-même donner lieu à une réduction par B, par exemple Ak/LT(B) B2 peut être divisible par LT(B). Le procédé de réduction doit toutefois s’arrêter, sinon on pourrait construire une suite décroissante infinie pour l’ordre > avec les Ak. On peut même diviser A par plusieurs polynômes B,C,.. en utilisant cet algorithme.

9.2  Idéaux

En dimension 1, les idéaux sont engendrés par un polynôme P et c’est la division euclidienne par P qui permet de savoir si on est dans l’idéal. En dimension plus grande, l’analogue est la base de Gröbner de l’idéal (relativement à un ordre monomial <) et on utilise la réduction par rapport aux polynômes de l’idéal pour savoir si on est dans l’idéal. On commence par montrer que les idéaux de monomes sont engendrés par les monômes minimaux, qui ne sont divisibles par aucun autre monôme de l’idéal. Supposons qu’ils soient en nombre infini. Considérons le premier indice des monomes, s’il est borné, on aura une infinité de monomes ayant le même indice, sinon on aura une suite infinie de monômes d’indice croissant, dans les deux cas on peut extraire une suite infinie dont la première composante est croissante au sens large. On fait le même raisonnement sur la suite extraite pour la 2ème composante, etc. et on aboutit à une suite infinie de monômes qui se divisent les uns les autres ce qui est absurde. Donc les monômes minimaux sont en nombre fini.

Une base de Gröbner s’obtient en prenant des polynômes correspondant aux monômes minimaux de l’idéal de monômes LT(I) des coefficients dominants de I. La réduction par rapport aux éléments de cette base donne alors 0 pour tous les éléments de I, ce qui montre que I est engendré par cette base.

On appelle “s-polynôme” d’une paire de polynômes A et B le polynôme obtenu en calculant le monôme PPCM L de LT(A) et LT(B) et en créant la différence qui annule ce monôme PPCM L/LT(A)AL/LT(B)B.

On peut montrer que la base de Gröbner peut se calculer à partir d’une famille génératrice en effectuant la boucle suivante : on calcule tous les s-polynômes de la famille génératrice courante, on les réduit par rapport à la famille génératrice courante, si tous les s-polynomes sont nuls la famille courante est la base cherchée, sinon on garde les s-polynômes réduits non nuls, on réduit la famille génératrice courante par rapport à ces s-polynômes réduits non nuls et on fusionne les polynômes non nuls en la famille génératrice courante pour l’itération suivante de la boucle.

Le problème est que cela devient très vite très long. Il existe des méthodes permettant d’accélérer l’algorithme, par exemple on peut savoir à l’avance qu’un s-polynôme se réduit à 0 (règles de Gebauer-Möller) il est donc inutile de le calculer. On peut aussi précalculer tous les multiples des polynômes par rapport auxquels on réduit et réduire simultanément tous les polynômes à réduire en ramenant la réduction à un algorithme de pivot de Gauß (c’est la partie algèbre linéaire de l’algorithme F4). L’ordre choisi est aussi très important pour l’efficacité. Enfin, pour le cas des coefficients entiers, des méthodes modulaires permettent d’accélérer les calculs. Xcas implémente un algorithme modulaire très compétitif pour l’ordre revlex, présenté dans l’article en anglais qui suit.

Les instructions Xcas correspondantes sont gbasis, greduce.

9.3  Introduction

During the last decades, considerable improvements have been made in CAS like Maple or specialized systems like Magma, Singular, Cocoa, Macaulay... to compute Groebner basis. They were driven by implementations of new algorithms speeding up the original Buchberger ([]) algorithm: Gebauer and Möller criterion ([]), F4 and F5 algorithms from J.-C. Faugère ([], []), and are widely described in the literature if the base field is a finite field. Much less was said about computing over . It seems that implementers are using the same algorithm as for finite fields, this time working with coefficients in or in (sometimes with fast integer linear algebra), despite the fact that an efficient p-adic or Chinese remaindering algorithm were described as soon as in year 2000 by E. Arnold ([]). The reason might well be that these modular algorithms suffer from a time-consuming step at the end: checking that the reconstructed Groebner basis is indeed the correct Groebner basis.

Section 9.4 will show that if one accepts a small error probability, this check may be fast, so we can let the user choose between a fast conjectural Groebner basis to make his own conjectures and a slower certified Groebner basis once he needs a mathematical proof.

Section 9.5 will explain learning, a process that can accelerate the computation of a Groebner basis modulo a prime pk once the same computation but modulo another prime p has already been done ; learning is an alternative to the F5 algorithm in order to avoid computing useless critical pairs that reduce to 0. The idea is similar to F4remake by Joux-Vitse ([]) used in the context of computing Groebner basis in large finite fields.

Section 9.6 will show in more details how the gbasis algorithm is implemented in Giac/Xcas ([]) and show that - at least for the classical academic benchmarks Cyclic and Katsura - the deterministic modular algorithm is competitive or faster than the best open-source implementations and the modular probabilistic algorithm is comparable to Maple and slower than Magma on one processor (at least for moderate integer coefficient size) and may be faster than Magma on multi-processors, while computation modulo p are faster for characteristics in the 24-31 bits range. Moreover the modular algorithm memory usage is essentially twice the memory required to store the basis on , sometimes much less than the memory required by other algorithms.

9.4  Checking a reconstructed Groebner basis

Let f1,..,fm be polynomials in [x1,..,xn], I=<f1,...,fm> be the ideal generated by f1,...,fn. Without loss of generality, we may assume that the fi have coefficients in by multiplying by the least common multiple of the denominators of the coefficients of fi. We may also assume that the fi are primitive by dividing by their content.

Let < be a total monomial ordering (for example revlex the total degree reverse lexicographic ordering). We want to compute the Groebner basis G of I over (and more precisely the inter-reduced Groebner basis, sorted with respect to <). Now consider the ideal Ip generated by the same fi but with coefficients in /p for a prime p. Let Gp be the Groebner basis of Ip (also assumed to be inter-reduced, sorted with respect to <, and with all leading coefficients equal to 1).

Assume we compute G by the Buchberger algorithm with Gebauer and Möller criterion, and we reduce in (by multiplying the s-poly to be reduced by appropriate leading coefficients), if no leading coefficient in the polynomials are divisible by p, we will get by the same process but computing modulo p the Gp Groebner basis. Therefore the computation can be done in parallel in and in /p except for a finite set of unlucky primes (since the number of intermediate polynomials generated in the algorithm is finite). If we are choosing our primes sufficiently large (e.g. about 231), the probability to fall on an unlucky prime is very small (less than the number of generated polynomials divided by about 231, even for really large examples like Cyclic9 where there are a few 104 polynomials involved, it would be about 1e-5).

The Chinese remaindering algorithm is as follow: compute Gp for several primes, for all primes that have the same leading monomials in Gp (i.e. if coefficient values are ignored), reconstruct Gpj by Chinese remaindering, then reconstruct a candidate Groebner basis Gc in by Farey reconstruction. Once it stabilizes, do the checking step described below, and return Gc on success.

Checking step : check that the original fi polynomials reduce to 0 with respect to Gc and check that Gc is a Groebner basis.

Théorème 14   (Arnold, Greuel and Pfister) If the checking step succeeds, then Gc is the Groebner basis of I.

This is a consequence of ideal inclusions (first check) and dimensions (second check), for a complete proof, see [] and theorem 7.5.1 in Greuel,G.-M., Pfister,G., 2007, A Singular Introduction to Commutative Algebra, Springer.

Probabilistic checking algorithm: instead of checking that s-polys of critical pairs of Gc reduce to 0, check that the s-polys reduce to 0 modulo several primes that do not divide the leading coefficients of Gc and stop as soon as the inverse of the product of these primes is less than a fixed ε>0.

Deterministic checking algorithm: check that all s-polys reduce to 0 over . This can be done either by integer computations (or even by rational computations, I have not tried that), or by reconstruction of the quotients using modular reduction to 0 over /p for sufficiently many primes. Once the reconstructed quotients stabilize, we can check the 0-reduction identity, and this can be done without computing the products quotients by elements of Gc if we have enough primes (with appropriate bounds on the coefficients of Gc and the lcm of the denominators of the reconstructed quotients).

9.5  Speeding up by learning from previous primes

Once we have computed a Groebner basis modulo an initial prime p, if p is not an unlucky prime, then we can speedup computing Groebner basis modulo other lucky primes. Indeed, if one s-poly reduce to 0 modulo p, then it reduces most certainly to 0 on (non zero s-poly have in general several terms, cancellation of one term mod p has probability 1/p, simultaneous cancellation of several terms of a non-zero s-poly modulo p is highly improbable), and we discard this s-poly in the next primes computations. We name this speedup process learning. It can also be applied on other parts of the Groebner basis computation, like the symbolic preprocessing of the F4 algorithm, where we can reuse the same collection of monomials that were used for the first prime p to build matrices for next primes (see Buchberger Algorithm with F4 linear algebra in the next section).

If we use learning, we have no certification that the computation ends up with a Groebner basis modulo the new primes. But this is not a problem, since it is not required by the checking correctness proof, the only requirement is that the new generated ideal is contained in the initial ideal modulo all primes (which is still true) and that the reconstructed Gc is a Groebner basis.

9.6  Giac/Xcas implementation and experimentation

We describe here briefly some details of the Giac/Xcas gbasis implementation and give a few benchmarks.

The optimized algorithm runs with revlex as < ordering if the polynomials have at most 15 variables (it’s easy to modify for more variables, adding multiples of 4, but this will increase a little memory required and slow down a little). Partial and total degrees are coded as 16 bits integers (hence the 15 variables limit, since 1 slot of 16 bits is kept for total degree). Modular coefficients are coded as 31 bit integers (or 24).

The Buchberger algorithm with linear algebra from the F4 algorithm is implemented modulo primes smaller than 231 using total degree as selection criterion for critical pairs.
Buchberger algorithm with F4 linear algebra modulo a prime

  1. Initialize the basis to the empty list, and a list of critical pairs to empty
  2. Add one by one all the fi to the basis and update the list of critical pairs with Gebauer and Möller criterion, by calling the gbasis update procedure (described below step 9)
  3. Begin of a new iteration:
    All pairs of minimal total degree are collected to be reduced simultaneously, they are removed from the list of critical pairs.
  4. The symbolic preprocessing step begins by creating a list of monomials, gluing together all monomials of the corresponding s-polys (this is done with a heap data structure).
  5. The list of monomials is “reduced” by division with respect to the current basis, using heap division (like Monagan-Pearce []) without taking care of the real value of coefficients. This gives a list of all possible remainder monomials and a list of all possible quotient monomials and a list of all quotient times corresponding basis element monomial products. This last list together with the remainder monomial list is the list of all possible monomials that may be generated reducing the list of critical pairs of maximal total degree, it is ordered with respect to <. We record these lists for further primes during the first prime computation.
  6. The list of quotient monomials is multiplied by the corresponding elements of the current basis, this time doing the coefficient arithmetic. The result is recorded in a sparse matrix, each row has a pointer to a list of coefficients (the list of coefficients is in general shared by many rows, the rows have the same reductor with a different monomial shift), and a list of monomial indices (where the index is relative to the ordered list of possible monomials). We sort the matrix by decreasing order of leading monomial.
  7. Each s-polynomial is written as a dense vector with respect to the list of all possible monomials, and reduced with respect to the sparse matrix, by decreasing order with respect to <. (To avoid reducing modulo p each time, we are using a dense vector of 128 bits integers on 64 bits architectures, and we reduce mod p only at the end of the reduction. If we work on 24 bit signed integers, we can use a dense vector of 63 bits signed integer and reduce the vector if the number of rows is greater than 215).
  8. Then inter-reduction happens on all the dense vectors representing the reduced s-polynomials, this is dense row reduction to echelon form (0 columns are removed first). Care must be taken at this step to keep row ordering when learning is active.
  9. gbasis update procedure
    Each non zero row will bring a new entry in the current basis (we record zero reducing pairs during the first prime iteration, this information will be used during later iterations with other primes to avoid computing and reducing useless critical pairs). New critical pairs are created with this new entry (discarding useless pairs by applying Gebauer-Möller criterion). An old entry in the basis may be removed if it’s leading monomial has all partial degrees greater or equal to the leading monomial corresponding degree of the new entry. Old entries may also be reduced with respect to the new entries at this step or at the end of the main loop.
  10. If there are new critical pairs remaining start a new iteration (step 3). Otherwise the current basis is the Groebner basis.

Modular algorithm

  1. Set a list of reconstructed basis to empty.
  2. Learning prime: Take a prime number of 31 bits or 29 bits for pseudo division, run the Buchberger algorithm modulo this prime recording symbolic preprocessing data and the list of critical pairs reducing to 0.
  3. Loop begin: Take a prime of 29 bits size or a list of n primes if n processors are available. Run the Buchberger algorithm. Check if the output has the same leading terms than one of the chinese remainder reconstructed outputs from previous primes, if so combine them by Chinese remaindering and go to step 4, otherwise add a new entry in the list of reconstructed basis and continue with next prime at step 3 (clearing all learning data is probably a good idea here).
  4. If the Farey -reconstructed basis is not identical to the previous one, go to the loop iteration step 3 (a fast way to check that is to reconstruct with all primes but the last one, and check the value modulo the last prime). If they are identical, run the final check : the initial polynomials fi must reduce to 0 modulo the reconstructed basis and the reconstructed basis s-polys must reduce to 0 (this is done on either directly or by modular reconstruction for the deterministic algorithm, or checked modulo several primes for the probabilistic algorithm). On success output the Groebner basis, otherwise continue with next prime at step 3.

Benchmarks
Comparison of giac (1.1.0-26) with Singular 3.1 from sage 5.108 on Mac OS X.6, Dual Core i5 2.3Ghz, RAM 2*2Go:

 giac mod pgiacsingulargiac prob.giac singular
 24, 31 bitsrun2mod p1e-7, 1e-16 certifiedstd
Cyclic70.5, 0.580.12.03.5, 4.221, 29.3>2700
Cyclic87.2, 8.91.852.5103, 106258, 679>>
Cyclic9633, 1340200?1 day>>>>
Kat80.063, 0.0740.0090.20.33, 0.536.55, 4.354.9
Kat90.29, 0.390.051.372.1, 3.254, 3641
Kat101.53, 2.270.311.6514, 20.7441, 335480
Kat1110.4, 13.82.886.8170, 2104610?
Kat1276, 103278851950, RAMRAM>>
alea60.83, 1.08.264.18202, 204738, >>>1h


This leads to the following observations :

Example of Giac/Xcas code:

alea6 := [5*x^2*t+37*y*t*u+32*y*t*v+21*t*v+55*u*v,
39*x*y*v+23*y^2*u+57*y*z*u+56*y*u^2+10*z^2+52*t*u*v,
33*x^2*t+51*x^2+42*x*t*v+51*y^2*u+32*y*t^2+v^3,
44*x*t^2+42*y*t+47*y*u^2+12*z*t+2*z*u*v+43*t*u^2,
49*x^2*z+11*x*y*z+39*x*t*u+44*x*t*u+54*x*t+45*y^2*u,
48*x*z*t+2*z^2*t+59*z^2*v+17*z+36*t^3+45*u];
l:=[x,y,z,t,u,v];
p1:=prevprime(2^24); p2:=prevprime(2^29);
time(G1:=gbasis(alea6 % p1,l,revlex));
time(G2:=gbasis(alea6 % p2,l,revlex));
threads:=2; // set the number of threads you want to use
// debug_infolevel(1); // uncomment to show intermediate steps
proba_epsilon:=1e-7; // probabilistic algorithm.
time(H0:=gbasis(alea6,indets(cyclic5),revlex));
proba_epsilon:=0; // deterministic
time(H1:=gbasis(alea6,indets(cyclic5),revlex));
time(H2:=gbasis(alea6,indets(cyclic5),revlex,modular_check));
size(G1),size(G2),size(H0),size(H1),size(H2);
write("Halea6",H0);

Note that for small examples (like Cyclic5), the system performs always the deterministic check (this is the case if the number of elements of the reconstructed basis to 50).

9.7  Conclusion

I have described some enhancements to a modular algorithm to compute Groebner basis over which, combined to linear algebra from F4, gives a sometimes much faster open-source implementation than state-of-the-art open-source implementations for the deterministic algorithm. The probabilistic algorithm is also not ridiculous compared to the best publicly available closed-source implementations, while being much easier to implement (about 10K lines of code, while Fgb is said to be 200K lines of code, no need to have highly optimized sparse linear algebra).

This should speed up conjectures with the probabilistic algorithm and automated proofs using the deterministic algorithm (e.g. for the Geogebra theorem prover []), either using Giac/Xcas (or one of it’s interfaces to java and python) or adapting it’s implementation to other open-source systems. With fast closed-source implementations (like maple or magma), there is no certification that the result is a Groebner basis : there might be some hidden probabilistic step somewhere, in integer linear system reduction for example. I have no indication that it’s the case but one can never know if the code is not public, and at least for my implementation, certification might take a lot more time than computation.

There is still room for additions and improvements

Acknowledgements
Thanks to Frédéric Han for interfacing giac with Python. Thanks to Vanessa Vitse for insightfull discussions.

9.8  Représentation rationnelle univariée (rur).

Lorsqu’on résoud un système polynomial, on a (en général) autant d’équations que d’inconnues et en principe un nombre fini de solutions. On peut utiliser une base de Groebner dans l’ordre lexicographique, résoudre par rapport à la dernière variable, puis remonter, mais d’une part le calcul d’une base de Groebner dans l’ordre lexicographique est significativement plus long que dans l’ordre revlex, et d’autre part il faut calculer des PGCD et factoriser des polynômes sur des extensions algébriques dont la taille peut augmenter au fur et à mesure que l’on remonte (ou faire des calculs approchés...). Il serait plus intéressant de calculer d’un seul coup une extension algébrique de qui permette d’exprimer toutes les variables. Ceci peut se faire si on arrive à trouver une forme linéaire en les variables qui sépare les solutions (la valeur de la forme est distincte si les points solutions sont distincts). On rajoute cette variable et on résoud l’équation obtenue en cette variable, pour chaque solution on aura une unique solution en remontant les autres variables. La représentation univariée rationnelle fait précisément cela, et donne même les autres variables comme polynôme en la forme linéaire séparante.

La présentation classique de la représentation univariée rationnelle utilise des calculs de trace (cf. par exemple le rapport de l’Inria 1998 de Fabrice Rouillier), l’algorithme implémenté dans Giac/Xcas (versions 1.1.1 et ultérieures) est un algorithme modulaire. On commence par se ramener au cas d’un idéal radical (c’est-à-dire que les points solutions du système sont de multiplicité 1) en ajoutant aux générateurs de l’idéal les parties squarefree des polynômes minimaux de toutes les variables. Pour un idéal radical, on montre qu’il existe une forme linéaire séparante, le degré du polynôme minimal de cette forme linéaire séparante est exactement égal à la dimension du quotient de l’anneau de polynômes par l’idéal radical. On peut donc tester si une forme linéaire est séparante en calculant son polynôme minimal. En pratique, on commence par calculer une base de Groebner pour l’ordre revlex (le plus efficace). On génère la liste des monomes du quotient en commençant par majorer les degrés en chacune des variables, puis on élimine parmi les monomes possibles ceux qui sont divisibles par le monome dominant d’un élément de la base de Groebner. On calcule ensuite la classe d’un polynôme dans le quotient en effectuant une réduction par la base de Groebner, on obtient un vecteur de coordonnées dans cette base de monome. Le calcul du polynôme minimal d’une forme linéaire devient ainsi un problème d’algèbre linéaire. Le calcul de chaque variable en fonction des puissances d’une forme linéaire séparante est également un problème d’algèbre linéaire (on le fait simultanément pour toutes les variables, si on veut optimiser on peut même faire une décomposition LU lors du calcul du polynôme minimal et la réutiliser). Pour éviter les problèmes de croissance de coefficients dans les calculs intermédiaires, ce calcul est effectué modulo plusieurs nombres premiers dans giac, jusqu’à pouvoir reconstruire par les restes chinois le polynôme minimal de la forme séparante sur et les expressions des variables comme polynôme de la forme séparante (on n’a alors pas besoin de reconstruire la base de Groebner sur ). Bien entendu, il faut traiter le cas des mauvaises réductions, pour cela on regarde si les monomes de la base du quotient de l’anneau par l’idéal sont indépendants du nombre premier choisi, en cas de différence, il faut conserver le nombre premier correspondant à la liste de monômes la plus grande (l’autre nombre premier est de mauvaise réduction), ou rejeter les deux nombres premiers si aucune des deux listes de monomes ne contient l’autre.

Les fonctions solve, fsolve et cfsolve utilisent cet algorithme pour des systèmes polynomiaux qui s’y prêtent (en cherchant une forme séparante d’abord parmi les variables puis avec des combinaisons linéaires aléatoires à petits coefficients entiers), solve essaie de renvoyer des solutions exactes si le polynome minimal de la forme linéaire séparante se factorise sur , fsolve (en mode réel) localise les racines réelles par la méthode d’Akritas, cfsolve localise les racines complexes par factorisation de Schur de la matrice companion. La fonction gbasis(eqs,vars,rur) avec comme paramètre optionnel rur effectue le calcul de la représentation univariée rationnelle et renvoie une liste contenant le polynôme minimal P (exprimée arbitrairement en fonction de la 1ère variable du système), sa dérivee P′ et les P1,...,Pn qui permettent d’exprimer la i-ème variable d’une solution comme étant Pi(r)/P′(r) avec r racine de P. On peut alors vérifier que l’on a bien une solution en remplaçant la variable xi par Pi/P′ dans les équations, le reste de la division euclidienne du numérateur de la fraction obtenue par le polynome minimal P doit donner 0.

Exemple :

Calcul en mode pas à pas :
l’élément d’indice 1 est la forme linéaire séparante, en indice 2, le polynôme minimal de l’élément séparant exprimé en fonction de x, en indice 3 sa dérivée qui sera le dénominateur commun de la solution, en indices de 4 à la fin le numérateur de la solution en fonction de x
On vérifie :
On trouve les solutions approchées

La représentation rationnelle univariée a des applications au-delà de la seule résolution de systèmes polynomiaux. On peut s’en servir pour trouver une extension algébrique unique de permettant de calculer toutes les racines d’un polynôme, il suffit de poser le système formé par les relations racines-coefficients de ce polynôme et d’en chercher la représentation rationnelle univariée, cf. la section 18.6. On peut également s’en servir pour trouver une extension algébrique unique contenant plusieurs extensions de dont on a le polynôme minimal. Par exemple pour travailler dans [√2,√3,√5], on pose
G:=gbasis([a^2-2,b^2-3,c^2-5],[a,b,c],rur),
on a alors ± √2=rootof(G[4],G[2])/rootof(G[3],G[2]),
± √3=rootof(G[5],G[2])/rootof(G[3],G[2]),
± √5=rootof(G[6],G[2])/rootof(G[3],G[2])
(on peut utiliser normal ou evalf pour décider du signe).

10  Courbes paramétriques et polaires

10.1  Introduction

Le graphe d’une fonction f: I ↦ (I un intervalle) est un exemple de courbe du plan, mais il n’est pas assez général pour représenter tous les types de courbe du plan, par exemple un segment de droite vertical, ou un cercle, car deux points distincts d’un graphe doivent avoir des abscisses différentes. D’autre part, il apparait naturellement d’autres types de courbes que les graphes de fonction, par exemple la trajectoire d’un mobile dans le plan dont les coordonnées x,y dépendent du temps (selon une équation différentielle ou un système différentiel), ce sont les courbes paramétriques, ou des courbes vérifiant une équation cartésienne (par exemple en géométrie le cercle x2+y2=1, ou en cinématique des courbes de niveau de l’énergie totale dans le plan position-impulsion) ce sont les courbes implicites.

Dans cette section, on va étudier les courbes en paramétriques, donnée par un couple de fonctions (x(t),y(t)) définies pour t dans un sous-ensemble des réels et à valeurs dans . (Ceci ne restreint pas trop la généralité, on peut montrer sous des hypothèses assez générales que l’allure locale d’une courbe implicite est identique à celle d’une courbe paramétrique, sauf en certains points dits singuliers, c’est le théorème des fonctions implicites).

Exemples :

10.2  Représentation graphique

La plupart des calculatrices graphiques et de nombreux logiciels de maths permettent de représenter graphiquement un arc de courbe en donnant des valeurs extrêmes t et t+ (souvent notées tmin et tmax) et un pas Δ t (tstep). Le logiciel évalue la valeur de x(t) et y(t) en t, tt, t+2Δ t, ... puis relie les points de la courbe obtenue par des segments (parfois avec des autres arcs de courbes). La plupart du temps cela donne une bonne idée de la courbe, mais parfois on peut manquer un détail intéressant (valeur de Δ t trop grande), ou un morceau de courbe (mauvaises valeurs de t et t+). Par exemple

Il peut être nécessaire d’ajuster le cadrage graphique à l’affichage (xmin, xmax, ymin, ymax) ou de l’affiner avec un menu de zoom. Sur les calculatrices les opérations de changement de cadrage graphique provoquent un nouveau calcul complet qui peut durer une dizaine de secondes.

Mise en oeuvre :

Exemples : essayez de tracer quelques courbes en paramétriques

(2cos(t),3sin(t)),    (cos(2t),sin(3t)),    (t2,t3),    (t+1/tt2+2/t),    (
t−1
,
2−t

10.3  Paramétrage

On adoptera souvent la convention d’appeler temps le paramétre t. Mais cela ne signifie pas que le paramétrage est réellement le temps mesuré en secondes. On peut très bien paramétrer une courbe avec un paramètre autre, qui peut être un multiple constant ou variable du temps (c’est d’ailleurs conforme au principe de la relativité). Le paramétrage n’est jamais unique, on peut changer de paramétrage pourvu que la fonction donnant le nouveau en fonction de l’ancien paramétrage soit une bijection (qui peut même renverser le sens de déroulement du temps c’est-à-dire le sens de parcours de la courbe). On utilisera d’ailleurs plus loin un paramétrage par la longueur, où la courbe est parcourue à vitesse constante égale à 1.

Le choix d’un paramétrage est ce qui fait la différence entre la cinématique (on prend le temps comme paramètre) et la géométrie (où on cherche à décrire les propriétés intrinséques de la courbe indépendamment du paramétrage). Ainsi, l’équation cartésienne d’une courbe est une propriété géométrique, indépendante du choix de paramétrage choisi pour l’obtenir.

On observe aussi que l’opération inverse, trouver un paramétrage à partir d’une équation cartésienne de courbe n’est pas possible de manière explicite, sauf dans quelques cas particuliers. C’est pour cette raison qu’il est beaucoup plus difficile (et couteux en temps) d’obtenir une représentation graphique d’une courbe donnée par son équation cartésienne.

10.4  Étude analytique d’une courbe en paramétrique

10.4.1  Rappel sur les graphes de fonction

Pour tracer le graphe d’une fonction f, on commence par déterminer le domaine de définition, on restreint éventuellement l’intervalle d’étude (parité, périodicité). Puis on calcule les limites aux bornes du domaine de définition :

On calcule la dérivée première de f pour déterminer le sens de variations, les points d’annulation correspondent à des tangentes horizontales. On peut étudier la convexité de f (signe de f′′), les points d’inflexion de la courbe se produisent lorsque f′′ s’annule. On trace alors la courbe en faisant apparaitre les points particuliers et les asymptotes.

Pour une courbe en paramétrique, le plan général est analogue, mais l’étude est un peu plus compliquée puisqu’on a deux fonctions tx(t) et ty(t) au lieu d’une seule xy(x).

10.4.2  Domaine et périodicité

On supposera dans toute la suite que les fonctions x(t) et y(t) sont continument dérivables au moins 2 fois, sauf peut-être en un nombre fini de réels d’un intervalle I de .

On commence par déterminer le domaine de définition de x(t) et de y(t), et on essaie de le réduire si possible, soit par périodicité (par exemple pour le cercle ci-dessus, t ∈ [0,2 π]) soit par l’existence de symétries si les fonctions x(t) et y(t) sont paires ou impaires. Par exemple, si x et y sont paires, alors on parcourt deux fois le même arc de courbe sur + et , on peut restreindre le domaine d’étude à t≥ 0. Si x est pair et y impair, alors (x(−t),y(−t))=(x(t),−y(t)), il y a une symétrie par rapport à l’axe des x, on se restreint à tR+. Dans le cas périodique, on peut tester des symétries correspondant à des demi (voire quart) de période. Exemple : (3cos(t)+2cos(3t),3sin(t)−2sin(3t))

10.4.3  Branches infinies

On s’intéresse ensuite aux bornes du domaine de définition et aux points où x ou/et y ne sont pas définis. Si x et y admettent une limite finie, on peut prolonger la courbe. Si les limites existent mais ne sont pas finies, on a une branche infinie (x ou y). Si l’une des deux valeurs tend vers l’infini, l’autre restant finie, on a une asymptote (horizontale si x tend vers l’infini, verticale si y tend vers l’infini), on peut déterminer la position de l’arc de courbe par rapport à l’asymptote en cherchant le signe de yl ou xl lorsque t tend vers la valeur particulière (limite à droite et limite à gauche). Enfin si x et y tendent vers l’infini tous les deux, on cherche la limite de y/x, Si y/xa ≠ 0, on a une branche parabolique de direction asymptotique y=ax, on cherche alors la limite yax, si cette limite est finie et vaut b on a une asymptote oblique y=ax+b (on peut déterminer la position en cherchant le signe de y−(ax+b).

Exemples :

(
t2
t+1
,t+
1
t2+1
),  (t2,t3),  (
t3
t2+1
,t+
1
t2+2
),  (
1
t2−1
,
t
t+1
), 

On peut utiliser la commande limit dans Xcas pour étudier une asymptote, par exemple dans le premier cas, pour étudier la branche infinie pour t → +∞9
On définit les fonctions

puis on calcule les limites lorsque t → +∞


Si a est fini non nul et b fini, on en déduit que y=ax+b est asymptote à la courbe. Il y a une autre asymptote pour t=−1 (Y fini, X tend vers l’infini)

Autre exemple :

10.4.4  Étude locale

On se place en une valeur de t0x et y sont continument dérivables au moins deux fois. On notera la dérivation par rapport au paramètre par le signe ’ (en physique on utilise aussi le point). On a alors un développement de Taylor à l’ordre 2 du vecteur


M(t0)M(t0+h)
 
=(x(t0+h)−x(t0),y(t0+h)−y(t0)) 
 =
h (x′(t0),y′(t0))+
h2
2
(x′′(tx),y′′(ty))

tx et ty sont compris entre t0 et t0+h. Si le vecteur vitesse v=(x′(t0),y′(t0)) est non nul, on en déduit un équivalent


M(t0)M(t0+h)
 
≈ h (x′(t0),y′(t0))

Lorsque h est proche de 0, le vecteur M(t0)M(t0+h) est équivalent à un vecteur colinéaire à v=(x′(t0),y′(t0)) (supposé non nul), qui est donc vecteur tangent à la courbe en (x(t0),y(t0)).

Définition 15   On appelle point régulier d’une courbe paramétrique un point où la vitesse v(t)=(x′(t),y′(t)) est non nulle. En un point régulier, la courbe est tangente au vecteur vitesse (la direction du vecteur vitesse est donc une propriété géométrique, alors que le vecteur vitesse est une propriété cinématique). On notera en particulier que la tangente est horizontale si y′=0 et verticale si x′=0.

On appelle point singulier un point où la vitesse s’annulle.

On verra dans la suite comment étudier la tangente en un point singulier d’une courbe. Génériquement, une courbe n’a pas de points singuliers, car il faut annuler simultanément les deux dérivées, or on n’a qu’un seul paramètre libre t. Par contre une famille de courbes (xm(t),ym(t)) dépendant d’un paramètre m (par exemple xm(t)=t2mt, ym(t)=m/(1+t2)+t) possède en général un nombre discret de valeurs du paramètre pour lesquelles la courbe admet un point singulier. Dans l’exemple, xm′=2tm, ym′=−2mt/(1+t2)2+1, les deux dérivées s’annulent si m=−2 (en t=−1, x=−1, y=−2) ou m=2 (en t=1). Commandes Xcas :
x:=t^2-m*t; y:=m/(1+t^2)+t;
solve([diff(x,t),diff(y,t)],[m,t]);
supposons(m=[-2.0,-5,5,0.1]);
plotparam([x,y],t=((-3) .. 3));

=
Not evaled

Remarque : en cinématique, si la vitesse et l’accélération sont nulles en un point et que les équations ne dépendent pas explicitement du temps, on reste indéfiniment en ce point qui est un point d’équilibre, la notion de tangente à la courbe n’a alors pas de sens. On peut aussi suivre une trajectoire qui se rapproche de plus en plus d’un point d’équilibre (la limite de (x(t),y(t)) est alors ce point, pour t → +∞ si l’équilibre est stable ou t → − ∞ si l’équilibre est instable).

Pour faire une étude locale plus précise dans le cas d’un point régulier, ou pour déterminer la tangente en un point singulier, il faut poursuivre le développement de Taylor à un ordre plus grand. Á l’ordre 2, si x et y sont 3 fois continument dérivables, il existe tx,ty∈ [t0,t0+h] tels que :


M(t0)M(t0+h)
 
h (x′(t0),y′(t0))+
h2
2
(x′′(t0),y′′(t0)) +
h3
6
(x′′′(tx),y′′′(ty))

Si les vecteurs vitesse v=(x′(t0),y′(t0)) et accélération a=(x′′(t0),y′′(t0)) ne sont pas colinéaires, alors {v,a} forme une base, et dans cette base M(t0)M(t0+h) a pour coordonnées (h,h2/2)+un terme d’ordre 3 en puissances de h, l’arc de courbe est à l’ordre 2 identique à un arc de parabole. On parle de point birégulier. Si {v,a} est une base directe, l’arc est convexe (la vitesse “tourne” dans le sens trigonométrique), sinon il est concave. On peut tester cela en calculant le déterminant des coordonnées de {v,a} ou le sens de variations de m, la pente de la tangente

m=
y
x
,    m′=
xy′′−x′′y
x2
 
Théorème 16   Si det(v,a)=xy′′−x′′y′>0 [resp <0] sur un intervalle du domaine de définition, la courbe n’a que des points réguliers, la direction de la tangente en un point est donnée par le vecteur vitesse, et la courbe est convexe [resp. concave]. Si xy′′−x′′y′=0, on parle de point d’inflexion analytique.

Remarque : pour un graphe de fonction, x=t, on retrouve le critère usuel y′′>0.

Exemple : point d’inflexion en t=0 de

(
t3
t2+1
,t+
1
t2+2

La courbe admet deux autres points d’inflexion (t=−3.16... et t=1.31...) qu’on peut déterminer avec les commandes Xcas suivantes :






Note : on utilise comme paramètre x au lieu de t pour pouvoir utiliser la notation ' pour dériver (si on utilise t comme paramètre, il faut utiliser diff(.,t) pour calculer la dérivée par rapport à t). L’instruction fsolve effectue une résolution numérique, pour tenter une résolution exacte, utiliser solve, mais on risque alors de manquer certaines solutions.

On observe que la convexité est (presque) une propriété géométrique, en effet si on change de paramétrage

x′=
dx
dt
 =
dx
ds
 s′ 

on dérive par rapport à t

x′′ = (
dx
ds
 s′)′=
d2x
ds2
 s2 + 
dx
ds
s′′ 

puis :

xy′′− y′ x′′ = 
dx
ds
 s′ (
d2y
ds2
 s2 +
dy
ds
 s′′ ) − 
dy
ds
 s′ (
d2x
ds2
 s2 +
dx
ds
 s′′ )  = s3 (
dx
ds
 
d2y
ds2
 − 
dy
ds
 
d2x
ds2
 )

on retrouve en facteur s3 qui est positif si on parcourt la courbe dans le même sens ou négatif sinon.

La convexité décrit qualitativement la géométrie de la courbe (orientée) à l’ordre 1. On verra plus loin que le rayon de courbure décrit quantitativement la géométrie de la courbe à l’ordre 2 (comme la tangente décrit la géométrie de la courbe à l’ordre 1).

Dans le cas d’un point singulier (v=0), si l’accélération a≠ 0, alors la tangente est portée par a. L’étude complète de la nature d’un point singulier ou de la convexité d’un point régulier tel que a est colinéaire à v nécessite de faire un développement de Taylor en t=t0 jusqu’au premier ordre q, s’il existe, tel que :

Dans la base { T,A}, les composantes de M(t0)M(t) sont alors respectivement équivalentes à hp/p! et hq/q! où h=tt0. On en déduit que la tangente à la courbe est portée par T.

Exemples de points singuliers en t=0 avec dans l’ordre rebroussement de 1ère puis 2ième espèce, méplat et inflexion :

(t2,t3),  (t2+t4,t4+t5),  (t3,t4),  (t3,t5






Les deux derniers cas peuvent être reparamétrés (au prix de la perte de dérivabilité seconde) en posant t′=t1/3.

Pour faire l’étude d’un point singulier avec Xcas, on peut utiliser la fonction series sur x(t) et y(t) (ici c’est inutile, le développement de Taylor est déjà fait).

Remarque : il peut arriver dans des cas pathologiques que toutes les dérivées de (x,y) s’annulent en un point sans que la fonction (x,y) soit nulle (par exemple si x et y contiennent un facteur exp(−1/t2) en t=0, on parle de fonction plate). Il peut aussi arriver que toutes les dérivées soit colinéaires à la première dérivée non nulle, si on se déplace sur une droite ou si la tangeance est plate.

10.5  Plan d’étude d’une courbe

  1. On détermine et on restreint le domaine de définition (périodicité, symétries).
  2. On étudie les branches infinies (point exclus du domaine, ± ∞) : asymptotes horizontales, verticales, directions asymptotiques, asymptotes obliques.
  3. Recherche de x′ et y′, on étudie l’annulation conjointe des deux (points singuliers).
  4. Signe de x′ et y′, double tableau de variations faisant apparaitre x,x′,y,y′ et mise en évidence des tangentes horizontales et verticales
  5. Pour préciser le tracé, on peut chercher la convexité en étudiant le signe de xy′′−x′′y′.
  6. Tracé des points remarquables et des asymptotes et on les relie entre eux en suivant les sens de variations du tableau de variations.

Exemple

10.6  Courbes en polaires

Une courbe en polaire est essentiellement donnée par la distance au centre O d’un point M de la courbe en fonction de l’angle θ entre la direction Ox et le vecteur OM :

OM = r(θ)

On s’autorise toutefois des valeurs négatives pour r, si c’est le cas, on prend alors le symétrique par rapport à l’origine du point situé à distance −r et d’angle θ.

Représentation graphique : avec Xcas, on utilise la commande plotpolar, sur calculatrices graphiques, sélectionner le mode de tracé en polaire (touche Mode sur TI89/92/V200) ou l’application Polaire ou Géométrie sur les HP Prime. Par exemple





Remarque : une courbe en polaires est toujours parcourue dans le sens trigonométrique.

C’est un cas particulier de courbe en paramétriques puisque

(x,y)=(r(θ) cos(θ), r(θ) sin(θ))

mais on préfère souvent faire l’étude directement sur la fonction r. Le plan d’étude est calqué sur celui d’une courbe en paramétrique, mais on n’a qu’une seule fonction r à étudier.

  1. domaine de définition de r, recherche de périodicités et symétries (θ → −θ ou ajout d’une demi ou d’un quart de période). Si la période n’est pas un multiple de 2π, cela correspond à obtenir un arc de la courbe par rotation à partir d’un autre arc de la courbe.
  2. branches infinies pour θ0 (non infini) où r n’est pas défini. La branche a pour direction asymptotique la droite faisant un angle θ0 avec l’axe des x. On calcule alors la limite si elle existe de r sin(θ−θ0)), c’est l’ordonnée dans le repère obtenu par rotation d’angle θ0, si la limite est finie et vaut l on a une asymptote (d’équation Y=l dans le repère tourné).
    Exemple r=1/(1+2cos(θ)). r n’est pas défini pour cos(θ)=−1/2, donc θ=± 2π/3. Pour θ0=2π/3, on calcule limθ → 2π/3rsin(θ−2π/3)



    La tangente est donc l’image par la rotation de centre O et d’angle 2π/3 de la droite Y=−√3/3
  3. si la fonction n’est pas périodique, il y a lieu d’étudier l’existence de limites de r en ± ∞, si la limite est nulle on s’approche en spiralant de l’origine, si elle est finie, il y a un cercle asymptote, si elle est infinie une spirale.
  4. comme OM=r er, er=(cos(θ),sin(θ)), la vitesse (si le temps est θ) est donnée par

    v
     
    r′ er + r eθ
    où { er,eθ} est une base orthonormée directe.
    Donc si r≠ 0 ou r′ ≠ 0, le point est régulier et l’angle V de la tangente avec er vérifie
    tan(V)=
    r
    r
     ∈  ⋃ { ± ∞ } 
    (si r ≠ 0 et r′=0, la tangente est portée par eθ). Si r=0, la tangente est portée par er.10
  5. On ne peut avoir un point singulier que pour r=0. On ne fait pas leur étude comme en paramétriques, en effet la tangente est toujours portée par er, si r change de signe la courbe a la même allure que pour un point régulier, si r ne change pas de signe on a un rebroussement de première espèce (puisqu’on traverse la tangente lorsque θ augmente)
    Exemple :
  6. Convexité : pour avoir un point d’inflexion, il faut que
    1
    r
     + 


    1
    r



    ′′=0 ⇔ r2+2r2rr′′=0
    On peut le montrer de différentes manières :
    • En calculant le déterminant de {vitesse,accélération } par rapport à θ dans le repère er,eθ, on a
      v=rer+reθ   a=r′′er+2reθrer
      ⇒ det(v,a)=

      rr′′−r 
      r2r


      =2r2r r′′+r2
    • En calculant la dérivée de l’angle fait avec l’axe Ox θ+arctan(r/r′)
    • avec Xcas en se ramenant en paramétriques
      où on a noté x l’angle θ pour pouvoir dériver avec ' et X et Y les deux coordonnées.
  7. de même on calcule la courbure définie en section 11.2
    κ = 
    r2+2r2rr′′
    r2+r2
    3

10.7  Coniques

Les coniques sont des courbes implicites dont l’équation cartésienne est du second degré 

ax2+cy2+bxy+dx+ey+f=0

Exemples:


On va voir qu’elles sont de trois types : ellipses, hyperbole, parabole11 et on va les paramétriser, à partir de leur équation cartésienne ou à partir de leurs éléments géométriques (le calcul des éléments géométrique à partir de l’équation cartésienne fait intervenir l’étude des formes quadratiques, il ne sera pas abordé dans ce cours). Les coniques sont des courbes importantes en géométrie, ce qui a un intérêt en optique (parabole), mais aussi en cinématique (première loi de Kepler : l’orbite décrite par une planète est une ellipse dont le Soleil occupe un foyer).

10.7.1  Ellipse

Définition 17   L’ellipse E de foyers F1 et F2 de demi-grand axe a est l’ensemble des points M du plan tels que
MF1+MF2=2a

Exemple : ouvrir un niveau de géométrie 2d dans Xcas, choisir le mode ellipse cliquer 2 points (ce sont les foyers) puis un 3ème point (point de l’ellipse), passer en mode pointeur et faire bouger l’un des points, observer la forme de l’ellipse qui en résulte. Ou dans une ligne de commande normale taper la commande ellipse() avec en arguments les 2 points foyers et un point de l’ellipse ou l’équation cartésienne de l’ellipse, par exemple ellipse(-1,1,3+i) trace l’ellipse de foyers (−1,0), (1,0) et passant par le point (3,1).

On note 2c=F1F2 la distance entre les deux foyers, qui doit être plus petite que 2a pour que l’ellipse soit non vide. L’excentricité de l’ellipse est définie par e=c/a < 1. Si e=0, on obtient un cercle de centre F1=F2 et de rayon a. Si e≠ 0, on va voir qu’il s’agit d’un cercle contracté selon l’axe perpendiculaire à F1F2 dans un rapport de √1−e2. On va également calculer l’équation en coordonnées polaires de E (c’est sous cette forme que l’on montre que la Terre décrit une ellipse dont le Soleil occupe un foyer).

Soit O le milieu de F1 et F2, on se place dans le repère orthonormé dont le premier axe Ox contient F1 et F2 donc les coordonnées de F1 sont (c,0) et celles de F2 sont (−c,0). Soit M(x,y) un point de l’ellipse, on a d’une part :

MF12 − MF22 = (xc)2−(x+c)2 = −4cx 

et d’autre part :

MF12 − MF22 = (MF1 + MF2)(MF1 − MF2 ) = 2a (MF1 − MF2 )

donc :

MF1 − MF2 = 
−2cx
a
 

en additionnant avec MF1+MF2=2a et en appliquant c=ea, on en déduit :

  MF1 = a − 
cx
a
 = aex      (15)

En prenant le carré, on a :

(xea)2 + y2 = (aex)2

d’où :

y2 + x2 (1−e2) = a2(1−e2

finalement :

x2 + 
y2
1−e2
 = a2 

qui est bien la contraction selon Oy de rapport √1−e2 du cercle de centre O et de rayon a (appelé grand cercle de l’ellipse).

En coordonnées paramétriques, on peut utiliser le paramétrage suivant :

(x,y)=(acos(t),bsin(t))

En coordonnées polaires, on note ρ la distance de F1 à M, et θ l’angle entre l’axe Ox et F1M. L’abscisse de M est donc :

xea + ρ cos(θ)

que l’on combine avec (15) pour obtenir :

ρ = aex =a(1−e2) − e ρ cos(θ) 

donc :

ρ = 
a(1−e2)
1+ecos(θ)
 

Remarques :

Exemple : faites varier la valeur de l’excentricité ci-dessous, que voit-on pour E=0.0, E un peu inférieur à 1 (par exemple 0.8) et un peu supérieur à 1 (par exemple 1.3)

10.7.2  Parabole

Si F est un point et D une droite ne passant pas par F, la parabole de foyer F et directrice D est l’ensemble des points équidistants de F et D. En choisissant un repère tel que la droite D ait pour équation y=0 et en prenant F(0,1), M(x,y) appartient à la parabole si

|y|=d(M,D)=d(M,F)=
(y−1)2+x2
 

donc en passant au carré :

y2=(y−1)2+x2 ⇒ y=
x2+1
2

La parabole est donc (ici) un graphe de fonction, donc un cas particulier de courbe paramétrique. On peut trouver son équation en polaire, en prenant F comme origine et la directrice verticale (donc l’équation de la droite devient par exemple y=−1) sous la forme

ρ=
1
1+sin(θ)

cf. l’exercice sur les coniques données par foyer et directrice, qui traite aussi le cas des hyperboles. On peut aussi faire à titre d’exercice l’étude de la courbe en polaire :

ρ = 
A
1+ecos(θ)

lorsque e=1 et e>1.

Un intérêt majeur de la parabole en optique est que les rayons incidents perpendiculaires à la directrice se réfléchissent en passant par le foyer (on peut même montrer que cela caractérise une parabole). Illustration-démonstration avec Xcas dans un niveau de géométrie taper les commandes

P:=plotfunc(x^2/2+1/2,x=-5..5);
supposons(a=[-1.4,-5,5,0.1]);
D:=line(x=a,color=red);
M:=inter_unique(P,D);
T:=tangent(P,M);
R:=symetrie(T,D,color=red);
trace(R);

puis faire varier a en cliquant sur les flèches. Pour tester en ligne, commencez par initialiser la trace en exécutant
puis faites varier a en cliquant sur le bouton + ou - :

=
Not evaled
Noter la valeur
inter_unique(R,line(x=0))
elle est indépendante de a et est le foyer. On peut montrer qu’une courbe ayant cette propriété est une parabole.

10.7.3  Hyperbole

Une hyperbole de foyers F et F′ est définie comme l’ensemble des points M tels que :

|MFMF′|=2a

a est une constante telle que 2a>2c=FF′, avec une excentricité e=c/a>1.

En physique, les hyperboles interviennent dans les trajectoires non périodiques en mécanique céleste, mais aussi comme courbes de déphasage constant entre deux sources situées aux deux foyers (les figures d’interférence font apparaitre des hyperboles).

On peut faire un calcul analogue à celui de l’ellipse,

MFMF′=± 2a,  MF+MF′=
MF2MF2
MFMF
=−± 2ex

on en déduit que

MF=± (aex)

l’équation cartésienne de l’hyperbole dans le repère centré au milieu des foyers, d’axe Ox l’axe des foyers est donc :

x2
a2
y2
a2(e2−1)
=1

On peut paramétrer les deux branches de l’hyperbole par

x(t)=± acosh(t), y(t)=a
e2−1
 sinh(t)

et en polaires

ρ=
a(1−e2)
1+ecos(θ)

Exercice : faire l’étude de la courbe paramétrée et montrer que l’hyperbole admet deux asymptotes d’équation y = ± b/a x.

10.7.4  Paramétrisation rationnelle

Si on connait un point d’une conique, on peut effectuer un changement d’origine en ce point, l’équation cartésienne devient

P(x,y)=ax2+bxy+cy2+dx+ey=0

On suppose que (d,e)≠(0,0)13. On cherche alors l’intersection de la conique avec la droite y=tx (de pente t), on va voir que la droite coupe en général la conique en deux points, l’origine et un autre point dont on calcule les coordonnées en fonction de t14. Graphiquement, par exemple


=
Not evaled
puis faire varier la valeur de t ou d’un des coefficients de l’équation. En effet on obtient une équation du second degré en x, qui se factorise par x, l’autre solution donne alors x comme fraction rationnelle en t, puis y=tx.

(ax+btx+ct2x+d+et)x=0 ⇒ x=0, x=
det
ct2+bt+a

Comme dans le premier exemple sur le cercle trigonométrique, on n’obtient pas toujours toute la conique (s’il existe un autre point d’abscisse x=0).

Si on cherche les points où le dénominateur en t s’annule, on doit calculer (pour c≠ 0 et en supposant que la fraction −det/ct2+bt+a est irréductible15) le discriminant16 de l’équation du second degré

Δ= b2−4ac

Il y a trois cas possibles:

Exercice : paramétrer et faire l’étude des coniques :
x2+4y2+2xy=4, x2−3y2+2xy=4

Remarque : on a vu que les ellipses, paraboles, hyperboles admettent une équation réduite du second degré. On en déduit facilement que leur équation dans un repère quelconque est toujours du second degré. Réciproquement, pour une équation cartésienne on a calculé une paramétrisation rationnelle, mais pas démontré que c’était forcément une conique. Pour faire cela, l’outil adapté est l’étude des formes quadratiques. On peut toutefois le faire à la main en dimension 2, en faisant une rotation x,yX,Y pour annuler le coefficient de XY. Par exemple


on voit que l’angle de la rotation à effectuer vérifie

(ca)sin(2d)+bcos(2d)=0    ⇒    tan(2d)=
b
ac

11  Propriétés métriques des courbes.

11.1  Longueur d’arc

La longueur ds d’un morceau de courbe régulier parcouru pendant un petit intervalle de temps dt est égal au premier ordre à la longueur du segment tangent parcouru, ou encore au produit de la norme de la vitesse instantanée par dt

ds=
x2+y2
 dt

On remarque que cette quantité est invariante par changement de paramétrage, si t=t(τ) alors

ds=
dx
dt
2+
dy
dt
2
 dt 
 =
 


dx
dτ
2+
dy
dτ
2





dτ
dt



2



 
 |
dt
dτ
dτ 
 =
 
dx
dτ
2+
dy
dτ
2
 dτ

On en déduit

Proposition 18   La longueur d’un arc de courbe entre les points de paramètre t0 et t1 vaut
t1


t0
 
x2+y2
 dt

En coordonnées polaires :
θ1


θ0
 
r2+r2
 dθ

Remarque : il est très rare que l’on puisse effectuer le calcul explicite d’une primitive de √x2+y2, il faut alors se contenter d’une valeur approchée de l’intégrale lorsque t0 et t1 ont des valeurs numériques, calculée par des méthodes numériques qui généralisent la méthode des rectangles (cf. le cours de mat249). Ce calcul se fait avec Xcas (ou une calculatrice formelle) en donnant une valeur approchée à l’une des bornes. Il y a quelques exceptions  par exemple la longueur d’un arc de parabole se calcule avec une formule explicite (essayez la commande int(sqrt(1+4t^2),t,t0,t1) ou



La cycloïde17

x(t)=R(t−sin(t)), y(t)=R(1−cos(t))


admet aussi une formule simple pour sa longueur

Par contre, la longueur d’un arc d’ellipse ne se calcule pas avec les fonctions usuelles (pour pouvoir le faire, il faut introduire des fonctions spéciales adaptées, appelées intégrales elliptiques) :

11.2  Courbure, repère de Frenet, accélération normale et tangentielle.

Si on choisit s, la longueur d’arc, comme nouveau paramètre de temps, la longueur parcourue est égale au temps, donc la vitesse instantannée par rapport à s est de norme 1. On peut aussi le voir en notant M(t)=(x,y) :

dM
dt
 =
dM
ds
 
ds
dt
   ⇒  || 
dM
dt
 || = || 
dM
ds
 || |
ds
dt
|  ⇒  v = || 
dM
ds
 || v

v est la norme de la vitesse avec t comme paramètre, donc || dM/ds || est bien égal à 1.

Calculons maintenant l’accélération avec ce nouveau paramètre s. Comme la vitesse est de norme constante égale à 1, donc de carré 1, en dérivant (dM/ds)2 par rapport à s, on vérifie que l’accélération est perpendiculaire à la vitesse pour ce paramétrage par la longueur d’arc s. L’accélération par rapport à s est donc portée par la normale à la trajectoire, et sa mesure algébrique est appelé courbure (signée), notée κ, la valeur absolue de l’inverse de κ est appelé le rayon de courbure (la direction de l’accélération pointe vers le centre de courbure).

d2M
ds2
 ⊥ 
v
 
,    || 
d2M
ds2
|| = |κ| = 
1
R
 

Si on se déplace sur un cercle de centre O et de rayon R à vitesse 1, alors x(t)+iy(t)=Reit/R, la vitesse est donnée par x′+iy′=ieit/R donc de norme 1, et l’accélération par x″+iy″=−1/R eit/R, sa norme vaut 1/R et sa direction pointe vers le centre du cercle. Donc la courbe est, à l’ordre 2 au point considéré, identique à un cercle de rayon R.

Revenons au paramètrage initial t. Dérivons par rapport à t la vitesse dM/dt = v dM/ds, on obtient :

 
a
 
 =
d2M
dt2
=
dv
dt
 
dM
ds
 + v 
d
dt
 


dM
ds
 


 =
dv
dt
 
dM
ds
 + v 
ds
dt
 
d2M
ds2
 =
dv
dt
 
dM
ds
 + v2 
d2M
ds2

L’accélération se décompose donc en deux parties 

Autre formule de calcul du rayon de courbure : l’accélération normale an vaut v2/R donc

|| 
a
 
 ∧ 
v
 
 ||= an ||
v
 
||=
v3
R
   ⇒  R =v3/||
a
 
 ∧ 
v
 
|| =
x2+y2
3
|xy′′−yx′′|
Proposition 19   On appelle repère de Frenet en un point M régulier d’une courbe, le repère orthonormé direct formé par le point de la courbe, le vecteur tangent T et le vecteur normal N. On a alors
v
 
=v 
T
 
=
ds
dt
T
 
,    
d
ds
T
 
=κ 
N
 
d
ds
N
 
=−κ 
T
 
,    R
1
κ
(l’avant-dernière formule vient du fait que { T ,N } est une base orthonormée directe, le signe ± est déterminé par la convexité de la courbe), et :
a
 
=
d
dt
v
 
=
dv
dt
 
T
 
 ±
v2
R
 
N
 
,    R=
x2+y2
3
|xy′′−yx′′|
On appelle centre de courbure le point Ω=M+1/κN. Le cercle de centre Ω passant par M (de rayon R) est appelé cercle osculateur en M à la courbe.

Exemple : calcul du cercle osculateur en un point d’une parabole (t,t2).

x′=1, y′=2t,  
T
 
=(
1
1+4t2
,
2t
1+4t2
),    y′′=2   R=
1+4t2
3
2

=
Not evaled
Avec Xcas version 1.1.1-18 ou supérieure, on peut taper directement :
C:=cercle_osculateur(G,M)

Remarques :

12  Représentation des courbes implicites.

Certaines représentations graphiques nécessitent peu d’outillage mathématique, ainsi les fonctions, les courbes paramétrique et polaires peuvent être représentées en échantillonant une ou plusieurs expressions selon une discrétisation donnée explicitement par l’utilisateur ou par des paramètres par défaut, les points obtenus étant ensuite reliés par des segments. On pourrait bien sur automatiser avec le calcul formel l’étude de la courbe (tableaux de variations, asymptotes, points singuliers, etc.).

Par contre les courbes données par une équation implicite font intervenir des algorithmes et des mathématiques plus intéressantes. En dimension 2, on se donne donc une équation f(x,y)=0 et on suppose f suffisamment régulière. Supposons la courbe non vide, soit (x0,y0) un point de cette courbe, si (∂x f,∂y f)(x0,y0) ≠ 0 on peut appliquer le théorème des fonctions implicites et la courbe est localement comme une courbe de fonction (en x ou en y). On en calcule la tangente et on peut suivre cette tangente un pas de discrétisation puis utiliser une méthode numérique de recherche de solution près de la tangente. Ces points sont appelés points réguliers

Les points où (∂x f,∂y f)=0 sont les points singuliers. Génériquement, il n’y en a pas puisque cela donne 3 équations à 2 inconnues, par contre si on s’intéresse à une famille de courbes dépendant d’un paramètre, il en apparait. En ces points, on calcule le développement de Taylor et on recherche le premier terme homogène non nul (homogène après translation bien sur), par exemple P2=x2y2 pour x3+x2y2 en (0,0). Supposons que le polynôme correspondant Pm est sans racines multiples, et (quitte à faire une rotation) que le coefficient de ym est non nul. Pm est un polynôme homogène donc se factorise au moins numériquement (en remplaçant une des variables par 1, on est ramené en dimension 1), et on montre qu’il y a m arcs de courbe complexes tangents aux droites d’équations ces m facteurs (et au plus m arcs de courbe réels si on ne garde que les racines réelles). En effet, on pose y=xY et on est amené à résoudre

f(x,xY)=0=xmPm(1,Y) + xm+1 g(x,Y)

g est un polynôme si f est un polynôme (plus généralement a la même régularité que f). Après simplification par xm, on peut appliquer le théorème des fonctions implicites pour déterminer Y en fonction de x au voisinage de x=0 et de chacune des racines de Pm(1,Y) en Y (puisque les racines sont simples). Le point est dit singulier-régulier ou singulier ordinaire. C’est ce que fait la commande implicitplot de Xcas (affichage des informations intermédiaires).

Si le point singulier n’est pas ordinaire, l’équation devient

(Yt)k 
 
i
 (Yti) + xg(x,Y)=0,    k>1

et il faut faire intervenir des puissances fractionnaires en x (dépendant de termes supérieurs du développement de Taylor de f en (0,0)) pour désingulariser les k arcs de courbes ayant même tangente y=tx en (0,0). Par exemple si g(0,t) ≠ 0, on pose X=x1/k, Y=t+XZ qui donne

Zk 
 
i
 (tti+XZ) + g(Xk,t+XZ)=0

pour X=0 on a alors k solutions non nulles Z qui se prolongent au voisinage de X=0 par le théorème des fonctions implicites.

Certains cas particuliers peuvent être traités en transformant la courbe implicite en courbe paramétrique, c’est le cas des courbes algébriques de degré 2, qui sont des coniques. On peut les paramétrer rationnellement si on en connait un point (en prenant la droite passant par ce point de pente m et en cherchant l’autre point d’intersection avec la conique (il y en a forcément un et un seul autre, parce que l’équation correspondant aux points d’intersection est de degré 2 et on connait déjà une solution), cette paramétrisation est intéressante pour faire du calcul formel, mais moins pour des représentations graphiques, on lui préferera une paramétrisation trigonométrique pour une conique ou exponentielle pour une hypebole, par exemple (cos(t),sin(t)) plutot que 1+it/1−it pour le cercle unité, paramétrisation obtenue en calculant les éléments propres de la conique (conique_reduite). Pour les courbes algébriques de degré plus grand, on commence par factoriser le polynôme, c’est une factorisation absolue (section 18.7) qui est nécessaire (ou au moins numérique dans [x,y]). Pour le moment, Xcas fait simplement une factorisation sur le corps des coefficients, et repère les équations de coniques.

13  Formes différentielles et intégrales curvilignes

Il s’agit dans cette section de calculer des intégrales le long de l’arc. Cela intervient par exemple pour calculer le travail d’une force au cours d’un déplacement le long d’une courbe ou la quantité de chaleur/travail pendant un cycle en thermodynamique (le long d’une courbe dans le plan défini par deux coordonnées indépendantes comme par exemple pression-température ou pression-volume). Dans les cas favorables, on a un analogue des primitives, on peut calculer un potentiel et faire la différence de potentiel entre les deux extrémités du chemin pour calculer l’intégrale curviligne. On va d’abord définir ce qu’on peut intégrer le long d’une courbe, à savoir une forme différentielle (aussi appelée 1-forme), puis on donnera quelques résultats sur les formes fermées et exactes (c’est le cas favorable, il correspond aux forces conservatives en mécanique ou aux différentielles totales de fonctions d’état en thermodynamique).

13.1  Forme différentielle

Soit V(x,y) une fonction de deux variables continument dérivable. On s’intéresse aux variations de V lorsqu’on se déplace dans le plan depuis le point M(x,y) dans une direction donnée à la vitesse w. On a alors une formule équivalente à celle de la dérivée d’une fonction d’une variable :

Proposition 20   Pour tout vecteur w=(w1,w2), la dérivée de V en (x,y) dans la direction w est donnée par :
 
lim
h→ 0
 
V((x,y)+wh)−V(x,y)
h
= ∂xVw1+∂yV w2
On appelle différentielle de V et on note dV l’application qui en un point (x,y) associe au vecteur w la valeur de la dérivée directionnelle de V en (x,y) selon w
dV(w)=∂xV w1+∂yV w2
Cette application est linéaire par rapport à w.

En effet :

 V(x+w1h,y+w2h)=V(x+w1h,y)+∂yV(x+w1h,yw2h + o(h)
 =V(x,y)+∂xV(x,yw1 h + ∂yV(x+w1h,yw2h + o(h)

donc

 
V(x+w1h,y+w2h)−V(x,y)
h
=xV(x,yw1 + ∂yV(x+w1h,yw2 +o(1) 
 h→ 0xV(x,yw1 + ∂yV(x,yw2 

Exemples :

Remarque : Différentielle et gradient
La différentielle dV a les mêmes composantes que le gradient de V (gradient(V,[x,y]) avec Xcas), mais ce ne sont pas les mêmes objets : en un point donné dV est une application linéaire (qui a un sens indépendamment de la définition d’un produit scalaire) alors que ∇ V est un vecteur (dont la relation avec la dérivée directionnelle dépend du produit scalaire), on a pour tout vecteur w la relation

dV(w)=∇ Vw 

On a la même relation entre le travail d’une force (qui est une forme linéaire qui s’applique sur les vecteurs déplacement) et la force correspondante (qui est un vecteur défini à l’aide du produit scalaire). On parle parfois de vecteur covariant pour la différentielle (et vecteur contravariant pour le gradient).

Applications :

On note donc dx [resp. dy] la différentielle de V(x,y)=x [resp. V(x,y)=y]21 on a :

dV=∂xV dx + ∂yV dy

Une forme différentielle ω est la généralisation de la différentielle d’une fonction, elle s’écrit sous la forme

ω=M(x,ydx + N(x,ydy

M et N sont des fonctions des deux variables x,y, mais pas forcément les dérivées partielles d’une fonction V.

La définition géométrique d’une forme différentielle ω est la donnée en tout point du plan (ou d’un domaine ouvert du plan) d’une application linéaire de 2 à valeur dans 22 (ou en tout point de l’espace d’une application linéraire de 3 à valeurs dans pour une courbe de 3). Si on prend la base canonique de 2, une application linéaire de 2 dans est caractérisée par sa matrice qui possède une ligne et deux colonnes et a donc deux coefficients M et N, une forme différentielle équivaut donc bien à la donnée d’un couple de fonction M(x,y),N(x,y).

13.2  Intégrale curviligne

Ceci permet de donner la :

Définition 21   Pour calculer l’intégrale curviligne d’une forme différentielle le long d’un arc de courbe orienté, on choisit un paramétrage de l’arc continument dérivable par morceaux (on suppose qu’il en existe un), et on calcule l’intégrale usuelle par rapport au paramètre de la forme différentielle appliquée au vecteur tangent entre les deux valeurs du paramètre correspondant à l’origine et extrémité de l’arc de courbe :
 


γ
ω = 
t1


t0
ω


dγ(t)
dt



dt 
En coordonnées,
 


γ
ω =
t1


t0
  (M(x(t),y(t)) 
dx
dt
 + N(x(t),y(t
dy
dt
)  dt

Exemple: on prend ω=ydx et on calcule l’intégrale curviligne le long de l’arc de parabole (t,t2) pour t∈[0,1], on obtient

1


0
 t2  dt =
1
3

En paramétrant par (u2,u4) avec u∈[0,1]

1


0
 u4 (2u  du) = 


2
u6
6



1



0
=
1
3

on retrouve le même résultat.

La valeur de l’intégrale est bien définie indépendamment du paramétrage, en effet si on change de paramétrage avec une bijection tu(t) envoyant [t0,t1] sur [u0,u1], on a (en utilisant la linéarité de ω à la deuxième ligne) :

 
u1


u0
 ω


dγ(u)
du



du
=
t1


t0
 ω


dt
du
 
dγ(t)
dt
 


du
dt
 dt 
 =
t1


t0
  
dt
du
 ω


dγ(t)
dt
 


du
dt
 dt 
 =
t1


t0
  ω


dγ(t)
dt
  


dt

Attention à l’orientation, si on change d’orientation, on change le signe, par exemple si on parcourt l’arc de parabole de (1,1) vers (0,0), en utilisant le paramétrage (1−t,(1−t)2), t ∈ [0,1], on obtient l’opposé :

1


0
 (1−t) (−dt) = 


(t−1)2
3



1



0
 = −
1
3

Remarque : le travail d’une force F=(Fx,Fy) le long d’un arc de courbe est donné par l’intégrale curviligne de la forme différentielle Fx dx+Fydy.

13.3  Forme différentielle exacte

Voyons maintenant à quelle condition il existe un analogue du calcul avec une primitive. On a:

 


γ
dV=V(γ(t1))−V(γ(t0)), 

En effet, si on est sur un morceau d’arc où on peut paramétrer par x alors

 


γ
dV=
x1


x0
 ∂x V  dx = V(γ(t1))−V(γ(t0))

De même si on peut paramétrer par y. On recolle alors les morceaux d’arcs (on peut paramétrer par x ou par y en tout point régulier de γ).

Pour une force qui dérive d’un potentiel, on a donc montré que le travail de la force se calcule en faisant la différence de potentiel entre les deux extrémités. Cette propriété, analogue au calcul d’intégrale classique en utilisant une primitive n’est pas automatique, car elle implique que l’intégrale curviligne ne dépend pas du chemin choisi pour relier les deux points. Or en thermodynamique, la chaleur est modélisée par une forme différentielle, mais la chaleur échangée dépend du chemin suivi (c’est vrai aussi en mécanique pour le travail de forces non conservatives comme les forces de frottement). En mathématiques, on parle de forme différentielle exacte ou non exacte.

Définition 22   Une forme différentielle ω est exacte s’il existe une fonction V telle que sur tout arc de courbe γ d’origine A et extrémité B
 


γ
ω = V(B)−V(A)
Attention la convention de signe est opposée à celle utilisée pour le potentiel d’une force en physique.

Si on choisit comme chemin un segment entre deux points A et B d’ordonnées identiques y et d’abscisses x et x+h, alors

x+h


x
 M dx+Ndy = V(x+h,y)−V(x,y

en faisant tendre h vers 0, on a

M=
 
lim
h→ 0
 
V(x+h,y)−V(x)
h
 = ∂x V

De même N=∂y V. Réciproquement, si M=∂x V et N=∂y V alors ω=dV donc ∫γω=V(B)−V(A)

Proposition 23   Une forme différentielle ω est exacte si et seulement si il existe une fonction V telle que :
ω=∂x V dx + ∂y V dy=dV

Si V est deux fois continument différentiable alors ∂yx V = ∂xy V. D’où une condition nécessaire pour que ω soit exacte :

y M = ∂yx V = ∂xy V = ∂x N
Définition 24   On appelle forme différentielle fermée une forme différentielle ω=Mdx+Ndy telle que y M=∂x N

Une forme exacte est toujours fermée, mais la réciproque n’est pas toujours vraie, une forme fermée n’est pas forcément exacte, cela dépend où elle est définie. Si elle est définie dans un domaine ouvert de 2 sans trou (2 tout entier, un rectangle, un disque, etc.), on peut montrer qu’une forme fermée est une forme exacte, en appliquant le théorème de Stokes (voir section suivante). Sinon, il existe des contre-exemples, comme sur le cercle unité

ω=
ydxxdy
x2+y2

La forme est fermée :
simplify(diff(y/(x^2+y^2),y)-diff(-x/(x^2+y^2),x))
mais elle n’est pas exacte :
x:=cos(t); y:=sin(t);
int((y*diff(x,t)-x*diff(y,t))/(x^2+y^2),t,0,2*pi)

Pour trouver le potentiel V dont une forme différentielle fermée ω=M dx+Ndy est la différentielle, on résoud d’abord par exemple M = ∂x V en intégrant M par rapport à x, y étant considéré comme un paramètre, on obtient V à une constante d’intégration près, cette constante d’intégration en x peut dépendre de y, c’est donc une fonction C(y), on remplace dans N=∂y V et on intègre en y pour trouver la valeur de C(y) (à une constante près). Cette opération est executée par la commande potential() de Xcas.

Si une forme n’est pas fermée, elle n’est pas exacte, et on ne peut pas calculer une intégrale curviligne par différence de potentiel. Il peut arriver qu’en multipliant la forme par une fonction, on trouve une nouvelle forme qui elle est fermée, on parle alors de facteur intégrant. Par exemple en thermodynamique, la forme chaleur n’est pas fermée, mais en divisant par la température on obtient une forme fermée dont le potentiel est l’entropie. Cela peut aussi servir à trouver des constantes du mouvement pour certaines équations différentielles. En effet, si on se déplace le long d’une courbe de niveau du potentiel d’une forme exacte, alors le long de cette courbe le potentiel est constant, donc la forme appliquée au vecteur tangent est nulle, on dit que la courbe de niveau est une courbe intégrale de la forme différentielle (exacte).

13.4  Intégrale curviligne et intégrales doubles.

Terminons ce chapitre par le lien entre intégrale curviligne sur un lacet (chemin fermé) et intégrale double à l’intérieur du lacet. C’est évidemment surtout intéressant pour les formes non exactes, car si γ est un lacet et ω une forme exacte, alors ∫γω=0. On a le théorème de Stokes, aussi appelé en dimension 2 formule de Green-Riemann :

Théorème 25   Si U est un domaine de frontière orientée γ continument dérivable (γ est donc un chemin fermé parcouru une fois que l’on oriente dans le sens trigonométrique), et si ω=Mdx + N dy est une forme différentielle continument dérivable alors :
 


γ
ω = 
 


U
 dω := 
 


U
 (∂x N −∂y M)  dx dy

Idée de la preuve : on commence par le cas où U est un rectangle [a,b] × [α,β], on peut alors calculer

 


U
 ∂x N   dx dy  = 
β


α
(
b


a
 ∂x N   dxdy
β


α
(N(b,y)−N(a,y)) dy 

on compare avec les intégrales curvilignes sur les segments verticaux {(a,y), y ∈ [α,β]} et {(b,y), y ∈ [β,α]}. De même pour M et les segments horizontaux.

Pour le cas d’un domaine d’intégration U plus général, on approche U par une réunion disjointe de petits rectangles.

Application : pour calculer l’aire d’un domaine U de frontière γ, il suffit de calculer l’une des intégrales curvilignes :

 


γ
y dx
 


γ
x dy = 
 


γ
y dx − x dy
2

Par exemple, l’aire à l’intérieur de l’ellipse x=acos(t), y=bsin(t) vaut



0
 
b sin(td(acos(t)) − acos(t) d(bsin(t))
2
 = abπ 

On peut aussi calculer des moments d’inertie ou la position d’un centre de gravité en se ramenant à une intégrale curviligne.
Exemple : Calculer la position du centre d’inertie d’un quart de cercle C={(cos(t),sin(t)), t ∈ [0,π/2]}.
On a donc U délimité par γ, réunion de {(x,0), x ∈ [0,1]} , C et {(0,y), y ∈ [1,0]}. Pour trouver la position du centre d’inertie en x (en y c’est identique), on doit calculer

 


U
 x  dx dy = 
 


γ
1
2
 x2  dy = 0 + 
1
2
 
π
2


0
 cos(t)2 cos(t)  dt + 0= 
1
3

et on divise par π/4 l’aire du quart de cercle, on trouve donc (4/3π,4/3π), on peut visualiser avec la commande cercle(0,1); G:=point(4/(3*pi),4/(3*pi))

14  Équations et systèmes différentiels.

14.1  Introduction et représentation graphique.

On s’intéresse à l’équation différentielle

  y′=
dy
dt
=f(y,t)     (16)

y(t) ∈ n et f: n × → n. Si n=1, c’est une équation différentielle, si n>1 c’est un système différentiel.

Exemple : en dimension n=1, y′=f(y,t)=ay. On sait résoudre cette équation, les solutions sont de la forme y(t)=Ceat. Si on trace la courbe représentative de ces solutions (appelée courbe intégrale), on observe que par tout point du plan, il passe une solution unique. La tangente à une courbe intégrale a pour pente y′=ay donc pour vecteur directeur le vecteur de composantes (1,ay).

C’est vrai de manière plus générale, le vecteur directeur de la tangente à une courbe intégrale est (1,f(y,t)). Si on représente dans le plan selon un quadrillage régulier les vecteurs (1,f(y,t)), une courbe intégrale doit être tangente à ces vecteurs chaque fois qu’elle passe en un point du quadrillage, (et à peu près tangente si elle passe à proximité). Un tel quadrillage est appelé champ des tangentes (commande plotfield en Xcas, mode également disponible sur certaines calculatrices).

Exercice : tracer le champ des tangentes et quelques solutions pour quelques exemples de fonction f(y,t), avec Xcas créer une figure 2d, puis choisir le mode Champ des tangentes du menu Geo, Graphe, entrer la fonction, puis cliquer en quelques points pour faire tracer la solution passant par ces points.

Par exemple pour y′=−y+cos(t)

=
Not evaled

L’équation (16) est d’ordre 1, or certaines équations différentielles se présentent naturellement comme des équations d’ordre 2, par exemple l’équation fondementale de la dynamique (accélération=somme des forces divisée par la masse). Mais on peut facilement se ramener à un système différentiel d’ordre 1, en augmentant la dimension de y. Par exemple, si on pose y=(x(t),v(t)), où x(t) est la position et v(t) la vitesse, alors l’équation devient un système d’ordre 1

d
dt
 

x(t
v(t







v(t
F
m
 
 




F est la force, qui dépend de la position x(t) (champ électrique, gravitation...) et éventuellement de la vitesse (force de frottement, champ magnétique...). On utilise aussi assez fréquemment y=(q(t),p(t)) où q(t) est la position, et p(t) la quantité de mouvement (qui dépend de la vitesse, linéairement en mécanique classique).

Représentation graphique : comme précédemment, on peut se placer dans l’espace des (t,x,v) (si x est en dimension 1), mais il est souvent plus difficile d’observer des phénomènes sur un graphe en 3-d que dans le plan, on préfère ne pas représenter explicitement le temps t, mais uniquement (x,v), on est donc naturellement ramené à représenter une solution (une courbe intégrale) par une courbe paramétrique en (x,v) (ou en position impulsion). On a encore la notion de champ des tangentes si f(y,t)=f(y) ne dépend pas explicitement du temps (on dit que le système est autonome), dans ce cas une courbe intégrale a pour tangente en y2 de direction portée par le vecteur f(y) ∈ 2.
Exemple : (x,v)′=5(−v,x). La commande

permet d’en représenter le champ des tangentes et d’avoir une idée approximative de l’allure des solutions (normalize affiche des vecteur tangents de norme 1, si on n’utilise pas cette option, la taille des vecteurs tangents donne la “vitesse” de déplacement). On sait résoudre ce système différentiel, soit en appliquant une technique matricielle présentée ci-dessous, soit en se ramenant à une équation linéaire d’ordre 2 à coefficients constants:

x′′=−5v′=−25x

donc x(t)=Acos(5t)+Bsin(5t), A, B étant déterminés par les conditions initiales sur (x,v).

Une équation donnée sous la forme (16) est appelée une équation résolue en y, car on a exprimé la dérivée en fonction de y et de t. Il existe (plus fréquemment en mathématiques) d’autres formes d’équations différentielles (non résolues) où le premier travail de résolution peut consister à exprimer y′ en fonction de y et t (ce qui n’est pas toujours possible explicitement).

Exemple : en dimension 1, ty′=y, on sait résoudre exactement cette équation à variables séparables, les solutions sont de la forme Ct.

On observe que contrairement à y′=ay où passe une solution et une seule par chaque point du plan, ici toutes les solutions valent 0 en t=0 : il passe une infinité de solutions par le point (0,0) et il n’en passe aucune par (0,a), a ≠ 0. Ce phénomène de non unicité/non existence vient de la mise sous forme résolue y′=y/t qui fait apparaitre une singularité de f(y,t) en t=0.

On présente dans la suite de cette section des résultats qualitatifs sur les équations sous forme résolue lorsqu’on ne sait pas les résoudre, ainsi que quelques méthodes explicites pour certaines équations différentielles que l’on sait résoudre.

14.2  Existence et unicité

Il s’agit ici de préciser dans quelles conditions le résultat intuitif suivant est vrai : étant donné une condition initiale y(t0)=y0, il y a une et une seule évolution possible, donc une solution unique y(t) de l’équation ou du système (16).

On a le :

Théorème 26   (Cauchy-Lipschitz) Si f est continument dérivable en y et t sur n × ou sur un domaine ouvert D inclus dans n × , alors l’équation (ou le système) résolu (16) admet pour toute condition initiale y(t0)=y0 une solution unique sur un intervalle maximal ouvert en temps contenant t0.

Remarques

On admettra ce théorème, voici quelques idées heuristiques de la preuve. L’équation y′=f(y,t) peut se réécrire sous la forme intégrale équivalente

y(t)=y(t0)+
t


t0
 y′(u)  du = y(t0)+
t


t0
 f(y(u),udu 

Si t est assez proche de t0, on peut approcher l’intégrale par

y(t) = y(t0) + (tt0f(y(t0),t0) + petite erreur

C’est exactement ce qu’on fait en suivant le champ des tangentes pour approcher une courbe intégrale graphiquement, et si on discrétise le temps avec un pas petit, cette méthode d’approximation est appelée méthode d’Euler. On peut bien sur utiliser d’autres approximations (meilleures) de l’intégrale pour avoir une meilleure approximation de la solution, et les méthodes dites de Runge-Kutta utilisent cette idée. D’un point de vue théorique, la preuve repose plutôt sur ce qu’on appelle le théorème du point fixe, on met la valeur approchée de y(t) trouvée dans l’équation intégrale pour avoir une nouvelle valeur approchée de y(t), on recommence, ainsi de suite, et on montre que le processus converge (il s’agit mathématiquement parlant d’une suite récurrente de fonctions, la preuve rigoureuse de la convergence nécessite des outils mathématiques de niveau L3-M1 de maths, c’est l’analogue des suites récurrentes de réels qui permettent de résoudre numériquement des équations comme x=cos(x) abordées en mat249).

Conséquence du théorème 26 : deux courbes intégrales de la même équation différentielle ne peuvent se couper dans D. Donc si on connait une courbe intégrale C de D et qu’on prend une condition initiale en-dehors de cette courbe, la courbe intégrale unique passant par cette condition initiale restera du même coté de D. Si on connait deux courbes intégrales de D, une courbe intégrale passant par une condition initiale entre les deux courbes restera entre les deux courbes.

Exemple : y′=y(1−y) (équation logistique).

=
Not evaled
Cette équation autonome admet deux solutions évidentes y=0 et y=1. Donc pour toute condition initiale y(t0) ∈ ]0,1[, on a y(t) ∈ ]0,1[23. On en déduit que y′=y(1−y)>0 donc la solution y est strictement croissante, comme elle est bornée par 0 et 1, elle admet une limite pour t → ± ∞, donc y′ tend vers 0 pour t → ± ∞, donc y tend vers 0 ou 1, et comme y croit, y → 0 en t=−∞ et y → 1 en t=+∞. Le comportement à l’infini est donc indépendant de la valeur précise de la condition initiale, pourvu qu’elle soit dans ]0,1[.

Exercice : toujours pour y′=y(1−y) que se passe-t-il pour une condition initiale y(t0)>1 ?

14.3  Quelques méthodes de résolution explicite.

14.3.1  Équations à variables séparables

Si on peut factoriser f(y,t) en termes ne dépendant que de y ou ne dépendant que de t, on dit que l’équation est à variable séparable

y′=f(y,t)=g(t)h(y)

Cette équation admet des solutions constantes y=y0 lorsque h(y0)=0. Si h(y(t0)) ≠ 0, par le théorème de Cauchy-Lipschitz h(y(t)) ne s’annule nulle part sur son domaine de définition. On peut donc diviser par h(y) et intégrer :

⇒ 
dy
h(y)
 = g(tdt

On obtient une équation implicite de la forme H(y)=G(t)+CG est une primitive de g, H de 1/h et C une constante arbitraire. Dans les cas favorables, on peut exprimer y en fonction de t (par exemple si l’équation est linéaire sans second membre, on a h(y)=y donc H est le log que l’on sait inverser). Dans les cas moins favorables, on peut exprimer y et t en fonction d’un paramètre u : la courbe intégrale est une courbe paramétrée. Dans les cas défavorables, on reste sous forme implicite.

Exercice : résoudre explicitement l’équation y′=y(1−y) et retrouver les résultats qualitatifs de la section précédente.

14.3.2  Équations linéaires

On commence par résoudre l’équation sans second membre (aussi appelée homogène)

an(ty[n] +...+a1(t)y′+a0(t)y=0

sur un intervalle ouvert sur lequel an(t) ≠ 0. L’ensemble des solutions est un espace vectoriel (car l’équation est linéaire) et de dimension l’ordre de l’équation : pour le prouver on peut appliquer le théorème de Cauchy-Lipschitz au système d’ordre 1 équivalent, ce système est tel que y est un vecteur de n, on a ensuite un isomorphisme entre les solutions et la condition initiale.

Si l’ordre est 1, on a une équation à variables séparables y′/y=−a0(t)/a1(t) et la solution est une exponentielle :

y(t)=Ce
a0
a1
  dt
 

Exemple : y′−ty=0, on a y(t)=Cet dt=Cet2/2

Si l’ordre est plus grand que 1, on n’a en général pas de solution explicitable avec les fonctions usuelles et des primitives, pour certaines équations importantes en physique, des fonctions spéciales ont été créées pour exprimer les solutions, par exemple les fonctions de Bessel. Il existe quelques cas particuliers où le calcul explicite est possible, dont le cas où les coefficients sont constants (section suivante). Si on connait une solution w d’une équation linéaire, alors en posant y=wz, la fonction z′ vérifie une équation linéaire d’ordre un de moins, ainsi si on connait une solution d’une équation linéaire d’ordre 2, on peut la résoudre complètement.

Le calcul d’une solution particulière d’une équation linéaire avec second membre se fait en faisant varier les constantes d’intégration : on prend la forme générale de la solution de l’équation homogène, on remplace les constantes d’intégration par des fonctions inconnues, on remplace dans l’équation avec second membre et on résoud en les fonctions inconnues, la méthode détaillée dans le cas des coefficients constants s’applique à l’identique. La solution générale est la somme d’une solution particulière et de la solution générale de l’équation sans second membre.

Exemple : y′−ty=−t, solution générale de l’équation homogène y(t)=Cet2/2, variation de la constante on remplace y(t)=C(t)et2/2 dans y′−ty=−t et on obtient Cet2/2=−t, donc C′=−tet2/2 et C=et2/2+K, d’où la solution générale y(t)=(et2/2+K)et2/2=1+Ket2/2.

14.3.3  Équations linéaires à coefficients constants

On peut chercher des solutions de l’équation sans second membre sous la forme d’exponentielles ert, r doit alors vérifier une équation polynomiale P(r)=0 appelée équation caractéristique, de degré le degré de l’équation différentielle. Plus précisément, si on remplace ert dans

an y[n]+...+a1 y′+a0y=0

alors

an rn +...+a1r +a0=P(r)=0

Si P n’a que des racines simples r1,...,rn ∈ , l’ensemble des solutions est alors l’espace vectoriel engendré par { er1t, ... , ernt }. En effet, on a le bon nombre d’éléments (n), il suffit donc de montrer qu’il s’agit d’une famille libre. Cela se fait par récurrence. Au rang n=1 c’est évident. Si n>1 et si (λ1,...,λn) vérifient :

n
j=1
 λj erjt = 0

on factorise ern t et on dérive, on a

n−1
j=1
 λj (rjrne(rjrn)t =0 

on est ramené à l’identité précédente au rang n−1 donc par récurrence, λj (rjrn)=0 et λj=0 si jn, puis λn=0 avec la relation du départ.

Si P a des racines multiples, on peut montrer que pour chaque racine rk de multiplicité m>1, il faut rajouter { terkt, ..., tm−1 erkt } pour former une base de solutions. En effet

(ty)[j] = t y[j] + j y[j−1]

donc si y est solution de l’équation alors ty est encore solution si :

nan y[n−1] + (n−1)an−1 y[n−2]+...+a1=0

on reconnait l’équation différentielle linéaire à coefficients constants dont l’équation caractéristique est P′. La preuve de l’indépendance linéaire de ces fonctions est plus technique, elle peut se faire par récurrence sur le nombre de racines. S’il y en a une, on a un polynôme. Sinon, on factorise ernt, et on dérive la multiplicité de rn pour appliquer le résultat au cran n−1, on a alors un système triangulaire sur le groupe d’inconnues de la même exponentielle.

Si P est à coefficients réels et admet une racine non réelle z alors z est encore racine, on peut réécrire avec des fonctions trigonométriques les combinaisons linéaires de ezt et ezt :

(α + i β) e(a+ib)t + (α − i β)  e(aib)t =  eat ( 2 α cos(bt) − 2 β sin(bt) )

Exemples :

On peut trouver une solution particulière de l’équation avec second membre comme dans le cas général (méthode de variation des constantes). Si la solution générale est engendrée par y1,...,yn, on pose :

y=
n
i=1
 λi yi

On pose

n
i=1
 λi′ yi=0  ⇒   y′=
n
i=1
 λi yi

et ainsi de suite jusqu’à la dérivée d’ordre n de y, ces n−1 équations et léquation différentielle donnent alors un système linéaire n,n en les λi′.

Pour des second membre combinaison linéaire de termes b(t)ert avec b polynôme, on peut chercher une solution particulière combinaison linéaire de a(t)erta est de même degré que b si r n’est pas racine de P, ou de degré le degré de b plus la multiplicité de r comme racine de P. On peut aussi utiliser la transformation de Laplace et son inverse.

14.3.4  Systèmes différentiels linéaires à coefficients constants d’ordre 1.

Il s’agit donc de systèmes de la forme

y′=Ay+b(t)

y(t)∈ n, A est une matrice carrée de taille n indépendante du temps, et b(t) ∈ n.

On commence par résoudre l’équation homogène y′=Ay. Si la matrice A est diagonalisable, alors A=PDP−1D=diag(d1,...,dn) est diagonale et P inversible, le système devient :

y′=PDP−1 y

donc en posant y=Pz, on a (puisque P est indépendant du temps) :

z′=Dz    ⇔    zk′=dkzk,  k=1..n

donc zk=ck edkt, puis la solution générale

y(t)=P


c1 ed1t 
... 
cn ednt
 


Le calcul avec Xcas se fait en utilisant la commande desolve, par exemple
desolve(y'=[[1,2],[2,1]]*y)


ou avec conditions initiales
desolve(y'=[[1,2],[2,1]]*y and y(0)=[1,2])

On peut aussi utiliser la fonction exp avec comme argument At (on généralise ainsi la notation eat de la dimension 1), multiplié par la condition initiale :
exp([[1,2],[2,1]]*t)*[1,2]
Les calculs intermédiaires pour diagonaliser la matrice A sont exécutés par les commandes eigenvals, eigenvects, jordan.

On peut ensuite calculer une solution particulière par la méthode de variation des constantes, ou encore en résolvant z′=Dz+P−1b(t) composante par composante (ou par transformation de Laplace). Avec Xcas, il suffit d’ajouter le second membre dans la commande desolve
desolve(y'=[[1,2],[2,1]]*y+[x,x+1])

Si la matrice A n’est pas diagonalisable (ce qui entraine qu’elle a au moins une valeur propre de multiplicité plus grande que 1), on peut alors la trigonaliser, on se ramene à résoudre un système triangulaire, ce qui revient à résoudre pour chaque composante une équation différentielle linéaire d’ordre 1 avec un éventuel second membre.

14.3.5  Systèmes et équations

Il y a un lien entre systèmes différentiels linéaires et équations linéaires. En effet une équation d’ordre n peut s’écrire comme un système différentiel d’ordre 1, on peut calculer le polynôme caractéristique de la matrice on retrouve alors l’équation caractéristique. Inversement, toute matrice A admet un polynôme P annulateur tel que P(A)=024, le polynôme caractéristique de A est un polynôme annulateur (théorème de Cayley-Hamilton). Les composantes des solutions du système différentiel sont des solutions de l’équation différentielle dont l’équation caractéristique est P(x)=0. En effet :

0=P(A)y=
n
k=0
 pk Ak y = 
n
k=0
 pk y[k]

Exemple en dimension 2. Soit

A=

ab 
cd


Si b=0 alors y1′=ay1 on en déduit y1 puis y2. Supposons donc b≠ 0, alors

P(x)=x2 − x (a+d) +a db c

(on peut vérifier que P(A)=0) donc si y′=Ay alors

y1′′−(a+d)y1′+adbc=0

et y2 s’en déduit avec y1′−ay1=by2 (on peut du reste partir de cette relation pour établir l’équation d’ordre 2 vérifiée par y1). On peut ainsi résoudre tous les systèmes de dimension 2, même si la matrice A n’est pas diagonalisable.

Exercice : Résoudre de cette manière le système
desolve(y'=[[1,2],[2,1]]*y and y(0)=[1,2])

Autre exemple : système d’ordre 2 se ramenant à une équation d’ordre 2 à coefficients complexes. Les équations pour une particule chargée dans un champ magnétique constant porté par l’axe Oz et un champ électrique constant perpendiculaire (donc dans le plan Oxy), avec vitesse initiale nulle ou contenue dans le plan Oxy donnent une trajectoire plane



mx=qBẏ +qEx
my=qBẋ +qEy
 

Si on pose z=x+iy alors z vérifie l’équation

z=−i
qB
m
 ż+
qE
m
,    E=Ex+iEy

Le polynôme caractéristique de cette équation

r2=−i
qB
m
r

possède deux racines distinctes 0 et −iqB/m (mais pas le conjugué, l’équation n’est pas à coefficients réels!) donc la solution homogène est

z=α + β e
i
qB
m
t
 
,    α,β ∈ 

Le champ électrique joue ici le rôle de second membre, comme 0 est solution de l’équation caractéristique, la forme de la solution particulière est z=At, en remplaçant on obtient A=iE/B donc

z=α + β e
i
qB
m
t
 
+i
E
B
 

La forme générale des solutions est un cercle si E=0 parcouru une infinité de fois, qui se déforme sous l’effet du champ électrique en une sorte de spirale de ressort, pour une valeur seuil de E on obtient une cycloïde.

14.3.6  Allure des courbes en dimension 2.

Si on se place dans le repère propre (en prenant les vecteurs propres comme vecteurs de base), et si A a deux valeurs propres distinctes (A est alors diagonalisable), alors chaque coordonnée suit une exponentielle, dans ce repère y(t)=(α eat, β ebt) avec ab. Si a et b sont réels, l’une des exponentielles domine l’autre lorsque t→ +∞ et c’est l’inverse lorsque t→ −∞, la courbe est donc asymptote aux directions propres. Si a et b sont complexes conjugués de partie réelle non nulle, on a une spirale qui tend vers 0 d’un coté et vers l’infini de l’autre (selon le signe de la partie réelle). Si A est symétrique, alors a et b sont réels, ce cas ne peut pas se produire, de plus on peut choisir un repère propre orthonormé, les courbes ressemblent à des hyperboles. Ce sont des hyperboles si trace(A)=0 (la somme des valeurs propres vaut 0 donc le produit des coordonnées dans le repère propre vaut une constante), ces hyperboles sont équilatères si A est symétrique. Quelques exemples :




Remarque :pour un système différentiel à coefficients non constants, il n’existe pas de méthode générale de résolution. Il arrive que dans certains cas particuliers, on puisse résoudre le système, par exemple si on trouve une matrice de passage indépendante du temps ramenant le système à un système diagonal ou triangulaire : un exemple avec

A=

1+tt 
t1+t


Ou si ∫A(t) dt commute avec A, on peut prendre exp(∫A(t)) comme solution.

14.3.7  Systèmes d’ordre plus grand que 1

On se ramène à un système d’ordre 1. Par exemple deux ressorts couplés



x1=−2ω2 x12 x2
x2=ω2 x1−2ω2x2

on pose Y=(x1,x2,x_1,x_2), on a

Ẏ=



x_1
x_2
x1
x2




=



001
000
−2· ω2ω20
ω2−2· ω200








x1
x2
x_1
x_2




On délègue le calcul des valeurs propres à la machine :

Les valeurs propres sont ± i ω, ± i3 ω imaginaires pures, donc les solutions du système sont périodiques de fréquence ω et √3ω, qui sont des fréquences intrinsèques du système. Si on ajoute un second membre périodique de période Ω, lorsque Ω ≠ ω et Ω ≠ √3ω, il y a une solution particulière de fréquence Ω et les solutions sont bornées (3 fréquences), par contre si Ω=ω ou Ω=√3ω, il y a résonance.

14.3.8  Intégrales premières.

Lorsqu’on ne sait pas résoudre explicitement une équation ou un système différentiel, il peut arriver qu’on connaisse une ou des constantes du mouvement en cinématique, appelées aussi intégrales premières.

C’est le cas par exemple de l’énergie totale (mécanique plus cinétique) pour des forces conservatives. En dimension un, la connaissance de l’intégrale première énergie totale permet de ramener l’équation fondamentale de la dynamique d’ordre 2 à une équation du premier ordre à variables séparables :

1
2
 m x2V(x) = E 

soit

dx
dt
 = 
2(EV(x))
m

donc

dx
2(EV(x))
m
=dt

on peut ainsi calculer le temps en fonction de x et tracer le graphe de t en fonction de x puis le graphe de x en fonction de t par symétrie par rapport à la première bissectrice.

Exemple : calcul de la période d’un pendule, on repère une masse reliée à un fil de longueur l à un point fixe par l’angle θ formé avec la verticale (orienté vers le bas), de sorte que l’énergie potentielle de la masse est −mglcos(θ) on a donc

1
2
 m l2θ2mglcos(θ)=E

puis

θ=
2(E+mglcos(θ))
ml2

Si on lache sans vitesse initiale la masse avec un angle θ0 ∈ ]−π,π[ alors E=−mglcos(θ0) donc

θ=
2
g
l
(cos(θ)−cos(θ0))

puis

dθ
2
g
l
(cos(θ)−cos(θ0))
dt 

Pour des raisons de symétrie, la période du pendule est donc

T=4
t/θ(t)=θ0


t/θ(t)=0
 dt  =4 
θ0


0
 
dθ
2
g
l
(cos(θ)−cos(θ0))
 

L’expression à intégrer n’admet pas de primitive avec les fonctions usuelles, on l’appelle intégrale elliptique (il y a un lien avec la longueur d’un arc d’ellipse). On peut calculer une valeur numérique approchée de cette intégrale si θ0 est donné.

Pour de petites valeurs de θ0, on peut approcher cos(θ) par 1−θ2/2 et calculer l’intégrale

qui ne dépend pas de θ0. On observe que cette approximation est encore assez bonne pour θ0<π/4 (erreur<4%).

En dimension plus grande, l’existence d’intégrales premières peut permettre de connaitre la forme de la courbe intégrale et même parfois de résoudre complètement l’équation (cas du problème à deux corps ci-dessous).

Autre exemple, la découverte d’un facteur intégrant pour la forme différentielle Mdx+Ndy donne une intégrale première pour l’équation dy/dx=M/N, en effet ω=φ(Mdx+Ndy)=dV(x,y) est nul sur une courbe intégrale, donc V(x,y) est constant, les courbes intégrales sont donc les courbes de niveau de V(x,y). Une équation à variables séparables est un cas particulier, avec M ne dépendant que de x et N de y.

Pour un système autonome, E est une intégrale première si grad(E).f=0, en effet

d
dt
 E(y(t))= 
n
j=1
 
∂ E
∂ yj
 fj

Problème à deux corps  Cas d’un point de 3 soumis à une force centrale comme la gravité ou la force coulombienne :

d2 r
dt2
=−µ 
r
r3

on montre

Si on prend l’axe des x porté par E, en faisant le produit scalaire avec r :

rE cos(θ)=r.E
1
µ
  (
dr
dt
 ∧ L) . r − r

on obtient en appliquant les propriétés du produit mixte et la définition de L :

r = 
L2
µ(1+E cos(θ))

la courbe intégrale est donc une conique d’excentricité E ayant l’origine pour foyer et parcourue selon la loi des aires (l’aire balayée par le segment origine-point mobile est proportionnelle au temps).

14.3.9  Le modèle proie-prédateur

C’est un système autonome en dimension 2 pour lequel on sait calculer une intégrale première. Il se présente sous la forme

=x(aby)