1. Du sinus
2. De la racine carrée
3. Chaotique bon, loyal mauvais : l'alignement

Assembleur : calculs en virgule flottante

Tout bouffis de l'inénarrable orgueil que le précédent épisode nous a offert, continuons notre périple dans les joies des calculs vectoriels en virgule flottante. On est des oufs. L'objectif est ambitieux : faire plus rapide que le code C naïf optimisé par GCC. Je vous dis, des oufs.

De quoi avons-nous besoin ?

De la même chose que l'épisode précédent. A peu près. Presque. Peut-être un petit

sudo apt-get install yasm

Est-ce qu'à part frimer devant la machine à café, ça a une utilité ?

La racine carrée d'un vecteur de double en 5 fois plus vite, le sinus en 2 fois plus vite, ça ne vous intéresse pas ? Non mais vous faites comme vous voulez, hein. Moi je vous propose juste des fonctions qui prennent moins de place et plus rapides, après vous en faites ce que vous voulez, hein. Je ne suis pas là pour juger. Si vous aimez faire lent et compliqué, qui suis-je pour critiquer cela ?

Mais pourquoi, comment ?

Deux secrets ici : la vectorisation et les instructions, le tout en version matérielle.

Nos processeurs modernes incluent dans leur petit bidou de silicium d'autres petits processeurs, qui eux-même incluent d'autres petits processeurs. C'est vrai, un quad core par exemple contient 4 processeurs, et chacun contient à son tour l'Unité Arithmétique et Logique, le coprocesseur arithmétique dit x87, et encore d'autres. Et bien dans ces petits processeurs, il y en a un qui fait du calcul vectorisé en virgule flottante. Il fait d'autres genres de calcul, mais il fait celui-là. Si bien qu'au lieu de faire une somme, puis une somme, puis une somme... Il en fait deux d'un coup.

N'importe quelle fonction codée en quelque langage que ce soit ira toujours moins vite que sa version gravée dans le semi-conducteur. Il est donc particulièrement rentable d'utiliser les instructions fournies par le processeur. A tel point que sans hésiter, il faut abandonner le calcul vectorisé pour privilégier les instructions natives.

Et il se trouve que si gcc utilise bien la vectorisation, il n'utilise pas les instructions natives.

Des annoncées fonctions

Du sinus

Il n'y a pas d'instruction SSE qui calcule le sinus de quoi que ce soit. Alors j'avais ressorti mes leçons de dérivées, hop hop hop, développement limité on est partis. J'étais arrivé à quelque chose d'un peu plus rapide que GCC, quand soudain, avec une brusquerie de Hun dans une cristallerie d'art, le wikipédia anglophone me crache à la figure que depuis le coprocesseur arithmétique, le sinus est une fonction native. Pan dans les dents. Non. Mais. Depuis-le-coprocesseur-arithmétique i387. Sérieux les gars. Depuis 1987, il y a là, dans nos petites machines, une fonction sinus codée en dur ET PERSONNE NE L'UTILISE !!!!!!!!!!!!!!!

Votre serviteur
Pourquoi ? Est-elle nulle ? Lente ?
Est-elle molle et ballante, monsieur, comme une lente ?
Qu'a-t-elle d'hétéroclite ?
Est-ce un phénomène ?

math.lib
De l'inclure dans mon code j'avais su me garder !

Votre serviteur
Et pourquoi, s'il vous plaît, ne pas l'utiliser ?

math.lib
J'avais...

Votre serviteur
Elle vous dégoûte alors ?

math.lib
Monsieur...

Votre serviteur
Malsain
Semble son opcode ?

math.lib
Mais !

Votre serviteur
Son mnémonique, crétin ?

math.lib
Mais pas du tout !

Votre serviteur
Pourquoi coder en la dénigrant ?
Peut-être que monsieur veut faire du code beaucoup plus lent ?
Pour les trois du fond qui dormaient au collège : Cyrano de Bergerac

Voici le code pour l'essayer, c'te FSIN.

sinDV:
	mov rcx, rdi					;Nombre d'éléments à traiter mis dans RCX parce que c'est le registre compteur
.sin:
	fld qword [rsi - 8 + 8 * rcx]	;On met dans un registre du processeur à virgule flottante un double
	fsin
	fstp qword [rsi - 8 + 8 * rcx]	;On stocke le résultat
	loop .sin						;Et on boucle sur RCX
	ret

Avant de discuter des instructions proprement dites, il faut exposer comment ce vieux coprocesseur travaille. En une phrase, il travaille avec une pile. Plus précisément, il dispose de huit registres qui forment une pile. On n'accède qu'au premier registre, les opérations prennent leurs paramètres dans la pile. Je ne sais pas si vous connaissez les calculatrices de poche HP (Hewlett-Packard) : elles utilisent le même principe.

Bon, le MOV, vous le connaissez quand même ?

Ensuite, FLD, première instruction de toute notre vie du coprocesseur arithmétique, ho là là, quelle émotion ! F comme virgule Flottante et LD comme LoaD : chargement de la donnée sise à l'adresse indiquée par le paramètre dans la pile de registres du coprocesseur.

FSIN, comme SINus en virgule Flottante, pif-pouf le registre du hut de la pile, je crois que tu as attrapé une sinusite !

FSTP comme STocke Pop virgule Flottante, ce qui fait le contraire de FLD : stocke dans la mémoire le registre du coprocesseur et le sort de la pile (pop!, comme un bouchon). Au passage, FST existe mais je vous la déconseille : elle est sacrément plus lente.

LOOP comme on décrément RCX et s'il n'est pas nul, on retourne à .sin

Chers amis, dans votre grande sagacité, saurez-vous faire la fonction cosDV ?

Remarques essentielles : la fonction FSIN n'est pas juste. Elle diverge de la fonction de la bibliothèque mathématique dès qu'on sort de [-2pi, 2pi].

De la racine carrée

Ici, on travaille en SSE, parce qu'en SSE, on a une instruction SQRTPD, comme SQuare RooT Packed Double.

Voyons le code :

sqrtDV:
	xor r8, r8		;Une MàZ d'un registre qui va se souvenir de la parité du nombre d'éléments à traiter.
	mov rcx, rdi	;Nombre d'éléments à sommer mis dans RCX parce que c'est le registre compteur
	shr rcx, 1		;qu'on divise par 2, parce qu'on va les traiter par 2
	jnc .racine		;Et si la retenue (carry flag) est posée, la taille est impaire
	inc r8			;On se souvient que c'est impair pour plus tard
	test cl, cl
	jz .dernier
.racine:
	movapd xmm0, oword [rsi];On met dans un registre du processeur à virgule flottante 2 double
	sqrtpd xmm0, xmm0		;L'instruction de la racine carrée
	movntpd oword [rsi], xmm0;On stocke le résultat dans l'accumulateur
	add rsi, 16		;On passe à la paire de double suivante
	loop .racine		;Et on boucle sur RCX

.dernier:
	test r8, r8		;Si la taille des vecteurs est paire
	jz .retour		;Merci, c'est fini. Sinon, c'est impair.
	movsd xmm0, qword [rsi]	;On met dans un registre du processeur à virgule flottante le dernier double
	sqrtsd xmm0, xmm0	;On fait une seule somme
	movsd qword [rsi], xmm0	;Et donc une seule sauvegarde. Et fin.
	
.retour:
	ret

XOR, le shérif de l'espace, X-OR... Si, je vous jure. Bien. Bien bienbien. Découvrons ensemble ce R8. R8 ? Pourquoi pas R2D2 ? R8 comme Registre n°8. Et l'intelligence supérieure de mes lecteurs de s'emballer : où sont les 7 premiers ? Puis, avec la réflexion : pardon, les HUIT premiers ? Je vous les donne : RAX = R0, RBX = R1, RCX = R2, RDX = R3, RBP = R4, RSI = R5, RDI = R6 et RSP = R7. Sur architecture 64 bits, on a rajouté huit registres de 64 bits qu'on a aussi bêtement que d'habitude appelés registre numéro tant, Rx. L'absence de poésie en électronique ne cessera jamais de m'étonner. Et donc, j'ai besoin d'un registre pour stocker un malheureux bit, je prends R8, je le mets à 0 et c'est tout. Oui, XOR x, x met x à 0 plus vite que MOV x, 0. Essayez, sortez vos tables de vérité, on ne fait pas, jamais, MOV x, 0, on fait XOR x, x; le shérif annule et prend la pose.

MOV, le cultissime MOV, récupère le nombre d'éléments du vecteur pour le traiter nativement en tant que variable compteur.

SHR comme SHift Right, décalage à droite. De un bit, quantité passée en paramètre. Ca équivaut à diviser RCX par deux, binaire oblige. Parce qu'en fait, on ne va faire que la moitié des éléments, chaque instruction de racine carrée en traitant deux d'un coup. SSE power.

JNC comme Jump if Not Carry, bondit si la retenue n'est pas posée. Il faut savoir que SHx (décalage vers la droite ou vers la gauche), transfère le dernier bit décalé dans le bit de retenue. Vous savez qu'il y a un bit de retenue, tout de même ? Il y a un bit qu'on dit "de retenue" ou "carry" qui est mis à 1 par certaines opérations et qu'on peut tester avec, précisément, JNC. Donc, si le dernier bit de RCX est à 1, RCX est... impair ! C'est important comme information : ça veut dire qu'il y a un dernier nombre à traiter séparément. Et donc, justement, si RCX est impair :

INC, comme INCrémente. On ajoute 1, ici à R8 qui vaudra donc... sacrebleu, suivez ! 1 !

TEST, comme heu... C'est pas clair ? L'explication est un peu ardue. Ca teste la présence de bits à 1 dans l'opérande 1 selon ce qui est donné en opérande 2. En vrai, ça fait un ET logique entre les 2 opérandes, et ça met les drapeaux correspondants. Les drapeaux ? Non, les drapeaux, vous savez ? Bon, ben le bit de retenue est un drapeau. Il en existe un autre, le bit de zéro. A chaque fois qu'une opération a la tête à toto comme résultat, il met à 1 ce bit. Il est à 0 dans tous les autres cas. Et bien, le seul cas où un ET logique entre un nombre et lui-même est à 0, c'est quand ce nombre est lui-même 0 ! C'est fou. Donc, tout comme on ne fait pas MOV x, 0, on ne fait pas CMP x, 0. On fait TEST x, x.

JZ comme Jump if Zero, bondit si le bit de zéro est posé. On est dans le cas où on aura un dernier nombre à traiter séparément parce que le nombre d'éléments est impair. Si en plus on n'a aucun double élément à traiter, on saute directement au cas impair.

MOVAPD, comme MOVe Aligned Packed Double, déplace un paquet de nombres à virgule flottante en double précision alignés sur des adresses mémoires multiples de 16 octets.

Votre serviteur
Tant de choses en un mnémonique ?

Covielle
Oui, la langue assembleur est comme cela, elle dit beaucoup en peu de paroles. Allez vite où il souhaite.
Le bourgeois gentilhomme

SQRTPD comme SQuare RooT Packed Double, racine carrée d'un paquet de double. Et oui, pas d'algorithme compliqué, juste une fonction codée en dur.

MOVNTPD comme MOVe Non Temporal Packed Double, déplace un paquet de double non temporel. Non temporel ? Spirituel, alors ? En quelque sorte, oui. La temporalité des données a à voir avec le cache. Vous savez ce qu'est un cache ? Un cache, c'est un bout de mémoire rapide. Le cache est à la RAM ce que le post-it est au livre. C'est petit, rapide mais pas pertinent sur le long terme. Le processeur a un cache qu'il gère tout seul comme un grand, dans le but de limiter l'accès à la RAM. Nous programmeurs n'y avons pas accès, même en assembleur. Tout au plus peut-on dire qu'il est inutile de stocker une valeur en cache. Ca sert à prévenir le processeur qu'on ne l'utilisera plus. Ici, en l'occurence, la racine carrée est le résultat, on ne l'utilisera plus, donc on peut faire fi du cache et écrire directement en mémoire, ce que fait MOVNTPD. Si la RAM est le paradis, la non temporalité permet donc de passer outre la case incarnation dans le cache pour atteindre directement la béatification : si c'est pas spirituel, ça.

ADD comme ADDitionne. Une paire de double, c'est 16 octets. La prochaine adresse à lire est donc RSI + 16.

LOOP, la boucle, RCX <- RCX - 1, tout ça...

TEST comme test, si 0, pair et on a fini, sinon on racinecarrise le dernier.

JZ comme qu'est-ce que je viens de dire.

MOVSD comme MOVe Single Double, déplace un double célibataire (et sans enfant). Cette instruction met dans les 64 bits de poids faible d'un registre XMM le double sis dans l'opérande source, ici la mémoire, le dernier élément du vecteur.

SQRTSD comme - vous pouvez deviner, là, quand même - SQuare Root Single Double, la racine carrée du célibataire qu'on vient de charger dans ce même registre.

MOVNTSD n'existe pas, donc on refait MOVSD pour remettre le contenu du XMM en mémoire.

RET comme RETour à l'appelant ! Fi-ni !

Chaotique bon, loyal mauvais : l'alignement

Au paragraphe précédent, j'ai utilisé le terme "alignés". Loin de moi de faire du militarisme, ça a plus à voir avec les tailles de trucs et les positions des bidules. On va regarder les adresses mémoires des choses, et plus précisément la fin de ces adresses. Mettons que nous ayons trois entiers de 4 octets. Les adresses référencent des octets. Si les adresses desdits octets se terminent par 11, 17 et 43, elles sont aussi alignées que peuvent l'être autant de hippies dans une danse chamanique. Par contre, si elles se terminent par, au hasard vraiment, 0, 4 et 8, là on parle d'alignement.

C'est souvent pratique, l'alignement. Par exemple, pour parcourir une liste de variables, si elles sont toutes espacées pareillement, on fera toujours la même opération pour passer à la suivante, on n'a pas à tester la taille de celle qu'on vient de lire. En informatique, on aligne sur des puissances de 2. Nos trois entiers, là, pourraient être alignés sur 1, 5, 9 théoriquement. Sachez que Mme de la Fressange l'a formellement interdit. Ca ne se fait pas. C'est 0, 4, 8 ou rien. Même pas 2, 6, 10, non. On aligne sur des multiples de la taille d'alignement. Ainsi, les double s'alignent sur... Alors, 64 bits divisés par 8... 8. Ca fait des adresses en 0, 8, 10 (Arg ! En hexa, le salaud !), 18... 0, 8, 0, 8, trop drôle ! Du coup, nos registres XMM de 128 bits, ils prennent leurs données sur 16 octets, soient des adresses qui se terminent uniquement par 0 !

Du coup, tant qu'à faire, autant aligner les données sur 16 octets. Surtout quand on en a plusieurs millions, c'est pas 15 octets maximum (l'espace perdu pour aligner 128 bits dans le plus mauvais cas) qui valent la peine de s'énerver.

Ta. Dam.

Du code

Il se trouve ici.