\documentclass{article}
\usepackage[francais]{babel}
%\usepackage[OT1]{fontenc}
\usepackage[T1]{fontenc}
\usepackage{amsmath,amsfonts,amssymb}
\usepackage{xspace}
\newcommand{\R}{{\mathbb{R}}}
\newcommand{\C}{{\mathbb{C}}}
\newcommand{\Z}{{\mathbb{Z}}}
\newcommand{\N}{{\mathbb{N}}}

\begin{document}
Dans cette s\'erie d'articles, je vais présenter quelques aspects
de l'implémentation d'un système de calcul formel.

Comme tout programme, un système de calcul formel manipule des données. 
Une implémentation doit donc d'abord définir des 
structures de données, ensuite y appliquer des algorithmes efficaces. 
Les deux aspects sont d'ailleurs liés et il est 
nécessaire d'introduire des structures de données ad hoc pour traiter 
certains problèmes. On verra aussi que les 
problèmes inhérents au calcul formel sont souvent difficiles, leur 
résolution nécessitant de conjuguer des connaissances 
en mathématiques et en informatique. C'est d'ailleurs un champ de 
recherches très actif aujourd'hui. 

Ce premier article pr\'esente les structures de donn\'ees utilis\'ees 
dans les syt\`emes de calcul formel ; il aura donc 
une coloration informatique plus marqu\'ee que les articles suivants 
qui pr\'esenteront uniquement des algorithmes. 

\section{Problèmes spécifiques au calcul formel}

\subsection{Calcul exact et approché, types, évaluation.}
Dans les langages de programmation traditionnel (C, Pascal,...), il existe 
déjà des types permettant une représentation 
exacte des données (type entier) ou une représentation approchée 
(type flottant). Mais ces types de donnée de base 
occupent une taille fixe en mémoire, le type entier est donc
limité à un intervalle d'entiers (par exemple $[0,2^{32}-1]$ 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\`etres dont la valeur n'est 
pas connue d'autre part ; il faut donc s'affranchir de ces limites~: 
\begin{itemize}
\item pour les entiers relatifs, on utilise des entiers de 
{\em précision arbitraire}
dont la taille en mémoire est dynamique (déterminée pendant l'exécution et non
à la compilation),
\item pour les nombres complexes, on utilise un couple de nombres réels,
\item pour les rationnels, on utilise un couple d'entiers relatifs,
\item pour les irrationnels algébriques (par exemple $\sqrt{2}$), 
on utilise un polyn\^ome irréductible dont ils sont racines,
\item pour les param\`etres ($x,y,z,t...$), on utilise un type 
structuré contenant un champ de type chaine de caract\`eres pour 
repr\'esenter le nom du param\`etre et
un champ pour attribuer une valeur à (ou une hypoth\`ese sur) ce param\`etre,
\item pour les nombres transcendants (par exemple $\pi$), on est obligé
d'introduire un paramètre auquel on attribue une valeur numérique, 
qui ne sera utilisée qu'au moment où on veut une 
approximation numérique d'une expression contenant ce nombre transcendant,
on parle de constante,
\item lorsqu'on a besoin d'une approximation numérique d'un nombre,
on peut utiliser des conversions de ces types en un type flottant. On peut 
aussi pour lutter contre les erreurs 
d'arrondi utiliser des nombres flottants étendus dont la précision est 
dynamique ou même des intervalles de flottants étendus,
\item il faut aussi
un nouveau type, appelé expression ou symbolique, permettant d'appliquer
une fonction qu'on ne peut évaluer directement sur les objets pr\'ec\'edents,
par exemple $\sin(x)$. Il
doit s'agir d'une op\'eration de clôture, au sens où appliquer une fonction \`a
un objet symbolique ne nécessite pas la création d'un nouveau type
(en général on renvoie un objet symbolique).
\end{itemize}

Enfin, il faut pouvoir \'evaluer un objet (en particulier symbolique)~:
par exemple évaluer $\sin(x)$ lorsqu'on assigne une valeur \`a $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 \verb|'| 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.


\subsection{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\^eme nombre ou, ce qui revient au 
m\^eme, par soustraction décider quand un nombre est nul. 
Par exemple $4/2$ et 2 représentent le m\^eme 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~: 
\begin{itemize}
\item on ne connaît pas toujours le statut de certaines constantes
(par exemple la constante d'Euler),
\item il n'existe pas d'algorithmes permettant de déterminer
s'il existe des relations algébriques entre constantes,
\item il n'existe pas forcément une seule forme plus simple, par exemple~:
\[ \frac{(\sqrt{2}+1)x+1}{x+\sqrt{2}+1}=\frac{x+\sqrt{2}-1}{(\sqrt{2}-1)x+1} \]
Ce cas se présente fréquemment avec les extensions algébriques.
\item en pratique il peut être trop coûteux d'utiliser une forme
normale, par exemple le polynôme $\frac{x^{1000}-1}{x-1}$ possède 1000 monômes
\end{itemize}
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).
% Ici je rajouterai qq mots sur le problème de la reconnaissance du 0.

On démontre que le problème de la reconnaissance du zéro pour une classe 
d'expressions "assez générale" est un problème
pour lequel il ne peut exister d'algorithme solution. 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\^ole central dans tout syst\`eme de calcul formel
puisque sauf dans les cas les plus simples (fractions d'entiers par exemple), 
la simplification d'expressions
fait appel \`a un moment ou \`a un autre \`a des calculs
de PGCD de polyn\^omes. Le PGCD de polynômes est un algorithme 
très sollicité auquel nous consacrerons le prochain article. 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 combinatoire
(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\'eriv\'ee (les deux polyn\^omes sont premiers entre eux)~:
\begin{eqnarray*}
7\* (x+1)^{6}-6\* (x-1)^{5} &&\\
\frac{162}{49} \* x^{5}+\frac{-390}{49} \* x^{4}+\frac{1060}{49} \* x^{3}+\frac{-780}{49} \* x^{2}+\frac{474}{49} \* x+\frac{-78}{49}&&\\
\frac{157780}{729} \* x^{4}+\frac{-507640}{2187} \* x^{3}+\frac{290864}{729} \* x^{2}+\frac{-101528}{729} \* x+\frac{28028}{729}&&\\
\frac{1}{49} \* (\frac{1400328}{2645} \* x^{3}+\frac{-732888}{2645} \* x^{2}+\frac{1133352}{3703} \* x+\frac{-732888}{18515})&&\\
\frac{1}{2187} \* (\frac{2161816376832}{4669921} \* x^{2}+\frac{-555436846944}{4669921} \* x+\frac{301917024864}{4669921})&&\\
\frac{1}{907235} \* (\frac{469345063045455}{129411872} \* x+\frac{-47641670106615}{129411872})&&\\
\frac{5497465490623352995840}{209648836272383412129}
\end{eqnarray*}
Le lecteur voulant tester d'autres exemples pourra utiliser le programme 
\verb|giac| (cf. l'appendice) suivant~:
\begin{verbatim}
pgcd(a):={
  local b,r,res;
  b:=diff(a,x);
  res:=[];
  for (;b!=0;){
    res:=append(res,b);
    r:=rem(a,b);
    a:=b;
    b:=r;
  }
  return(res);
}
\end{verbatim}

\subsection{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 combinatoire évoqués dans la section précédente). 
%? je ne vois pas où sont évoqués ces pbs d'explosion combinatoire.
La plupart des systèmes de calcul formel contournent la difficulté en 
supposant que le paramètre possède une valeur 
générique (par exemple la solution de $(t^2-1)x=t-1$ sera $x=1/(t+1)$) ou 
choisissent une branche pour les fonctions 
possédant un point de branchement (par exemple pour résoudre $x^2=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é.

\section{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. Il faut un type générique
recouvrant ces divers types de scalaire.

On peut 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ée\footnote{Les
systèmes de calcul formel embarqués (calculatrices) ou devant 
gérer des calculs intensifs utilisent d'ailleurs des
méthodes spécifiques pour gérer le problème de la fragmentation de
la mémoire, appelés ``garbage collector''. Ce type de méthode
est intégré dans des langages comme Lisp ou Java, en C/C++ on trouve
des libraries pour cela, par exemple GC de Boehm.}). 
En programmation orientée objet, on utiliserait plutôt un
type abstrait dont dérivent ces différents scalaires. 

Il faut aussi un type pour les vecteurs, les matrices et les
listes. 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
\verb|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\'e, les syst\`emes de calcul formel
utilisent souvent des repr\'esentations particuli\`eres pour les polyn\^omes
dont on a dit qu'ils jouaient un r\^ole central.
Pour les polyn\^omes \`a une variable,
on peut utiliser la liste des coefficients du polyn\^ome, on parle
alors de repr\'esentation dense. On peut aussi d\'ecider de ne stocker
que les coefficients non nuls, on parle alors de repr\'esentation creuse
(on stocke alors un couple form\'e par le coefficient et le degr\'e
du mon\^ome correspondant). Pour les polyn\^omes \`a plusieurs variables,
on peut les consid\'erer comme des polyn\^omes \`a une variable \`a
coefficients polynomiaux, on parle alors de repr\'esentation r\'ecursive.
On peut aussi d\'ecider de ne pas briser la sym\'etrie entre les
variables (pas de variable principale), on parle alors de repr\'esentation
distribu\'ee, le plus souvent les repr\'esentation distribu\'ees
sont creuses car les repr\'esentations
denses n\'ecessitent tr\`es vite beaucoup de coefficients. Les m\'ethodes
de repr\'esentation creuses sont parfois aussi utilis\'ees pour les
matrices ayant beaucoup de coefficients nuls.

Voyons maintenant plus précisément sur quelques exemples de logiciels
de calcul formel r\'epandus quelles structures de données sont
utilisées. Plusieurs \'el\'ements entrent en compte dans les choix faits~:
\begin{itemize}
\item le(s) profil(s) d'utilisation (enseignement, ing\'eni\'erie,
calcul intensif, recherche)
\item les ressources disponibles (m\'emoire, puissance du processeur...)
\item la facilit\'e d'impl\'ementation (choix du langage, outils
disponibles en particulier d\'ebuggueurs, ...)
\item l'histoire du syst\`eme (un syst\`eme conçu avec les outils
disponibles aujourd'hui est forc\'ement diff\'erent d'un syst\`eme 
conçu il y a 20 ans)
\end{itemize}
Nous allons d'abord parler des calculatrices formelles HP
et TI (le lecteur pourra facilement les tester grâce aux \'emulateurs gratuits
pour PC). Ce sont des systèmes plutôt destinés à l'enseignement, soumis 
\`a de fortes contraintes en termes de taille m\'emoire, et destinés
à traiter des petits problèmes.
Puis nous pr\'esenterons des syst\`emes pour ordinateur où les ressources
(par exemple m\'emoire) sont moins limit\'ees ce qui permet 
d'utiliser des langages de programmation de plus haut niveau.

\subsection{Calculatrices formelles HP}
Les langages utilis\'es pour programmer ces calculateurs sont l'assembleur
et le RPL (Reverse Polish Lisp) adapt\'e \`a l'\'ecriture de code
en m\'emoire morte tr\`es 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 quartets\footnote{un quartet=un demi octet}) 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 et HP40 utilisent 
le type programme pour représenter les entiers de Gau\ss\ (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 \verb|SEMI|). Ces types sont eux-m\^emes des 
objets et peuvent donc \^etre utilis\'es de mani\`ere
r\'ecursive. 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\^omes denses (avec 
représentation récursive pour les polyn\^omes à plusieurs variables). 

Les calculatrices HP4xG utilisent une pile\footnote{Plus pr\'ecis\'ement deux
piles, la pile de donn\'ee et la pile g\'erant le flux d'ex\'ecution. Cette
derni\`ere n'est pas visible par l'utilisateur}, 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 \verb|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 \verb|{a b +}|.
Les objets symboliques sont donc représenté par une suite d'objets écrit
en syntaxe polonaise inversée. L'\'evaluation d'un objet symbolique se fait
dans l'ordre polonaise invers\'e~: les arguments sont \'evalu\'es
puis les fonctions leur sont appliqu\'es. Pour des raisons d'efficacit\'e, 
on repr\'esente souvent les objets composites (listes, symboliques) par 
leurs composants plac\'es sur la pile (appel\'e meta-objets).

Une rigidité de la syntaxe polonaise est
que les fonctions ont toujours un nombre fixe d'arguments\footnote{Sauf si
on utilise comme dernier argument le nombre d'arguments de la fonction ou 
si on utilise (cf. infra) un tag de début de liste d'arguments}, 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 \verb|a b + c +| ou \verb|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. \`A noter l'existence d'une commande 
\verb|QUOTE| permettant \`a  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 \verb|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 \verb|INTEGERASSUME| à cet effet.
Il n'y a pas de possibilité de faire des hypothèses ayant une portée
locale.

\subsection{Calculatrices formelles TI}
Le langage utilis\'e pour programmer ces calculatrices est le langage C
(on peut aussi bien sur \'ecrire du code en assembleur).
On retrouve ici les diff\'erents types de donn\'ees regroup\'e en un
type g\'en\'erique qui est un tableau d'octets (aussi appelé quantum). 
Le champ type
est appel\'e tag dans la documentation TI. Contrairement \`a ce qui
pr\'ec\`ede, ce champ type est plac\'e en m\'emoire \`a la fin de l'objet,
ce qui est possible car la longueur d'un objet est toujours indiqu\'ee
au d\'ebut de l'objet. Ceci est fait afin de faciliter l'évaluation (cf.
infra).

Les entiers en précision arbitraire sont codés par deux tags possibles (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
\footnote{Toutefois une adaptation du logiciel utilisant comme
quantum de base par exemple 32 bits porterait cette limite
à $65536^{65535}-1$} $(256^{255}-1)$.
Il existe un tag spécifique pour les rationnels, pour les constantes 
réelles et entières qui apparaissent par exemple en r\'esolvant une \'equation.
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\'esenter
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 sur les HP, les TI utilisent une pile (non visible par
l'utilisateur) appelée expression stack
afin de traduire un expression math\'ematique sous forme d'un texte
en un objet symbolique cod\'e exactement comme ci-dessus en syntaxe
polonaise. Toutefois, la pr\'esence du champ longueur
permet d'\'evaluer un objet symbolique sans perdre en efficacit\'e
en partant de l'op\'erateur
final et en redescendant ensuite sur ces arguments, c'est la stratégie
adoptée par ces calculatrices. 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\^echer 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.

\subsection{Maple, MuPAD, Mathematica, ...}
Ces syst\`emes ont un noyau ferm\'e, au sens o\`u l'utilisateur n'a pas
acc\`es du tout, ou en tout cas pas facilement, aux structures de donn\'ees
de base. Je ne dispose donc pas d'information sur les structures de donn\'ees
utilis\'ees par le noyau (pour MuPAD, on pourrait sans doute en savoir
plus en achetant de la documentation sur la programmation des
modules dynamiques).

L'interaction système-utilisateur se fait quasiment toujours en utilisant le
langage de programmation propre au syst\`eme, langage interpr\'et\'e
par le noyau du syst\`eme (ce qui ralentit l'exécution). Ces langages 
utilisateurs sont essentiellement
non typ\'es~: on travaille avec des variables du type générique sans pouvoir
accéder aux types sous-jacents. On ne bénéficie donc pas des
v\'erifications 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.

\subsection{Giac/xcas}
Il s'agit du système de calcul formel que j'implémente actuellement sous 
forme d'une biblioth\`eque C++ (ce qui
permettra aux programmes tiers d'utiliser beaucoup plus facilement du 
calcul formel qu'avec les syst\`emes pr\'ec\'edents). L'objectif est 
d'avoir un syst\`eme facile \`a programmer directement en C++, proche 
du langage utilisateur, lui-m\^eme compatible avec Maple ou MuPAD, 
tout cela sans trop perdre en performances comparativement aux
librairies sp\'ecialis\'ees \'ecrites en C/C++. Ce qui explique un choix 
de type g\'en\'erique (\verb=gen=) non orient\'e objet, avec un champ type 
et soit une donn\'ee imm\'ediate pour les nombres flottants par exemple, 
soit un pointeur vers un objet du type correspondant au champ type pour 
les donn\'ees de taille non fixe (on pourrait donc se
contenter du langage C, mais le langage C++ permet de red\'efinir 
les op\'erateurs sur des types utilisateurs ce qui
am\'eliore consid\'erablement la lisibilit\'e du code source). 
Les donn\'ees dynamiques ne sont pas dupliqu\'ees, Giac
utilise un pointeur sur un compteur de r\'ef\'erence pour d\'etruire 
ces donn\'ees lorsqu'elles ne sont plus r\'ef\'erenc\'ees.

Les entiers en pr\'ecision arbitraire sont h\'erit\'es de la biblioth\`que
GMP (\'ecrite en C) du projet GNU. Les flottants en pr\'ecision arbitraire
utiliseront aussi GMP (plus précisément MPFR).
Il y a un type fraction, structure C compos\'e d'un champ num\'erateur
et d'un champ d\'enominateur, et un type nombre complexe.

Les listes, vecteurs, matrices utilisent le type paramétré \verb|vector<>|
de la librairie standard C++ (Standard Template Library).
Les objets symboliques sont des structures compos\'es d'un champ sommet
qui est une fonction prenant un argument de type \verb|gen|
et renvoyant un r\'esultat
de type \verb|gen|, et d'un champ feuille qui est de type \verb|gen|.
Lorsqu'une fonction poss\`ede plusieurs arguments, ils sont rassembl\'es
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\'ees pour repr\'esenter vecteurs, matrices
et polyn\^omes en une variable en repr\'esentation dense. Les polyn\^omes
en repr\'esentation creuse ou en plusieurs ind\'etermin\'ees sont \'egalement
disponibles.

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 simplement une valeur spéciale
affectée au paramètre, valeur ignorée par la routine d'évaluation.

\section{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. Le lecteur trouvera en appendice une brève présentation
de certains de ces algorithmes, mes références en la matière sont le livre
de Henri Cohen, et les livres de Donald Knuth (cf. appendice).

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).
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\'eduire la
complexit\'e (donc le temps de calcul) par des adaptations au
probl\`eme d'une m\^eme id\'ee \`a condition de faire des
hypoth\`eses sur les donn\'ees (autrement dit en abandonnant la volont\'e
d'impl\'ementer un algorithme g\'en\'erique).

Par exemple lorsqu'on travaille
avec des entiers (ou des polyn\^omes \`a coefficients entiers, ou
des matrices \`a coefficients entiers...) on utilise souvent des algorithmes
modulaires et $p$-adiques. Comme le calcul exact n\'ecessite
presque toujours de calculer avec des entiers, ces m\'ethodes
ont un r\^ole central en calcul formel, nous les pr\'esentons donc
maintenant bri\`evement. Dans les prochains articles, nous utiliserons
ce type de m\'ethode, par exemple pour le calcul de PGCD ou la factorisation
de polyn\^omes \`a coefficients entiers.

Les m\'ethodes modulaires consistent \`a r\'eduire un probl\`eme dans 
$\Z$ \`a son \'equivalent dans $Z/n\Z$ pour une ou 
plusieurs valeurs de $n$, nombre premier. Le calcul dans $\Z/n\Z$
a l'avantage de se faire avec des entiers dont la taille est bornée.
Ensuite \`a l'aide d'estimations 
\`a priori sur la taille des solutions 
\'eventuelles du probl\`eme initial, on reconstruit la solution au problème
initial avec le th\'eor\`eme des restes chinois. 

Par exemple, on peut calculer un d\'eterminant d'une matrice
\`a coefficients entiers en cherchant ce d\'eterminant dans $\Z/n\Z$
pour plusieurs premiers $n$, dont le produit est plus grand qu'une 
estimation \`a priori de la taille du d\'eterminant 
(donnée par exemple par l'inégalité d'Hadamard, cf. Cohen, p. 50). 

Les m\'ethodes $p$-adiques commencent de mani\`ere identique par un 
calcul dans $\Z/n\Z$, on augmente ensuite la
pr\'ecision de la solution en la \og liftant\fg de $\Z/n^k \Z$ vers 
$\Z/n^{k+1}\Z$ ou vers $\Z/n^{2k}\Z$ (lift
lin\'eaire ou lift quadratique), on s'arr\^ete lorsque $k$ est assez grand 
(\`a l'aide d'estimations \`a priori) et on
reconstruit alors la solution initiale. L'\'etape de \og lift\fg est en 
g\'en\'eral un lemme de Hensel dont on verra quelques exemples dans les
prochains articles. L'algorithme
commun au lemme de Hensel et au th\'eor\`eme des restes chinois est 
l'identit\'e de B\'ezout, 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\^ome $P(X)=a_d X^d + \cdots + a_0$ 
\`a coefficients entiers ou polynomiaux, avec $a_d$ et $a_0$ non nuls. 
L'algorithme g\'en\'erique (assez connu) consiste 
\`a chercher les diviseurs de $a_0$ et de $a_d$ et \`a tester toutes 
les fractions de ces diviseurs, on montre en effet 
ais\'ement que si $X=p/q$ fraction irr\'eductible est racine de $P$ 
alors $q$ divise $a_d$ et $p$ divise $a_0$. Cet 
algorithme est tr\`es inefficace si $a_d$ ou $a_0$ est un grand entier 
(car on ne sait pas forc\'ement le factoriser) ou 
s'il a beaucoup de facteurs premiers (la liste des diviseurs \`a tester 
est alors tr\`es grande). 

Lorsque les coefficients de $P$ sont entiers, la recherche pr\'ec\'edente 
revient \`a trouver un facteur \`a
coefficients entiers $qX-p$ de $P$, on peut donc r\'eduire le probl\`eme 
modulo un entier premier $n$ qui ne divise pas $a_d$~: si un tel facteur 
existe dans $\Z$ alors ce facteur (r\'eduit modulo $n$) est un facteur 
de $P$ dans $\Z/n\Z$
donc $P$ admet une racine dans $\Z/n\Z$ (puisque $q$ est inversible 
modulo $n$ car on a choisi $n$ premier ne divisant pas $a_d$). On
\'evalue maintenant $P$ en les $n$ \'el\'ements de $\Z/n\Z$. 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 $\Z/n^k\Z$ dans $\Z/n^{2k}\Z$.

On suppose donc que pour $k\geq 1$, il existe un entier $p_k$ tel que
\[ P(p_k)=0 \pmod{n^k} \]
Il s'agit de trouver un entier $x$ tel que $p_{k+1}=p_k+n^k \* x$
vérifie
\[ P(p_{k+1})=0 \pmod{n^{2k}} \]
On applique la formule de Taylor \`a l'ordre 1 pour $P$ en $p_k$, le
reste est nul modulo $n^{2k}$, donc~:
\[ P(p_k)+ n^k \* x P'(p_k)=0 \pmod{n^{2k}} \]
soit finalement~:
\[ x=-\frac{P(p_k)}{n^k} \* ( P'(p_k) \pmod{n^k}) ^{-1} \]
On reconnaît au passage la méthode de Newton, pour qu'elle fonctionne 
il suffit que $P'(p_k) \neq 0 \pmod n$ ce qui
permet de l'inverser modulo $n^k$ (et c'est ici qu'intervient 
l'identit\'e de B\'ezout). En pratique quand on factorise
un polyn\^ome, on commence par retirer les multiplicités, 
on peut donc supposer que $P$ est sans facteur multiple dans
$\Z$. Ceci n'entraîne pas forcément qu'il le reste dans $\Z/n\Z$ 
ce qui crée une contrainte supplémentaire sur le choix
de $n$, à savoir que $P$ et $P'$ restent premier entre eux dans $\Z/n\Z$ 
(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 $\Z$).

Reste donc à revenir dans $\Z$ à partir d'une racine $p_k$ dans $\Z/(n^k \Z)$
(où on peut choisir $k$). 
On va maintenant utiliser la repr\'esentation modulaire sym\'etrique~:
on prend comme représentant modulaire d'un entier $z$ dans $\Z/n^k\Z$
l'unique entier congru \`a $z$ modulo $n$ qui est strictement compris entre
$-n^k/2$ et $n^k/2$ (si $n$ est pair, la deuxi\`eme in\'egalit\'e
est choisie large).

Si $qX-p$ est un facteur de $P$, alors $a_dX-\frac{a_d}{q}p$ est encore 
un facteur de $P$ (le quotient de $P$ par $a_dX-\frac{a_d}{q}p$
est \`a coefficients rationnels mais le facteur est \`a coefficients entiers). 
Si on a choisi $k$ tel que $n^k>2|a_d a_0|$, l'\'ecriture en représentation
modulaire symétrique de $a_dX-\frac{a_d}{q}p$ est inchang\'ee,
en effet on a des estimations à priori sur les entiers $p$ et $q$~: 
$|q|\leq |a_d|$ et $|p| \leq |a_0|$ puisque $q$ 
divise $a_d$ et $p$ divise $a_0$. 
Comme $a_dX-\frac{a_d}{q}p$ est \'egal \`a $a_d(X-p_k)$ dans $\Z/(n^k \Z)$,
il nous suffit d'écrire en repr\'esentation modulaire 
sym\'etrique $a_d(X-p_k)=a_d X-p'$.
Pour conclure, on sait que $a_d X-p'$ est un multiple entier de $qX-p$.
On divise donc le facteur $a_d X-p'$ par le pgcd de $a_d$ et $p'$ et on
teste la divisibilité de $P$ par ce facteur réduit.

{\bf Exemple}\\
Consid\'erons le polyn\^ome $2 X^3-X^2-X-3$ qui est sans facteur carr\'e.
On ne peut pas choisir $n=2$ car on r\'eduirait le degr\'e, pour $n=3$,
on a $P'=X-1$ qui est facteur de $P$, pour $n=5$, $P'=6X^2-2X-1$,
on v\'erifie que $P$ et $P'$ sont premiers entre eux (par exemple
avec \verb|GCDMOD| sur une HP49 o\`u on aura fix\'e la variable \verb|MODULO|
\`a 5).

On teste ensuite les entiers de -2 \`a 2 sur $P$. Seul -1 est racine
modulo 5 ($P(-1)=-5$), on va maintenant lifter $p_1=-1$. 

L'estimation \`a priori est $2|a_d||a_0|=12$ donc $k=2$ ($5^2=25>12$), 
une it\'eration suffira. On a $P'(-1)=7$, l'inverse de $P'(-1) \pmod 5$
est -2 donc:
\[ x= -\frac{P(-1)}{5} (-2) = -(-1) \* (-2)=-2 \]
et $p_2=-1+5\times(-2)=-11$ est racine de $P$ dans $\Z/25\Z$.
On calcule ensuite $a_d(X-p_k)=2(X+11)=2X+22=2X-3$ en repr\'esentation
sym\'etrique, 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\'e 1 qui est $2X-3$.

%\pagebreak


\appendix

\section{Quelques algorithmes fondamentaux en arithmétique.}
\begin{itemize}
\item Les algorithmes de multiplication et division dit rapides
des entiers et polyn\^omes (Karatsuba, FFT, ...). Cf. par exemple Knuth.
ou pour les entiers la documentation de GMP.
\item Au lieu de la division euclidienne, on utilise très souvent la
pseudo-division pour les polynômes~: étant donné deux polyn\^omes $A$
et $B$ de degrés $a$ et $b$ à coefficients dans un anneau contenu dans un corps
(par exemple $\Z$), on multiplie $A$ par une puissance du coefficient
dominant $B_b$ de $B$, plus précisément par $B_b^{a-b+1}$, ce qui permet 
d'effectuer la division par $B$ sans que
les coefficients sortent de l'anneau.
\[ B_b^{a-b+1} A= B Q + R \]
On utilise cette méthode lorsqu'on peut multiplier les polyn\^omes par
des constantes sans changer le problème (par exemple pour l'algorithme
d'Euclide).
\item L'algorithme d'Euclide est un algorithme \og générique\fg de calcul
de PGCD. Il n'est en général pas utilisé tel quel. Pour les entiers 
on utilise une variation adaptée à la
représentation binaire des entiers (cf. Cohen ou le manuel de GMP version 4 
pour plus de détails). Nous décrirons des
algorithmes de PGCD plus efficaces pour les polynômes dans le prochain article.
\item l'identité de Bézout, aussi appelée PGCD étendu. \'Etant donné
deux entiers ou deux polyn\^omes $a$ et $b$ on calcule $u$, $v$ et
$d$ tels que $au+bv=d$. On écrit la matrice~:
\[ \left( \begin{array}{lll}
a & 1 & 0 \\
b & 0 & 1
\end{array} \right) \]
où on remarque que pour chaque ligne le coefficient de la 1ère colonne 
est égal à $a$ multiplié par le coefficient de la
2ème colonne additionné à $b$ multiplié par le coefficient de la 
3ème colonne. Ce qui reste vrai si on effectue des
combinaisons linéaires de lignes (type réduction de Gau\ss). 
Comme on travaille dans les entiers ou les polyn\^omes, on remplace la
réduction de Gau\ss\ des matrices à coefficients réels par une combinaison 
linéaire utilisant le quotient {\em euclidien\/} $q$
de $a$ par $b$. On obtient alors le reste $r$ en 1ère colonne~:
\[ L_3=L_1-qL_2 \quad \left( \begin{array}{lll}
a & 1 & 0 \\
b & 0 & 1 \\
r & 1 & -q
\end{array} \right) \]
et on recommence jusqu'à obtenir 0 en 1ère colonne.
L'avant-dernière ligne obtenue est l'identité de Bézout (la dernière
ligne donne le PPCM de $a$ et $b$). Si l'on veut l'inverse de $a$ modulo
$b$ on remarque qu'il n'est pas utile de calculer les coefficients
appartenant à la 3ème colonne. Enfin, les lignes intermédiaires
peuvent servir à reconstruire une fraction d'entier représentée
par un entier de $\Z/n\Z$ lorsque le numérateur et le dénominateur
sont de valeur absolue inférieure à $\sqrt{n/2}$.
\item Le théorème des restes chinois. Si on connaît $x=a \pmod m$
et $x= b \pmod n $ avec $m$ et $n$ premiers entre eux,
on détermine $c$ tel que
$x=c \pmod{m\times n}$ ($c=a+mu=b+nv$ et on applique
Bézout pour trouver $u$ et $v$, on en d\'eduit $c$).
\item Les tests de pseudo-primalité. Il est essentiel d'avoir une
méthode rapide permettant de générer des nombres premiers pour appliquer
des méthodes modulaires et $p$-adiques. On utilise par exemple le
test de Miller-Rabin, qui prolonge le petit théorème de Fermat
(si $p$ est premier, alors $a^p=a \pmod p$).
\end{itemize}

\section{Pour en savoir plus.}
Sur des aspects plus th\'eoriques~:
\begin{itemize}
\item Knuth: TAOCP (The Art of Computer Programming), volumes 1 et suivants
\item Henri Cohen: A Course in Computational Algebraic Number Theory
\item Davenport, Siret, Tournier: Calcul formel: Syst\`emes et algorithmes 
de manipulations  alg\'ebriques
%\item quelques articles en ligne dont un article sur l'\'evaluation:\\
%\verb|http://www.cs.berkeley.edu/~fateman/algebra.html|
\end{itemize}

Sur des aspects plus pratiques, quelques r\'ef\'erences en ligne, 
la plupart sont accessibles  gratuitement~:
\begin{itemize}
\item le code source de Giac disponible à l'URL~:\\
\verb|http://www-fourier.ujf-grenoble.fr/~parisse/giac.html|
\item le code source de GiNaC, cf.~:
\verb|http://www.ginac.de|
\item le site \verb|http://www.hpcalc.org| pour les calculatrices HP,
on y trouve tout, de la documentation, des \'emulateurs de calculatrices HP, des outils de d\'eveloppement pour Windows
et Unix/Linux, ... Pour ce qui concerne cet article, je conseille de lire\\
\verb|http://www.hpcalc.org/hp48/docs/programming/rplman.zip|
\item le site \verb|http://www.ticalc.org|, on y trouve le portage
tigcc du compilateur C de GNU, des \'emulateurs, etc. Des informations de 
cet article ont leur source dans le guide du
d\'eveloppeur TI89/92\\
\verb|http://education.ti.com/developer/8992/download/download.html|

\item la librairie du système \verb|MuPAD| (archivée dans le fichier
\verb|share/lib.tar| des distributions Unix), on peut se procurer
MuPAD par exemple sur~:\\ \verb|ftp://ftp.inria.fr/lang/MuPAD|\\ cf. \verb|www.sciface.com| pour obtenir une licence
d'utilisation.
\item pour les personnes ayant acc\`es \`a Maple, il est possible de
décompiler une instruction \verb|Maple| avec la commande\\
\verb|eval(instruction);|\\
après avoir tapé\\
\verb|interface(verboseproc=2);|
\item le source du plus ancien système de calcul formel \verb|maxima|
(devenu logiciel libre) pour les personnes famili\`eres du langage Lisp\\
\verb|http://www.math.utexas.edu/users/wfs|
\item le source de librairies plus spécialisées (GMP, GP-PARI, Singular,
NTL, Zen, ALP, Cocoa, ...), cf. le site Scientific Application for Linux\\
\verb|http://www-sor.inria.fr/mirrors/sal/index.shtml|
\end{itemize}


\end{document}