<-Précédent | Retour à l'accueil | Contact : etienne"at"les-sauvages.fr | suivant-> |
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.
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.
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.
Pourquoi, je viens de le dire, et comment : nous allons nous inspirer très fortement de cette page.
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.
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.
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.
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); }
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.
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
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
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.