👩🏽‍💻 🥢 🤦🏾 Genaue und schnelle Berechnungen für Gleitkommazahlen am Beispiel der Sinusfunktion. Einleitung und Teil 1 🍽️ 🚵 🚧

Lesen Sie sorgfältig sehr gute Artikel aus ArtemKaravaevüber die Hinzufügung von Gleitkommazahlen. Das Thema ist sehr interessant und ich möchte es fortsetzen und anhand von Beispielen zeigen, wie man in der Praxis mit Gleitkommazahlen arbeitet. Nehmen Sie die GNU glibc library (libm) als Referenz. Und damit der Artikel nicht zu langweilig wird, werden wir eine wettbewerbsfähige Komponente hinzufügen: Wir werden versuchen, den Bibliothekscode nicht nur zu wiederholen, sondern auch zu verbessern, um ihn schneller / genauer zu machen.

Ich habe als Beispiel die trigonometrische Sinusfunktion gewählt. Dies ist eine weit verbreitete Funktion, deren Mathematik aus Schule und Universität bekannt ist. Gleichzeitig wird es mit seiner Implementierung viele anschauliche Beispiele für "korrektes" Arbeiten mit Zahlen geben. Ich werde double als Gleitkommazahl verwenden.

In dieser Artikelserie ist viel geplant, von Mathematik über Maschinencodes bis hin zu Compileroptionen. Die Sprache zum Schreiben eines Artikels ist C ++, jedoch ohne "Schnickschnack". Im Gegensatz zur C-Sprache sind Arbeitsbeispiele auch für Personen besser lesbar, die mit der Sprache nicht vertraut sind und weniger Zeilen benötigen.

Die Artikel werden nach der Immersionsmethode geschrieben. Es werden Unteraufgaben besprochen, die dann zu einer einzigen Lösung des Problems zusammengefasst werden.

Zerlegung des Sinus in eine Taylor-Reihe.

Die Sinusfunktion erweitert sich zu einer unendlichen Taylor-Reihe.

\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$

Es ist klar, dass wir keine unendliche Reihe berechnen können, außer in Fällen, in denen es eine analytische Formel für eine unendliche Summe gibt. Dies ist jedoch nicht unser Fall.))) Angenommen, wir möchten den Sinus im Intervall berechnen

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

$[0, \frac{\pi}{2}]$ ... Wir werden die Arbeit mit Intervallen in Teil 3 genauer besprechen

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

$\sin(\frac{\pi}{2})=1$ Schätzen Sie, finden Sie den ersten Begriff, der verworfen werden kann, basierend auf der Bedingung, dass

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

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

e

$e$ Es ist der Unterschied zwischen der Zahl 1 und der kleinsten Zahl, der größer als 1 ist. Grob gesagt ist dies das letzte Bit der Mantisse ( Wiki ). Es ist einfacher, diese Gleichung mit roher Gewalt zu lösen. Zum

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

$e \approx 2.22\times10^{-16}$ ... Ich schaffte es

n = 23

$n=23$ kann bereits fallen gelassen werden. Die richtige Wahl der Anzahl der Begriffe wird in einem der nächsten Teile besprochen, daher werden wir für heute die Begriffe "rückversichern" und übernehmen

n = 25

$n=25$ inklusive.

Die letzte Amtszeit ist ungefähr 10.000 Mal kürzer als

e

$e$ ...

Einfachste Lösung

Hände jucken schon, wir schreiben:

Volltext des Programms zum Testen

#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;
}

Wie man das Programm beschleunigt, denke ich, dass viele Leute sofort herausgefunden haben. Wie oft können Ihre Änderungen Ihrer Meinung nach das Programm beschleunigen? Die optimierte Version und die Antwort auf die Frage unter dem Spoiler.

Optimierte Version des Programms.

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)

Verbesserung der Genauigkeit

Methodik

Die Berechnungsgenauigkeit der Funktion wird durch 2 Standardparameter bestimmt.

Die Standardabweichung von der wahren Sünde (float128) und der Mittelwert dieser Abweichung. Der letzte Parameter kann wichtige Informationen über das Verhalten unserer Funktion geben. Sie kann das Ergebnis systematisch unterschätzen oder überschätzen.

Zusätzlich zu diesen Parametern werden wir zwei weitere einführen. Zusammen mit unserer Funktion nennen wir die in die Bibliothek integrierte sin (double) -Funktion. Wenn die Ergebnisse zweier Funktionen: unserer und der eingebauten nicht übereinstimmen (Stück für Stück), fügen wir der Statistik hinzu, welche der beiden Funktionen weiter vom wahren Wert entfernt ist.

Summationsreihenfolge

Kehren wir zum ursprünglichen Beispiel zurück. Wie können Sie die Genauigkeit "schnell" erhöhen? Diejenigen, die den Artikel sorgfältig gelesen haben Ist es möglich, N Zahlen vom Doppeltyp am genauesten hinzuzufügen? wird höchstwahrscheinlich sofort eine Antwort geben. Es ist notwendig, den Zyklus in die entgegengesetzte Richtung zu drehen. Vom kleinsten zum größten Modulo hinzufügen.

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;
}

Die Ergebnisse sind in der Tabelle gezeigt.

Funktion	Durchschnittlicher Fehler	STD	Unsere ist besser	Bessere 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%

Es mag den Anschein haben, als würde die Verwendung der intelligenten Summierungsalgorithmen den Fehler auf fast 0 beseitigen, aber dies ist nicht der Fall. Natürlich erhöhen diese Algorithmen die Genauigkeit, aber um Fehler vollständig zu beseitigen, sind auch intelligente Multiplikationsalgorithmen erforderlich. Sie existieren, sind aber sehr teuer: Es gibt viele unnötige Operationen. Ihre Verwendung ist hier nicht gerechtfertigt. Wir werden jedoch später in einem anderen Kontext darauf zurückkommen.

Es ist sehr wenig übrig. Kombinieren Sie schnelle und genaue Algorithmen. Kehren wir dazu noch einmal zur Taylor-Serie zurück. Beschränken wir es zum Beispiel auf 4 Mitglieder und nehmen die folgende Transformation vor.

\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!))))$

Sie können die Klammern erweitern und überprüfen, ob der ursprüngliche Ausdruck erhalten wird. Diese Ansicht lässt sich sehr einfach in eine Schleife einfügen.

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;
}

Funktioniert schnell, verliert aber im Vergleich zu e3 an Genauigkeit. Auch hier ist das Runden ein Problem. Schauen wir uns den letzten Schritt der Schleife an und transformieren den ursprünglichen Ausdruck ein wenig.

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

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

Und der entsprechende Code.

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;
}

Die Genauigkeit hat sich im Vergleich zu libm verdoppelt. Wenn Sie erraten können, warum die Genauigkeit gestiegen ist, schreiben Sie in die Kommentare. Darüber hinaus gibt es noch eine weitere, viel unangenehmere Sache an sin_e4, die in sin_e5 in Bezug auf die Genauigkeit fehlt. Versuchen Sie zu erraten, wo das Problem liegt. Im nächsten Teil werde ich auf jeden Fall ausführlich darauf eingehen.

Wenn Ihnen der Artikel gefällt, werde ich Ihnen im nächsten Artikel erklären, wie die GNU libc einen Sinus mit einem maximalen ULP von 0,548 berechnet.

Genaue und schnelle Berechnungen für Gleitkommazahlen am Beispiel der Sinusfunktion. Einleitung und Teil 1

Zerlegung des Sinus in eine Taylor-Reihe.

Einfachste Lösung

Verbesserung der Genauigkeit

Methodik

Summationsreihenfolge

More articles: