Assembleur : des algorithmes de matheux, M. Horner

Au chapitre précédent, nous étions restés sur notre faim en terme de performances, ce qui a irrité plus que de raison le tenancier. Lors, puisque le coprocesseur arithmétique ne voulait plus travailler à une vitesse décente (sûrement un gréviste), il m'a fallu me tourner vers les algorithmes d'approximation des fonctions trigonométriques.

De quoi avons-nous besoin ?

Donc on garde tout. On se prend une verveine ou une camomille, quelque chose qui détend et qui ouvre les chakras. On s'ouvre une petite page wikipedia, celle-ci en l'occurence.

Aaaaaah ! Des maths ! Des vraies, avec des exposants, des indices et des symboles chelous ! On va décéder dans d'atroces souffrances.

Ce à quoi je répondrai que très certainement, mais que ces choses n'ont aucun rapport entre elles.

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

Oui. On va pouvoir calculer des fonctions trigonométriques sur des vecteurs à une fraction de la vitesse rendue possible par le calculateur. Et c'est très formateur pour la programmation, car pour gagner du temps d'accès mémoire, nous allons devoir modifier un algorithme jusque-là bien tenu.

Mais pourquoi, comment ?

Pourquoi, je viens de le dire, et comment : nous allons nous inspirer très fortement de cette page.

Du cosinus par les approximants de Padé

Incalculable

Vous l'avez sans doute déjà remarqué, mais la fonction cosinus n'est pas calculable par les quatre opérations standards que sont la multiplication, la division, la soustraction et l'addition. Or, et c'est dommage, ce sont les seules dont nous disposons nativement hors du coprocesseur arithmétique non vectoriel. Comme nous souhaitons faire du vectoriel, on l'a dans le baba.

Mais les mathématiciens sont des gens pleins de ressources, et ils ont trouvé des astuces pour calculer une grande partie des fonctions qui ne se calculent pas avec des +, -, *, /. Voyons cela.

De l'utilité de M. Padé.

M. Padé, Henri de son prénom et mathématicien de son état, s'est penché sur le problème de l'approximation de fonctions. Il a inventé une méthode qui aboutit à des polynômes, lesquels doivent être divisés l'un par l'autre pour aboutir à l'approximation de la fonction souhaitée. Ne me demandez pas comment il calcule cela, mais je peux vous donner une approximation de la fonction cosinus : (127*x**8-16632*x**6+753312*x**4-10926720*x**2+22619520)/(2352*x**4+383040*x**2+22619520). La commande dans Xcas est : "pade(cos(x), x, 12, 9)" pour les obtenir.

De l'utilité de Mr Horner.

Mr Horner, William Georges de son prénom, quoique mathématicien de métier, était anglais de son état. Nul n'est parfait. Il a donné son nom à une méthode de calcul des polynômes en un point quelconque. Nous allons appliquer sa méthode sur le numérateur de la fonction obtenue précédemment, puis sur le dénominateur correspondant. Il ne nous restera plus qu'à diviser l'un par l'autre et nous aurons notre approximation du cosinus (et du sinus, ça marchera pareil). L'avantage de la méthode Horner tient à sa bonne résistance aux débordements de représentation, mais allez lire la page Wikipédia, ils y expliquent cela mieux que moi.

Moi, ce que je sais faire, c'est transformer un algorithme en langage vernaculaire en assembleur. Parfois en passant par le C++, ce que je fis ici :

PUBLIC horner
horner PROC
;	for (auto i = 0; i < nb; i++)
;	{
;		double x0 = coeff[0];
;		double x = v[i];
;		v[i] = x0;
;		for (int j = 1; j < 9; j++)
;		{
;			v[i] = v[i] * x + coeff[j];
 ;		}
;	}
    shr rdx, 2
    mov rax, rcx
    mov r10, r8
main:
    vmovsd xmm1, QWORD PTR [r8];coeff[0]
    vpermpd ymm1, ymm1, 0
    vmovapd ymm0, YMMWORD PTR [r9];v[0]
    vmovapd ymm2, ymm0 ;
    vmovapd ymm0, ymm1
    add r8, 8
    dec rcx
boucle:
    vmovsd xmm1, QWORD PTR [r8];coeff[j]
    vpermpd ymm1, ymm1, 0
    vfmadd132pd ymm0, ymm1, ymm2
    add r8, 8
    loop boucle
    vmovntpd YMMWORD PTR [r9], ymm0
    add r9, 32
    mov rcx, rax
    mov r8, r10
    dec rdx
    jnz main
    ret
horner ENDP

La fonction s'appelle avec en premier argument, le nombre de coefficients du polynôme. En second, le nombre de points à traiter. En troisième, le pointeur sur la liste des coefficients, en commençant par le n-ième. En quatrième, l'adresse de la liste des valeurs à traiter, qui doit contenir un nombre divisible par 4 de points.

Du calcul du cosinus en un seul passage.

Tout cela est magnifique, c'est même superbe, mais on doit passer par la fonction horner une fois pour le numérateur, une fois pour le dénominateur, et on doit encore une fois parcourir l'ensemble des points pour en faire la division. Cela fait trois fois le parcours des valeurs dont on doit calculer le cosinus, ce qui fait deux fois de trop. Nous allons condenser cela en une seule fonction :

PUBLIC cosPade
cosPade PROC
    push r12

    shr rcx, 2
    mov r11, 9
    mov r10, 5
    mov r9, offset coeffNumCosPade
    mov r12, offset coeffDenomCosPade
main:
    vmovsd xmm1, QWORD PTR [r9];coeff[0]
    vpermpd ymm1, ymm1, 0
    vmovsd xmm3, QWORD PTR [r12];coeff[0]
    vpermpd ymm3, ymm3, 0
    vmovapd ymm0, YMMWORD PTR [rdx];v[0]
    vmovapd ymm2, ymm0
    vmovapd ymm4, ymm0
    vmovapd ymm0, ymm1
    vmovapd ymm5, ymm3
    add r9, 8
    add r12, 8
    dec r11
    dec r10
boucleNum:
    vmovsd xmm1, QWORD PTR [r9];coeff[j]
    vpermpd ymm1, ymm1, 0
    vfmadd132pd ymm0, ymm1, ymm2
    add r9, 8
    dec r11
    jnz boucleNum
boucleDenom:
    vmovsd xmm3, QWORD PTR [r12];coeff[j]
    vpermpd ymm3, ymm3, 0
    vfmadd132pd ymm5, ymm3, ymm4
    add r12, 8
    dec r10
    jnz boucleDenom
fin:
    vdivpd ymm0, ymm0, ymm5    ;L'instruction de la division
    vmovntpd YMMWORD PTR [rdx], ymm0
    add rdx, 32
    add r8, 32
    mov r11, 9
    mov r10, 5
    mov r9, offset coeffNumCosPade
    mov r12, offset coeffDenomCosPade
    dec rcx
    jnz main
    pop r12
    ret
cosPade ENDP

Cette fois, les coefficients des polynômes sont en dur dans le code. Plus précisément dans la section TEXT :

coeffNumCosPade dq 127., 0., -16632., 0., 753312., 0., -10926720., 0., 22619520.;
coeffDenomCosPade dq 2352., 0., 383040., 0., 22619520.;
coeffNumSinPade dq 551., 0., -22260., 0., 166320., 0.;
coeffDenomSinPade dq 75., 0., 5460., 0., 166320.;

Pour ceusses et celles qui veulent un exemple d'appel de cette fonction, voici un bout de code C++ :

void testCos(int nb)
{
	const int nbReps = 5000;
	double* a = (double*)_aligned_malloc(nb * sizeof(double), 32);
	double* c = (double*)_aligned_malloc(nb * sizeof(double), 32);
	for (auto i = 0; i < nb; i++)
	{
		a[i] = rand() * (2 * M_PI) / (double)RAND_MAX - M_PI;
		c[i] = a[i];
	}
	float temps;
	clock_t t1 = clock();
	for (auto i = 0; i < nbReps; i++)
	{
		for (auto i = 0; i < nb; i++)
		{
			a[i] = c[i];
		}
		cosPade(nb, a);
	}
	clock_t t2 = clock();
	temps = (float)(t2 - t1) / CLOCKS_PER_SEC;
	for (auto i = 0; i < 4; i++)
	{
		printf("cos %f = %f\n", c[i], a[i]);
	}
	for (auto i = 1; i < 5; i++)
	{
		printf("cos %f = %f\n", c[nb-i], a[nb-i]);
	}
	std::cout << "temps de cosinus : " << temps << std::endl;
	t1 = clock();
	for (auto i = 0; i < nbReps; i++)
	{
		for (auto i = 0; i < nb; i++)
		{
			a[i] = c[i];
		}
		cosVect2(nb, a);
	}
	t2 = clock();
	temps = (float)(t2 - t1) / CLOCKS_PER_SEC;
	for (auto i = 0; i < 4; i++)
	{
		printf("cos %f = %f\n", c[i], a[i]);
	}
	for (auto i = 1; i < 5; i++)
	{
		printf("cos %f = %f\n", c[nb - i], a[nb - i]);
	}
	std::cout << "temps de cosinus : " << temps << std::endl;
	_aligned_free(c);
	_aligned_free(a);
}

De la norme

Quand on parle de vecteurs, on arrive rapidement à la norme des vecteurs, c'est-à-dire, peu ou prou, à leur longueur. La formule est à peu près celle_ci : √∑(xb - xa)². Or, on voit que le premier élément est une racine carrée, chose coûteuse à calculer ! Donc, on essaie autant que possible de travailler avec la norme au carré, qui n'utilise que des additions, des soustractions et des multiplications.

De la norme au carré.

Ci-dessous vous trouverez une implémentation du calcul de la norme pour n vecteurs de dimension m. La déclaration de la fonction est celle-ci : void norme(int nbDim, int nbVect, double* v, double* results);

PUBLIC normeSquare_asm
normeSquare_asm PROC
    mov rax, rcx
    shl rdx, 3 ;Conversion du nombre de vecteurs en nombre d'octets
    mov r10, r8
    mov r11, r8
    add r11, rdx    ;Calcul de la fin de la première ligne
depart:
    vmovapd ymm0, YMMWORD PTR [r8]
    vmulpd ymm0, ymm0, ymm0
    dec rcx
boucle:
    add r8, rdx
    vmovapd ymm1, YMMWORD PTR [r8]
    vfmadd231pd ymm0, ymm1, ymm1
    loop boucle
    vmovntpd YMMWORD PTR [r9], ymm0
    add r9, 32
    add r10, 32
    mov r8, r10
    mov rcx, rax
    cmp r8, r11
    jb depart
    ret
normeSquare_asm ENDP

De la vraie norme.

Pour la vraie norme, il suffit de prendre la racine carrée (on a une instruction pour) du résultat précédent. L'ordre des paramètres ne change pas.

PUBLIC norme_asm
norme_asm PROC
    mov rax, rcx
    shl rdx, 3 ;Conversion du nombre de vecteurs en nombre d'octets
    mov r10, r8
    mov r11, r8
    add r11, rdx    ;Calcul de la fin de la première ligne
depart:
    vmovapd ymm0, YMMWORD PTR [r8]
    vmulpd ymm0, ymm0, ymm0
    dec rcx
boucle:
    add r8, rdx
    vmovapd ymm1, YMMWORD PTR [r8]
    vfmadd231pd ymm0, ymm1, ymm1
    loop boucle
    vsqrtpd ymm0, ymm0
    vmovntpd YMMWORD PTR [r9], ymm0
    add r9, 32
    add r10, 32
    mov r8, r10
    mov rcx, rax
    cmp r8, r11
    jb depart
    ret
norme_asm ENDP

Résultats.

Sur ma machine, on gagne un facteur 2 en temps d'exécution des sinus et cosinus, aux dépens d'un peu de précision. On gagne ce même facteur 2 pour le calcul des normes sans perte de précision. Je reste donc sur mon idée qu'un code assembleur humainement écrit est plus rapide qu'un code généré par un compilateur.

Assembleur : des algorithmes de matheux, M. Horner

De quoi avons-nous besoin ?

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

Mais pourquoi, comment ?

Du cosinus par les approximants de Padé

Incalculable

De l'utilité de M. Padé.

De l'utilité de Mr Horner.

Du calcul du cosinus en un seul passage.

De la norme

De la norme au carré.

De la vraie norme.

Résultats.