Genaue und schnelle Berechnungen für Gleitkommazahlen am Beispiel der Sinusfunktion. Teil 2: libm

Ich setze die Artikelserie über die Arbeit mit Gleitkomma fort. Im ersten Artikel gab ich eine kleine mathematische Einführung und zeigte die einfachste und naheliegendste Methode zur Berechnung des Sinus anhand von Beispielen für Programme mit unterschiedlichen "Fallstricken". Der heutige Artikel wird etwas anders aussehen. Hier wird es keine Übung geben, aber wir werden tiefer in die Mathematik eintauchen und uns mit dem Allerheiligsten befassen - dem Standard-Bibliothekscode. Ich werde auch am Ende des ersten Artikels eine Antwort auf die Frage geben. So lass uns gehen.

Wieder etwas Mathe

Je schneller der Absolutwert der Taylor-Reihe abnimmt, desto weniger Terme werden offensichtlich benötigt, um die erforderliche Genauigkeit zu erreichen. Und so scheint es, dass das Ergebnis genauer sein wird (im Folgenden wird dies ausführlicher besprochen). Nehmen Sie zum Vergleich beispielsweise einen Term des siebten Grades der Taylor-Reihe ($$$\frac{x^7}{7!}$) beim $$$x=1.0$ und $$$x=0.1$... Die Ausdruckswerte sind$$$1.198\cdot 10^{-4}$ und $$$1.198\cdot 10^{-11}$beziehungsweise. Ein großer Unterschied, nicht wahr? Versuchen wir also, einen Weg zu finden, um die Obergrenze des Berechnungsintervalls für die Sinusfunktion zu verringern.

Serienerweiterung um vorgegebene Werte

Um diese Methode zu verstehen, müssen wir zum ersten Jahr des Instituts zurückkehren und uns an die Definitionen der Taylor-Reihe ( Wiki ) erinnern . Kurz gesagt: Wenn Sie die Funktion und ihre Ableitungen irgendwann kennen, können Sie die Werte der Funktion in der Nähe dieses Punkts finden, indem Sie sie in eine Taylor-Reihe erweitern. Für die Sinusfunktion bedeutet dies Folgendes

$$

$$\sin(x_0+\Delta x) \approx \sin(x_0)+\sin'(x_0)\frac{\Delta x}{1!}+\sin''(x_0)\frac{\Delta x^2}{2!}+\sin'''(x_0)\frac{\Delta x^3}{3!}+\cdots$$

Was gibt uns dieser Ansatz aus praktischer Sicht? Stellen Sie sich vor, wir haben ein Intervall von$$$0$ Vor $$$\pi/2$... Wählen wir 10 linear verteilte Punkte in diesem Intervall (die Auswahl ist nicht optimal):$$$x_0=0$, $$$x_1=\pi/20$, $$$x_2=2\pi/20$, $$$x_i=i\cdot \pi/20$... Berechnen Sie für jeden Punkt die Platte mit dem Sinus und seinen Ableitungen an diesem Punkt. Jetzt können Sie die Funktion so ändern, dass beim Abrufen des Wertes$$$x$ Die Funktion nimmt den nächstgelegenen Wert an $$$x_i$ und liegt in einer Reihe um den Punkt herum $$$x_i$, nicht um Null ($$$\Delta x = x - x_i$).

Trigonometrische Transformationen verwenden

Wenn wir noch weiter zurückgehen, zu den höheren Klassen der Schule, dann können wir uns an eine sehr wichtige Formel erinnern:

$$

$$\sin(x_0+\Delta x)=\sin(x_0)\cos(\Delta x)+\cos(x_0)\sin(\Delta x)$$

Und dann ist alles das gleiche wie im vorherigen Absatz. Wir wählen Punkte innerhalb des Intervalls aus, berechnen den Sinus und den Cosinus für sie, und wenn wir die Sinusfunktion aufrufen, suchen wir nach dem nächsten und berechnen den Sinus anhand der obigen Formel mit einem kleinen Wert$$$\Delta x$...

Überlegen Sie, welche dieser beiden Methoden besser zu wählen ist, aber jetzt werden wir von der Mathematik zu praktischen Berechnungen übergehen.

Verteilungseigenschaft der Multiplikation in der Gleitkommawelt

Ich musste das Internet um Rat fragen, wie es heißt $$$a\cdot(b+c)=a\cdot b + a\cdot c$... Es stellt sich heraus, dass es sich um eine Verteilungseigenschaft handelt. Kehren wir zu der Frage zurück, die ich am Ende des ersten Teils gestellt habe. Nämlich warum mathematisch äquivalente Ausdrücke$$$a\cdot(b+c)$ und $$$a\cdot b + a\cdot c$kann bei Gleitkommaberechnungen unterschiedliche Ergebnisse liefern? Dies lässt sich am einfachsten anhand eines Beispiels veranschaulichen. Nehmen wir ein hypothetisches System, das mit Gleitkommazahlen im Dezimalformat mit 4 Stellen Genauigkeit arbeitet. Stellen wir uns das vor$$$b=1.000\text{E}0$, $$$c=1.234\text{E-}2$, und $$$a=5.678\text{E-}2$... Nehmen wir zunächst einen Ausdruck in Klammern und berechnen ihn Schritt für Schritt. Denken Sie daran, bei jedem Schritt abzurunden:

1)$$$b+c=1.000\text{E}0 + 1.234\text{E-}2=1.012\text{E}0$

2) $$$a\cdot(b+c)=5.678\text{E-}2\times 1.012\text{E}0=5.746\text{E-}2$

Antwort erhalten $$$y=5.746\text{E-}2$

Berechnen wir nun Schritt für Schritt den zweiten Ausdruck auf die gleiche Weise:

$$$a\cdot b + a\cdot c=5.678\text{E-}2\times 1.000\text{E}0 + 5.678\text{E-}2\times 1.234\text{E-}2=5.678\text{E-}2 + 7.007\text{E-}4=5.748\text{E-}2$

Antwort erhalten $$$y=5.748\text{E-}2$

Die wahre Antwort lautet 0.0574806652.



Wie Sie sehen können, ist die Antwort im zweiten Fall der wahren viel näher als im ersten. Wenn wir dies an den Fingern erklären, stellen Sie sich vor, dass wir im ersten Fall die Zahl zu 1.0 addieren$$$c=1.234\text{E-}2$Wir werfen nur die letzten beiden Ziffern weg. Sie sind nicht mehr. Im zweiten Fall erfolgt das Verwerfen ganz am Ende nach der Multiplikation. Jene. im zweiten Fall ist die Multiplikationsoperation (en) genauer.



Es scheint, dass Sie damit fertig werden können, aber schauen Sie sich die erste Methode genauer an und sagen Sie mir, was das Ergebnis der Berechnung sein wird$$$b+c-b$... Und ... wir haben eine Möglichkeit, Gleitkommazahlen zu runden! Verpassen Sie dieses Beispiel nicht. Geben Sie sich Zeit, um es herauszufinden. Das Runden von Zahlen wird später in diesem und den folgenden Artikeln von uns sehr intensiv genutzt.



Beachten wir noch ein Merkmal dieses Ausdrucks. Stellen Sie sich vor, die 4-stellige Genauigkeit einer Variablen reicht uns nicht aus. Was zu tun ist? Und hier haben wir bereits die Antwort - die Zahl in der Form darzustellen$$$b+c$und speichern Sie es im Speicher als die Summe von zwei Ziffern. Führen Sie dementsprechend Operationen (z. B. Multiplikation) für beide Begriffe getrennt aus. Diese Technik wird im Artikel Hinzufügen von zwei Gleitkommazahlen ohne Genauigkeitsverlust ausführlicher beschrieben .



Im vorigen Artikel habe ich auch geschrieben, dass die Methode$$$a\cdot(b+c)$Es gibt eine unangenehme Eigenschaft. Und es ist wie folgt. Nummer$$$c$ immer an der letzten signifikanten Stelle einer Zahl abgeschnitten $$$b$... Dies bedeutet, dass unabhängig von der Anzahl$$$c$, wenn ein $$$b+c \neq b$, dann ist ein Fehler im letzten Vorzeichen auch bei kleinen immer möglich $$$c$... Dies ist im Ansatz im nächsten Kapitel nicht zulässig.

Wie es am Beispiel der GNU-Bibliothek funktioniert

Und wie? Haben Sie gewählt, welche der beiden am Anfang des Artikels beschriebenen Methoden Sie für die genaue Berechnung des Sinus gewählt haben? Welche Methode Sie auch wählen, beide sind korrekt. Darüber hinaus sind sie absolut identisch. Glauben Sie mir, probieren Sie es aus. Im Folgenden werde ich Schulformeln verwenden. Sie sind leichter zu erklären.



Mit den im vorherigen Artikel und in diesem Artikel gewonnenen Erkenntnissen können Sie den Code der Standardbibliothek leicht verstehen. Öffnen wir die Datei s_sin.c und suchen dort die Funktion __sin :



Der Code ist recht einfach. Es ist leicht zu verstehen, dass es abhängig von den Grenzen der Eingabevariablen unterschiedliche Funktionen aufruft. In diesem Artikel werden wir den Codeabschnitt 218-224 für die Winkel 2 ^ -26 <| x | <0,855469 diskutieren. Sie können sehen, dass in diesem Abschnitt des Codes die Funktion do_sin (x, 0) aufgerufen wird. Wir werden auf diese Funktion näher eingehen:

1. , dx=0 .
2. 129-130 , abs(x)<0.126, .. x , . , , , .
3. 136-137. , . x 2 . u x. , 0.345678. u=0.34, 0.005678.
4. 140-142. ( s ) ( c ) x . , cos(x)=1-c, 1.0, (. ), .
5. 143. u. , u=0.34 34. sin(u)=sn+ssn, cos(u)=cs+ccs. sn cs — «» u, ssn ccs — .
6. 144-145. sin(u+x)=(sn+ssn)*(1-c)+(cs+ccs)*s. , , 144-145. — .

Tatsächlich habe ich nur den einfachsten Teil der Berechnung des Sinus auf diese Weise beschrieben. Es bleibt viel Mathematik zurück. Wie berechnen Sie beispielsweise die Größe einer Tabelle und der darin enthaltenen Elemente? Woher kommen die magischen Zahlen 0.126 und 0.855469? Wann sollte die Berechnung anhand der Taylor-Zahl abgebrochen werden? Korrekturen an den Koeffizienten der Taylor-Reihe, um das Ergebnis zu verfeinern.



All dies ist natürlich interessant, aber objektiv hat die vorgestellte Methode viele Nachteile: Es ist notwendig, Sinus (e) und Cosinus (c) gleichzeitig zu berechnen, was doppelt so viele Berechnungen der Taylor-Reihe ¹ erfordert . Wie wir sehen können, ist die Multiplikation mit Tabellenwerten ebenfalls nicht kostenlos. Das Speichern einer Tabelle mit 3520 Bytes im RAM ist natürlich kein Problem, aber der Zugriff darauf (auch im Cache) kann teuer sein.



Daher werden wir im nächsten Teil versuchen, die Platte loszuwerden und den Sinus im Intervall [0.126, 0.855469] direkt, aber genauer als im ersten Kapitel zu berechnen.



Vor dem Ende - eine Frage des schnellen Verstandes. Die große Zahl in diesem Beispiel ist 52776558133248 = 3 * 2 ⁴⁴ . Woher kam eine solche Zahl, nicht zum Beispiel 2 ⁴⁵ ? Ich werde die Frage genauer formulieren. Warum ist die Zahl 3 * 2 ^N optimal, wenn Zahlen gerundet werden , und nicht beispielsweise 2 ^{N + 1} ? Eine andere Frage, welches N sollten Sie wählen, um eine Zahl auf eine ganze Zahl zu runden?



¹ Es ist zu beachten, dass ein wesentlicher Vorteil dieses Ansatzes auftreten kann, wenn Sinus und Cosinus gleichzeitig aus demselben Winkel berechnet werden. Die zweite Funktion kann fast kostenlos berechnet werden.