Si vous vous êtes frotté au calcul numérique, c'est à dire le calcul, par des méthodes informatiques, des valeurs approchées des solutions d'un problème trop complexe ou trop volumineux, alors peut-être avez-vous remarqué que la précision interne des nombres à virgule flottante est parfois insuffisante. Pour ne citer que quelques exemples: manipuler des matrices de grande taille, afficher des fractals avec un facteur de zoom extrêmement élevé, utiliser des algorithmes dans lesquels les valeurs intermédiaires peuvent être très grandes ou nécessitent une précision accrue (des algorithmes instables, très sensibles aux conditions initiales voire carrément chaotiques), calculer des orbites planétaires fiables sur plusieurs millions d'années, etc...

Il est alors nécessaire de réaliser les calculs numériques avec des représentations de nombres en grande précision. Lorsque cette précision n'est pas connue à priori ni même majorée au moment de l'écriture du code, on doit alors programmer des algorithmes fonctionnant avec une précision quelconque (dans la limite de la mémoire disponible biensûr). On parle alors de calculs en précision arbitraire. Il existe des librairies dites en précision accrue capables de manipuler des représentations "double-double" ou "quad-doubles" constituées respectivement de deux ou quatre nombres à virgule flottante de type double precision et fournissant une précision d'environ respectivement 30 ou 60 décimales. Mais ces librairies ne permettent pas d'aller plus loin (en tout cas, pas facilement).

Il existe bien entendu plusieurs librairies de manipulation des nombres en précision arbitraire, la plus connue étant probablement GMP (pour GNU Multi-Precision). Cette librairie est portable sur une multitude de plateformes et très performante, mais elle est aussi assez complexe dans son implémentation et son utilisation. C'est le prix à payer pour ses performances et sa portabilité (comme bien souvent). Elle convient donc aux développeurs qui ont besoin de l'outil sans se préoccuper de son fonctionnement "interne".

En revanche, ici nous sommes intéressés par le fonctionnement interne. Vous allez voir que s'il est simple de programmer des additions, des soustractions et des multiplications sur des représentations de quelques centaines de décimales, il est autrement plus complexe de programmer la division, la racine carrée et les fonctions transcendantes (exponentielle, logarithme, sinus, cosinus, etc...) sur ces mêmes représentations. Il est encore bien plus complexe de gérer des précisions de plusieurs milliers (voire centaines de milliers) de décimales avec une bonne efficacité.

C'est pourquoi nous allons nous fixer l'objectif de programmer une librairie capable de manipuler des nombres jusqu'à 1 million de décimales (et même plus...) et d'effectuer toutes les opérations et les fonctions habituelles avec des performances "proches" des meilleures implémentations. Ne croyez pas que nous allons faire "mieux" que GMP. Nous n'allons programmer notre librairie que pour processeurs x86 (en mode 64 bits) et sous Windows (un portage linux est possible sans trop de difficultés mais nous ne le ferons pas ici). Malgré cette restriction notre librairie restera plus lente que GMP. Mais nous serons entre \(1.2\) et \(1.8\) fois plus lents avec (probablement) \(100\) fois moins de code que GMP.

Usage de l'assembleur

Comme la portabilité n'est pas d'une importance capitale et que nous voulons de bonne performances, l'usage de l'assembleur est tout indiqué. Cependant, que les possesseurs d'anciens microprocesseurs ne s'irritent pas, je fournirai l'équivalent C au code assembleur essentiellement SSE4. Des routines AVX sont en cours d'écriture et, au moment où j'écris ces lignes, ne sont pas encore "opérationnelles".

Sous Visual Studio

Sous Visual Studio 2015, la première chose à savoir, c'est qu'en mode 64 bits (celui qui nous intéresse), il n'est plus possible de "mixer" directement du C++ avec de l'assembleur (contrairement au mode 32 bits avec la directive _asm). Il faut donc impérativement passer par un assembleur externe. Heureusement, Microsoft fournit son assembleur (MASM) avec la version "community" de Visual Studio (ce n'est pas le cas avec toutes les version de l'IDE). Il est biensûr possible d'utiliser YASM, NASM, FASM ou d'autres assembleurs (par exemple, le site de YASM propose un plugin Visual Studio).

Dans le menu de l'IDE, après avoir sélectionné le projet, allez dans le menu "Projet" puis "personnalisations de la build...". Là, vous devriez voir une ligne "MASM" . Cochez la case correspondante puis "OK". En théorie, les fichiers assembleur devraient maintenant "compiler" pourvu qu'on aille dans leurs propriétés et que dans "type d'élément" on choisisse "Microsoft Macro Assembler".

Sous Visual Studio Code

Avec le plugin CMAKE, il est trés simple d'intégrer de l'assembleur. Dans la partie 'LANGUAGES' de la directive 'project', il suffit de mettre: 'ASM-MASM' (ou ASM-NASM, ou autre cf. la documentation de CMAKE). Ensuite, vous pouvez mettre des fichiers .asm dans vos fichiers source.

Conventions sous Windows 64 bits

Quelques règles sont à respecter pour correctement interfacer de l'assembleur x64 avec du C sous Windows 64 bits (n'oubliez pas que ces règles sont différentes sous Linux).

Premièrement, le passage de paramètres entiers et pointeurs se fait par les registres RCX, RDX, R8 et R9. C'est à dire que, de gauche à droite, le premier paramètre entier ou pointeur est placé dans RCX, le second dans RDX, le troisième dans R8 et la quatrième dans R9. S'il y a d'autres paramètres entiers ou pointeurs, il sont poussés sur la pile dans l'ordre de droite à gauche. Il est important de noter que, bien que passés par des registres, la place nécessaire à ces paramètres DOIT tout de même être réservée en pile, afin de permettre à l'appelé de les sauver en cas de besoin. De plus, les bits inutilisés lors du passage de paramètre de taille inférieure à 8 octets ne sont PAS mis à zéro. Si, par exemple, le premier paramètre est un entier 32 bits, il sera placé dans les 32 bits de poids faible de RCX et les 32 bits de poids fort ne seront pas nécessairement mis à zero, bref, ça sera du flan.

Deuxièmement, le passage de paramètres à virgule flottante se fait via les registres XMM0, XMM1, XMM2 et XMM3, dans cet ordre. Comme pour les entiers, si plus de quatres paramètres à virgule flottante doivent être passés, les autres sont poussés sur la pile. L'espace nécessaire aux premiers paramètres doit également être réservé en pile. A noter que, bien que les registres XMM puissent contenir deux flottants, dans le cas du passage de paramètre, un seul flottant est stocké dans chaque registre XMM (en partie basse).

Les paramètres de type structure ou union dont la taille ne rentre pas dans 64 bits sont passés par pointeur. Il est de la responsabilité de l'appelé de ne pas y écrire si la sémantique de la fonction le spécifie.

C'est à l'appelant de réserver puis de libérer l'espace en pile. De plus, même si la fonction appelée a moins de quatre paramètres, un espace de 32 octets DOIT être réservé.

La pile doit TOUJOURS être alignée à 16 octets. Donc, lors d'un appel de fonction, l'adresse de retour (8 octets) étant poussée par le CALL, il faudra "ajuster" le pointeur de pile (en faisant un nombre impair de PUSH, ou en réservant 16n+8 octets par un SUB RSP,16n+8) avant tout appel à une sous fonction.

Un exemple étant plus parlant qu'un long discours, une dernière subtilité est à prendre en compte. Ainsi, si la fonction appelée a un entier en premier paramètre, un flottant en deuxième et troisième puis un entier en quatrième, ils seront placés respectivement dans RCX, XMM1, XMM2 et R9 (les registres RDX, R8 et XMM0 ne sont pas utilisés).

Enfin, dans une fonction assembleur, les registres RAX, RCX, RDX, R8, R9, R10 et R11 peuvent être modifiés sans devoir être sauvegardés auparavant. En revanche, les registres RBX, RDI, RSI, R12, R13, R14 et R15 sont considérés comme non modifiés par le compilateur. Ainsi, s'ils sont utilisés dans la fonction, ils doivent être sauvegardés puis restaurés avant le retour. De la même manière, les registres XMM0, XMM1, XMM2, XMM3, XMM4 et XMM5 peuvent être détruits, mais XMM6 à XMM15 doivent être conservés.

Bien entendu, tout manquement à l'une de ces règles se soldera la plupart du temps par un crash, et parfois par une corruption de données. De plus, je n'ai pas énoncé toutes les règles. Dans le cas de retour de type structure ou union, de gestion d'exception ou du mécanisme de "stack unwinding", il faudra se référer à une documentation plus précise. Nous n'en aurons pas besoin ici.

Si vous n'avez pas tout compris aux quelques lignes précédentes, ne vous faites pas trop de soucis, en voyant le code cela paraîtra bien plus simple.

struct st_mpInt
{
	long long				sgn;
	long long				size;
	unsigned long long		m[0];
};
typedef struct st_mpInt MpInt;

///////////////////////////////////////

struct st_mpFloat
{
	long long				sgn;
	long long				size;
	long long				exp;
	long long				prec;
	unsigned long long		m[0];
};
typedef struct st_mpFloat MpFloat;