sin() und arccos() : Assembler für die FPU

[helgew]

Zitat von himitsu:

Deine Codeformatierung und vorallem die Einrückung ist grauenhaft.

Da hast du leider recht und von dir höre ich mir das auch gerne an. Als ich das programmiert habe, war ich zudem mit Maple und gnuplot beschäftigt. Da sollte das Produkt nicht drunter leiden, tut es aber, zumindest bei mir und wenn man es aus Spaß macht. Aufräumen kommt dann später

Da ich mir aber zumindest einen kleinen Lernwert bei den vorgestellten Prozeduren verspreche, sollte ich mein posting auch leichter erschließbar und nach allgemein anerkannten Gesichtspunkten hin akzeptabel formatieren. Ich habe die gröbsten Unklarheiten nun umformatiert. für die globalen Variablen, die nur in einer Prozedur verwendet werden und eigentlich unter das implementation - Schüsselwort gehören würden, fällt mir aber nichts saubereres ein.

Zitat von jfheins:

Mich würde das Verfahren dahinter interessieren.

Taylorreihen vom Rand des Intervalls ab sind a-priori eine der schlechtesten Entwicklungsansätze. Ganz gleich, aus welchem Bedürfnis heraus man eine Funktion in eine Reihe entwickelt (sei es stetige Differenzierbarkeit der Approximation, Invertierbarkeit oder Kürzbarkeit mit anderen Termen) ist es immer nützlich, das Gültigkeitsintervall vorher zu kennen. Ein kleines Beispiel, wo es unablässlich ist:

In der komplexen Ebene wird eine Approximation innerhalb eines Quadrates 1 auf i benötigt, die Reihentwicklung der betrachteten Funktion ist aber nur gültig innerhalb eines Kreises mit Radius 1 (man spricht von Konvergenzradius). Legt man den Kreismittepunkt auf eine Seitenmitte, hat man schon ungültige Bereiche, legt man sie auf eine Ecke, sind nurnoch ~78% des benötigten Intervalls abgedeckt.

Ist der Fehler der Approximation in der Zahlenebene etwa ein Paraboloid, würde man den Kreismittelpunkt in die Mitte des Quadrats legen. Ist er asymmetrisch, so schiebt man den steileren Teil weiter zum Rand hin.


Genaugenommen habe den Fehler des Polynoms der arccos-Funktion gegen arccos(x) / sqrt(1-x) (der Hinweis kam von hier) über das Intervall [0;1] integriert und dieses Funktional minimiert. So kommt die 0.42 als Entwicklungspunkt zustande, wobei der Fehler am unteren Ende kritischer ist. Da auch die einseitige Ableitung der arccos-funktion an den Intervallgrenzen divergiert, nimmt man zunächst eine sqrt()-Funktion. Der Quotient ist dann ein sehr langsam veränderliches Polynom, fast eine Gerade.

Eine ebenfalls beliebte Vorgehensweise ist das fitting der Parameter. Für sin_6() habe ich mir 1000 Sinuswerte von 0..pi berechnet und diese um das Maximum bei pi/2 herum polynomial approximiert. Dort fängt man an mit etwas wie "1-x^2". Durch das gerade Verhalten fallen alle ungeraden Potenzen raus, das berücksichtigt man schon im Ansatz.

Eine dritte Methode, die hier keine Anwendung findet, die ich aber auch schon praktiziert habe, ist das Optimieren nach einzelnen Parametern hin, auf das Fehlerfunktional oder den dem Betrage nach größten, auftretenden Fehler (globales Extremum). So kann es schon den Fehler um eine Größenordnung verringern, aus einem Koeffizienten von 1 eine 0.94 zu machen und bei der Berechnung ist das natürlich vollkommen nebensächlich, welchen Wert man dann einsetzt. An dem Punkt ist man auch frei, beliebige Funktionen zur Approximation des Ausgangsproblems heranzuziehen, zum Beispiel Gausskurven für bestimmte, gebrochenrationale Funktionen oder eben diese für das Integral einer Gausskurve, welches man so nicht bestimmen kann.


function RndSingleParity(...) Wenn ich schon beim Erklären bin, noch ein paar Worte zu der wohl kryptischsten Funktion hier:
32bit floats werden IEEE -konform dargestellt als

• 1 Signumsbit (MSB)
• 8 Exponentbits
• 23 bit Mantisse (bis LSB)

00: 0 00000000 00000000000000000000000
01: 0 01111111 00000000000000000000000 ( wenn ich mir das so ansehe, könnte man sich das einfügen der eins sogar sparen... )
02: 0 10000000 00000000000000000000000
03: 0 10000000 10000000000000000000000
04: 0 10000001 00000000000000000000000
05: 0 10000001 01000000000000000000000
06: 0 10000001 10000000000000000000000
07: 0 10000001 11000000000000000000000
08: 0 10000010 00000000000000000000000
09: 0 10000010 00100000000000000000000


Je nach Architektur/Betriebssystem (wer weiß dazu näheres?) kann es u.U. noch Unterschiede geben, ob little / big endian gespeichert wird.
Glücklicherweise ist nach FRNDINT der Single-Wert nach vorzeichenbefreitem Integer übersetzbar und die Parität (gerade Zahlen haben binär eine 0 in der "Einser"-Stelle). Wo man diese Null oder Eins finden kann, sagt der Exponent, nur ist dieser bei 0 : $00 und bei 1 : $7F, bei 2 dann $80 und danach aufsteigend. Das hat rührt von der Interpretation der Kommastelle her und dass das MSB des Oktetts als Exponentenvorzeichen interpretiert wird. Da es maximal 24 Stellen gibt, an denen das Paritätsbit stehen kann, habe ich den Exponenten auf 5 bit getrimmt und wegen der Problematik mit 0, +/-1 um eins erhöht. Da für eins an einer Stelle gelesen wird, die vorher nicht definiert war, muss nach Extrahieren des Exponenten eine 1 gesetzt werden. Das shifting ist zwar nicht gerade elegant und die Schleife wäre besser mit einer schnellen Sprungtabelle zu ersetzen, doch ich weiß noch nicht, wie.
shl / shr sieht weiters nur Konstanten vor, da müsste man schon zur Laufzeit im Maschinencode herumschreiben und das lasse ich an der Stelle.