Bibliothek – 61 Mathematische Funktionen für FPU iX87

Auflistung · Beschreibung · Zubehör

Links auf zugehörige Dateien + historisches Bild

Library: C-Quelle mit Inline-ASM-Code

Testprogramm als Beigabe

In dem Datenheft wird auch tabellarisch dargelegt, daß ein Koprozessor 150-fach schneller sein kann als eine Emulation in Software. Im Heft wird der i8087 aus den 1970ern beschrieben, der sich bis heute kaum geändert hat, bis auf eine beträchtlich gesteigerte Geschwindigkeit. Das Bild ist vom Browser klein skaliert worden. Es hat real 1037×1344 pixel.

Historie

Es war im Jahr 1991, als eine umfangreiche Funktions-Bibliothek unter dem BS SCO Unix/386 in Assembler für den Gleitkomma-Koprozessor iX87 und für den Assembler MASM von Microsoft (Xenix) entwickelt wurde. Die Assembler-Dateien wurden zu datei.o assembliert, und in ein ar-Archiv aufgenommen, das als Library asm87.a diente. Diese Lösung wurde etwa fünf Jahre lang erfolgreich benutzt, bis sich herausstellte, daß sie kaum portabel ist:

;       sc/24.1.91
        TITLE atan87
        .386
        .387
        .MODEL small
PUBLIC  _atan87
        .DATA
COMM _deg_87:DWORD
        .DATA?
        .CONST
$radtodeg    DT    57.295779513082320876798
ALIGN 4
        .CODE
_atan87 PROC
        fld  QWORD PTR [esp+4]
        fld1
        fpatan
        mov  eax, _deg_87
        cmp  eax, 0
        jg   SHORT $deg
        ret
        ALIGN 4
$deg:
        fld  $radtodeg
        fmul
        ret
        ALIGN 4
_atan87 ENDP
        END

Neufassung

Diese Library von damals wurde vor einigen Wochen wieder aufgegriffen und übersetzt in in C eingebetteten Erweiterten Inline-Assembler, für clang, gcc und bcc32x.exe. Der letztgenannte Compiler läuft unter Windows. Es wurden Erweiterungen und Verbesserungen hinzugefügt. Es sind zum Start 57 Funktionen daraus geworden.

STaTIc long double atan87l(long double x)
{
   __asm__ ("\n\t"
            "fld1 \n\t"
            "fpatan \n\t"
            "cmpl  $0, %[deg] \n\t"
            "jle  ATANe \n\t"
            "fldt  %[rtd] \n\t"
            "fmulp \n\t"
            "ATANe: \n\t"
            : [x]"+t"(x)
            : [rtd]"m"(radtodeg), [deg]"m"(deg_87)
            :
           );
   return x;
}

Diese neue Lösung ist zwar spezifisch für X87 auf X86-Plattform, jedoch herrscht ansonsten volle Portabiltät. Compiler unter allen gängigen Betriebssystemen können die Quelle übersetzen. Der 32bit- und 64bit-Modus werden mit allen Änderungen automatisch berücksichtigt. Die Anbindung an die Funktions-Parameter und den return-Wert erfolgt automatisch. Unterschiede zwischen long double, double und float werden vom Compiler berücksichtigt. Die angegebene Zugriffsgröße bei einigen Instruktionen läßt sich per Makro automatisieren.

Im 64bit-Modus: Der Typ long double wird im Stack-Speicher (16 Byte) an eine Funktion übergeben. Retourniert wird im obersten Gleitkomma-Register st(0). Für den Typ double wird beides in xmm0- (und xmm1)-Register(n) abgewickelt. Im 32bit-Modus werden 12 Byte im Stack übergeben, Return in st(0); etc.

In der vorstehenden Tabelle wurden die meisten Prozessoren ausgelassen. Es wurden nur wichtige Meilensteine aufgeführt.
Den weitaus größten Performance-Sprung brachte der erste Pentium 1993 (mit 60/66 MHz). Die 3 Takte für eine Multiplikation wurden damals fast als Sensation herausgehoben. In neuerer Zeit fallen die Steigerungen der Geschwindigkeit bei der Division und der Quadratwurzel auf. Es kann nicht gesagt werden, daß sich bei der Fortentwicklung der FPU nichts mehr tut. Die FPU in manchen Prozessoren von AMD scheint bemerkbar langsamer zu sein als die in den Prozessoren von Intel. Intel gewinnt noch bei der SingleThread-Performance gegenüber AMD.
Der i8087 enthält noch nicht die Instruktionen fsin und fcos. Sie müssen mit Hilfe der Instruktion fptan und zusätzlichen Formeln, die relativ umfangreich sind, gebildet werden. Das bedeutet ungefähr 1200 Takte für je sin() und cos().
Daten von 1994: Messwerte und aus Datenbüchern

Es wird empfohlen, die Library-C-Datei in eine C-Datei zu inkludieren:

• Die übersetzte Code-Größe beträgt komplett lediglich etwa 2900 Byte. Das ist winzig.
• Wird static (STaTIc, s.u.) gegeben, werden unbenutzte Funktionen gar nicht kompiliert.
• Der benötigte Maschinen-Code ist stets in der Executable vorhanden. Egal, ob diese dynamisch oder statisch gelinkt wurde.
• Es kann nicht passieren, daß bei Ausführung der (mitgenommenen) Executable irgend etwas fehlt.
• Im benutzten Betriebssystem muß nicht das für die Kompilierung verwendete Entwicklungssystem installiert sein.
• Bis auf deg_87 sind diese Funktionen87l thread-sicher.

Eine Überraschung kam auf, die ein paar Tage für Verwirrung sorgte: Es mußte bei diversen Instruktionen die jeweils andere Variante verwendet werden, damit korrekte Funktion vorliegt, gegenüber dem Code von 1991. Und zwar fsub ↔ fsubr, fsubp ↔ fsubrp, fdiv ↔ fdivr und fdivp ↔ fdivrp betreffend, sofern das Ziel nicht st(0) ist. Vor langer Zeit hatte man aus irgendwelchen historischen Gründen die AT&T-Assembler diesbezüglich abweichend von den Definitionen von Intel und AMD verdreht. Es arbeiten nun seitdem alle Compiler (auch die drei oben genannten) und andere Tools wie objdump falsch, und es kann nichts geändert werden - man hängt da drin und traut sich nicht!  Bug-Report

Inhalt der Bibliothek

Beginn Dokument
Diverse Funktionen
Kombinatorik
Koordinaten-Systeme
Trigonometrie
Hyperbolische Funktionen
Konstanten-Funktionen
Kontroll-Funktionen
Potenzreihen
BCD
Portabilität
Zeitstempel-Funktion

STaTIc int deg_87;
static long double const radtodeg= 57.295779513082320876798L;
static long double const degtorad= 1.7453292519943295769e-2L;
static float const zwei= 2.0f, vier= 4.0f;

Vorstehend Datei-globale Variable und Konstanten. Siehe auch weiter unten.
Vor den Funktionen steht das Makro STaTIc, welches leer oder mit static definiert werden kann. Compiler können bei static Objekte weglassen, wenn sie in der sichtbaren Quelle nicht benutzt werden. Der Typ long double wurde durchgehend gewählt, weil die FPU in ihren Registern grundsätzlich mit 80 Bit Breite (10 Byte) rechnet.

Hinter den Funktions-Prototypen steht der Zeitbedarf der jeweiligen Funktion in Takten. Gemessen wurde mit der Instruktion rdtsc auf einem Intel Core 2 Duo 3333 E8600. Ein Takt dauert hier ⅓ ns. Es ist der vollständige Zeitbedarf angegeben, inklusive des Speicherns aller zurückgegebenen Werte. Inlining ist abgeschaltet.

Der verwendete Prozessor (Architektur Wolfdale, 10.08.2008) hat die erste Architektur, bei der Division und Quadratwurzel stark beschleunigt sind, gegenüber früheren Prozessoren. Offensichtlich war hier der Schritt von 65 nm zu 45 nm ausschlaggebend. Auch der Cache wurde vergrößert.

Diverse Funktionen ^

Übergebene Null-Pointer (*xyz) werden von den Funktionen beachtet. Die Zugriffe werden übersprungen oder es wird auf dummy-Zugriffe umgelenkt.

long double log87l(long double b, long double x);  //133-149
long double lg87l(long double x);                  //66
long double ln87l(long double x);                  //66
long double ld87l(long double x);                  //66
long double pow87l(long double x, long double y);  //336-373
long double ilg87l(long double x);                 //197-254
long double iln87l(long double x);                 //197-254
long double ild87l(long double x);                 //190-244

Es sind Logarithmus-Funktionen mit den Basen 10, e und 2 vorhanden, einschließlich eine universelle Funktion mit angebbarer Basis. Zu allen diesen sind die inversen Funktionen vorhanden.
Eine Besonderheit hinsichtlich ihrer Ausgestaltung ist die Funktion pow87l() x^y. Diese reagiert auf die Vorzeichen von x sowie von y. Weiterhin darauf, ob der Exponent y ganzzahlig ist, und falls dies der Fall ist, ob y ungerade ist, wodurch der retournierte Wert im Wechsel negativ oder positiv gesetzt werden kann. Bei negativem y wird der Kehrwert als Resultat errechnet. Bei negativem x wird -π als Fehler anzeigender Wert retourniert, falls y nicht ganzzahlig ist.
Bei x=0.0 wird 0.0 und bei y=0.0 wird 1.0 retourniert.  x^0.5 ist die Quadratwurzel von x (0.5=½). Diese Funktion ist die größte der 61 und hat etwa 230 Byte Maschinen-Code. Sie verwendet die Funktionen rem87l() und irnd87l() inline, bei der Untersuchung des Exponenten y.

Bei den letzten vier Funktionen fällt deren vergleichsweise hoher Zeitbedarf auf. Das ist der Instruktion f2xm1 geschuldet, deren Exponent den Bereich -1.0 … +1.0 erfüllen muß.

long double sqrt87l(long double x);                //21
long double rem87l(long double q, long double d);  //40
long double irnd87l(long double x);                //44

Die Quadratwurzel.
Der Rest einer Division r = q / d  (remainder).
Quasi die Inversion einer Rundung zu einem Integer: Es wird der Wert vor dem Komma entfernt: aus 468.1234 wird 0.1234.

long double rnd87l(long double x);   //43
long double dwn87l(long double x);   //41
long double up87l(long double x);    //41
long double chop87l(long double x);  //43

Es wird x mittels vier verschiedener Rundungsarten gerundet:

Rundung zur nächsten Ganzzahl
Rundung -∞ ←
Rundung → +∞
Rundung → 0 ←

long double abs87l(long double x);  //13
long double chs87l(long double x);  //13

Absolutwert |x|, ein positives Vorzeichen wird gesetzt.
Das Vorzeichen wird komplementiert.

long double scale87l(long double x, int e);    //63
long double xtract87l(long double x, int *e);  //27

y = x · 2^e
Es ist zu beachten, daß bei negativem e der Kehrwert der Potenz gebildet wird.
s;e = x
Der Signifikant s und der Exponent e werden aus dem Parameter x extrahiert.
Der Signifikant wird retourniert, der Exponent auf *e gespeichert.

long double f2pi87l(long double f);                    //15
long double sincos87l(long double x, long double *c);  //163

ω = 2 π f – Die Kreisfrequenz aus der Elektrotechnik.
Retourniert wird der Sinus, der Cosinus wird auf *c gespeichert. Es liegen hier Geschwindigkeitsvorteile vor. Beide Funktionen (fsincos) brauchen insgesamt nur wenig mehr Zeit als eine einzelne.

Kombinatorik ^

long double fak87l(unsigned n);           //20-8835
long double nk87l(int32_t n, int32_t k);  //25-225086

Die Fakultät  6! = 1 · 2 · 3 · … · 6 = 720 ; 0! = 1
max = 1754! = 1.97926189010501006e+4930
Fakultäten werden oft in Reihen verwendet.

n über k = (ⁿ_k) =  n!/[k!·(n-k)!] ; (⁴⁹₆) = 13983816
Zuletzt dargestellt ist die Lotto-Wahrscheinlichkeit »6 aus 49« (//97).
(³²₅) = 201376  ;  2 (⁸₂) (⁴₃) (⁴₂) = 1344 Full House ; Leerkombinationen = 53040
Vorstehend die Situation beim Pokerspiel »5 aus 32« Karten.
Anzahl der verschiedenen Hände, Anzahl der FullHouse-Hände.

Die Berechnungsformel wurde zerlegt, so daß ein Limit für Fakultäten nicht erreicht werden kann.
Resultatwerte, die keine Ganzzahl sind, sind von eingeschränktem Nutzen.

Der Zeitbedarf der Funktionen ist bei normal hohen Argument-Werten ungefähr im Rahmen der trigonometrischen Funktionen.
Bei nk87l(16000,8000) = 1.9046005613914896e+4814 kann das nicht mehr so sein.

Polare und Kartesische Koordinaten ^

long double ptoc87l(int dor, long double grad, long double r, long double *x);
long double ctop87l(long double x, long double y,                      //168
                    long double *rad, long double *deg);               //198

Die Funktionen wandeln polare Koordinaten in kartesische um und umgekehrt.
Das Argument dor zeigt an, ob der Wert des Arguments grad in Grad (deg, degree) oder in Radiant (rad) angegeben ist.
Wenn nicht 'r' und nicht 'R' angegeben wird, gilt deg, andernfalls rad.
Koordinate x wird auf *x geschrieben, retourniert wird y.
Mit einem Null-Pointer x wird nicht zugegriffen.

Die zweite Funktion speichert den Winkel auf *rad und *deg, sie retourniert r.
Mit einem Null-Pointer rad oder deg wird nicht zugegriffen.
Bei diesen Funktionen konnten die Instruktionen fsincos und fpatan in idealer Weise verwendet werden. Bei fpatan sogar die dort eingebaute Division y⁄x.

Polardiagramm
Polare und kartesische Sicht

Trigonometrische Funktionen ^

Die folgenden 12 Funktionen agieren mit dem Wert 2 π oder 360° (rad oder deg).
Um mit Grad zu arbeiten, muß die globale Variable deg_87 auf einen Wert >0 gesetzt werden.

long double sin87l(long double x);  //113
long double cos87l(long double x);  //101
long double tan87l(long double x);  //160
long double cot87l(long double x);  //156
long double sec87l(long double x);  //101
long double csc87l(long double x);  //114

Vorstehend die sechs trigonometrischen Grundfunktionen.
Diese Funktionen sind periodisch, mit wiederkehrenden Werten im Abstand von 2 π.
Weniger bekannt sind die Funktionen sec und csc (Sekans und Kosekans).
 Kurvendiagramm, 6 Kurven
Nachfolgend deren Umkehrfunktionen (invers).
 Kurvendiagramm, 2 Kurven
 Kurvendiagramm, 2 Kurven
 Kurvendiagramm, 2 Kurven

long double asin87l(long double x);  //202
long double acos87l(long double x);  //196
long double atan87l(long double x);  //143
long double acot87l(long double x);  //178
long double asec87l(long double x);  //191
long double acsc87l(long double x);  //203

Hyperbolische Funktionen ^

long double sinh87l(long double x);  //210
long double cosh87l(long double x);  //210
long double tanh87l(long double x);  //208
long double coth87l(long double x);  //215
long double sech87l(long double x);  //219
long double csch87l(long double x);  //219

Die Berechnungsformeln der vorstehenden Funktionen enthalten Potenzen der Zahl e.
 Kurvendiagramm, 3 Kurven
 Kurvendiagramm, 3 Kurven
Nachfolgend deren Umkehrfunktionen (invers).
 Kurvendiagramm, 6 Kurven

long double asinh87l(long double x);  //85
long double acosh87l(long double x);  //85
long double atanh87l(long double x);  //85
long double acoth87l(long double x);  //85
long double asech87l(long double x);  //121
long double acsch87l(long double x);  //122

Die Berechnungsformeln der vorstehenden Funktionen enthalten den Logarithmus zur Basis e.

Konstanten ^

long double pi87l(void);   //7
long double l2t87l(void);  //7
long double l2e87l(void);  //7
long double lg287l(void);  //7
long double ln287l(void);  //7

pi          3.141592653589793239
log₂(10)    3.321928094887362348
log₂(e)     1.442695040888963407
log₁₀(2)    0.3010299956639811952
log_e(2)     0.6931471805599453094
Die Zahl e = iln(1) = 2.718281828459045235
Die zweite Konstante ist als Multiplikant und Divisor geeignet, um zwischen Bit-Stellen und Dezimal-Stellen umzurechnen:  64 ⁄ 3.322 = 19.27 : 19 Dezimalstellen aus 64 Bit.
Die Konstanten haben intern signifikante 66 Bit und werden gerundet in st(0) geschrieben.

Kontroll-Funktionen ^

int xam87l(long double x, char *buf);  //29|60
int top87(void);                       //7
int rndctl87(int rnd);                 //41|9
int prectl87(int prec);                //41|9

Die Prüffunktion (examine) untersucht die Eigenschaften von x und gibt eine Ergebniszahl zurück. Es sind 7 Ergebniszahlen definiert: 0, 1, 2, 3, 4, 5, 6. Bei negativem Vorzeichen von x sind die Ergebniszahlen -1..-6 negativ. Ist der Pointer buf nicht 0, wird eine beschreibende Zeichenkette auf *buf geschrieben. Es müssen mindestens 12 Byte Platz durch *buf zur Verfügung stehen. Die assoziierten Zeichenketten lauten: Unsupported NaN Normal Infinity Zero Empty Denormal.

Die Position des obersten Gleitkomma-Registers st(0) im Stack (TOP) wird retourniert. Eine Position von 0 zeigt einen geleerten Stack an. Wenn eine Funktion, die die FPU benutzt, zurückgekehrt ist und der Rückgabewert abgenommen wurde, muß dieser Wert 0 sein. Ist der Wert beispielsweise 7, befindet sich (noch) ein Wert im Stack. Der Wertebereich beträgt 0..7 mit Wrap-around.

Einer von 4 Rundungsmodi ist mittels der Werte 0,1,2,3 einstellbar.
Ein negativer Wert rnd wird nicht gesetzt.
Der vorherige Wert wird retourniert.
0→round 1→down 2→up 3→chop   (Siehe oben: die vier Rundungs-Funktionen.)

Eine von 3 Berechnungs-Präzisionen ist mittels der Werte 0,2,3 einstellbar.
Ein negativer Wert prec wird nicht gesetzt.
Der vorherige Wert wird retourniert.
Die Präzisionen:  0→24 2→53 3→64  signifikante Bits (float, double, long double).
Diese Einstellung wirkt nur auf die Instruktionen der Grundrechenarten und auf fsqrt.

Potenzreihen ^

long double lnx87l(long double x, unsigned n);                   //n·23+32
long double aex87l(long double lna, long double x, unsigned n);  //n·22+35
ax = aex87l( lnx87l(a, na), x, nx );
ex = aex87l( 1.0L      , 1.0L, 18 ) = 2.718281…;                 //398
ax = aex87l( lnx87l(10.0L, 100), 3.0L, 60 ) = 1000.0;            //3675
     aex87l( lnx87l( 2.0L,  18), 0.5L, 14 ) = 1.4142135623…;     //766

Es wurden zwei Funktionen gemäß entsprechender Potenzreihen entwickelt, die mit den vier Grundrechenarten auskommen. Die Zielfunktionalität ist grundlegend die von pow87l(). Die zweite Funktion benötigt die erste Funktion, um eine Potenz mit beliebiger Basis zu berechnen. Für volle Genauigkeit von 19 signifikanten Stellen sind die Berechnungen von etwa 100 bzw. 60 Gliedern der Reihen notwendig. Dadurch ist der Zeitbedarf etwa 10-fach höher als der von pow87l(). Die erste Funktion konvergiert bei kleineren Werten x, die in der Nähe von 1.0 bleiben, wesentlich besser. Dies ist anhand der beiden letzten Beispiele deutlich erkennbar.

Die Funktionen wurden natürlich für höchstmögliche Effizienz entwickelt. So wird beispielsweise der Gleitkomma-Stack maximal ausgenutzt. Dies nützt aber nicht genügend, denn es ist nicht überraschend, daß Hardware stets wesentlich schneller ist als Software-Lösungen. Und in der FPU werden zudem besondere Verfahren zur Erzielung von hoher Geschwindigkeit genutzt, wie der CORDIC-Algorithmus und Lookup-Tabellen, worauf in dieser Bibliothek verzichtet wurde, weil es ebenso nicht genug Nutzen hätte. Diese beiden Funktionen wurden zwecks Demonstration und für Versuche und Analysen implementiert.

//Takte-Berechnung:
c = nl+k ; d = ml+k ; l = (c-d)/(n-m)
Formeln für erste und zweite Messung, Dauer l für 1×Schleife.
Die Dauer k mußte per Messung mit n=1 festgestellt werden, da die Formel k = c-nl durch die Meßwerte viel zu ungenau ist.

1.414213562373095049 = sqrtl(2.0L);                         //21,054
1.414213562373095049 = sqrt87l(2.0L);                       //21,052
1.414213562373095049 = pow87l(2.0L, 0.5L);                  //346
1.414213562373095049 = aex87l(lnx87l(2.0L, 18), 0.5L, 14);  //766
1.414213562373095048801688724209698 = √2
2.718281828459045235360287471352662 = e

Es ist sehr wahrscheinlich, daß die Funktion sqrtl() aus der C-Lib die FPU per fsqrt benutzt.
Eine Berechnung per Reihe hat Referenzcharakter, da eine mathematische Konvergenz stattfindet und die Resultate bei genügend hohem n korrekt sind, im Rahmen der zur Verfügung stehenden Menge der signifikanten Bits. Es wurden schon Werte √2 gesehen, die in den letzten 4 Stellen ungleich den vorstehenden Werten waren und auch noch mehrfach mehr Zeit für dieses falsche Resultat brauchten. Das liegt wohl daran, daß es scheinbar mancherorts und bei manchen Personen Mode geworden ist, die FPU X87 per Behauptung schlechtzumachen.

4.605170185988091368 = lnx87l(100.0L, 900);  //20785
4.605170185988091369 = lnx87l(100.0L, 950);  //21934
4.605170185988091368 =  ln87l(100.0L);       //66
4.605170185988091368035982909368728
4.61512051684125945  = lnx87l(101.0L, 950);  //21917
4.615120516841259451 =  ln87l(101.0L);       //58
4.615120516841259450884198266912989

Vorstehend ist ein Ergebnis der Potenzreihe um ein Zehntel in der letzten Stelle höher. Höher wird das Ergebnis nicht mehr, bei noch höherem n. Die Ursache für die ±Abweichungs-Zehntel ist nicht klar. Es kann sein, daß komplexe Funktionen der FPU mit weiteren Schutzstellen rechnen. Bei den Konstanten ist zu lesen, daß diese intern in 66 Bit vorliegen. Die Ausgabefunktion printf() könnte mit den End-Bits Probleme haben. Am wahrscheinlichsten ist aber, daß Divisionen in den Gliedern der Reihe zum Schluß fast keine Unterschiede zwischen Dividend und Divisor aufweisen. Dieser dritte Grund und der erste könnten die Ursache sein.

Vorstehend wurden einige besonders lange Konstanten gezeigt, wie sie ein Gleitkomma-Format mit insgesamt 128 Bit aufweist; 16 Byte statt 10 Byte (X87). Die Logarithmus-Reihe ist für Werte oberhalb von ungefähr x=30.0 eigentlich unbrauchbar.

1:  log(a·b) = log(a) + log(b)
2:  c+c+c+c+c = 5·c
3:  4.605170185988091368 =  ln87l(100.0L);                        //68
4:  4.605170185988091368 = lnx87l(100.0L, 900);                   //20784
5:  1.584893192461113485 = pow87l(100.0L, 0.1L);                  //345
6:  4.605170185988091367 = 10·lnx87l(1.584893192461113485L, 14);  //325
7:  1.5625·26            = xtract(100.0)                          //27
8:  4.605170185988091368 = 6·ln2 + lnx87l(1.5625L, 14);
9:  y= xtract87l(100.0L, &v);
    y= aex87l(v·ln287l() + lnx87l(y, 14), 3.0L, 60);              //1716

Das Problem mit der sehr langsamen Konvergenz der Reihe in lnx87l() bei größeren Werten für x ist durch Beachtung einer mathematischen Logarithmus-Regel (1) lösbar, wie vorstehend gezeigt. Eine weitere verwendete Regel zeigt (2). In (3) die normale Bibliotheks-Funktion. Der Zeitbedarf in (4) ist unerträglich. In (5) wird die 10-te Wurzel von 100.0 berechnet. Der Erfolg derer Verwendung ist in (6) zu sehen. Die beste Methode zeigt (7), deren Erfolg in (8) sichtbar ist. In (8) ist ln2 die Konstante der FPU fldln2 (Instruktion).

Die Rückkehr zur allgemeinen Potenz mit Gesamtausdruck in (9). Die Reihe in aex87l() war vor der immensen Verbesserung für lnx87l() deutlich besser als diese. Nun steht sie schlechter da. Sie konvergiert fortschreitend langsamer, je größer die Zahlen werden. Bei 10000³³,n=500 laufen //11468 Takte auf. Der Resultatwert sieht aber etwas besser aus als mit pow87l():
1.000000000000000009e+132 :: 9.999999999999999955e+131
Referenzwert:  10000³³ = 10^4·33 = 1·10¹³²
Eine Verbesserung für die aex-Reihe ist aktuell nicht vor dem Horizont.

Die aex-Reihe:  a^x = 1 + [x·ln(a)] + [x·ln(a)]²/2! + [x·ln(a)]³/3! + ...
Für e^x muß ln(a) entfernt werden.  Und: e = Σ x=1 ; 1 = Σ x=0
Die Potenz und die Fakultät kämpfen bei der Division gegeneinander. Die Glieder der Reihe müssen für Konvergenz im Wert fortlaufend kleiner werden. Wegen des Maximums der Fakultät ist n auf 1750 gedeckelt.

BCD (dezimal) ^

int bcd87(int ctl, void *val);  //149

Die FPU kann mit 18-stelligen BCD-Integer-Objekten umgehen.
Diese sind 80 Bit breit, wovon die unteren 9 Byte die 18 BCD-Nibbles enthalten.
Das oberste zehnte Byte enthält in Bit 7 das Vorzeichen.

Die bcd-Funktion nimmt die Zeichen 's', 'l' oder 'L' in ctl entgegen.
Die Zeichen bedeuten: BCD-Speichern, BCD-Laden, Gleitkomma-Laden.
Beim Speichern schreibt die Funktion auf *val, beim Laden liest sie von *val.
Beim Speichern 's' wird der oberste Wert im Gleitkomma-Stack entfernt.
Parameter *val muß auf mindestens 10 Byte zeigen.

Bei ungültigem Inhalt von ctl wird -2 retourniert.
Bei 's' und leerem Stack wird -1 retourniert.
Ohne Fehler wird 0 retourniert.
Falls val ein 0-Pointer ist, wird auf *val nicht zugegriffen.

static int cpbcd(char *s, long double bcd)
{
   char *b= (char*)&bcd;
   int i, k;
   s[0]= b[9]&0x80 ? '-' : '+';
   for (k=1,i=8;  i>=0;  --i,k+=2)  {
      s[k+0]= (b[i]>>4)+'0';
      s[k+1]= (b[i]&15)+'0';
   }
   s[k]= 0;
   return 0;
}

Die vorstehende Funktion kopiert ein BCD-Objekt in einen 0-terminierten String.
Diese Funktion ist in dem Test-C-Programm enthalten, das auch ganz oben verlinkt ist.

Es gibt eine Reihe von x86-Instruktionen, außerhalb der FPU-Instruktionen, die BCD-Operationen an einzelnen Bytes oder Nibbles vornehmen können.

Portabilität ^

Durch die Verwendung dieser Bibliothek erfolgt zwar eine Festlegung auf eine X86-Plattform mit X87 - mehr aber nicht. Das bedeutet u. a., es herrscht Little Endian. Durch die Verwendung des Erweiterten Inline-Assemblers entsteht eine hohe sonstige Portabilität.

• Compiler unter allen gängigen Betriebssystemen können die Quelle übersetzen.
• Der 32bit- und 64bit-Modus werden mit allen Änderungen automatisch berücksichtigt.
• Die Anbindung an die Funktions-Parameter und den return-Wert erfolgt automatisch.
• Unterschiede zwischen long double, double und float werden vom Compiler berücksichtigt.
• Die angegebene Zugriffsgröße bei einigen Instruktionen läßt sich per Makro automatisieren.

Eine Portabilität mit dem 16bit-Modus wurde nicht explizit hergestellt, weil das heutzutage einfach zu weit weg wäre. Prinzipiell wäre eine solche Portabilität möglich. Der damals so genannte NDP 8087 konnte nahezu das, was er heute kann. Heute als FPU nur wesentlich schneller.

^

uint64_t rdtsc(void)
{
   uint32_t d, a;
   __asm__ volatile ("\n\t"
            "lfence \n\t"
            "rdtsc \n\t"
            : "=a"(a), "=d"(d)
            :
            :
           );
   return (uint64_t)d<<32|(uint64_t)a;
}

;64-Bit-Modus:
	pushq	%rbp
	movq	%rsp, %rbp
	#APP
	lfence
	rdtsc
	#NO_APP
	shlq	$32, %rdx
	movl	%eax, %eax
	orq	%rdx, %rax
	popq	%rbp
	retq

;32-Bit-Modus:
	pushl	%ebp
	movl	%esp, %ebp
	#APP
	lfence
	rdtsc
	#NO_APP
	popl	%ebp
	retl

Vorstehend eine C-Quelle mit Inline-Assembler und zwei Compiler-Ausgaben dazu.
Die Instruktion rdtsc setzt die Register edx:eax. Dies wird dem Compiler nach : mitgeteilt. Im 32-Bit-Modus werden diese Register bei 64-Bit-return per Konvention benutzt, im 64-Bit-Modus wird rax benutzt, weshalb der Compiler rax entsprechend befüllt. Die mov-Instruktion ist nicht sinnlos, denn sie löscht die obere Hälfte von rax. Den Inhalt von __asm__(...); grenzt der Compiler in seiner Ausgabe durch #APP und #NO_APP ab.

Die folgenden Werte vom tsc wurden mit folgendem C-Code gewonnen:
  tsc= rdtsc(); sleep(1); tsc= rdtsc()-tsc;
3337607930 3343052710 3366005900 3350808460 3344029890 3362758630 3339670120
Der Mittelwert ist 3349133377. Die größte Abweichung beträgt etwa 0,5 %. Das ist sehr gering im Rahmen von Zeitdauern in Takten.

long double asin87l(long double x)
{
   __asm__ ("\n\t"
            "fld  %%st(0) \n\t"
            "fmul  %%st(0), %%st(0) \n\t"
            "fld1 \n\t"
             fsubrp "\n\t"
            "fsqrt \n\t"
             fdivp "\n\t"
            "fld1 \n\t"
            "fpatan \n\t"
            "cmpl  $0, %[deg] \n\t"
            "jle  ASINe \n\t"
            "fldt  %[rtd] \n\t"
            "fmulp \n\t"
            "ASINe: \n\t"
            : [x]"+t"(x)
            : [rtd]"m"(radtodeg), [deg]"m"(deg_87)
            :
           );
   return x;
}

Compiler-Ausgabe:

	fldt	16(%rbp)
        #APP
        ;...
	cmpl	$0, deg_87(%rip)
	jle	ASINe
	fldt	radtodeg(%rip)
        ;...
	#NO_APP
	popq	%rbp
	retq

für double statt long double:

	movsd	%xmm0, -16(%rbp)
	fldl	-16(%rbp)
        #APP
        ;...
	#NO_APP
	fstpl	-8(%rbp)
	movsd	-8(%rbp), %xmm0
	popq	%rbp
	retq

Der Compiler verbindet auch alle Funktionsparameter, den retournierten Wert und eventuell vorhandene lokale Objekte, jeweils mit passender Syntax. Er setzt für jede Sequenz %[name] den passenden Zugriff ein. Auch vor #APP und nach #NO_APP handelt er automatisch, so daß die Situation recht portabel gestaltet ist. Vor #APP lädt der Compiler wegen "+t"(x) den Parameter x auf den Gleitkomma-Stack. Und er beläßt ihn nach #NO_APP als return-Wert auf dem Gleitkomma-Stack.

Nach Änderung des Typs zu double reagiert der Compiler entsprechend den Konventionen. Er leert den Stack (fstp) und speichert den return-Wert in das xmm0-Register. Der Quellcode bleibt dazu vollkommen unverändert. Das alles ist eben portabel. Nichts davon war in der alten Ausgabe der Bibliothek von 1991 portabel.