Genaue und schnelle Berechnungen für Gleitkommazahlen am Beispiel der Sinusfunktion. Teil 3: Festpunkt

Wir setzen den Vorlesungszyklus fort ( Teil 1 und Teil 2 ). In Teil 2 haben wir uns angesehen, was sich in der libm-Bibliothek befindet, und in dieser Arbeit werden wir versuchen, die Funktion do_sin leicht zu ändern, um ihre Genauigkeit und Geschwindigkeit zu erhöhen. Ich werde diese Funktion noch einmal zitieren ( do_sin ):







Wie im vorherigen Artikel, Teil 132-145 gezeigt. Wird für x im Bereich [0,126, 0,855469] ausgeführt. Und was. Versuchen wir, eine Funktion zu schreiben, die innerhalb der vorgegebenen Grenzen genauer und möglicherweise schneller ist.



Die Art und Weise, wie wir es verwenden, ist ziemlich offensichtlich. Die Genauigkeit der Berechnungen muss um weitere Dezimalstellen erweitert werden. Die naheliegende Lösung wäre, den langen Doppeltyp auszuwählen, darin zu zählen und dann zurück zu konvertieren. In Bezug auf die Genauigkeit sollte die Lösung gut sein, aber in Bezug auf die Leistung kann es Probleme geben. Long Double ist jedoch eine ziemlich exotische Art von Daten, und ihre Unterstützung in modernen Prozessoren hat keine Priorität. Unter x86_64 funktionieren SSE / AVX-Anweisungen mit diesem Datentyp nicht. Der mathematische Coprozessor wird "weggeblasen".



Was solltest du dann wählen? Schauen wir uns die Argument- und Funktionsgrenzen genauer an.



Sie befinden sich in der Region 1.0. Jene. Tatsächlich brauchen wir keinen Gleitkomma. Verwenden wir bei der Berechnung der Funktion eine 64-Bit-Ganzzahl. Dies gibt uns zusätzliche 10-11 Bits zur ursprünglichen Genauigkeit. Lassen Sie uns herausfinden, wie man mit diesen Zahlen arbeitet. Eine Zahl in diesem Format wird als a / d dargestellt , wobei a eine Ganzzahl und d ein Divisor ist, den wir für alle Variablen als Konstante auswählen und "in unserem Speicher" und nicht im Speicher des Computers speichern. Nachfolgend sind einige Operationen für solche Nummern aufgeführt:

$$

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

Wie Sie sehen, ist daran nichts Kompliziertes. Die letzte Formel zeigt die Multiplikation mit einer beliebigen Ganzzahl. Beachten Sie auch eine ziemlich offensichtliche Sache, dass das Ergebnis der Multiplikation von zwei vorzeichenlosen ganzzahligen Variablen der Größe N häufiger eine Anzahl von Größen bis zu 2 * N einschließlich ist. Das Hinzufügen kann einen Überlauf von bis zu 1 zusätzlichen Bit verursachen.



Lassen Sie uns versuchen , den Divisor wählen d . In der binären Welt ist es natürlich am besten, sie als Zweierpotenz zu wählen, um nicht zu teilen, sondern nur das Register zu verschieben. Welche Zweierpotenz solltest du wählen? Den Hinweis finden Sie in den Anweisungen der Multiplikationsmaschine. Beispielsweise multipliziert der Standard-MUL-Befehl im x86-System 2 Register und schreibt das Ergebnis auch in 2 Register, wobei 1 der Register der "obere Teil" des Ergebnisses und das zweite der untere Teil ist.



Wenn wir beispielsweise zwei 64-Bit-Nummern haben, ist das Ergebnis eine 128-Bit-Nummer, die in zwei 64-Bit-Register geschrieben wird. Nennen wir RH - "Großbuchstaben" und RL - "Kleinbuchstaben" ¹ . Dann kann das Ergebnis mathematisch wie folgt geschrieben werden$$$R=R_H \cdot 2^{64}+R_L$ . Nun verwenden wir die obigen Formeln und schreiben die Multiplikation für$$$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}}$$

Und es stellt sich heraus, dass das Ergebnis der Multiplikation dieser beiden Festkommazahlen das Register ist $$$R=R_H$ .



Für das Aarch64-System ist es noch einfacher. Der Befehl "UMULH" multipliziert zwei Register und schreibt den "oberen" Teil der Multiplikation in das 3. Register.



Na dann. Wir haben eine Festkommazahl angegeben, aber es gibt immer noch ein Problem. Negative Zahlen. In der Taylor-Reihe geht die Erweiterung mit einem variablen Vorzeichen einher. Um dieses Problem zu lösen, transformieren wir die Formel zur Berechnung des Polynoms nach der Goner-Methode in die folgende Form:

$$

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

Überprüfen Sie, ob es mathematisch genau mit der ursprünglichen Formel übereinstimmt. Aber in jeder Klammer gibt es eine Nummer der Form$$$1/(2n + 1)! - x^2\cdot(\cdots)$ immer positiv. Jene. Durch diese Konvertierung kann der Ausdruck als vorzeichenlose Ganzzahlen ausgewertet werden.

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

Wo $$$tsx[i]=1/i!$im Festkommaformat. Dieses Mal habe ich der Einfachheithalberden gesamten Code auf dem Fast_Sine-Github veröffentlichtund Quadmath aus Gründen der Kompatibilität mit Clang und Arm entfernt. Und ich habe die Methode zur Berechnung des Fehlers ein wenig geändert.



Der Vergleich der Standard-Sinusfunktion und der Festkommafunktion ist in den beiden folgenden Tabellen angegeben. Die erste Tabelle zeigt die Berechnungsgenauigkeit (für x86_64 und ARM ist sie völlig gleich). Die zweite Tabelle ist ein Leistungsvergleich.







Während des Tests wurde der "wahre" Sinuswert unter Verwendung der MPFR- Bibliothek berechnet... Der maximale ULP-Wert wurde als maximale Abweichung vom "wahren" Wert angesehen. Prozentsatz der Fehler - Die Anzahl der Fälle, in denen der berechnete Wert der Sinusfunktion von uns oder von libm binary nicht mit dem auf den doppelten Sinus aufgerundeten Wert übereinstimmte. Der Mittelwert der Abweichung zeigt die "Richtung" des Berechnungsfehlers: Überschätzung oder Unterschätzung des Wertes. Wie Sie der Tabelle entnehmen können, überschätzt unsere Funktion die Sinuswerte. Dies kann behoben werden! Wer hat gesagt, dass die tsx-Werte genau den Koeffizienten der Taylor-Reihe entsprechen sollten? Eine ziemlich offensichtliche Idee bietet sich an, die Werte der Koeffizienten zu variieren, um die Genauigkeit der Approximation zu verbessern und die konstante Komponente des Fehlers zu entfernen. Es ist ziemlich schwierig, eine solche Variation korrekt vorzunehmen. Aber wir können es versuchen. Nehmen wir zum Beispiel4. Wert aus dem Array der tsx-Koeffizienten (tsx [3]) und ändern Sie die letzte Zahl a in f. Lassen Sie uns das Programm neu starten und die Genauigkeit (sin_e7a) sehen. Schau, es ist ein wenig, aber erhöht! Wir fügen diese Methode unserem Sparschwein hinzu.



Nun wollen wir sehen, was die Leistung ist. Zum Testen nahm ich das vorhandene i5 mobile und eine leicht übertaktete vierte Himbeere (Raspberry PI 4 8 GB), GCC10, aus der Ubuntu 20.04 x64-Distribution für beide Systeme.







Ich gebe nicht vor, bei diesen Messungen genauer zu sein. Abhängig von der Prozessorlast sind Abweichungen von mehreren zehn Prozent möglich. Die Hauptschlussfolgerung kann so gezogen werden. Das Umschalten auf Ganzzahlarithmetik führt bei modernen Prozessoren ² nicht zu einem Leistungsgewinn . Die unvorstellbare Anzahl von Transistoren in modernen Prozessoren ermöglicht die schnelle Durchführung komplexer Berechnungen. Ich denke jedoch, dass dieser Ansatz sowohl bei Prozessoren wie Intel Atom als auch bei schwachen Controllern zu einem erheblichen Leistungsgewinn führen kann. Was denken Sie?



Während dieser Ansatz zu einem Genauigkeitsgewinn geführt hat, scheint dieser Genauigkeitsgewinn eher interessant als nützlich zu sein. In Bezug auf die Leistung kann sich dieser Ansatz beispielsweise im IoT befinden. Für Hochleistungsrechner ist es jedoch kein Mainstream mehr. In der heutigen Welt bevorzugen SSE / AVX / CUDA die parallele Funktionsberechnung. Und in Gleitkomma-Arithmetik. Es gibt keine parallelen Analoga der MUL-Funktion. Die Funktion selbst ist eher eine Hommage an die Tradition.



Im nächsten Kapitel werde ich beschreiben, wie Sie AVX effektiv für Berechnungen verwenden können. Lassen Sie uns noch einmal in den libm-Code gehen und versuchen, ihn zu verbessern.



¹ Es gibt keine mir bekannten Register mit solchen Namen in Prozessoren. Die Namen wurden zum Beispiel gewählt.

²Hierbei ist zu beachten, dass mein ARM mit der neuesten Version des Mathe-Coprozessors ausgestattet ist. Wenn der Prozessor Gleitkommaberechnungen emuliert, können die Ergebnisse erheblich abweichen.