Juni 2022

Bibliothek – 59 Mathematische Funktionen mit FPU iX87

Auflistung · Beschreibung · Zubehör


Links auf zugehörige Dateien + historisches Bild

Library: C-Quelle mit Inline-ASM-Code

Testprogramm als Beigabe

intel
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 in datei.o assembliert, und in ein ar-Archiv aufgenommen, das als Library 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 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 iX87 auf Intel-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.

Es wird empfohlen, die Library-C-Datei in eine C-Datei zu inkludieren. Die übersetzte Code-Größe beträgt lediglich etwa 2700 Byte. Das ist winzig.

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

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.

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);
long double lg87l(long double x);
long double ln87l(long double x);
long double ld87l(long double x);
long double pow87l(long double x, long double y);
long double ilg87l(long double x);
long double iln87l(long double x);
long double ild87l(long double x);

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() xy. 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 y=0.0 wird 1.0 retourniert. Bei x=0.0 wird 0.0 retourniert.  x0.5 ist die Quadratwurzel von x (0.5= ½). Diese Funktion ist die größte der 59 und hat etwa 230 Byte Maschinen-Code. Sie verwendet die Funktionen rem87l() und irnd87l() inline, bei der Untersuchung des Exponenten y.

long double sqrt87l(long double x);
long double rem87l(long double q, long double d);
long double irnd87l(long double x);

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);
long double dwn87l(long double x);
long double up87l(long double x);
long double chop87l(long double x);

Es wird x mittels vier verschiedener Rundungsarten gerundet:

5/4-Rundung (zur nächsten Ganzzahl)
Rundung → +∞
Rundung -∞ ←
Rundung → 0 ←

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

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

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

y = x · 2e
Es ist zu beachten, daß bei negativem e der Kehrwert der Potenz gebildet wird.
s = x ; 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);
long double sincos87l(long double x, long double *c);

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

Kombinatorik

long double fak87l(unsigned n);
long double nk87l(int32_t n, int32_t k);

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 = (nk) =  n!/[k!·(n-k)!] ; (496) = 13983816
Zuletzt dargestellt ist die Lotto-Wahrscheinlichkeit »6 aus 49«.
(325) = 201376  ;  2 (82) (43) (42) = 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.

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,
                    long double *rad, long double *deg);

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);
long double cos87l(long double x);
long double tan87l(long double x);
long double cot87l(long double x);
long double sec87l(long double x);
long double csc87l(long double x);

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);
long double acos87l(long double x);
long double atan87l(long double x);
long double acot87l(long double x);
long double asec87l(long double x);
long double acsc87l(long double x);

Hyperbolische Funktionen

long double sinh87l(long double x);
long double cosh87l(long double x);
long double tanh87l(long double x);
long double coth87l(long double x);
long double sech87l(long double x);
long double csch87l(long double x);

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);
long double acosh87l(long double x);
long double atanh87l(long double x);
long double acoth87l(long double x);
long double asech87l(long double x);
long double acsch87l(long double x);

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

Konstanten

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

pi          3.14159265358979324
log2(10)    3.32192809488736235
log2(e)     1.44269504088896341
log10(2)    0.301029995663981195
loge(2)     0.693147180559945309
Die Zahl e = iln(1) = 2.71828182845904523
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.

Kontroll-Funktionen

int xam87l(long double x, char *buf);
int top87(void);
int rndctl87(int rnd);
int prectl87(int prec);

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:  024253364  signifikante Bits (float, double, long double).
Diese Einstellung wirkt nur auf die Instruktionen der Grundrechenarten und auf fsqrt.

BCD (dezimal)

int bcd87(int ctl, void *val);

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.



Copyright © Helmut Schellong - 2022

Coprocessor Koprozessor FPU NPU iX87 Library Mathematische Funktionen