\documentclass{article} \usepackage[latin1]{inputenc} \usepackage[T1]{fontenc} \usepackage{amsmath,amsfonts,amssymb} \usepackage[francais]{babel} \usepackage{xspace} \usepackage{times} \newcommand{\bs}{\symbol{92}} \newcommand{\Z}{\mathbb{Z}} \newcommand{\Q}{\mathbb{Q}} \newcommand{\N}{\mathbb{N}} \newcommand{\R}{\mathbb{R}} \newcommand{\C}{\mathbb{C}} \newtheorem{theorem}{Theorem} \newtheorem{defn}{Definition} \begin{document} \title{Représentation des objets mathématiques et algorithmes} \author{Préparation agrégation, option C} \date{Révision du 27/01/2011} \maketitle Les objets de base à représenter sur machine sont d'abord les scalaires. On distingue les objets sur lesquels le microprocesseur opère directement dont la taille mémoire est fixée, ce sont les entiers courts et les flottants, et les objets dont la taille n'est pas fixée et qui sont gérés par des programmes, dont les entiers longs, les flottants multiprécision, les polynômes à une ou plusieurs variables. Les opérations arithmétiques (addition, multiplication, division) sur les premiers sont rapides et en temps constant, ce qui n'est pas le cas des autres, le temps dépendant de la taille des objets, on utilise des algorithmes plus ou moins sophistiqués mathématiquement parlant pour réaliser les opérations arithmétiques de base le plus rapidement possible, dont voici quelques exemples. \section{Les entiers} On distingue les entiers machines (codés sur 32 ou 64 bits, signés ou non) et les entiers en précision arbitraire (dont les opérations arithmétique de base ne sont pas effectuées par le microprocesseur mais doivent être programmées). Pour représenter un entier en précision arbitraire, on utilise par exemple un mot de 32 bits pour représenter le signe et la taille de l'entier, puis l'écriture en base $2^{32}$ de la valeurs absolue de l'entier. La plupart des implémentations utilisent des algorithmes efficaces pour~: \begin{itemize} \item la multiplication~: Karatsuba, Toom-3, FFT, dont certains sont esquissés plus loin pour les polynômes (l'adaptation de ces algorithmes aux entiers nécessitant de tenir compte des retenues) \item la division~: la m\'ethode ``de base'' consiste \`a poser la division et n\'ecessite un temps proportionnel \`a la taille du diviseur par la taille du quotient. Il existe aussi des am\'eliorations. On dit qu'un nombre \'ecrit en base $\beta^N$ est normalis\'e si son \'ecriture commence par un digit $B\geq \beta^N/2$. Si on divise un nombre $A$ de taille $2N$ en base $\beta$ par un nombre $B$ normalis\'e de taille $N$, avec $N$ pair, l'algorithme ``diviser pour r\'egner'' consiste \`a faire une division de base dont les ``chiffres'' sont des entiers longs de taille $N/2$ en base $\beta$ (on divise alors un nombre de 4 chiffres \'ecrit en base $\beta^{N/2}$ par un nombre de 2 chiffres), les multiplications que l'on effectuera entre les ``chiffres'' du quotient et du diviseur peuvent alors utiliser un algorithme de multiplication rapide. Pour diviser un nombre $a$ de 4 chiffres par un nombre $b$ de 2 chiffres, on calcule le premier chiffre du quotient en utilisant les 2 premiers chiffres de $a$ et le premier chiffre de $b$ (en appelant r\'ecursivement l'algorithme de division rapide), on obtient un candidat quotient $q$ que l'on d\'ecale et multiplie par $b$ et retranche de $a$, on ajuste alors $q$ en fonction du signe du reste (on diminue $q$ si le reste est n\'egatif). Si $b$ est normalis\'e, donc au moins \'egal \`a la moiti\'e du plus petit nombre \'ecrit sur 3 ``chiffres'' (par exemple en base 10, si $b$ est sup\'erieur \`a 50), on montre qu'on doit diminuer $q$ de 0, 1 ou 2. On recommence pour calculer le deuxi\`eme chiffre du quotient. Par exemple pour calculer le quotient de 4513 par 68, on commence par diviser 45 par 6, il y va 7 fois, 7 fois 68 font 476, donc le reste est 451-476=-25, donc on diminue le premier chiffre du quotient de 7 \`a 6, on ajoute 68 au reste qui vaut 43 il reste \`a diviser 433 par 68, 43 divis\'e par 6 il y va 7 fois, puis 433 moins 7 fois 68 font -43, on diminue donc le deuxi\`eme chiffre du quotient (le quotient final est 66), et on ajoute 68 au reste qui devient 25. Les op\'erations longues sont les 2 multiplications de 1 chiffre par 2 chiffres. Si le diviseur n'est pas normalis\'e, on multiplie les 2 nombres \`a diviser par une puissance de 2 pour qu'il soit normalis\'e. Il existe aussi une m\'ethode th\'eorique de complexit\'e meilleure appel\'ee me\'ethode de Newton (cf. par exemple Gerhard et Von zur Gathen). \item l'extraction de racine carrée~: on peut distinguer deux algorithmes, savoir si le nombre est un carré parfait, extraire la partie entière de la racine carrée. Pour le premier, la documentation de GMP indique `mpz\_perfect\_square\_p' is able to quickly exclude most non-squares by checking whether the input is a quadratic residue modulo some small integers. The first test is modulo 256 which means simply examining the least significant byte. Only 44 different values occur as the low byte of a square, so 82.8\% of non-squares can be immediately excluded. Similar tests modulo primes from 3 to 29 exclude 99.5\% of those remaining, or if a limb is 64 bits then primes up to 53 are used, excluding 99.99\%. A single Nx1 remainder using `PP' from `gmp-impl.h' quickly gives all these remainders. Pour le deuxi\`eme, il existe aussi des algorithmes r\'ecursifs utilisant la possibilit\'e de multiplier rapidement deux entiers longs de m\^eme taille. \end{itemize} \section{Les flottants} On va d'abord d\'ecrire les flottants double pr\'ecision (double). La représentation d'un double en mémoire se compose de 3 parties~: le bit de signe $s=\pm 1$ sur 1 bit, la mantisse $M \in [0,2^{52}[$ sur 52 bits, et l'exposant $e \in [0, 2^{11}[$ sur 11 bits. Pour les nombres ``normaux'', l'exposant est en fait compris entre 1 et $2^{11}-2$, le nombre repr\'esent\'e est le rationnel \[ (1+\frac{M}{2^{52}}) 2^{e+1-2^{10}} \] L'exposant 0 est r\'eserv\'e aux nombres d\'enormalis\'es de la forme $M/2^{52} 2^{2-2^{10}}$, qui permettent de représenter des nombres plus petits que $2^{2-2^{10}}$ sous forme d\'enormalis\'ee, pour \'eviter que l'on puisse avoir simultan\'ement $x\neq y$ et $x-y=0.0$), l'exposant $2^{11}$ sert \`a coder les d\'epassements de capacit\'e (+ et - l'infini ainsi que NaN, not a number, obtenu par exemple par division de 0 par 0) Exemples~: \begin{itemize} \item -2 \\ Signe n\'egatif. Il faut diviser sa valeur absolue 2 par $2^1$ pour \^etre entre 1 et 2 dont $e+1-2^{10}=1$, l'exposant est $e=2^{10}$. On a alors $r=1$, $r-1=0$. Repr\'esentation \\ \verb|1 10000000000 00000000...0000| \item 1.5=3/2\\ Signe positif, compris entre 1 et 2 dont l'exposant v\'erifie $e+1-2^{10}=0$ soit $e=2^{10}-1=2^9+2^8+2^7+2^6+2^5+2^4+2^3+2^2+2^1+2^0$. On a $r-1=1/2=2^{-1}$. D'o\`u la repr\'esentation\\ \verb|0 01111111111 10000000...0000| \end{itemize} Lorsque le d\'eveloppement en base 2 est infini, on adopte une r\`egle d'arrondi du dernier bit au nombre sup\'erieur si le bit suivant est un 1. L'arrondi au plus proche entraine une erreur relative de repr\'esentation de $2^{-53}$. Les op\'erations sur les flottants sont tr\`es rapides, car impl\'ement\'ees directement par les micro-processeurs mais ne sont pas exactes. Ces erreurs se propagent selon la règle approximative~: l'erreur absolue d'une somme/différence est inférieure ou égale à la somme des erreurs absolues sur les arguments, l'erreur relative sur un produit ou un quotient est inférieure ou égale à la somme des erreurs relatives sur les arguments (en toute rigueur, il faut y ajouter l'erreur relative d'arrondi et par exemple le produit des erreurs relatives pour le produit). On peut diminuer les erreurs d'arrondis en utilisant des flottants multipr\'ecision, dont le principe est identique aux flottants machines, mais dont la mantisse est cod\'ee sur un nombre param\'etrable de bits. Les op\'erations sont alors plus lentes, car elles n\'ecessitent de faire des op\'erations sur des entiers (correspondant \`a la mantisse) de taille arbitraire. Certaines implémentations permettent de spécifier l'arrondi à effectuer pour chaque opération, ce qui permet d'implémenter une arithmétique des intervalles et donc d'obtenir un encadrement rigoureux du résultat obtenu. Mais on ne peut pas empêcher les problèmes des opérations en flottant comme par exemple~: \begin{itemize} \item la non associativité de l'addition (par exemple $(2^{-53}+1)-1 \neq 2^{-53}$ dans les flottants) \item le fait que par exemple $0.3-3 \times 0.1$ n'est pas nul si les flottants sont représentés en base 2 (certains microprocesseurs peuvent travailler en base 10, on parle de BCD, certains logiciels comme par exemple maple, utilisent par défaut des flottants codés en base 10). \end{itemize} On peut d'ailleurs utiliser ces arrondis pour d\'eterminer la base de repr\'esentation utilis\'ee par le programme suivant (\'ecrit en syntaxe maple) \begin{verbatim} A:=1.0; B:=1.0; while ((A+1.0)-A)-1.0=0.0 do A:=2*A; od; while ((A+B)-A)-B <> 0.0 do B:=B+1.0; od; B; \end{verbatim} \section{Les polynômes à une variable.} On peut représenter les polynômes à une variable par une liste de coefficients, classés par ordre croissant ou décroissant de degré. Cette représentation est appelée dense. La représentation creuse sert pour les polynômes ayant peu de coefficients non nuls (par exemple $x^{100}-1$). On suppose ici que les opérations de base sur les coefficients sont en temps constant (par exemple les coefficients sont dans un corps fini ou dans les flottants). L'opération d'addition de deux polynômes est proportionnelle en temps au nombre de coefficients représentés. Pour la multiplication de deux polynômes creux, le temps de calcul est proportionnel au produit du nombre de coefficients représentés. Pour les polynomes denses, il existe plusieurs méthodes permettant d'améliorer les performances. \subsection{L'algorithme de Karatsuba} Si $P$ et $Q$ sont de degrés $< 2n$, on peut écrire \[ P=A+Bx^n, Q=C+Dx^n \] avec $A,B,C,D$ de degré $1$ arbitraire (mais la constante devant $n^r$ est d'autant plus grande que $r$ est petit). \subsection{La FFT} Soient $P$ et $Q$ deux polynômes dont le degré du produit $PQ$ est inférieur strict à $N$ (par exemple $P$ et $Q$ de degrés strictement inférieur à $N/2$). Pour le calculer il suffit de connaitre sa valeur en $N$ points distincts $\xi_k$, ce qui peut se faire en faisant le produit terme à terme $P(\xi_k)Q(\xi_k)$. Pour que cette méthode soit intéressante, il faut pouvoir passer rapidement de la représentation polynomiale habituelle $P$ à la représentation de la liste des $P(\xi_k)$ et inversement. Ceci est possible lorsque $N$ est une puissance de 2, et qu'il existe des $\xi_k=\omega^k$ racines $N$-ième de l'unité (par exemple si les coefficients sont complexes flottants mais aussi dans certains corps finis), il s'agit de la FFT (Fast Fourier Transform) qui nécessite alors $N \ln(N)$ opérations. Soit donc $P$ un polyn\^ome de degr\'e strictement inf\'erieur \`a $N$~: \[ P(X)=\sum_{j=0}^{N-1} p_j X^j \] que l'on souhaite \'evaluer aux $N$ points $1$, $\omega$, ..., $\omega^{N-1}$. Si $N=2M$, on d\'ecoupe $P$ en 2 parties de m\^eme longueur par division euclidienne par $X^M$~: \[ P(X)=X^M Q(X) + R(X) \] on a alors \[ P(\omega^{2k})= (Q+R) ((\omega^2)^{k}), \quad P(\omega^{2k+1})= (-Q +R)_\omega ((\omega^2)^k) \quad S_\omega(X)=\sum s_k \omega^k X^k \] On est donc ramen\'e \`a deux additions de 2 polyn\^omes de degr\'e $M$, une multiplication coefficient par puissances ($4M$ op\'erations), et au calcul des transform\'ees de Fourier discr\`ete de $Q+R$ et $R-Q$. Si $N=2^n$, ces transform\'ees peuvent encore se calculer par division successives par 2, on v\'erifie que cela n\'ecessite $O(n2^n)$ op\'erations, donc $N\ln(N)$ opérations. \section{Les polynômes à plusieurs variables} Le problème principal des polynômes à plusieurs variables est d'avoir beaucoup de coefficients même pour des degré totaux petits, ainsi le nombre de monomes possible en degré total $d$ et nombre de variables $n$ est \verb|comb(d+n,n)|, en degr\'es partiels il est major\'e par $\prod_{i=1}^n (d_i+1)$ o\`u $d_i$ d\'esigne le degr\'e par rapport \`a la $i$-i\`eme variable. Il est donc souvent impraticable de travailler avec des polynomes denses. On utilise alors des représentations creuses, par exemple une table de coefficients indiciées par des monomes (des n-uplet d'entiers représentant les degrés partiels par rapport aux variables du polynôme). Pour faire le produit de deux polynômes, il faut donc faire le produit terme à terme et additionner les monomes obtenus de même degré. Pour accélerer le temps de calcul d'un produit de polynômes, on peut utiliser des astuces~: \begin{itemize} \item Si les polyn\^omes sont denses par rapport \`a chaque variable~:\\ On peut se ramener \`a une seule variable, par exemple pour deux polyn\^omes $P(x,y), Q(x,y)$ en deux variables $x$ et $y$ dont les degr\'es partiels sont $n_1,m_1$ par rapport \`a $x$ et $n_2,m_2$ par rapport \`a $y$, on pose $x=y^{m_1+m_2+1}$, et on multiplie $P'(y)=P(y^{m_1+m_2+1},y)$ par $Q'(y)=Q(y^{m_1+m_2+1},y)$ puis on d\'etermine $PQ(x,y)$ ``en \'ecrivant le produit $P'Q'$ en base $y^{m_1+m_2+1}$''. On peut alors utiliser l'une des m\'ethodes de multiplication rapide d\'ecrites pr\'ec\'edemment. Il y a toutefois un inconvenient majeur \`a utiliser cette m\'ethode, c'est que $P'$ et $Q'$ ont beaucoup de coefficients nuls, ils sont creux si on essaye de g\'en\'eraliser \`a une dimension plus grande que 2. Il est plus int\'eressant d\`es les petites dimensions d'interpoler le produit des deux polyn\^omes. On donne \`a l'une des variables $d_i+1$ valeurs ($d_i$ d\'esignant le degr\'e partiel par rapport \`a cette variable), on fait (r\'ecursivement) le produit des polyn\^omes \`a une variable de moins et on applique l'algorithme des diff\'erences divis\'ees. \item Si on multiplie terme \`a terme~:\\ Il faut un conteneur pour les monomes du polyn\^ome produit dans lequel on peut rapidement d\'eterminer si un monome identique existe d\'ej\`a (dans ce cas on ajoute le produit des coefficients des termes que l'on multiplie au coefficient du monome existant) ou n'existe pas (dans ce cas on ajoute ce monome, coefficient au conteneur). On peut utiliser des arbres binaires (tables dont l'indice est le monome et la valeur est le coefficient) ou des tables de hachage, qui permettent d'accéder plus rapidement à un élément de la table du polynôme produit (si la fonction de hachage est adapt\'ee). \item Dans les deux cas, en majorant chaque degré partiel du produit, on peut utiliser un entier au lieu d'un n-uplet d'entiers pour représenter un monome, ceci facilite les tests de comparaison entre monomes (l'ordre lexicographique devenant l'ordre habituel sur les entiers). Pour cela, on g\'en\'eralise l'\'ecriture en base $b$ d'un entier, en calculant les restes successifs par les $d_i+1$ ($d_i$ d\'esignant le degr\'e partiel du produit par rapport \`a la $i$-\`eme variable). \end{itemize} \section{Les corps finis} Ils sont isomorphes à $G=GF(p,n)$ où $p$ est leur caractéristique et $p^n$ leur cardinal, $GF(p,n)$ étant le quotient de $\Z/pZ[X]$ par un polynôme minimal irréductible de degré $n$. On peut donc représenter un élément de $G$ par un polynôme et effecter les opérations arithmétiques sur $G$ en faisant les opérations sur les polynômes suivi par une division euclidienne par le polynôme minimal. Si le polynôme minimal est primitif, alors $G$ est composé de 0 et des quotients des puissances de $x$ pour $i\in[0,p^n-2]$, sinon on cherche un élément primitif $a$. On peut alors représenter les éléments de $G$ par un entier, qu'on choisit égal à -1 pour 0 et compris entre 0 et $p^n-2$ pour les puissances correspondantes de $a$. La multiplication et l'inversion dans cette représentation est immédiate, on construit ensuite une table d'addition et d'opposés, ce qui permet d'effectuer les opérations de base sur le corps de manière optimale en temps. \section{Idées de développement} \begin{itemize} \item Justifier la complexité de l'algorithme de Karatsuba (ou/et des variantes). \item Illustrer la multiplication de polynômes à coefficients réels par FFT. \item Chercher un corps fini $\Z/q\Z$ permettant de faire la multiplication par FFT de 2 polyn\^omes de degr\'e $m=2^n$ \`a coefficients entiers de valeur absolue plus petite que $C$. Indication~: on cherchera $q$ premier sous la forme $q=t2^m+1, \quad t \in \N$ (existence garantie par le th\'eor\'eme de Dirichlet), tel que $q/2$ soit plus grand que le plus grand coefficient du produit, et on cherchera $\omega$ une racine $2^n$-i\`eme primitive de 1 dans $\Z/q\Z$ ($\omega=\beta^t$). \item Flottants et erreurs d'arrondi. \item Opérations sur les entiers multiprécision. \item Discuter la multiplication des polynômes à plusieurs variables, selon qu'ils sont denses ou creux. \item Opérations sur les corps finis. Comment trouver un polyn\^ome irr\'eductible, voire primitif. \end{itemize} \end{document}