Citer : Posté le 21/06/2021 22:17 | #

---

Calcul en point fixe et en point flottant (Programmation générale — ★☆☆)

Format des types en point flottant

En plus des types entiers (int et compagnie), le langage C spécifie deux types de nombres «à virgule» : float et double. Le premier est nommé comme ça parce que c'est un nombre à «virgule flottante», et le second parce que c'est un float avec double précision.

Quoi qu'on fasse, pour un ordinateur un nombre c'est une suite de chiffres, du style 37492. Si le nombre était entier ça s'arrêterait là, mais pour un nombre décimal on veut pouvoir indiquer la position de la virgule. Le principe du nombre en point flottant c'est de stocker, en plus de la suite de chiffres, la position de la virgule (sous la forme d'une puissance de 2, ou pour l'exemple ci-dessous, de 10).

Par exemple (conceptuellement), le nombre en point flottant «valeur=37492, puissance=10²» va représenter le nombre décimal 37492×10² = 3749200. Le nombre «valeur=37492, puissance=10⁻³» représenterait 37492×10⁻³ = 37.492 qui a 3 chiffres après la virgule. Généralement on peut mettre des puissances 10 assez fines, imaginez-vous 10⁻³⁸ pour une float ou 10⁻³⁰⁰ pour un double, et pareil dans les puissances positives. La plage est immense, on peut vraiment tout représenter.

À cela, on ajoute un bit pour spécifier si le nombre est positif ou négatif (le complément à 2 ne marchant pas vraiment), et donc on obtient un nombre qui a 3 champs :

• Le signe, qui indique si le nombre est positif ou négatif (1 bit).
• L'exposant, qui indique la puissance de 10 servant de base au nombre. (Exemple : -5 pour 2⁻⁵)
• La mantisse, qui est la séquence de chiffres. (Exemple : 37492)

Tout ça est compacté ensemble dans 32 bits (pour un float) ou 64 bits (pour un double) selon le format fixé par la norme IEEE 754. (Il existait d'autres formats avant, et c'était un enfer, mais aujourd'hui c'est fini et il n'existe plus que l'IEEE 754, ce qui est vraiment vraiment bien.)

La complexité de calculer en point flottant

Imaginons que vous voulez additionner deux nombres. Pour des entiers, c'est vraiment facile : il suffit de les jeter dans un circuit d'addition, qui fait la même chose qu'un humain qui pose l'addition au papier. Les entiers négatifs sont encodés d'une façon particulière (complément à 2) conçue pour qu'ajouter un nombre négatif ne demande aucune action particulière, le circuit est identique. En assembleur, ça prend 1 seule instruction qui occupe l'ALU (unité de calcul arithmétique) pendant 1 seul cycle d'horloge, on ne peut pas faire plus simple.

Pour des nombres en point flottant, c'est beaucoup plus compliqué.

• D'abord s'il y a un nombre négatif il faut faire une vraie soustraction qui est une opération différente ; il n'y a pas de complément à 2 donc il faut séparer les deux cas.
• Ensuite, les deux nombres n'ont pas forcément le même exposant, et donc il faut les «aligner». Par exemple, si on veut poser l'addition de 478×10² et 3921×10⁻¹ au papier, on les aligne pour que le 8 (qui représente 800) soit en face du 3 (qui représente 300). Ici, c'est pareil, sauf qu'électroniquement c'est assez casse-pieds, d'autant que si on a 478×10¹²⁷ et 3921×10⁻⁷⁵ on n'a pas envie de calculer une addition sur des nombres de 200 chiffres.
• Une fois qu'on a aligné les opérandes, et calculé l'addition, il faut encore trouver la puissance de 10 correspondant au résultat (on ne peut pas prendre n'importe laquelle, il y a des restrictions).
• Enfin, les nombres en points flottant ont des comportements un peu différents des entiers : si on dépasse la valeur maximale représentable le résultat doit être une valeur spéciale indiquant «l'infini» ; si l'opération est invalide il y a des exceptions à lever ; il y a des cas très particuliers avec des nombres «dénormalisés» qui ne font rien comme les autres ; et j'en passe.

Vous ne serez donc pas surpris quand je dirai qu'électroniquement, c'est beaucoup plus compliqué, et ça ne prend pas qu'1 cycle. C'est assez compliqué pour qu'un coprocesseur complet y soit dédié d'ailleurs ; on appelle ça un FPU (pour Floating-Point Unit ou unité de calcul en point flottant).

Grâce à des techniques électroniques (pipelining), un processeur moderne avec un FPU à qui on donne un gros paquet d'additions à faire va réussir à en compléter une par cycle en moyenne (il triche, il en exécute plusieurs en même temps). Donc on arrive à s'en sortir, mais c'est plus difficile de réaliser cette moyenne qu'avec des entiers.

Point flottant logiciel sur la calculatrice

Sur la calculatrice, vous serez peut-être surpris de savoir qu'il n'y a pas de FPU. Le processeur est un SuperH 4AL-DSP, qui est variante du SuperH 4A où le FPU est remplacé par un DSP (un autre type de coprocesseur plutôt orienté vers le traitement du signal). Pour une calculatrice, ça peut paraître surprenant étant donné la grande quantité de calcul, mais c'est surtout parce que CASIO préfère utiliser des nombres en point flottant en base 10, pour éviter les problèmes d'approximation de la base 2 (c'est un sujet à part entière, je ne détaille pas). Notez juste qu'en Basic CASIO toutes les variables sont des nombres en point flottant en base 10, sauf en mode BASE où c'est des entiers en binaire.

Je pense que vous voyez venir la tragédie. Toute la complexité des nombres en point flottant ne peut pas être cachée dans le circuit électronique d'un FPU s'il n'y a pas de FPU. Et donc tout doit être fait en logiciel, et on parle de point flottant logiciel (ou soft-fp pour les anglophones).

Voilà le code assembleur pour ajouter deux entiers (stockés dans les registres r4 et r5).


Et voilà le code assembleur pour ajouter deux float (pris dans _addsub_fs.o de libgcc ; vous pouvez ignorer les nombres en début de ligne).


Vous voyez la complexité ? Non ? Voilà le code pour ajouter deux double (pris dans _addsub_df.o de libgcc).


Ajouter deux entiers, ici, prend environ 8 ns (si on ignore le parallélisme bas-niveau qui est un sujet beaucoup plus compliqué). Ajouter deux float quelconques m'a pris autour de 2-3 µs (~300 fois plus), et ajouter deux double autour de 3-4 µs (~400 fois plus).

Et là ce n'est qu'une addition, ça ne s'améliore pas quand il s'agit de multiplier (ordre d'idée en double : dans les 5 µs), et vous pouvez oublier la division si vous tenez aux performances (ordre d'idée en double : dans les 10-11 µs). Diviser est lent aussi pour les entiers (et donne lieu à d'autres optimisations), mais ça n'a rien à voir avec ce niveau-là.

Soit dit en passant, en plus de prendre du temps ce code prend aussi de la place, et la lib mathématique en point flottant logiciel (OpenLibm pour nous) prend aussi beaucoup de place. Chaque fonction occupe vite 4 kio de code sinon plus, donc les add-ins grossissent vite (et puis ça ruine un peu le cache).

Les nombres en point fixe

Les nombres en point fixe sont une autre option pour représenter des nombres décimaux. Comme son nom l'indique, la virgule est fixée à une certaine position, ce qui élimine l'exposant. Par exemple, pour un point fixe «de puissance 10⁻³», on stockera la valeur 2874 pour représenter le nombre décimal 2.874. La valeur de l'exposant, -3, n'est pas représentée dans la mémoire mais est implicite.

L'avantage c'est qu'éliminer l'exposant signifie qu'on perd un champ sur les 3, et du coup on peut récupérer le complément à 2 et éliminer le bit de signe. Ultimement, le format complexe à trois champs redevient un entier.

Les nombres en point fixe sont vraiment aussi simples qu'écrire 1383 centimes dans un entier pour signifier 13.83 €. Il n'y a pas d'arnaque, on décide simplement qu'on utilise la valeur entière X pour représenter le nombre décimal X/100 (ou X/1000, ou autre puissance de 10), et ça marche exactement pareil avec les nombres négatifs.

Notez quand même que les puissances de 10 ça ne marche pas très bien sur un processeur qui compte en binaire, on préfère utiliser des puissances de 2. Du coup, dans la suite de cette partie, je vais utiliser du binaire.

Les entiers du processeur SuperH font 32 bits, et on a le choix sur la position de la virgule. Ici je vais prendre peut-ête le format le plus évident, où je mets la virgule au milieu, ce qui me laisse 16 chiffres pour la partie entière et 16 pour la partie décimale. On dit que j'utilise le «format de point fixe 16:16».

Comme il y a 16 bits de partie décimale, l'entier x représentera la valeur décimale x/2¹⁶, c'est-à-dire x/65536. Les entiers sur 32 bits vont de -2147483648 à 2147483647, et donc le format en point fixe va aller de -2147483648/2¹⁶ = -32768, à 2147483647/2¹⁶ = 32767.99998.

Comme on va le voir, les opérations sont pas tout à fait les mêmes que les entiers, mais ça se passe assez bien.

Opérations en point fixe

Pour faire les opérations en point fixe sans se planter, je vous conseille de poser les étapes suivantes :

1. Écrivez l'opération mathématique que vous voulez sur les nombres décimaux (z = ... un truc en fonction de x et y ...).
2. Remplacez les nombres décimaux (x, y, z) par leur représentation en point fixe (2¹⁶x, 2¹⁶y, 2¹⁶z).
3. Ensuite seulement, cherchez une séquence d'opération sur les entiers qui donne le résultat à partir des entrées.

Par exemple, pour additionner. Si j'ai deux nombre décimaux x et y, alors :

1. Je veux obtenir le nombre décimal x + y à partir de x et y.
2. Comme je calcule en point fixe 16:16, je veux donc obtenir l'entier 2¹⁶(x+y) à partir de 2¹⁶x et 2¹⁶y.
3. En observant que 2¹⁶(x+y) = 2¹⁶x + 2¹⁶y, je réalise mon opération ; en C, j'écris z = x + y.

La multiplication montre un exemple où il ne suffit pas d'utiliser l'opération équivalente sur les entiers :

1. Je veux obtenir le nombre décimal x · y à partir de x et y.
2. Comme je calcule en point fixe 16:16, je veux donc obtenir 2¹⁶(x · y) à partir de 2¹⁶x et 2¹⁶y.
3. J'observe que 2¹⁶(x · y) = (2¹⁶x · 2¹⁶y) / 2¹⁶, et donc en C j'écris z = (x * y) >> 16.

On remarque aussi que la multiplication de 2¹⁶x et 2¹⁶y est un nombre entier de 64 bits dont on récupère les 32 bits du milieu, donc il faut bien faire la multiplication 64 bits. On écrit donc, en réalité, z = ((int64_t)x * (int64_t)y) >> 16.

(Le compilateur est assez intelligent pour faire une multiplication 32×32→64 bits au lieu de 64×64→64 bits, et utilise même l'instruction SuperH xtrct pour sortir les 32 bits du milieu.)

Implémentation en pratique

Je ne voudrais pas trop m'étendre sur le code sinon il y en a pour des pages, mais voici les grandes idées. Je vous conseille de donner un nom différent au type des nombres en point fixe, parce que si vous utilisez int vous allez tout le temps confondre les vrais entiers avec les nombres en point fixe, et croyez-moi vous n'avez pas envie de debugger ça.

Par exemple, on peut appeler ce type fixed.


Notez que le compilateur ne vous empêchera pas d'additioner un int et un fixed. Le nommage ne sert qu'à vous pour bien vous souvenir de ce que représente chaque variable. En C++, si vous enrobez l'entier dans une classe, vous pourrez empêcher les interactions de ce genre, et vous ne le regretterez pas ; je vous le conseille. Vous pourriez enrober l'entier dans une structure en C et obtenir le même effet, mais ça rendrait le code beaucoup plus verbeux, donc c'est pas forcément intéressant.

Pour obtenir un fixed à partir d'un entier ou d'un nombre décimal, multipliez-le simplement par 2¹⁶ = 65536.


Vous pourriez être tentés de faire << 16 à la place de la multiplication, mais ça donne des warnings (en C) ou des erreurs (en C++) pour les nombres négatifs. C'est parce que le C/C++ n'impose pas que le complément à 2 soit utilisé pour les entiers, et si le complément à 2 n'est pas utilisé, un décalage de bits d'un nombre négatif vers la gauche serait illégal/insensé. Ne vous embêtez pas, écrivez *65536, le compilateur optimisera ça en un décalage de bits absolument toutes les fois.

Notez que le 1.625 * 65536 est calculé par le compilateur, vous n'allez pas vous retrouver avec une multiplication de double pour initialiser votre nombre en point fixe (ce serait fort !).

Pour additioner et multiplier, ça donne les opérations que j'ai mentionnées plus haut ; ça marche aussi pour soustraire, diviser, le négatif...


Je vous conseille d'enrober ça dans des fonctions pour ne pas avoir à le réécrire à chaque fois. Dans ce cas, il est de bon goût de les avoir en static inline dans un header, pour que le compilateur les insère sur place au lieu de faire un vrai appel de fonction (qui prendrait aussi longtemps que l'opération en elle-même).

Pendant que j'y suis, je vous invite très fortement à avoir une fonction fix() pour générer les constantes, histoire d'éviter d'écrire x+1 par erreur (ça m'est arrivé, c'est aussi chiant à debugger que ça en a l'air).


Vous ne pouvez pas écrire x*y à moins de faire du C++, il faut se contenter la notation mul_fixed(x,y) et du boost énorme de performances qui l'accompagne par rapport aux nombres en point flottant.

Avantages, inconvénients, conclusion

Voilà pour les grandes idées du calcul en point fixe. Les avantages sont assez clairs, vous allez avoir un code beaucoup plus compact et beacoup plus rapide, ce qui devient vite très important pour les programmes qui calculent beaucoup ou les jeux qui ont pas mal de physique.

Il y a plusieurs inconvénients qu'il convient de mesurer pour éviter les déconvenues.

• D'abord il y a la plage de valeurs. Généralement on se souvient assez bien que la précision c'est 1/2¹⁶ = 0.000015. Par contre on a tendance à oublier que la limite supérieure c'est ±32768, donc attention à ne pas dépasser sous peine d'overflow (et contrairement aux flottants il n'y a pas de valeurs infinies qui garantissent un semblant de résultat correct).
• Ensuite il y a le risque de confondre un int et un fixed. On peut écrire fixed y = x+1 mais ça ne donne pas ce qu'on veut (ça donne y = x + 1/2¹⁶ en termes décimaux). Beaucoup de rigueur est nécessaire ici (ou alors une définition à base de structure/classe qui l'empêche).
• Enfin, il y a la syntaxe. On ne peut pas écrire z = x * y, il faut utiliser une fonction. C'est un peu irritant par moments, mais croyez-moi ce n'est pas le pire sacrifice que vous ferez pour les performances. x3

Et dans un moment de publicité qui fera rager Yatis, j'observerai qu'on peut implémenter en C++ une classe de point fixe qui élimine les deux derniers points (ie. pas de risque de mélanger avec des entiers, et syntaxe naturelle avec les opérateurs normaux) sans compromis sur l'optimisation du code résultant. En fait, le C++ donne même des outils pour exiger que certains calculs soient faits à la compilation, ce qui est encore plus avantageux que la version C.

J'ai joué avec ça dans une bibliothèque numérique pour mon moteur de jeu en C++ (un truc plus gros avec différentes tailles de point fixes, et de l'arithmétique de vecteurs/matrices/etc), ça existera sûrement de façon réutilisable un jour. Honnêtement c'est assez agréable à utiliser.

---

Citer : Posté le 23/06/2021 09:58 | #

---

Hey ! Super, bravo !
J’ai juste une petite question :

Pourquoi dans 2874 la valeur -3 est implicite ?
Le nombre pourrait tout autant être 287,4, et à ce moment là ce ne sera pas -3 mais le nombre stocke dans la mémoire sera quand meme 2874…

---

Citer : Posté le 23/06/2021 10:05 | #

---

Mais, si tu veux 287,4 tu fais comment vu que c’est un chiffre avant la virgule ?

---

Citer : Posté le 23/06/2021 10:25 | #

---

Pour faire simple, et c'est ce que Lephe explique, ça revient exactement à compter en centimes.

Cette baguette coûte 110 centimes, ce croissant 150 centimes. Au total j'ai payé 260 centimes, mais vu que je sais que ce sont des centimes, je sais que ça représente 2,60 €

Et si je veux plus de précision, par exemple sur le prix de l'essence, je peux compter en millièmes d'euros. Le litre de gasoil coûte 1345 millièmes d'euros, si j'en prends 10 litres ça donne 13450 millièmes d'euro, et vu que je sais que je compte en millièmes, ça fait 13,45 €

Mais dans tous les cas il faut se fixer une convention et s'y tenir, parce que si je ne te dis pas dans quel unité est le nombre 4235, tu ne saura pas si c'est 4235 €, 42,35 € ou 4,235 €.

---

Citer : Posté le 23/06/2021 16:35 | #

---

Super tuto, bien construit et fourni

Je me demandais juste : l'objectif de tout ceci (la virgule fixe) est d'éviter toute opération sur des float/double qui soit gourmande non ?
Alors je ne comprends pas trop cette ligne-ci :


Si j'ai bien compris, il serait plus aisé (ça prendrait moins de temps) d'exécuter ce code si on calculait directement le résultat (pour éviter de faire une opération avec une virgule flottante), mais un problème qui surviendrait serait de se souvenir ce qu'on y a fait...

Encore jamais fait d'add-in, mais ce genre de tuto me dire super envie de commencer ! (mais bon, vu tout ce que je dois finir en ce moment...)

---

Citer : Posté le 23/06/2021 18:14 | #

---

Mais …



Attend, 1 x 65536 = 65536 non ?

Et pourquoi 2 puissance 16 ? Il y a un chiffre donc 1 bit ?!

---

Citer : Posté le 24/06/2021 09:48 | #

---

Merci beaucoup pour le travail derrière l'article. La lecture est fluide et le contenu aussi précis qu'intéressant.

Ce à quoi je ne m'attendais pas : "Notez juste qu'en Basic CASIO toutes les variables sont des nombres en point flottant en base 10".

Est-il possible d'avoir un microcontrolleur disposant d'un FPU qui ne sache traiter que les float et non les double ? L'ordre de grandeur des temps présentés dans l'article entre l'addition de deux entiers et celui d'un flottant pour une calculatrice qui n'a pas de FPU me fait penser à celui déjà rencontré sur certains cibles entre les float et les double émulés.

Dans le cas d'une target avec FPU, les nombres en point fixe font-ils concurrence aux nombre en virgule flottante ? La réponse est peut-être hardware specific...

Hâte de lire les prochains articles !

---

Citer : Posté le 24/06/2021 15:21 | #

---

Notons que dans le cas des GPU, les performances sont souvent exprimées sous la forme de trois valeurs en flops : half, single, double.
Et effectivement plus on réduit la taille, plus les performances sont élevées

Exemple ici avec la RX 5700XT : https://www.techpowerup.com/gpu-specs/radeon-rx-5700-xt.c3339

---

Citer : Posté le 27/06/2021 00:08 | #

---

Précalcul des divisions (Maths — ★★☆)

Introduction

• Motivation pour éviter les divisions
• Multiplication en point fixe

Exploration mathématique et preuve

• Pourquoi le format 0:32 n'est pas suffisant
• Une approximation correcte avec m = 32+k+1
• Le 33ème bit de l'approximation
• Vérifions que ça marche
• Le cas particulier où m = 32+k suffit

Implémentation

• Code de précalcul et code de division
• Programme pour générer les inverses et expérimenter avec les formules
• Comment le compilateur C++ peut tout calculer à votre place
• Comparaison du code assembleur

Cette technique est longue mais c'est très progressif. Ne soyez pas intimidé·e ! Si c'est trop dur regardez juste la partie «Implémentation», vous aurez le code sans les maths.

Ajout plus tard : Sur injonction de Dark Storm j'ai fait plus de tests, et on dirait que le compilateur fait parfois cette optimisation tout seul ! Regardez bien le code assembleur avant de vous lancer, ce serait dommage d'optimiser pour rien.


Motivation pour éviter les divisions

Les opérations arithmétiques sur les entiers ne coûtent pas très cher ; les additions, soustractions, négations, décalages de bits, et autres carry et overflow sont des opérations basiques qui occupent l'ALU pendant 1 seul cycle. Les multiplications prennent plus longtemps ; 1 cycle pour les multiplications 16 bits × 16 bits → 32 bits (très rares en C, on s'amuse rarement à ça) et 3 cycles pour les multiplications 32 bits × 32 bits → 64 bits. Le multiplieur est séparé de l'ALU donc on peut quand même travailler pendant les 2 cycles de délai tant qu'on n'a pas besoin du résultat de la multiplication, donc en général ça ne pose aucun problème de performances.

La division, par contre, est une autre affaire. Diviser est de loin l'opération la plus compliquée, et la méthode classique ressemble à ça.


C'est exactement comme on le fait à l'école primaire ; typiquement pour diviser 84934 par 23 on commence à regarder combien de fois 23000 passe dans 84934 (→ 3) ; on soustrait 3 fois 23000, il reste 15394 et maintenant on cherche combien de fois 2300 passe dedans (→ 6) ; et ainsi de suite jusqu'à trouver le résultat 3692.

En binaire, c'est un tout petit peu plus simple parce que chaque chiffre ne peut être que 0 ou 1, donc on ne demande pas « combien de fois » chaque multiple du diviseur passe, mais simplement si le multiple passe. C'est matérialisé par la condition « q × 2^i ≤ p ».

Dans le processeur SuperH, la méthode est similaire. On calcule un bit à la fois, et donc ça prend autant d'opérations qu'il y a de bits dans le résultat. Spécifiquement, la séquence d'instruction est la suivante :

• On commence par indiquer si on veut faire une division signée ou non signée, en utilisant l'instruction div0s ou div0u.
• Ensuite on exécute l'instruction div1 pour calculer le premier bit du résultat, qui est stocké dans T (le bit de test où est habituellement stocké le résultat des comparaisons).
• On récupère T pour le mettre dans un registre, soit à la main soit par effet de bord de div1.
• Et ensuite on recommence jusqu'à ce que tous les bits soient calculés.

Quand le résultat fait 8 ou 16 bits on peut récupérer T lors du div1 suivant (par un mécanisme bien foutu), mais GCC divise toujours les int de 32 bits et pour ceux-là il faut récupérer T à la main avec rotcl. La fonction _sdivsi3 de GCC qui calcule la division entre deux entiers signés ressemble donc à ça :


Le coût de cette procédure comprend le test pour la division par 0, l'appel de la fonction (pas négligeable, surtout que les conventions d'appel obligent à sauvegarder des registres), et bien sûr 32 coups de div1 et rotcl, le tout comptant probablement dans les ~80 cycles. C'est pas la mort, mais dès qu'on met ça dans une boucle ça commence à piquer pas mal.

Le principe de cette technique d'optimisation est de remplacer les divisions par des multiplication 32×32 → 64 bits, qui prend en tout autour de ~10 cycles, dans les cas où on connaît le diviseur à l'avance.

Multiplication en point fixe

Imaginons que je veux diviser un entier n par une constance c=10. Je peux diviser par 10... ou je peux multiplier par 1/10. Bien sûr 1/10 ce n'est pas un nombre entier, mais 2³²/10 c'est presque un nombre entier, si je peux multiplier n par ça j'obtiendrais à peu près 2³² × n/10, et ensuite c'est une affaire de diviser par 2³² (décalage de bits) pour récupérer le résultat.

Il faut faire attention quand même, parce que 2³²/10 n'est pas exactement un entier, et donc on n'aura qu'une approximation. L'erreur d'approximation se fait multiplier par n, et il faut qu'elle soit assez petite pour que le résultat reste correct.

Dans cette technique, je vais dérouler les calculs pour vous montrer comment on approche une question mathématique de ce genre, et comment on peut trouver la solution par le raisonnement plutôt que par le test.

Je présuppose que les signes sont traités séparément (donc les nombres sont positifs). Si c est une puissance de 2, on peut faire un décalage de bits, donc on ne s'embête pas avec cette méthode. Pour l'approximation, je vous demande d'admettre qu'il faut arrondir vers le haut sinon la méthode ne marche pas du tout.

Il y a pas mal d'entiers et de nombres décimaux dans les calculs, donc je vais utiliser les notations courtes pour les parties entières : ⌊x⌋ représente la partie entière inférieure de x, ⌈x⌉ sa partie entière supérieure, et {x} sa partie décimale (on a donc x = ⌊x⌋ + {x}, et aussi ⌊x⌋ ≤ x ≤ ⌈x⌉).

Et voilà les variables du problème :

• n est l'entier qu'on veut diviser (qu'on ne connaît pas).
• c est le diviseur constant connu à l'avance (par exemple c=10).
• x est notre approximation de 1/c, fois une certaine puissance 2^m, arrondi vers le haut : x = ⌈2^m / c⌉. Pour l'instant on ne sait pas quelle valeur de m va marcher (spoiler : pas 32), le but va être de la trouver.
• ϵ c'est l'erreur d'approximation : ϵ = x - 2^m/c (et donc 0 ≤ ϵ < 1).
• Et enfin k est le logarithme en base 2 de c, le nombre tel que 2^k < c < 2^(k+1).

Ce qu'on veut faire, c'est profiter qu'on connaisse la valeur de c pour calculer ⌊n/c⌋ en remplaçant la division par une multiplication par x. La valeur de x qu'on a choisi c'est ⌈2^m / c⌉ = 2^m/c + ϵ, et donc la multiplication va donner

nx = n(2^m/c + ϵ) = 2^m × n/c + nϵ = ⌊2^m × n/c⌋ + {2^m × n/c} + nϵ.
À la fin on va re-diviser par 2^m (avec un décalage de bits) et on perdra au passage la partie décimale de ce résultat. Autrement dit, on obtiendra ⌊nx/2^m⌋, et toute la technique consiste à trouver la bonne valeur de x (plus spécifiquement la bonne valeur de m) qui fait que ⌊nx/2^m⌋ = ⌊n/c⌋.

Pourquoi le format 0:32 n'est pas suffisant

L'exemple introductif prend x = ⌈2³²/c⌉, ce qui revient, si vous avez lu la technique sur le calcul en point fixe, à représenter 1/c en point fixe 0:32. Il est facile de voir que ce n'est pas assez précis ; l'erreur causée par la partie entière est dans l'intervalle [0:2⁻³²[, donc quand on multiplie par un n qui peut être n'importe quelle valeur entre 0 et 2³²-1, l'erreur se fait multiplier par 2³²-1 (on regarde le pire cas) et on obtient une erreur qui peut être à peu près n'importe quel nombre dans [0:1[. C'est trop, puisque ça peut parfaitement nous faire passer d'un entier à l'autre (par exemple une erreur de 0.3 transforme un résultat correct de 88/10 = 8.8 tronqué à 8, en un résultat incorrect de 88/10 + 0.3 = 9.1 tronqué à 9).

La raison pour laquelle ça ne marche pas, c'est que quand on calcule ⌈2³²/c⌉ on perd trop de précision. Par exemple si je prend 1 million comme constante, j'obtiens 4295, un nombre de 13 bits, alors que j'ai 32 bits à ma disposition. Vous pourriez être tenté·e de penser que quand on divise par 1 million le résultat est forcément petit, donc il y a peu de bits à calculer, donc c'est normal si on a besoin de peu de bits de précision. C'est raisonnable, mais malheureusement ce n'est pas vrai, et on va voir pourquoi.

Une approximation correcte avec m = 32+k+1

Lorsqu'on calcule le produit nx, il y a deux facteurs qui contribuent à la partie décimale du résultat, mis en valeur dans l'équation présentée dans l'introduction :

• L'erreur multipliée, nϵ ;
• Et la partie décimale du résultat du produit, {2^m × n/c}. (On se souvient que n/c n'est pas un entier, et 2^m × n/c non plus.)

À la fin on veut diviser par 2^m (par décalage de bits, ce qui arrondit vers le bas) et obtenir ⌊n/c⌋, ce qui signifie que ces deux termes doivent disparaître (le troisième terme ⌊2^m × n/c⌋ nous donnant précisément ⌊n/c⌋ après le décalage de bits). Pour ça, il faut que {2^m × n/c} + nϵ < 2^m.

La difficulté c'est que dans la division, {n/c} peut être très proche de 1 (du genre 0.99 si jamais n/c = 5.99) et il faut que l'erreur soit malgré tout assez petite pour éviter qu'on déborde sur le multiple de 2^m supérieur.

L'astuce pour s'en sortir est de remarquer que {n/c} ne peut pas être n'importe quelle valeur : c'est soit 1/c, soit 2/c, etc, jusqu'à (c-1)/c. Par exemple si c=10, {n/c} ne peut être que 0.1, 0.2, ... ou 0.9 (et donc pas 0.99 par exemple). Si on arrive à garantir que nϵ < 2^m/c alors on réussira à éviter les erreurs, puisque même dans le pire cas on aura {2^m × n/c} + nϵ ≤ 2^m(c-1)/c + nϵ < 2^m, et donc tout disparaîtra sans laisser de traces lors du décalage de bits final.

On veut donc que le résultat de la multiplication soit précis à nϵ < 2^m/c, ce qui montre bien que plus c est grand (par exemple 1 million) plus la précision nécessaire est grande. Ainsi, contrairement à ce que notre intuition nous dit, même quand le diviseur est grand on a besoin d'absolument tous les bits disponibles dans x pour maintenir une précision élevée.

Maintenant qu'on connaît la précision nécessaire, on peut déterminer la bonne valeur de m. Comme n peut aller jusqu'à 2³²-1, pour avoir nϵ < 2^m/c il faut que ϵ < 2^(m-32)/c. On sait de plus que c est coincé entre 2^k et 2^(k+1), et qu'il faudra donc une précision de 1/2^(k+1) pour atteindre 1/c (comme on est en binaire on n'a pas d'intermédiaire entre k bits et k+1 bits ; on est obligés de prendre k+1). Par conséquent, il faut que ϵ < 2^(m-32-k-1).

On se rappelle que 0 ≤ ϵ < 1 compte tenu de la façon dont on a formé notre calcul, et donc pour que ça marche à tous les coups on est obligés de prendre au moins m = 32+k+1. Et ça, contrairement à notre idée initiale de m = 32, ça va passer.

On peut revenir à l'exemple de c = 1 million pour comprendre l'importance de m = 32+k+1. Avec m=32, on avait x = 4295, un nombre pitoyable de 13 bits qui laissait 19 bits sur 32 inutilisés. Avec m=32+k+1, on magnifie suffisament 1/c pour exploiter toute la place disponible, et on obtient x = 4503599628, un nombre qui ne gâche aucun des 32 bits.

Le 33ème bit de l'approximation

Si vous avez le regard affûté, vous aurez remarqué que x = 4503599628 ça ne tient pas dans 32 bits. En fait, m = 32+k+1 ça nous oblige à calculer 33 bits de 1/c, et donc dans l'absolu il nous manque 1 bit pour représenter x.

Cependant, le 33ème bit vaut toujours 1. C'est le principe en fait ; on a choisi m = 32+k+1 spécifiquement pour amener le premier bit non nul de 1/c à la 33ème position et stocker toutes les décimales binaires importantes dans les 32 bits disponibles.

Du coup on peut facilement stocker non pas x mais x - 2³² et ensuite faire une addition pendant le calcul pour rajouter le 2³². Le code est le suivant :


Vous pouvez voir qu'on décale n*x de 32+k+1 bits, comme promis ; on s'arrête en chemin après 32 pour rajouter n, ce qui revient à ajouter 2³²×n au résultat (ce qui compense le 2³² qui n'était pas présent dans x lors de la multiplication).

Vérifions que ça marche

Maintenant qu'on a bien choisi notre valeur, on peut (re)vérifier que tout marche bien. Puisque m = 32+k+1, notre approximation de 1/c est



x ∈ [0..2³²[ comme promis. Et comme on a arrondi vers le haut, on a bien le ϵ qu'on veut.



Maintenant, on peut prendre n'importe quelle valeur de n entre 0 et 2³²-1, et calculer le produit avec l'addition spéciale.



De ce résultat, on veut éliminer q durant le décalage de 32+k+1 bits vers la droite, et garder uniquement ⌊n/c⌋. On veut donc que q < 2^(32+k+1), et comme {n/c} peut atteindre (c-1)/c, on veut avoir une garantie forte que



Et tout se goupille puisque ϵ < 1, c < 2^(k+1), et n < 2³². Autrement dit, une fois le décalage de 32+k+1 bits vers la droite effectué, on récupère exactement ⌊n/c⌋.

Et voilà !

Le cas particulier où m = 32+k suffit

Le 33ème bit est un peu frustrant quand même, il nous oblige à faire une addition en plus. Séparer le décalage de bits en deux est inévitable (on verra pourquoi quand on comparera le code assembleur), mais l'addition c'est toujours un cycle en plus dont on aimerait bien se passer. Heureusement, il se trouve qu'on peut parfois l'éviter.

Il n'y a pas de miracles en fait. On a choisi m = 32+k+1 pour représenter un maximum de bits de 1/c, à tel point que le bit de poids fort (la k+1ème décimale) est arrivé à la 33ème position pour laisser la place au 32+k+1ème. Vous allez me dire, ok, mais si le 32+k+1ème bit c'est 0, ça ne sert à rien. Et oui, c'est exact. Si le 32+k+1ème bit est 0, on peut prendre m = 32+k et éviter d'avoir un 33ème bit implicite (parfois ça marche aussi par pure chance).

Dans ce cas, le code ressemble à ça :


Il y a une autre situation très courante où c'est suffisant. Si n est un entier signé, alors la plage de valeurs possibles pour n n'est pas [0:2³²[, mais [2⁻³¹:2³¹[. Une fois les signes gérés, on se retrouve avec [0:2³¹[. Comme n est 2 fois plus petit qu'avant, l'erreur nϵ a 2 fois plus de marge, ce qui est précisément ce dont on a besoin pour réduire la précision de 1/c d'un bit. Et donc là aussi, on peut prendre m = 32 + k.

La preuve est encore plus simple ; on prend



Puis on considère n'importe quel n ∈ [0..2³¹[, ce qui nous donne le produit



et une fois encore q < 2^(32+k) puisque {n/c} ≤ (c-1)/c et nϵ < 2^(32+k)/c grâce à la borne plus stricte sur n.

Code de précalcul et code de division

D'abord voici la fonction qui calcule k (le logarithme en base 2 de c) ; elle divise simplement c par 2 jusqu'à atteindre 0.


Cas général — Dans le cas général (m = 32 + k + 1), on calcule x de la façon suivante :


Le fait d'ajouter (c-1) avant de diviser par c produit un arrondi vers le haut au lieu d'un arrondi vers le bas (c'est une technique classique).

Ensuite la multiplication se fait comme ceci :


Par exemple, avec c = 1729 :


Cas restreint — Dans le cas restreint (m = 32 + k), on calcule x comme ça :


La multiplication se fait en un seul coup :


Par exemple avec n=10, qui a la chance de marcher avec tous les n (même non signés) parce que son x normal est pair :


Avec ça, vous pouvez tout précalculer ; mais continuez à lire, il y a des options plus rapides.

Programme pour générer les inverses et expérimenter avec les formules

Le programme ci-dessous calcule les valeurs de k et x pour les deux méthodes (m = 32+k+1, et m = 32+k), et vérifie toutes les divisions possibles sur 32 bits pour déterminer si m = 32+k suffit. Les tests sont parallélisés sur 4 threads et prennent moins de 20 secondes sur mon modeste Core i3-8100.

Programme : fixed-point-inverse.c

Compilez le programme avec les optis et la bibliothèque de threads.


Le programme rapporte des erreurs même pour les méthodes qui échouent complètement, donc n'hésitez pas à jouer avec et tester des variations. Par exemple, m = 32+k+1 peut être plus rapide en assembleur si k+1 vaut 2, 8 ou 16 parce que le jeu d'instructions peut faire ces décalages de bits sans spécifier la valeur dans un registre.

Voilà ce que donne le programme sur un diviseur où m = 32+k suffit sur toutes les entrées :


Et voilà un exemple sur un diviseur où m = 32+k ne marche pas sur tout le monde (mais quand même sur tous les nombres signés, puisqu'on a prouvé que c'était toujours le cas) :


N'hésitez pas à noter ces résultats dans des fonctions dans un header, histoire pas avoir à copier les constantes partout :


Vous voilà paré·e pour la version C de cette technique.

Comment le compilateur C++ peut tout calculer à votre place

Le compilateur C++ est très fort. Déterminer la valeur de x est un truc très facile qu'il peut parfaitement faire à votre place. Déterminer si m = 32+k suffit est plus difficile, mais on peut se contenter de détecter le cas où x est pair et laisser l'addition dans les quelques autres même si elle n'est pas nécessaire. L'avantage de tout ça c'est qu'il n'y a rien à précalculer pour nous, le compilateur fait tout.

On a besoin de trois outils du C++ pour faire ça :

• D'abord un template : on va donner la valeur de c et le compilateur va générer la fonction qui divise par c (c'est le principe du template ; on donne un code générique et ensuite le compilateur génère des versions spécialisées).
• Ensuite des fonctions constexpr, qui sont une catégorie particulière de fonctions que le compilateur peut précalculer durant la compilation.
• Et enfin le if constexpr, une version spéciale de if que le compilateur évalue aussi à la compilation.

Voilà comment ça se passe. On commence par définir les fonctions utiles pour calculer ; le u32log2(), la fonction fixdiv_k1() qui calcule x dans le cas où m = 32+k+1, et la fonction fixdiv_k() qui calcule x dans le cas où m = 32+k.


J'ai mis constexpr sur les trois, ce qui indique au compilateur que les fonctions peuvent être évaluées (exécutées) à la compilation. Bien sûr, le compilateur ne me croit pas sur parole ; il y a des règles qu'on doit respecter, et évidemment pour que la fonction puisse être évaluée à la compilation il faut que tous les paramètres soient connus à la compilation.

Ensuite, vient l'outil vraiment magique, le template. La façon dont ça marche est que je crée une fonction fixdiv mais je dis qu'il y a un paramètre c qui est à déterminer à la compilation (le paramètre du template). Chaque fois qu'on appelera fixdiv, le compilateur générera une copie de la fonction où c sera hardcodé (naturellement il faut que c soit une constante).


Le corps de la fonction implémente la technique complète de précalcul. Le compilateur commence par calculer x avec fixdiv_k1 (c'est-à-dire m = 32+k+1). Si x est pair, on décide d'utiliser m = 32+k donc on recalcule x et on utilise la méthode directe. Sinon, on reste sur la méthode avec l'addition supplémentaire.

Le if constexpr indique que la condition peut (doit) être calculée à la compilation (ce qui possible puisque tout ne dépend que de la valeur de c), et donc le compilateur ne génère qu'une des deux branches. Grâce aux optimisations, tous les appels à fixdiv_k1 et fixdiv_k sont remplacés par les résultats, si x est pair le calcul initial de fixdiv_k1 est entièrement éliminé, et on obtient vraiment la même chose qu'en C, mais sans avoir à travailler.


Cerises sur le gâteau, comme la fonction générée est aussi constexpr, elle est inlinée par défaut, ce qui veut dire qu'on a exactement le même résultat que si on avait précalculé la constante à la main et copié le code de la méthode appropriée à l'endroit de la division ; et si n est une constante le compilateur optimise aussi le calcul comme avec la division normale !

Programme : fixed-point-inverse-template.cpp

Comparaison du code assembleur

On peut obtenir le code assembleur de la division en définissant une fonction quelconque qui l'utilise (ou en prenant la version C, les résultats sont les mêmes) :


Et on obtient le code avec objdump -d :


La fonction commence par charger la valeur x = 0xcccccccd (dans r1), multiplie n (r4) par x, récupère les 32 bits du haut (mach) puis décale de 3 bits en deux temps (il n'y a aucune instruction qui décale de 3 bits d'un coup, l'autre option aurait été de charger la valeur 3 dans un registre et d'utiliser shld).

Temps total, sauf erreur de ma part : 9 cycles (3 pour le mov.l avec la dépendance RAW, 3 pour dmul.u avec la dépendance sur mach, ensuite 1 par calcul ; peut-être que sts et shlr2 parallélisent par forwarding et descendent le total à 8).

C'est le bon moment pour expliquer où est passé le décalage de bits par 32. En fait lorsqu'on multiplie 32 bits × 32 bits → 64 bits, le résultat est séparé en deux registres : mach (32 bits du haut) et macl (32 bits du bas). Le SuperH n'a pas de registres de 64 bits, donc on ne peut manipuler qu'une moitié à la fois. Le compilateur transforme donc le décalage de 32 bits du résultat vers la droite en le fait de lire mach et d'ignorer macl. Si on voulait un mélange des bits des deux ce serait assez casse-pieds, heureusement le cas ne se présente pas avec cette méthode.

Et voilà, avec un peu d'effort de preuve, du précalcul si on fait du C ou un template constexpr si on fait du C++, on a réduit d'un facteur 8-10 le temps nécessaire pour diviser par des constantes. Diviser par 10 en particulier est super utile pour les fonctions d'affichage de nombres en décimal.

---

Citer : Posté le 30/06/2021 11:53 | #

---

Alors j'ai pas inventé la méthode non. Enfin je l'ai réinventée (si j'avais su...), mais ensuite j'ai retrouvé des traces assez facilement jusque dans des vieux papiers de recherche.

Il faut voir d'abord que la méthode n'est pas entièrement directe. D'une part là je n'ai présenté que la version non signée, sachant que la version signée nécessite test et inversion de signe en plus. Et d'autre part j'ai éludé la division par 0, qui a une sémantique différente selon les architectures (exception en point flottant sur x86 par exemple). Et du coup tu peux probablement pas faire l'opti directement sur l'IR parce que tu sais pas exactement quel comportement tu veux obtenir, il faudrait la faire dans le back-end (ce qui ouvre tout un paquet d'explications parce que tu te doutes que le backend SuperH est pas le plus perfectionné).

Ensuite, ce n'est pas forcément plus rapide sur toutes les architectures. Ici la division prend 80 cycles, mais sur un x86 par exemple tu peux payer du genre 3-6 cycles pour la multiplication et 12-40 cycles pour la division donc c'est pas évident qu'avec le branchement signé et le surcoût de charger une grosse constante (au lieu par exemple de 10) ce soit plus rapide.

Enfin, les compilateurs l'implémentent parfois, et même si je n'ai pas d'exemple hyper concret sous la main, le manuel d'optimisation x86 d'Intel en parle justement pour les divisions 128-bits (page 360) :


La section 20.2.6.2 référence ce paragraphe dans le contexte de toutes les divisions, et l'exemple 14-21 montre le cas exact de la division par 10 pour convertir les donnés numériques en texte (le "%d" essentiellement).

Après ça reste des trucs assez tordus, tu peux sans doute trouver facilement un compilo qui le fait mais tu ne peux pas raisonnablement espérer qu'un compilateur tiré au hasard le fasse.

Sur SuperH spécifiquement, GCC a une option de stratégie de division, mais ça ne concerne que les cas génériques :


Il utilise call-table de base, c'est pour ça que diviser par 0 produit un TLB miss (l'accès à la table donne n'importe quoi).

Donc pour répondre à ta question, le compilo ne le fait pas tout seul probablement pour une combinaison de (1) c'est pas toujours rentable selon les architectures, (2) la méthode s'alourdit si on veut exactement la même sémantique que la division, sans hypothèses, et (3) implémentation incomplètes, back-ends pas très développées, etc.

La morale c'est que si tu veux que ton code fasse exactement ce que tu lui dis, même que c'est quelque chose d'exotique, il vaut mieux être le plus explicite possible. Si l'opti est tellement dépendante du back-end que le compilateur n'en implémente rien tu as aussi à gagner à le faire en assembleur (et les parties assembleur de ce topic montreront qu'il y a pas mal à gagner, sur SH4 en tous cas).

---