🐲 👨🏼‍🤝‍👨🏻 💄 Calculs précis et rapides pour les nombres à virgule flottante en utilisant l'exemple de la fonction sinus. Introduction et partie 1 ⬅️ 🕒 🖨️

Lisez attentivement de très bons articles de ArtemKaravaevsur l'addition de nombres à virgule flottante. Le sujet est très intéressant et je voudrais continuer et montrer avec des exemples comment travailler avec des nombres à virgule flottante dans la pratique. Prenez la bibliothèque GNU glibc (libm) comme référence. Et pour que l'article ne soit pas trop ennuyeux, nous ajouterons un composant compétitif: nous essaierons non seulement de répéter, mais aussi d'améliorer le code de la bibliothèque, le rendant plus rapide / plus précis.

J'ai choisi la fonction sinus trigonométrique comme exemple. Il s'agit d'une fonction répandue dont les mathématiques sont bien connues de l'école et de l'université. En même temps, quand il sera mis en œuvre, il y aura de nombreux exemples frappants de travail «correct» avec des nombres. J'utiliserai le double comme nombre à virgule flottante.

Dans cette série d'articles, beaucoup est prévu, des mathématiques aux codes machine et aux options du compilateur. Le langage d'écriture d'un article est C ++, mais sans «fioritures». Contrairement au langage C, les exemples de travail seront plus lisibles même pour les personnes qui ne connaissent pas le langage et occuperont moins de lignes.

Les articles seront écrits par méthode d'immersion. Les sous-tâches seront discutées, qui seront ensuite réunies dans une solution unique au problème.

Décomposition du sinus en une série de Taylor.

La fonction sinus se développe en une série de Taylor infinie.

\sin (x) = x - \frac{x^{3}}{3!} + \frac{x^{5}}{5!} - \frac{x^{7}}{7!} + \frac{x^{9}}{9!} - \dots

$\sin(x)=x-{\frac {x^{3}}{3!}}+{\frac {x^{5}}{5!}}-{\frac {x^{7}}{7!}}+{\frac {x^{9}}{9!}}-\cdots$

Il est clair que nous ne pouvons pas calculer une série infinie, sauf pour les cas où il existe une formule analytique pour une somme infinie. Mais ce n'est pas notre cas))) Supposons que l'on veuille calculer le sinus dans l'intervalle

[0, \frac{π}{2}]

$[0, \frac{\pi}{2}]$ ... Nous discuterons plus en détail du travail avec intervalles dans la partie 3. Sachant que

\sin (\frac{π}{2}) = 1

$\sin(\frac{\pi}{2})=1$ estimation, trouvez le premier terme qui peut être écarté en fonction de la condition que

\frac{(π / 2)^{n}}{n!} < e

$\frac{(\pi/2)^n}{n!}<e$ où

e

$e$ c'est la différence entre le nombre 1 et le plus petit nombre qui est supérieur à 1. En gros, c'est le dernier bit de la mantisse ( wiki ). Il est plus facile de résoudre cette équation par la force brute. Pour

e \approx 2.22 \times 10^{- 16}

$e \approx 2.22\times10^{-16}$ ... je me suis débrouillé

n = 23

$n=23$ peut déjà être abandonné. Le choix correct du nombre de termes sera discuté dans l'une des parties suivantes, donc pour aujourd'hui, nous allons "réassurer" et prendre les termes jusqu'à

n = 25

$n=25$ compris.

Le dernier terme est environ 10 000 fois moins que

e

$e$ ...

La solution la plus simple

Les mains démangent déjà, nous écrivons:

Texte intégral du programme de test

#include <iostream>
#include <iomanip>
#include <cmath>
#include <array>
#include <bitset>
#include <quadmath.h>
//      clang
//#include "/usr/lib/gcc/x86_64-linux-gnu/10/include/quadmath.h"
#include <numeric>
#include <limits>
#include <vector>

#include <boost/timer/timer.hpp>
#include <boost/math/special_functions/factorials.hpp>

namespace bm = boost::math;

using namespace std;

typedef union { uint32_t i[2]; double x; } mynumber;

array<double, 25> fc;

double sin_e1(double x) {
  double result = 0;
  int sign = 1;
  for(int i = 1; i < 25; i += 2) {
    result += sign * pow(x, i) / bm::unchecked_factorial<double>(i);
    sign = -sign;
  }
  return result;
}

double sin_e2(double x) {
  double result = 0;
  int sign = 1;
  double xx = x * x;
  double pw = x;
  double fti = 1.0;
  for(int i = 1; i < 25; i += 2) {
    fti /= i;
    result += sign * pw * fti;
    fti /= ( i + 1 );
    sign = -sign;
    pw  *= xx;
  }
  return result;
}

double sin_e3(double x) {
  double result = 0;
  for(int i = 25; i >= 1; i -= 2) {
    result += (((i - 1) % 4 == 0) ? 1 : -1 ) * pow(x, i) / bm::unchecked_factorial<double>(i);
  }
  return result;
}

double sin_e4(double x) {
  double xx = x * x;
  double res = fc[25];
  for(int i = 23; i >= 1; i -= 2) {
    res = fc[i] + xx * res;
  }
  return x * res;
}

double sin_e5(double x) {
  double xx = x * x;
  double res = fc[25];
  for(int i = 23; i >= 3; i -= 2) {
    res = fc[i] + xx * res;
  }
  return x + x * xx * res;
}

inline
double fsin(double x) {
  double result;
  asm ("fsin" :"=t" (result) : "0" (x));
  return result;
}

#define SIN(a) fsin(a)
//#define SIN(a) sin(a)
//#define SIN(a) sin_e5(a)
// ^^     . ^^

int main() {

  __uint128_t ft = 1;
  fc[1] = 1.0; //3 * 5;
  for(int i = 2; i < fc.size(); i++) {
    ft *= i;
    // factorial with sign for Taylor series
    fc[i] = (1.0 / ft) * (( (i - 2) % 4 < 2) ? -1 : 1);
  }
  vector<double> xv;
  xv.resize(8 * 2000000);
  //      0  M_PI/2
  for (int i = 0; i < xv.size(); i++) {
    //      .
    xv[i] = (M_PI / 2) * i / double(xv.size());
  }

  double res = 0;
  {
    boost::timer::auto_cpu_timer at;
    for(int i = 0; i < xv.size(); i++) {
      res += SIN(xv[i]);
    }
  }

  int co = 0, cn = 0;
  //      .
  __float128 avg = 0.0, div = 0.0;
  for(int i = 0; i < xv.size(); i++) {
    mynumber dold, dnew;
    dold.x = sin(xv[i]);
    dnew.x = SIN(xv[i]);
    __float128 q = sinq(xv[i]); // <= sinq  .
    __float128 dd = __float128(dnew.x) - q;
    //     .
    div += dd * dd;
    avg += dd;
    //  ,         .
    //   ,         .
    if( dold.i[0] != dnew.i[0] || dold.i[1] != dnew.i[1] ) {
      if( fabsq(q - dold.x) <= fabsq(q - dnew.x) )
        co++;
      else
        cn++;
    }
  }
  avg /= xv.size();
  div /= xv.size();

  cout << res << endl;

  //  ,           .
  cout << "Better libm: " <<  co << " / " << xv.size() << "(" << 100.0 * co / xv.size() << "%)" << endl;

  //  ,  ""         .
  cout << "Better new: " <<  cn << " / " << xv.size() << "(" << 100.0 * cn / xv.size() << "%)" << endl;

  //         .
  cout << "  Avg / std new: " << double(avg) << " / " << double(sqrtq( div - avg * avg )) << endl;
  return 0;
}

double sin_e1(double x) {
  double result = 0;
  int sign = 1;
  for(int i = 1; i < 25; i += 2) {
    result += sign * pow(x, i) / bm::factorial<double>(i);
    sign = -sign;
  }
  return result;
}

Comment accélérer le programme, je pense que beaucoup de gens l'ont tout de suite compris. Combien de fois pensez-vous que vos modifications peuvent accélérer le programme? La version optimisée et la réponse à la question sous le spoiler.

Version optimisée du programme.

double sin_e2(double x) {
  double result = 0;
  int sign = 1;
  double xx = x * x;
  double pw = x;
  double fti = 1.0;
  for(int i = 1; i < 25; i += 2) {
    fti /= i;
    result += sign * pw * fti;
    fti /= ( i + 1 );
    sign = -sign;
    pw  *= xx;
  }
  return result;
}

10000 (GNU C++ v10; -O2)

Améliorer la précision

Méthodologie

La précision de calcul de la fonction sera déterminée par 2 paramètres standard.

L'écart type du vrai péché (float128) et la moyenne de cet écart. Le dernier paramètre peut donner des informations importantes sur le comportement de notre fonction. Elle peut systématiquement sous-estimer ou surestimer le résultat.

En plus de ces paramètres, nous en présenterons deux autres. Avec notre fonction, nous appelons la fonction sin (double) intégrée à la bibliothèque. Si les résultats de deux fonctions: la nôtre et la fonction intégrée ne correspondent pas (bit par bit), alors nous ajoutons aux statistiques laquelle des deux fonctions est la plus éloignée de la valeur vraie.

Ordre de sommation

Revenons à l'exemple d'origine. Comment pouvez-vous augmenter sa précision "rapidement"? Ceux qui ont lu attentivement l'article Est-il possible d'ajouter N nombres de type double le plus précisément? donnera probablement une réponse tout de suite. Il est nécessaire de tourner le cycle dans le sens opposé. Pour ajouter du plus petit modulo au plus grand.

double sin_e3(double x) {
  double result = 0;
  for(int i = 25; i >= 1; i -= 2) {
    result += (((i - 1) % 4 == 0) ? 1 : -1 ) * pow(x, i) / bm::unchecked_factorial<double>(i);
  }
  return result;
}

Les résultats sont présentés dans le tableau.

Fonction	Erreur moyenne	STD	Mieux le nôtre	Meilleure libm
sin_e1	-1,28562e-18	8.25717e-17	0,0588438%	53,5466%
sin_e3	-3.4074e-21	3.39727e-17	0,0423%	10,8049%
sin_e4	8.79046e-18	4.77326e-17	0,0686%	27,6594%
sin_e5	8.78307e-18	3,69995e-17	0,0477062%	13,5105%

Il peut sembler que l'utilisation des algorithmes de sommation intelligente supprimera l'erreur presque à 0, mais ce n'est pas le cas. Bien sûr, ces algorithmes augmenteront la précision, mais pour éliminer complètement les erreurs, des algorithmes de multiplication intelligents sont également nécessaires. Ils existent, mais ils sont très coûteux: il y a beaucoup d'opérations inutiles. Leur utilisation n'est pas justifiée ici. Cependant, nous y reviendrons plus tard dans un contexte différent.

Il en reste très peu. Combinez des algorithmes rapides et précis. Pour ce faire, revenons à nouveau sur la série Taylor. Limitons-le à 4 membres par exemple et faisons la transformation suivante.

\sin (x) \approx x (1 + x^{2} (- 1 / 3! + x^{2} (1 / 5! + x^{2} (- 1 / 7! + x^{2} \cdot 1 / 9!))))

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

Vous pouvez développer les parenthèses et vérifier que l'expression d'origine est obtenue. Cette vue est très facile à insérer dans une boucle.

double sin_e4(double x) {
  double xx = x * x;
  double res = fc[25];
  for(int i = 23; i >= 1; i -= 2) {
    res = fc[i] + xx * res;
  }
  return x * res;
}

Fonctionne rapidement, mais perd en précision par rapport à e3. Encore une fois, l'arrondi est un problème. Regardons la dernière étape de la boucle et transformons un peu l'expression originale.

\sin (x) \approx x + x \cdot x^{2} (- 1 / 3! + \dots))

$\sin(x)\approx x+x \cdot x^2(-1/3!+ \cdots))$

Et le code correspondant.

double sin_e5(double x) {
  double xx = x * x;
  double res = fc[25];
  for(int i = 23; i >= 3; i -= 2) {
    res = fc[i] + xx * res;
  }
  return x + x * xx * res;
}

La précision a doublé par rapport à libm. Si vous pouvez deviner pourquoi la précision a augmenté, écrivez dans les commentaires. De plus, il y a une autre chose, beaucoup plus désagréable à propos de sin_e4, qui manque à sin_e5, liée à la précision. Essayez de deviner quel est le problème. Dans la partie suivante, je vais certainement vous en parler en détail.

Si vous aimez l'article, alors dans le prochain je vous dirai comment la libc GNU calcule un sinus avec un ULP maximum de 0,548.

Calculs précis et rapides pour les nombres à virgule flottante en utilisant l'exemple de la fonction sinus. Introduction et partie 1

Décomposition du sinus en une série de Taylor.

La solution la plus simple

Améliorer la précision

Méthodologie

Ordre de sommation

More articles: