Calculs précis et rapides pour les nombres à virgule flottante en utilisant l'exemple de la fonction sinus. Partie 3: virgule fixe

Nous continuons le cycle de conférences ( partie 1 et partie 2 ). Dans la partie 2, nous avons regardé ce qu'il y a à l'intérieur de la bibliothèque libm et dans ce travail, nous allons essayer de modifier légèrement la fonction do_sin pour augmenter sa précision et sa vitesse. Je citerai à nouveau cette fonction do_sin ):







Comme indiqué dans l'article précédent, partie 132-145. Exécuté pour x dans la plage [0,126, 0,855469]. Bien. Essayons d'écrire une fonction qui, dans les limites données, sera plus précise et, éventuellement, plus rapide.



La façon dont nous l'utiliserons est assez évidente. Il est nécessaire d'étendre la précision des calculs pour inclure davantage de décimales. La solution évidente serait de choisir le type double long, de le compter, puis de le reconvertir. En termes de précision, la solution doit être bonne, mais en termes de performances, il peut y avoir des problèmes. Pourtant, le double long est un type de données plutôt exotique et sa prise en charge dans les processeurs modernes n'est pas une priorité. Sur les instructions x86_64 SSE / AVX avec ce type de données ne fonctionnent pas. Le coprocesseur mathématique sera "époustouflé".



Que devriez-vous alors choisir? Examinons de plus près l'argument et les limites de la fonction.



Ils sont dans la région 1.0. Ceux. en fait, nous n'avons pas besoin d'une virgule flottante. Utilisons un entier 64 bits lors du calcul de la fonction. Cela nous donnera 10 à 11 bits supplémentaires à la précision d'origine. Voyons comment travailler avec ces chiffres. Un nombre dans ce format est représenté par a / d , où a est un entier et d est un diviseur que nous choisissons constant pour toutes les variables et stockons "dans notre mémoire", et non dans la mémoire de l'ordinateur. Voici quelques opérations pour ces nombres:

$$

$$\frac{c}{d}=\frac{a}{d}\pm\frac{b}{d}=\frac{a\pm b}{d}\\ \frac{c}{d}=\frac{a}{d}\cdot\frac{b}{d}=\frac{a \cdot b}{d^2}\\ \frac{c}{d}=\frac{a}{d}\cdot x=\frac{a \cdot x}{d}$$

Comme vous pouvez le voir, il n'y a rien de compliqué à ce sujet. La dernière formule montre la multiplication par n'importe quel entier. Notez également une chose assez évidente que le résultat de la multiplication de deux variables entières non signées de taille N est plus souvent un nombre de taille jusqu'à 2 * N inclus. L'ajout peut provoquer un dépassement de 1 bit supplémentaire.



Essayons de choisir le diviseur d . Evidemment, dans le monde binaire, il vaut mieux la choisir comme puissance de deux, pour ne pas diviser, mais simplement déplacer le registre, par exemple. Quelle puissance de deux devriez-vous choisir? Trouvez l'indice dans les instructions de la machine de multiplication. Par exemple, l'instruction MUL standard du système x86 multiplie 2 registres et écrit également le résultat dans 2 registres, où 1 des registres est la "partie supérieure" du résultat et le second est la partie inférieure.



Par exemple, si nous avons 2 nombres de 64 bits, le résultat sera un nombre de 128 bits écrit dans deux registres de 64 bits. Appelons RH - casse "majuscule" et RL - casse "minuscule" ¹ . Ensuite, mathématiquement, le résultat peut être écrit comme$$$R=R_H \cdot 2^{64}+R_L$ . Maintenant, nous utilisons les formules ci-dessus et écrivons la multiplication pour$$$d=2^{-64}$

$$

$$\frac{c}{d}=\frac{a}{2^{64}}\cdot\frac{b}{2^{64}}=\frac{a \cdot b}{2^{128}}=\frac{R_H \cdot 2^{64} + R_L}{2^{128}}=\frac{R_H + R_L \cdot 2^{-64}}{2^{64}}$$

Et il s'avère que le résultat de la multiplication de ces deux nombres à virgule fixe est le registre $$$R=R_H$ .



C'est encore plus simple pour le système Aarch64. L'instruction "UMULH" multiplie deux registres et écrit la partie "supérieure" de la multiplication dans le 3ème registre.



Eh bien. Nous avons spécifié un nombre à virgule fixe, mais il reste un problème. Nombres négatifs. Dans la série Taylor, l'expansion s'accompagne d'un signe variable. Pour faire face à ce problème, nous transformons la formule de calcul du polynôme par la méthode Goner, sous la forme suivante:

$$

$$\sin(x)\approx x(1-x^2(1/3!-x^2(1/5!-x^2(1/7!-x^2\cdot1/9!))))$$

Vérifiez qu'elle est mathématiquement exactement la même que la formule originale. Mais dans chaque parenthèse, il y a un nombre comme$$$1/(2n + 1)! - x^2\cdot(\cdots)$ toujours positif. Ceux. cette conversion permet à l'expression d'être évaluée comme des entiers non signés.

constexpr mynumber toint    = {{0x00000000, 0x43F00000}};  /*  18446744073709551616 = 2^64     */
constexpr mynumber todouble = {{0x00000000, 0x3BF00000}};  /*  ~5.42101086242752217003726400434E-20 = 2^-64     */

double sin_e7(double xd) {
  uint64_t x = xd * toint.x;
  uint64_t xx = mul2(x, x);
  uint64_t res = tsx[19]; 
  for(int i = 17; i >= 3; i -= 2) {
    res = tsx[i] - mul2(res, xx);
  }
  res = mul2(res, xx);
  res = x - mul2(x, res);
  return res * todouble.x;
}

Où $$$tsx[i]=1/i!$au format virgule fixe. Cette fois, pour plus de commodité, j'ai posté tout le code sur legithubfast_sine, débarrassé de quadmath pour la compatibilité avec clang et arm. Et j'ai changé un peu la méthode de calcul de l'erreur.



La comparaison de la fonction sinusoïdale standard et de la fonction virgule fixe est donnée dans les deux tableaux ci-dessous. Le premier tableau montre la précision du calcul (c'est tout à fait la même chose pour x86_64 et ARM). Le deuxième tableau est une comparaison des performances.







Lors des tests, la valeur sinus "vraie" a été calculée à l'aide de la bibliothèque MPFR... La valeur ULP maximale a été considérée comme l'écart maximal par rapport à la valeur «vraie». Pourcentage d'erreurs - le nombre de cas où la valeur calculée de la fonction sinus par nous ou par libm ne correspondait pas à la valeur sinus arrondie au double. La valeur moyenne de l'écart indique le "sens" de l'erreur de calcul: surestimation ou sous-estimation de la valeur. Comme vous pouvez le voir dans le tableau, notre fonction a tendance à surestimer les valeurs sinusoïdales. Cela peut être corrigé! Qui a dit que les valeurs tsx devraient être exactement égales aux coefficients de la série de Taylor. Une idée assez évidente se suggère de faire varier les valeurs des coefficients afin d'améliorer la précision de l'approximation et de supprimer la composante constante de l'erreur. Il est assez difficile de faire correctement une telle variation. Mais nous pouvons essayer. Prenons par exemple4ème valeur du tableau des coefficients tsx (tsx [3]) et changez le dernier nombre a en f. Redémarrons le programme et voyons la précision (sin_e7a). Regardez, c'est un peu, mais augmenté! Nous ajoutons cette méthode à notre tirelire.



Voyons maintenant quelle est la performance. Pour les tests, j'ai pris ce qui était à portée de main i5 mobile et une quatrième framboise légèrement overclockée (Raspberry PI 4 8 Go), GCC10 de la distribution Ubuntu 20.04 x64 pour les deux systèmes.







Je ne prétends pas être plus précis dans ces mesures. Des variations de plusieurs dizaines de pour cent sont possibles en fonction de la charge du processeur. La principale conclusion peut être tirée de cette manière. Le passage à l'arithmétique entière ne donne pas de gain de performance sur les processeurs modernes ² . Le nombre inimaginable de transistors dans les processeurs modernes permet d'effectuer rapidement des calculs complexes. Mais, je pense que sur des processeurs tels que Intel Atom, ainsi que sur des contrôleurs faibles, cette approche peut donner un gain de performances significatif. Qu'est-ce que tu penses?



Si cette approche se traduit par un gain de précision, ce gain de précision semble plus intéressant qu'utile. En termes de performances, cette approche peut se retrouver, par exemple, dans l'IoT. Mais pour le calcul haute performance, il n'est plus courant. Dans le monde d'aujourd'hui, SSE / AVX / CUDA préfèrent utiliser le calcul de fonction parallèle. Et en arithmétique en virgule flottante. Il n'y a pas d'analogues parallèles de la fonction MUL. La fonction elle-même est plutôt un hommage à la tradition.



Dans le chapitre suivant, je décrirai comment vous pouvez utiliser efficacement AVX pour les calculs. Encore une fois, allons dans le code libm et essayons de l'améliorer.



¹ Il n'y a pas de registres avec de tels noms dans les processeurs que je connaisse. Les noms ont été choisis par exemple.

²Il faut noter ici que mon ARM est équipé de la dernière version du coprocesseur mathématique. Si le processeur émulait les calculs en virgule flottante, les résultats pourraient être radicalement différents.