\documentclass[XUPS,XML,SOM,Unicode,francais, NoFloatCountersInSection, NoEqCountersInSection]{cedram}
\usepackage{xups97-03}
\setcounter{tocdepth}{2}
%\XUPScorrections

\begin{document}
\frontmatter

\title[Calcul formel: ce qu'il y a dans la boîte]{Calcul formel:\\ ce qu'il y a dans la boîte}
\author[\initial{P.} \lastname{Zimmermann}]{\firstname{Paul} \lastname{Zimmermann}}

\address{INRIA Lorraine}
\email{Paul.Zimmermann@inria.fr}
\urladdr{https://members.loria.fr/PZimmermann/}

\thanks{Journées X-UPS 1997. Calcul formel. Prépublication du Centre de mathématique de l'École polytechnique, 1997}

\maketitle
\tableofcontents
\mainmatter

\section*{Introduction}

Le calcul formel est aujourd'hui en plein essor. Pour preuve, on dénombre
à l'heure actuelle pas moins de sept systèmes généraux (Axiom, Derive,
Macsyma, Maple, Mathematica, MuPAD, Reduce), une multitude de systèmes
spécialisés (GAP, GB, Macaulay, Magma, Pari/GP, \ldots),
et une abondante littérature, y compris en français
\cite{CoTe95,DaSiTo86,GoSaZi95,LePo95}.
Cet essor s'explique d'une part par l'accessibilité de ces systèmes
(ils fonctionnent pour la plupart sur un Mac ou un PC),
leur apprentissage aisé, lié principalement au fait que ce sont des
langages interprétés (et non compilés comme Pascal, C ou Fortran)
disposant d'une interface agréable et d'une bonne aide en ligne,
d'autre part par l'effort d'enseignement qui est fait, notamment en
France avec l'introduction du calcul formel dans les classes préparatoires
aux grandes écoles.

Si les systèmes de calcul formel commencent à être bien connus,
du moins dans le monde universitaire (il faudra encore quelques années
pour que les formations actuelles portent leurs fruits dans l'industrie),
on sait en général moins bien comment ces systèmes fonctionnent.
Comment fait Axiom pour factoriser le polynôme $x^5+x+1$?
Quelle méthode emploie Macsyma pour factoriser l'entier $2^{64}+1$?
Comment Mathematica calcule-t-il $\textstyle\int_0^2 \sqrt{x+1/x-2} \, dx$?
Quel principe utilise Maple pour déterminer la valeur de
$\sum_{k=0}^n \binom{n}{k}^2$?
Comment MuPAD résout-il la récurrence $n u_{n+2} = 5 u_{n+1} + (n+1) u_n$?
La méconnaissance des méthodes utilisées entraîne de fait
celle des classes de problèmes que les systèmes savent résoudre.
La conséquence est que les systèmes de calcul formel sont un peu
considérés comme des \og boîtes magiques\fg, et que l'on n'a aucune
idée a priori de la difficulté intrinsèque d'un problème donné
pour les systèmes de calcul formel.
Lorsque le système n'arrive pas à résoudre le problème
donné, est-ce parce qu'il n'est pas assez puissant ou parce qu'il
n'y a pas de réponse?
Certaines méthodes utilisées par les systèmes
sont des algorithmes de {\em décision}, c'est-à-dire que l'absence de
réponse positive signifie \og il n'existe pas de solution dans la
classe considérée\fg.
C'est le cas notamment de l'intégration indéfinie des fonctions
élémentaires (algorithme de Risch) et de la sommation indéfinie
des fonctions hypergéométriques (algorithme de Gosper).
La connaissance des méthodes utilisées permet alors de
déduire plus d'information que le système n'en donne apparemment.

À ce jour,
aucun ouvrage en français ne décrit l'ensemble des algo\-ri\-thmes
utilisés en calcul formel. Dans \cite{GoSaZi95}, nous avons essayé
de décrire les classes de problèmes que les systèmes savent traiter,
mais sans trop rentrer dans les détails car nous ne voulions pas
introduire trop de concepts mathématiques.
En anglais, il existe un ouvrage général
mais malheureusement cher \cite{GeCzLa92},
et deux ouvrages spécialisés que nous recommandons,
le premier sur la sommation et la preuve d'identités discrètes
\cite{PeWiZe96}, le second sur l'intégration \cite{Bronstein97}.
Citons également le livre en préparation {\em Modern Computer Algebra\/},
par Joachim von zur Gathen et Jürgen Gerhard,

Ce document, issu de notes de cours de DEA
{\em Quelques algorithmes de calcul symbolique\/} donné à Nancy 1,
indique de manière succincte le fonctionnement de cinq algorithmes
fondamentaux utilisés par les systèmes de calcul formel:
l'algorithme de Gosper pour le calcul de sommes indéfinies $\sum_n f_n$,
l'algorithme de Zeilberger pour le calcul de sommes définies
$\sum_{k=0}^n f_{n,k}$,
l'algorithme de Berlekamp pour la factorisation de polynômes à coefficients
dans un corps fini,
l'algorithme de Zassenhaus pour la factorisation de polynômes à
coefficients entiers, et
l'algorithme de Lenstra pour la factorisation d'entiers.
Chaque section peut être lue indépendamment des autres, sauf
l'algorithme de Zeilberger qui utilise celui de Gosper,
et l'algorithme de Zassenhaus qui recourt à celui de Berlekamp.
Des références bibliographiques sont fournies, où l'on pourra
trouver des informations complémentaires sur ces algorithmes,
notamment leur complexité, qui n'est pas abordée ici.

\section{L'algorithme de Gosper pour le calcul de \texorpdfstring{$\sum_n f_n$}{sumf}}

On trouvera une description plus détaillée de cet algorithme dans~\cite{PeWiZe96}, avec notamment une preuve des théorèmes.
\begin{df}[Fonction hypergéométrique]
Une suite $f_n$ est dite {\em hyper\-géo\-mé\-tri\-que\/} (en $n$) ssi
le rapport $f_{n+1}/f_n$ est rationnel (en $n$).
\end{df}
\begin{ex}
$2^n$ et $\binom{m}{n}$ sont hypergéométriques en $n$ (et $m$),
mais pas $\sin n$.
\end{ex}

\begin{df}[Somme indéfinie]
La somme indéfinie de $f_n$ par rapport à $n$, notée
$\sum_n f_n$, est une expression $g_n$ telle que $f_n=g_{n+1}-g_n$.
\end{df}
\begin{rem}
Comme pour l'intégration indéfinie, $g_n$ est définie à une
constante près.
On peut aussi choisir la définition $f_n=g_n-g_{n-1}$, qui mène à
des calculs différents.
\end{rem}
\begin{ex}
On a $\sum_n n = n(n-1)/2$ et $\sum_n n^2 = n (n-1)(2n-1)/6$.
\end{ex}

En 1978, Gosper a mis au point un algorithme permettant de décider si
une fonc\-tion hypergéométrique admet ou non une somme (indéfinie)
hypergéométrique \cite{Gosper78}.

\noindent
\begin{quote}\smaller
\noindent
{\bf Algorithme de {\sc Gosper}.} \\
Entrée: une suite hypergéométrique $f_n$. \\
Sortie: une suite hypergéométrique $g_n$ telle que $g_n = \sum_n f_n$
ou bien {\tt fail}. \\
1. Calculer $r(n) = f_{n+1}/f_n$ (fraction rationnelle en $n$). \\
2. Mettre $r(n)$ sous la forme $\frac{a(n)}{b(n)} \frac{c(n+1)}{c(n)}$
où $a(n), b(n), c(n)$ sont des polynômes en $n$ tels que
si $\alpha$ est racine de $a(n)$, alors $\alpha+h$ n'est pas
racine de $b(n)$, i.e.
\begin{equation} \label{abc}
\pgcd(a(n),b(n+h))=1 \ \mbox{pour tout entier $h \ge 0$.}
\end{equation}
3. Chercher un polynôme $x(n)$ solution de l'équation
\begin{equation} \label{gosper}
a(n) x(n+1) - b(n-1) x(n) = c(n).
\end{equation}
4. Si un tel $x(n)$ existe, renvoyer $g_n = \frac{b(n-1) x(n)}{c(n)} f_n$,
sinon renvoyer {\tt fail}.
\end{quote}
\begin{ex}
Soit $f_n = (4n+1) \frac{n!}{(2n+1)!}$. Alors $r(n)=\frac{4n+5}{2(2n+3)(4n+1)}$.
% f:=(4*n+1)*fact(n)/fact(2*n+1);
% r:=normal(expand(subs(f,n=n+1)/f));
Les polynômes $a(n)=1, b(n)=2(2n+3), c(n)=4n+1$ satisfont la condition 2,
et l'équation (\ref{gosper}) devient $x(n+1)-2(2n+1)x(n)=4n+1$, dont une
solution évidente est $x(n)=-1$. Donc $\sum_n f_n = -2 n!/(2n)!$.
% \begin{verbatim}
% >> f:=(4*n+1)*n!/(2*n+1)!: g:=sum(f, n);
%
%  - 2 fact(n) - 4 n fact(n)
%  -------------------------
%  fact(2 n + 1)
%
% >> normal(expand(subs(g,n=n+1)-g-f));
%
%   0
% \end{verbatim}
\end{ex}

\subsection{Mise sous la forme \texorpdfstring{$\frac{a(n)}{b(n)} \frac{c(n+1)}{c(n)}$}{abc}}

% Z:=1/2: t:=4*n+5: u:=(2*n+3)*(4*n+1):
% R:=resultant(t,subs(u,n=n+h),n); % 32 (4 h + 1) (h - 1)
Soit $r(n) = Z t(n)/u(n)$, où~$Z$ est une constante,
et $t(n), u(n)$ sont des polynômes premiers
entre eux et unitaires (de coefficient dominant $1$).
On calcule tout d'abord $R(h) = {\rm res}_n(t(n),u(n+h))$, puis
$S=\{ h_1, \ldots, h_N \}$, $h_1 \le \cdots \le h_N$,
l'ensemble des racines entières positives ou nulles de $R(h)$.
\begin{quote}\smaller
\noindent
$p_0(n) \leftarrow t(n), q_0(n) \leftarrow u(n)$ \\
{\tt for} $j=1$ {\tt to} $N$ {\tt do} \\
\hspace*{1cm} $s_j(n) \leftarrow \pgcd(p_{j-1}(n),q_{j-1}(n+h_j))$ \\
\hspace*{1cm} $p_j(n) \leftarrow p_{j-1}(n)/s_j(n)$ \\
\hspace*{1cm} $q_j(n) \leftarrow q_{j-1}(n)/s_j(n-h_j)$ \\
$a(n) \leftarrow Z p_N(n)$ \\
$b(n) \leftarrow q_N(n)$ \\
$c(n) \leftarrow \prod_{i=1}^N \prod_{j=1}^{h_i} s_i(n-j)$
\end{quote}

\begin{thm} Soit $K$ un corps de caractéristique zéro, et $r \in K[n]$
une fraction rationnelle non nulle. Alors il existe un unique triplet $a,b,c$
de polynômes de $K[n]$
vérifiant {\rm(\ref{abc})} tels que $b, c$ sont unitaires,
$r(n) = \frac{a(n)}{b(n)} \frac{c(n+1)}{c(n)}$, $\pgcd(a(n),c(n))=1$
et $\pgcd(b(n),c(n+1))=1$.
\end{thm}
Les deux dernières conditions assurent l'unicité, qui n'est pas vraiment
requise pour la bonne marche de l'algorithme. Si on prend $r(n)=\sfrac{n+3}
{n (n+1)}$, alors
$$
a(n)=1,\quad b(n)=n,\quad c(n)=(n+1)(n+2)
$$
et
$$
a(n)=n-1,\quad b(n)=n^2,\quad c(n)=(n+1)(n+2)(n-1)
$$
vérifient (\ref{abc}) et
$$r(n) = \frac{a(n)}{b(n)} \frac{c(n+1)}{c(n)},
$$
mais le second triplet ne
vérifie pas $\pgcd(a(n),c(n))=1$.

\Subsection{Recherche des solutions \texorpdfstring{$x(n)$}{xn} de l'équation de Gosper}

\begin{thm}[\cite{Gosper78}]
Soient $a(n), b(n), c(n)$ des polynômes en $n$ vérifiant {\rm(\ref{abc})}.
Alors toute fraction rationnelle $x(n)$ solution de {\rm(\ref{gosper})}
est nécessairement un polynôme.
\end{thm}
Pour déterminer le degré $d$ d'une éventuelle solution $x(n)$ de
(\ref{gosper}),
on distingue deux cas suivant les degrés et coefficients dominants de
$a(n)$ et $b(n)$.
\subsubsection*{Cas 1: ${\rm deg}\,a(n) \neq {\rm deg}\,b(n)$ ou
${\rm lc}\,a(n) \neq {\rm lc}\,b(n)$}
Dans ce cas le degré du membre gauche de (\ref{gosper}) est $d +
{\rm max}\{ {\rm deg}\,a(n), {\rm deg}\,b(n) \}$, et donc on~a:
\[ d = {\rm deg}\,c(n) - {\rm max}\{ {\rm deg}\,a(n), {\rm deg}\,b(n) \}. \]

\subsubsection*{Cas 2: ${\rm deg}\, a(n) = {\rm deg}\, b(n)$ et
${\rm lc}\, a(n) = {\rm lc}\, b(n) = \lambda$}
Les termes dominants de (\ref{gosper}) s'annulent. On distingue à nouveau
deux sous-cas (qui ne sont pas incompatibles, par conséquent il faut
prendre le plus grand degré obtenu, cf.\ $\sfrac{1}{n(n+1)}$, où (2a) donne
$d=0$ et (2b) donne $d=1$):
\begin{quote}\smaller
\noindent
{\bf (2a)} Si les termes de second ordre du membre gauche de~(\ref{gosper}) ne
s'annulent pas, alors:
\[ d = {\rm deg}\, c(n) - {\rm deg}\, a(n) + 1. \]
{\bf (2b)} Si les termes de second ordre du membre gauche de~(\ref{gosper})
s'annulent, soient $a(n) = \lambda n^k + A n^{k-1} + \cdots$,
$b(n-1) = \lambda n^k + B n^{k-1} + \cdots$, alors nécessairement
$d = \psfrac{B-A}{\lambda}$.
\end{quote}

Voici deux exemples avec MuPAD 1.3: le premier (cas 1) indique que la somme
indéfinie de $(4n+1) \frac{n!}{(2n+1)!}$ est $-2 n!/(2n)!$ (après
simplification), le second (cas 1) que $n!$ n'a pas de somme indéfinie
hypergéométrique.
\begin{small}
\begin{verbatim}
>> sum((4*n+1)*n!/(2*n+1)!, n);

                         - 2 fact(n) - 4 n fact(n)
                         -------------------------
                               fact(2 n + 1)

>> sum(n!, n);       

                              sum(fact(n), n)
\end{verbatim}
\end{small}

\section{L'algorithme de Zeilberger pour le calcul de \texorpdfstring{$\sum_{k=0}^n F(n,k)$}{sumF}}

L'algorithme de Zeilberger \cite{Zeilberger91b}
permet de calculer des sommes définies
$$f(n) = \sum_{k=0}^n F(n,k)$$ même lorsque l'expression $F(n,k)$ n'admet pas
de somme indéfinie $\sum_k F(n,k)$, i.e. lorsque l'algorithme de Gosper
échoue. C'est un peu l'analogue de la méthode de Laplace pour calculer
des intégrales défi\-nies. L'exemple type est $F(n,k) =\binom{n}{k} x^k$:
\bgroup
\smaller
\begin{verbatim}
>> sum(binomial(n,k)*x^k, k), sum(binomial(n,k)*x^k, k=0..n);

                         k                        n
                    sum(x  binomial(n, k), k), (x + 1)
\end{verbatim}
\egroup
L'idée de Zeilberger est d'essayer de trouver une formule de récur\-rence
pour $f(n)$ en fonction de $n$, que l'on peut ensuite résoudre pour trouver
une forme close. On peut aussi utiliser directement la récurrence et le
calcul des premiers termes pour montrer des identités du genre:
\[ \sum_{s=0}^{2m} (-1)^s \binom{2m}{s}^3 = (-1)^m \frac{(3m)!}{(m!)^3}. \]
\begin{minipage}{\textwidth}
\smaller \smaller \smaller
\begin{verbatim}
>> sum((-1)^s*binomial(2*m,s)^3,s=0..2*m);

     /    /               3 u1(m) (3 m + 1) (3 m + 2)                     \ \
     |    | u1(m + 1) = - ---------------------------, u1(m), {u1(0) = 1} | |
solve| rec|                               2                               | |
     \    \                        (m + 1)                                / /
\end{verbatim}
\end{minipage}

\medskip
L'astuce de génie de Zeilberger est d'utiliser une version étendue de
l'algorithme de Gosper
(qui normalement calcule des sommes indé\-fi\-nies) pour justement calculer
des sommes {\em définies}, en ajoutant des paramètres qui seront
déterminés au cours de l'algorithme.

\Subsection{Fonctions proprement hypergéométriques}

\begin{df}[{\cite[p.~64]{PeWiZe96}}] \label{prhy}
Une fonction $F(n,k)$ est {\em proprement hypergéométrique\/}
si elle peut se mettre sous la forme
\[ F(n,k) = P(n,k) \frac{\prod_{i=1}^u (a_i n + b_i k + c_i)!}
{\prod_{i=1}^v (a'_i n + b'_i k + c'_i)!} \;x^k \]
où $P$ est un polynôme, les $a_i, b_i, a'_i, b'_i$ sont des entiers,
et $u, v$ sont des entiers positifs ou nuls.
\end{df}
\begin{ex}
L'expression $\binom{n}{k} 2^k=\frac{n!}{k! (n-k)!} 2^k$ est
proprement hypergéométrique, ainsi que $\sfrac{1}{n+3k+1} =
\sfrac{(n+3k)!}{(n+3k+1)!}$, mais pas $\sfrac{1}{n^2+k^2+1}$.
\end{ex}

\begin{thm}[Wilf, Zeilberger] \label{wz}
Soit $F(n,k)$ une fonction proprement hypergéométrique.
Alors $F$ satisfait une récurrence de la forme
\begin{equation} \label{telescope}
\sum_{j=0}^J a_j(n) F(n+j,k) = G(n,k+1)-G(n,k)
\end{equation}
où de plus $G(n,k)/F(n,k)$ est une fraction rationnelle en $n$ et $k$.
\end{thm}
Non seulement une telle récurrence existe donc pour toute fonction
proprement hypergéométrique, mais en plus Wilf et Zeilberger ont
prouvé une borne sur l'ordre $J$: $J \le |b_1| + \cdots +|b_u| +
|b'_1| + \cdots + |b'_v|$ où les $b_i, b'_i$ sont ceux de la
définition~\ref{prhy}.

\noindent
\begin{quote}\smaller
\noindent
{\bf Algorithme de {\sc Zeilberger}.} \\
Entrée: une expression $F(n,k)$ hypergéométrique en $n$ et~$k$
et un entier $J$. \\
Sortie: une récurrence linéaire d'ordre au plus $J$ pour $f(n) =
\sum_{k=0}^n F(n,k)$ ou bien {\tt fail}. \\
1. Calculer $r_{n,k} = F(n,k)/F(n-1,k)$. \\
2. Calculer les fractions rationnelles $R_0=1$ et $R_j=R_{j-1} r_{n+j,k}$
pour $1 \le j \le J$. \\
3. Appliquer l'algorithme de Gosper étendu (cf ci-dessous)
par rapport à $k$ à $(R_0 + C_1 R_1 + \cdots
+ C_J R_J) F(n,k)$ où les~$C_j$ sont des indéterminées
ne dépendant pas de $k$. \\
4. Si Gosper retourne {\tt fail}, faire de même. Sinon, soit $G(n,k)$
la somme indéfinie retournée, et les $C_j$ sont des fractions
rationnelles en $n$. \\
5. Retourner la récurrence
\begin{equation} \label{rec}
f(n) + C_1 f(n+1) + \cdots + C_J f(n+J) = G(n,n+1)-G(n,0).
\end{equation}
\end{quote}

\subsection{Version étendue de l'algorithme de Gosper}

Les paramètres $C_j$ ajou\-tés par Zeilberger interviennent de façon
linéaire dans le polynôme $a(k)$ de
la forme de Gosper-Petkov\v{s}ek
\[
\frac{a(k)}{b(k)}\, \frac{c(k+1)}{c(k)},
\]
et donc aussi
dans l'équation de Gosper
\[ a(k) x(k+1) - b(k-1) x(k) = c(k). \]
On applique donc le raisonnement usuel de l'algorithme de Gosper (cas 1,
2a et 2b) sauf qu'en plus des coefficients indéterminés de $x(k)$, on
a aussi les coefficients $C_j(n)$.

\subsection{Certificat}

Outre la récurrence (\ref{rec}), l'algorithme de Zeilberger fournit
un {\em certificat\/} $R(n,k)=G(n,k)/F(n,k)$,
qui est une fraction rationnelle d'après le théorème~\ref{wz},
et
qui permet de vérifier très facilement que cette récurrence est correcte.
En effet, l'algorithme de Gosper étendu calcule $R(n,k)$ et des $C_j$ tels
que
\[ \sum_k F(n,k) + C_1 F(n+1,k) + \cdots + C_J F(n+J,k) = G(n,k) \]
soit par définition de la sommation indéfinie par rapport à $k$,
\begin{multline}
F(n,k) + C_1 F(n+1,k) + \cdots + C_J F(n+J,k)\\
= R(n,k+1) F(n,k+1) - R(n,k) F(n,k).
\end{multline}
Pour vérifier une identité $\sum_{k=0}^n F(n,k) = C \neq 0$ à partir
de son certificat $R(n,k)$, il suffit de former l'expression
\[ F(n,k) - F(n+1,k) - R(n,k+1) F(n,k+1) + R(n,k) F(n,k), \]
de diviser tout par $F(n,k)$, ce qui donne une expression rationnelle
dont on peut facilement tester la nullité.

\section{L'algorithme de Berlekamp pour la factorisation dans~\texorpdfstring{$\F_p[x]$}{Fp}}

Berlekamp fut le premier à proposer une méthode systématique
pour factoriser des polynômes sur un corps fini \cite{Berlekamp67}:\enlargethispage{-2\baselineskip}
\begin{quote}\smaller
\noindent
{\bf Algorithme de Berlekamp.} \\
Entrée: un polynôme $f \in \F_p[x]$. \\
Sortie: l'ensemble des facteurs irréductibles de $f$ et leur
multiplicité. \\
1. [SQF] Effectuer une décomposition sans carré de $f$. \\
2. [DDF] Séparer les facteurs de degrés différents. \\
3. [EDF] Séparer les facteurs de même degré.
\end{quote}

\subsubsection*{Factorisation sans carré}

C'est l'étape la plus facile, car elle ne nécessite que des pgcds:
\begin{thm}
Soit $f(x)$ un polynôme unitaire de $R[x]$, où $R$ est un domaine à
factorisation unique de caractéristique zéro.
Alors $f(x)$ est
sans facteur carré si et seulement si $\pgcd(f(x),f'(x))=1$.
\end{thm}

\begin{small}
\begin{verbatim}
>> f:=poly(x^11+x+1,IntMod(6449)): gcd(f,diff(f,x));      

                        poly(1, [x], IntMod(6449))
\end{verbatim}
\end{small}
On en déduit donc que $f=x^{11}+x+1 \bmod 6449$ est sans facteur carré.
Il existe d'autres méthodes plus efficaces pour calculer la décomposition
sans carré, notamment l'algorithme de Yun \cite{Bronstein97}.
Dans le cas de $R=\F_p$, une difficulé apparaît car $f'(x)$ peut
être nul sans que $f(x)$ soit constant, par exemple $f(x)=x^{13}+1$ dans
$\F_{13}[x]$. Mais dans ce cas $f(x)$ est nécessairement de la forme
$g(x)^p$.

\subsubsection*{Séparation des facteurs de degrés différents}

Le théorème suivant peut servir à
trouver le produit des facteurs de même degré.
\begin{thm} \label{berlekamp}
Soit $f \in \F_p[x]$ sans facteur carré. Alors le produit des facteurs
de $f$ de degré divisant $k$ est $\pgcd(x^{p^k}-x,f) \bmod p$.
\end{thm}
Si $p^k$ est petit, on peut calculer directement ce pgcd:
\begin{small}
\begin{verbatim}
>> d1:=gcd(poly(x^6449-x,IntMod(6449)),f);

                               2
                poly(2391 x + x  + 2155, [x], IntMod(6449))
\end{verbatim}
\end{small}
Mais lorsque
$p^k$ devient grand, il vaut mieux calculer $x^{p^k} \bmod f$
par carrés et réduc\-tions successives modulo $f$
avant d'en calculer le pgcd avec $f$.
\bgroup
\smaller\smaller\smaller\begin{verbatim}
>> f2:=divide(f,d1,Quo): xp2:=powermod(poly(x,IntMod(6449)),6449^2,f2):
>> d2:=gcd(xp2-poly(x,IntMod(6449)),f2);

                    2         3         4         5    6
poly(-614x - 167x - 2321x - 2101x - 2548x + x - 395,[x],IntMod(6449))
\end{verbatim}
\egroup
Voici donc le produit des facteurs de degré $2$ (il y en a trois).
Dans cet exemple on a fini puisqu'il ne reste plus qu'un polynôme de
degré trois, qui est forcément irréductible de par le théorème
\ref{berlekamp}.
\begin{small}
\begin{verbatim}
>> f3:=divide(f2,d2,Quo);

                               2    3
            poly(3211x + 157x + x + 150, [x], IntMod(6449))
\end{verbatim}
\end{small}
Il nous reste donc à séparer les deux facteurs linéaires et les
trois facteurs de degré deux.\enlargethispage{\baselineskip}

\subsubsection*{Séparation des facteurs de même degré}
Cantor et Zassenhaus \cite{CaZa81} ont trouvé un algorithme probabiliste
pour séparer les facteurs de même degré, basé sur l'identité
$x^p-x = x (x^{(p-1)/2}-1)(x^{(p-1)/2}+1)$ pour $p$ impair,
qui impose aux facteurs de $x^p-x$ de se scinder entre $x^{(p-1)/2}-\nobreak 1$
et $x^{(p-1)/2}+1$.
\begin{quote}\smaller
\noindent
Entrée: un polynôme $f$ sans carré, produit de \hbox{facteurs de
degré~$i$.}\\
Sortie: les facteurs irréductibles de $f$. \\
0. Si ${\rm deg}\, f = i$, retourner $f$. \\
1. Générer un polynôme unitaire aléatoire $r \in \F_p[x]$ de degré
$2i-1$. \\
2. Calculer $q = r^{p^0} + r^{p^1} + \cdots + r^{p^{i-1}}$. \\
3. Calculer $g=\pgcd(q^{(p-1)/2}-1,f)$. Si $g \neq 1$ et $g \neq f$,
recommencer en $1$ avec $g$ et $f/g$. Sinon recommencer avec~$f$.
\end{quote}

\begin{thm}[\cite{CaZa81}]
La probabilité que $\pgcd(q^{(p-1)/2}-1,f)$ soit non trivial est
\[ 1-\left(\frac{p-1}{2p}\right)^k-\left(\frac{p+1}{2p}\right)^k
\ge \frac{4}{9}. \]
\end{thm}
Donc on scinde le polynôme $f$ à l'étape 3 avec une probabilité
au moins $4/9$, ce qui donne une espérance finie pour le nombre total
d'étapes.
\bgroup
\smaller\smaller
\begin{verbatim}
>> z:=random(6449): r:=poly(z()+z()*x+z()*x^2+x^3,IntMod(6449));

                                2    3
          poly(- 1711 x - 3055 x  + x  + 174, [x], IntMod(6449))

>> q:=r+powermod(r,6449,f2): 
>> gcd(powermod(q,(6449-1)/2,f2)-q^0,f2);

                                2
               poly(- 1372 x + x  + 2487, [x], IntMod(6449))
\end{verbatim}
\egroup
On a ainsi trouvé un facteur de degré $2$, il suffit de recommencer avec
le cofacteur de degré $4$.

Le livre \cite{GeCzLa92} donne une description plus détaillée.
Il est à noter que l'étape de séparation des facteurs de degrés
différents (DDF) a été notablement améliorée par Shoup
\cite{Shoup95}. Cette amélioration est implantée en MuPAD 1.3.

\section{L'algorithme de Zassenhaus pour la factorisation dans~\texorpdfstring{$\Z[x]$}{Zx}}

Factoriser un polynôme de $\Q[x]$ revient à factoriser un polynôme de
$\Z[x]$, modulo la factorisation d'entiers.
Nous supposons donc que le polynôme est à coefficient entiers.
À titre d'exemple, nous considérons
$f=x^8+4x^7-2x^6-20x^5+3x^4+44x^3+22x^2-4x+34$ tiré de \cite{GeCzLa92}.

\subsubsection*{Factorisation sans carré}

Comme pour l'algorithme de Berlekamp pour la factorisation dans $\F_p[x]$,
on n'a besoin que de pgcds. De plus, comme on est en caractéristique zéro,
on n'a pas de problème d'annulation de $f'(x)$.
% \begin{small}
% \begin{verbatim}
% >> f:=x^8+4*x^7-2*x^6-20*x^5+3*x^4+44*x^3+22*x^2-4*x+34:
% >> gcd(f,diff(f,x));
%
%   1
% \end{verbatim}
% \end{small}

\subsubsection*{Analyse des degrés possibles des facteurs}
Contrairement à la factorisation dans $\F_p[x]$ (algorithme de Berlekamp),
il n'est ici pas possible de séparer les facteurs de même degré.
Par contre, on peut utiliser le fait que si $a(x)$ divise $f(x)$ dans
$\Z[x]$, alors l'image de $a(x)$ divise l'image de $f(x)$ dans $\F_p[x]$ pour
n'importe quel entier premier $p$.
En particulier, si on trouve un $p$ tel que $f(x)$ est irréductible dans
$\F_p[x]$, alors $f(x)$ l'est aussi dans $\Z[x]$, comme $p=2$ pour
$x^9+x+1$.
% \begin{small}
% \begin{verbatim}
% >> Factor(poly(x^9+x+1,IntMod(2)));
%
%  9
%  poly(x + x + 1, [x], IntMod(2))
% \end{verbatim}
% \end{small}
\begin{quote}\smaller
\noindent
{\bf Analyse des degrés possibles des facteurs.} \\
Entrée: un polynôme $f \in \Z[x]$ de degré $n$ et un entier premier~$p$.\\
Sortie: les degrés possibles des facteurs de $f$. \\
%0. $s \leftarrow \{ 1, \ldots, n-1\}$ \\
% et tel que $f$ soit sans facteur carré dans $\F_p[x]$. \\
1. Séparer les facteurs de degrés différents (DDF) de $f$ dans $\F_p[x]$.
\\
% Soient $d_1, \ldots, d_k$ (avec multiplicités)
% les degrés des facteurs dans $\F_p[x]$.\\
% 3. $s \leftarrow s \cap \left( \bigcup d_{i_1} + \cdots + d_{i_j}
% \right)$.
% 4. Retourner $s$
2. Si $d_1, \ldots, d_k$ sont les degrés obtenus, retourner
$\bigcup d_{i_1} + \cdots + d_{i_j}$.
\end{quote}
\begin{rem}
On peut aussi imposer que $f$ soit sans facteur carré dans $\F_p[x]$,
ce qui revient à dire que $\pgcd(f,f')=1$ dans $\F_p[x]$,
ou encore que $p$ ne divise pas ${\rm res}(f,f')$.
En fait, on peut toujours trouver un tel $p$, car d'après l'inégalité
d'Hadamard, le résultant de $f$ et $f'$ est borné par $n^{2n} 2^n
B^{2n-1}$ où $n={\rm deg} \, f$ et $B = {\rm max} \left| f_i \right|$.
\end{rem}

Par exemple, pour $x^{10}+x+1$, la factorisation modulo $2$ donne
$(x^3+x+1)(x^7+x^5+x^4+x^3+1)$, donc
comme
degrés possibles $d_2 = \{ 3, 7 \}$.
Modulo $3$, on trouve
\[
x^{10}+x+1 \equiv (x-1)(x^3-x^2-x-1)
(x^6-x^5+x^4-x^3+x+1),
\]
donc $d_3 = \{ 1, 3, 4, 6, 7, 9 \}$,
et l'intersection des degrés possibles
est $d_{2,3} = d_2 \cap d_3 = \{ 3 \}$; la factorisation modulo $5$ donne
\[
(x^2-x+2)(x^8+x^7-x^6+2x^5-x^4+2x^2+2x-2),
\]
à savoir
$d_5 = \{ 2, 8 \}$, donc $d_{2,3,5} = \emptyset$, dont on déduit
que $x^{10}+x+1$ est irréductible.

Même lorsque $f$ n'est pas irréductible, les informations sur les degrés
possibles seront utiles par la suite.
Pour notre polynôme $f$, la factorisation modulo $5$,
à savoir
\[
f = (x^2-x+1)(x^2-2x-2)(x^2-x+2)(x^2-2x-1),
\]
montre que les facteurs possibles ont degré dans $d_5 = \{ 2, 4, 6 \}$.

\subsubsection*{Remontée de Hensel}

La remontée de Hensel ({\em Hensel-lifting}) consiste à \og remonter\fg\
une factorisation dans $\F_p[x]$ en factorisation dans $\F_{p^N}[x]$.
C'est un peu comme le théorème des restes chinois, sauf qu'ici on a
un seul nombre premier.
\begin{quote}\smaller
\noindent
{\bf Remontée de Hensel} \\
Entrée: une factorisation $f = g h$ modulo $p^k$ avec $g, h$ premiers
entre eux et $h$ primitif. \\
Sortie: une factorisation $f = g^{*} h^{*}$ modulo $p^{k+1}$ avec
$g^{*} = g \bmod p^k$ et $h^{*} = h \bmod p^k$. \\
0. Calculer $s$ et $t$ tels que $s g + t h = 1$ par l'algorithme d'Euclide
étendu. \\
1. Calculer $e = f - g h \bmod p^{k+1}$. \\
2. Effectuer la division euclidienne de $s e$ par $h$:
$s e = q h + r \bmod p^{k+1}$ avec ${\rm deg}\, r < {\rm deg}\, h$.\\
3. Retourner $g^{*} \!=\! g \!+\! t e \!+\! q g \bmod p^{k+1}$ et
\hbox{$h^{*} \!=\! h\! +\! r \bmod p^{k+1}$}.
\end{quote}
Supposons que l'on prenne $h=x^2-x+1$ et
\[
g=(x^2-2x-2)(x^2-x+2)(x^2-2x-1)
= x^6+2x^4+2x^3-2x^2-1 \bmod 5.
\]
La remontée de $5^1$ vers $5^2$ donne
$g^{*} = x^6+12x^4+7x^3-7x^2-5x-6$ et $h^{*} = x^2+4x+11$.
\begin{thm}[Unicité de la remontée de Hensel]
Soit $p$ un entier pre\-mi\-er, $N$ un entier positif,
$f, g_1, h_1, g_2, h_2 \in \Z[x]$
avec $h_1, h_2$ primitifs et de même degré, $g_1 \!\equiv\! g_2 \bmod p$,
$h_1 \!\equiv\! h_2 \bmod p$, $\pgcd(g_1,h_1)\!=\!1 \bmod p$, et
$f \!\equiv\! g_1 h_1 \!\equiv\! g_2 h_2 \bmod p^N$. Alors $g_1 \!\equiv\! g_2 \bmod p^N$
et $h_1 \equiv h_2 \bmod p^N$.
\end{thm}

\begin{quote}\smaller
\noindent
{\bf Algorithme de Zassenhaus \cite{Zassenhaus69}.} \\
Entrée: un polynôme $f \in \Z[x]$ primitif et sans carré de degré
$n$.\\
Sortie: deux polynômes $g, h \in \Z[x]$ primitifs et non constants
tels que $f = g h$ ou bien \og $f$ est irréductible\fg. \\
1. Choisir un entier premier $p$ ne divisant pas ${\rm res}(f,f')$, et
$N = \lceil \log_p(2^n \left| f \right|) \rceil$. \\
2. Calculer une factorisation $f = w_1 \cdots w_k$ de $f$ dans $\F_p[x]$.\\
\hspace*{1cm} Si $k=1$, retourner \og $f$ est irréductible\fg. \\
3. Effectuer les étapes 4 et 5 pour chaque sous-ensemble $S$ de
$\{ 1, \ldots, k \}$ ayant entre $1$ et $k/2$ éléments. \\
4. Soit $g_1 = \prod_{i \in S} w_i$ et $h_1 = \prod_{i \notin S} w_i$,
calculer $s, t$ tels que $s g_1 + t h_1 = 1 \bmod p$. \\
5. Calculer $g_N, h_N \in \Z[x]$ tels que $g_N \equiv g_1 \bmod p$,
$h_N \equiv h_1 \bmod p$ et $f \equiv g_N h_N \bmod p^N$. \\
\hspace*{1cm} Si $f = g_N h_N$ retourner $g_N, h_N$. \\
6. Retourner \og $f$ est irréductible\fg.
\end{quote}
À l'étape 4, si le degré de $g_N$ ne fait pas partie des degrés
possibles issus de l'analyse préliminaire, alors on peut passer
directement au sous-ensemble suivant.
La borne de Mignotte donnée par le théorème suivant garantit que
la valeur de $N$ choisie à l'étape 1 est assez grande pour garantir
l'exactitude du résultat \og $f$ est irréductible\fg\ en 6.
\begin{thm}[Mignotte 1989]
Soit $g = \sum g_i x^i \in \C[x]$ un facteur de degré $m$
de $f = \sum f_i x^i \in \C[x]$ de degré $n$. Alors
\[ \| g \| \le \left| \frac{g_m}{f_n} \right| 2^n \| f \| \]
où $\| f \|$ représente la norme euclidienne de $f$.
\end{thm}
L'algorithme de Zassenhaus est de complexité exponentielle dans le cas le
pire, par exemple avec les polynômes de Swinnerton-Dyer
\[ f_k = \prod (x \pm \sqrt{2} \pm \sqrt{3} \pm \cdots \pm \sqrt{p_k}) \]
où $p_k$ est le $k$e nombre premier. Le polynôme $f_k$ est à coefficients
entiers, irréductible sur $\Z$, mais se décompose en facteurs de degré
au plus deux sur $\F_p$ pour tout nombre premier $p$~!
Cependant Collins a montré que l'algorithme de Zassenhaus se comporte
\og en moyenne\fg\ de façon polynomiale; c'est encore aujourd'hui l'algorithme
utilisé en pratique par les systèmes de calcul formel.

\section{L'algorithme ECM de Lenstra pour la factorisation d'entiers}

\begin{df}
Soit $F$ un corps de caractéristique différente de $2$ et $3$ et
$x^3+ax+b$ un polynôme de $F[x]$ sans racines multiples. Alors
$E = \{ (u,v) \in F^2, v^2 = u^3 + a v + b \} \cup \{ \mathcal{O} \}$ est
une courbe elliptique sur $F$, où $\mathcal{O}$ est le \og point à
l'infini\fg\ sur $E$.
\end{df}
\begin{ex}
Le polynôme $p=x^3-x$ n'a pas de racines multiples si ${\rm car}\, F
\neq 2$, car ${\rm res}(p,p')=-4$. Donc $y^2 = x^3 - x$ définit une courbe
elliptique sur $\R$ ou $\C$.
\end{ex}
\begin{figure}
\centerline{\includegraphics[scale=.35]{ec}}
\caption{La courbe elliptique $y^2 = x^3 - x$}
\end{figure}
Le polynôme $p=x^3+ax+b$ n'a que des racines simples si le résultant
de $p$ et de $p'$ (i.e.~son discriminant) est non nul, ce qui revient à
dire que $4 a^3 + 27 b^2$ ne s'annule pas.
\begin{small}
\begin{verbatim}
>> p:=x^3+a*x+b: resultant(p,diff(p,x),x);

                                3       2
                               4 a  + 27 b
\end{verbatim}
\end{small}
L'équation $y^2 \!=\! x^3 \!+\! a x \!+\! b$ est l'équation de Weierstrass de la
courbe~$E$, et $a, b$ sont ses coefficients de Weierstrass.

Le point à l'infini se comprend mieux sur l'équation projective, où
on ajoute une troisième coordonnée $z=1$ pour les points \og standards\fg:
\[ y^2 z = x^3 + a x z^2 + b z^3 \]
et $x=0, y=1, z=0$ pour le point à l'infini $\mathcal{O}$.

\subsubsection*{Addition sur $E$}
On définit une structure de groupe sur une courbe elliptique de la
façon suivante:
\begin{df}
Soient $P_1$ et $P_2$ deux points sur une courbe ellip\-tique~$E$, et $R$ le
troisième point d'intersection de la droite $(P_1,P_2)$ avec $E$.
Alors si $R=(x,y)$, on définit $P_1 + P_2 = (x,-y)$.
Dans le cas particulier où $P_1=P_2$, on considère la tangente en $P_1$.
Si $P_1 = \mathcal{O}$, alors on prend la droite verticale passant par $P_2$.
\end{df}
Dans le cas $P_1 = \mathcal{O}$, $R$ est le symétrique de $P_2$ par rapport
à l'axe horizontal, donc $\mathcal{O} + P_2 = P_2$: $\mathcal{O}$ est l'élément
neutre du groupe.
L'opposé de $P=(x,y)$ pour l'addition dans $E$ est donc $(x,-y)$.

Les coodonnées $(x_3,y_3)$ de $P_1 + P_2$ en fonction de celles de
$P_1$ et de $P_2$ sont:
\[ x_3 = \left( \frac{y_2-y_1}{x_2-x_1} \right)^2-x_1-x_2, \qquad
y_3 = -y_1 + \frac{y_2-y_1}{x_2-x_1} (x_1 - x_3).
\]

\begin{ex}
Sur $F_5$, la courbe $y^2 = x^3 - 2$ a exactement six points
$a=\mathcal{O}$, $b=(3,0)$, $c=(2,1)$, $d=(1,2)$, $e=(1,3)$, $f=(2,4)$.
% et la table d'addition est la suivante:
% \[
% \begin{array}{|c|cccccc|} \hline
% + & a & b & c & d & e & f \\
% a & a & b & c & d & e & f \\
% b & b & b & e & f & c & d \\
% c & c & e & c & b & d & f \\
% \end{array} \]
\end{ex}
\begin{thm}[Hasse]
Le nombre $n$ de points d'une courbe elliptique sur $\F_q$ vérifie
\[ \left| n - (q+1) \right| \le 2 \sqrt{q}. \]
\end{thm}

Dans l'algorithme de Lenstra, on utilise des pseudo courbes elliptiques
dans $\Z/N\Z$ (qui n'est pas un corps).

\noindent
\begin{quote}\smaller
\noindent
{\bf Algorithme de {\sc Lenstra}.} \cite{Lenstra87} \\
Entrée: un entier composite $N$ non divisible par $2$ et $3$,
ni une puissante parfaite. \\
Sortie: un facteur non trivial de $N$ ou \og échec\fg. \\
1. Choisir des entiers positifs $B$ (borne sur les nombres premiers de la base)
et $C$ (borne sur un facteur de $N$). Soit
\[ k = \prod_{\substack{r \ {\rm premier} \\ r \le B}} r^{\alpha_r} \]
où $\alpha_r \in \N$ est déterminé par $r^{\alpha_r} \le
C + 2 \sqrt{C} + 1 < r^{\alpha_r+1}$. \\
2. Choisir au hasard une \og courbe elliptique\fg\ $E$ sur $\Z_N$ et un point $P$
sur $E$, et calculer $k P$, en utilisant les formules d'addition sur la
\og courbe elliptique\fg:
\[ k P = r \cdots r \cdots 3 \cdots 3 \cdot 2 \cdots 2 \cdot P \]
et en calculant de droite à gauche. Pour cela, on a besoin de faire des
divisions $u/v \bmod N$; si $g = \pgcd(v,N) = 1$, alors on peut
inverser $v$ et continuer, sinon on a trouvé un diviseur non trivial de $N$
que l'on retourne. \\
3. Si on arrive à calculer $k P = \mathcal{O}_E$, alors la factorisation est
sans succès.
\end{quote}

En notant $L(x) = e^{\sqrt{\log x \log\log x}}$, la valeur optimale
de $B$ est $B \simeq L(p)$ où $p$ est le plus petit facteur premier de $N$,
et le temps de calcul de l'algorithme de Lenstra est alors de l'ordre de
$L(p)^2$.

L'algorithme ECM de Lenstra n'est pas implanté dans MuPAD 1.3,
où par contre est implanté l'algorithme ECPP --- utilisant aussi les
courbes elliptiques --- de preuve de primalité dû à Atkin et Morain
\cite{AtMo90b} (fonction \verb|numlib::proveprime|).

\backmatter
\bibliographystyle{jepplain+eid}
\bibliography{xups97-03}
\end{document}