1. Einführung
Gesucht ist eine Funktion zur Berechnung einer Exponentialfunktion in Ganzzahl-Arithmetik.
Als Basis bietet statt der Eulerkonstante e der Wert 2 an, um mit der Exponentatialfunktion mit Basis 2 und Exponenten n ∊ ℕ₀ kompatibel zu sein, die trivial durch den Shift-Operator implementiert ist. Wir können deshalb den Argumentbereich auf 0…1 ohne die 1 beschränken, also auf das oben offene Intervall x ∊ [0…1[.
Wenn Eingabe, alle Zwischenergebnisse und das Ergebnis in diesem offenen Intervall liegen,
können wir dieses auf den C
-Datentyp uint32_t
, abbilden,
indem wir die Werte mit 232 multiplizieren.
Einem mathematischen Wert x entspricht dann der C
-Wert
.
Leider liegt Wertebereich unserer Funktion außerhalb des erwünschten Bereiches. Deshalb berechnen wir stattderer die Funktion , deren Wertebereich genau unser gewünschtes Intervall ist.
2. Mathematik
Aufgabe ist die Berechnung von für Eingaben , wobei auch alle Zwischenergebnisse im Intervall liegen müssen.
2.1. Vorbereitung für die Reihenentwicklung
Um zur Berechnung eine Potenzreihe nutzen zu können, wechseln wir zur Basis e:
Die Konstante liegt im erlaubten Bereich, und es gilt .
Damit die Potenzreihe schnell konvergiert, zerlegen wir z:
und nutzen das Potenzgesetz der Addition:
y kann nur die drei Werte 0, ¼ und ½ annehmen; deshalb können die konstanten Werte von e0, e¼ und e½ vorab berechnet und in einer Tabelle bereitgestellt werden.
Weil y und x nicht negativ sind, sind die Potenzen ey und x größer oder gleich 1 und liegen damit außerhalb des erlaubten Bereiches. Deshalb definieren wir:
Wie man leicht nachrechnen kann, liegen die Werte von g1(y) und g2(x) zwischen 0 und 1. Damit können wir g(z) innerhalb des erlaubten Zahlenbereiches berechnen:
Die Exponentialfunktion zur Basis e kann als Potenzreihe berechnet werden:
Dabei berücksichtigen wir aber nur die Terme bis zu einem Grad N:
2.2. Horner-Schema
Die Potenzreihe lässt sich nach dem Horner-Schema umformen, um Rechenoperationen zu sparen und Rundungsfehler zu minimieren. Wir benutzen im Beispiel N = 4:
Wir klammern nacheinander x, x/2, x/3 usw. aus:
Leider liegen wegen der Additionen der 1 die Zwischenergebnisse außerhalb des uns erlaubten Bereiches. Deshalb nutzen wir ein modifiziertes Horner-Schema.
2.3. Modifiziertes Horner-Schema
Wir klammern nacheinander x/2, x/3 usw. aus, lassen aber jeweils das führende x stehen:
Damit ergibt sich diese Berechnungsfolge, bei der für x<loge2 alle Zwischenergebnisse hj kleiner als 1 sind:
2.4. Vollständiger Algorithmus
Einige Konstanten werden vorab berechnet:
Berechnung von f(a):
3. Implementierung
Wir erinnern uns daran, dass die in C
genutzte skalierte Zahl um den Faktor 232
größer ist als die wirkliche
mathematische Zahl:
Die Ganzzahl-Addition auf den skalierten Zahlen ist identisch zur mathematischen Addition:
Die Ganzzahl-Multiplikation auf den skalierten Zahlen liefert ein Ergebnis, das um den Faktor 232 zu groß ist:
Deshalb muss:
- eine 32-Bit-Multiplikation mit 64-Bit-Ergebnis genutzt werden, und
- das 64-Bit-Ergebnis durch 232 geteilt werden, was durch einen Rechts-Shift um 32 Bits realisiert werden kann.
Da es in C
keine 32-Bit-Multiplikation mit 64-Bit-Ergebnis gibt,
müssen wir eine 64-Bit-Multiplikation und einen 64-Bit-Shift um 32 Stellen codieren (u32mulhi
).
Zum Glück sind aktuelle Compiler schlau genug, aus der Sequenz von 64-Bit-Multiplikation und
64-Bit-Shift eine 32-Bit-Multiplikation zu generieren, deren oberes Ergebniswort genutzt wird.
//-------------------------------------------------------------------------------------- // Auxiliary multiplication function with impicite division by 2^32 //-------------------------------------------------------------------------------------- static inline constexpr uint32_t u32mulhi( const uint32_t a, const uint32_t b ){ return (uint32_t)( ((uint64_t) a * b) >> 32 ); } //-------------------------------------------------------------------------------------- // Constants //-------------------------------------------------------------------------------------- // 2**32 * log(2) - unfortunately log(const) is not a constexpr static constexpr uint32_t LOG_2 = 0xb17217f7u; // two top bits are handled by a lookup table, // resulting in series expansion input to be less then 0.25 static constexpr unsigned LookupBits = 2; static constexpr unsigned WordLength = 32; static constexpr unsigned ShiftValue = WordLength - LookupBits; // ( exp( 0/4 ) - 1 ) * 2**32 = 0x.0000.0000 | unfortunately exp() is not a constexpr static constexpr uint32_t EXP_0_BY_4_MINUS_1 = 0; // ( exp( 1/4 ) - 1 ) * 2**32 = 0x.48b5.e3c3 | unfortunately exp() is not a constexpr static constexpr uint32_t EXP_1_BY_4_MINUS_1 = 0x48b5e3c3u; // ( exp( 2/4 ) - 1 ) * 2**32 = 0x.a612.98e1 | unfortunately exp() is not a constexpr static constexpr uint32_t EXP_2_BY_4_MINUS_1 = 0xa61298e1u; // precomputed reciprocal values static constexpr uint32_t ONE_BY_3 = (uint32_t)( ( 1ull << 32 ) / 3 ); static constexpr uint32_t ONE_BY_5 = (uint32_t)( ( 1ull << 32 ) / 5 ); static constexpr uint32_t ONE_BY_6 = (uint32_t)( ( 1ull << 32 ) / 6 ); static constexpr uint32_t ONE_BY_7 = (uint32_t)( ( 1ull << 32 ) / 7 ); //-------------------------------------------------------------------------------------- // main function //-------------------------------------------------------------------------------------- // Integer computation of u32bexpm1(x) := 2**32 * ( 2 ** (x/2**32) - 1 ) //-------------------------------------------------------------------------------------- uint32_t u32p32exp2m1( uint32_t arg ){ //------------------------------------------------------------------------------ // exp2( arg ) = exp( arg * log(2) ) //------------------------------------------------------------------------------ // exp2( arg ) = exp( z ) z := arg * log(2) //------------------------------------------------------------------------------ const uint32_t z = u32mulhi( arg, LOG_2 ); //------------------------------------------------------------------------------ // extract leading two bits (used by table lookup) //------------------------------------------------------------------------------ // exp( z ) = exp( y ) * exp( x ) //------------------------------------------------------------------------------ const uint32_t i = z >> ShiftValue; const uint32_t x = z - (i << ShiftValue); //------------------------------------------------------------------------------ // compute series expansion using modified horner scheme //------------------------------------------------------------------------------ //------------------------------------------------------------------------------ // h := x //------------------------------------------------------------------------------ uint32_t h = x; //------------------------------------------------------------------------------ // h := x + x/7 * x //------------------------------------------------------------------------------ h = u32mulhi( u32mulhi( h, x ), ONE_BY_7 ) + x; //------------------------------------------------------------------------------ // h := x + x/6 * (x + x/7 * x) //------------------------------------------------------------------------------ h = u32mulhi( u32mulhi( h, x ), ONE_BY_6 ) + x; //------------------------------------------------------------------------------ // h := x + x/5 * (x + x/6 * (x + x/7 * x)) //------------------------------------------------------------------------------ h = u32mulhi( u32mulhi( h, x ), ONE_BY_5 ) + x; //------------------------------------------------------------------------------ // h := x + x/4 * (x + x/5 * (x + x/6 * (x + x/7 * x))) //------------------------------------------------------------------------------ h = ( u32mulhi( h, x ) >> 2) + x; //------------------------------------------------------------------------------ // h := x + x/3 * (x + x/4 * (x + x/5 * (x + x/6 * (x + x/7 * x)))) //------------------------------------------------------------------------------ h = u32mulhi( u32mulhi( h, x ), ONE_BY_3 ) + x; //------------------------------------------------------------------------------ // h := x + x/2 * (x + x/3 * (x + x/4 * (x + x/5 * (x + x/6 * (x + x/7 * x))))) //------------------------------------------------------------------------------ h = ( u32mulhi( h, x ) >> 1) + x; //------------------------------------------------------------------------------ // Now: h = exp(x)-1 = x + x²/2! + x³/3! + x⁴/4! + x⁵/5! + x⁶/6! + x⁷/7! //------------------------------------------------------------------------------ // note the "-1" resulting from the missing "1" in the series expansion //------------------------------------------------------------------------------ //------------------------------------------------------------------------------ // Compute: f(a) = exp(z) - 1 = exp(y) * exp(x) - 1 //------------------------------------------------------------------------------ // f(a) = exp(z) - 1 // f(a) = exp(y) * exp(x) - 1 // f(a) = ( 1 + EXP_i_BY_4_MINUS_1 ) * (1 + h ) - 1 //------------------------------------------------------------------------------ // f(a) = EXP_i_BY_4_MINUS_1 * h + EXP_i_BY_4_MINUS_1 + h //------------------------------------------------------------------------------ switch( i ){ //---------------------------------------------------------------------- // i=0: multiply by exp( 0/4 ) == 1 //---------------------------------------------------------------------- // exp( 0/4 ) = 1 + EXP_0_BY_4_MINUS_1 = 1 //---------------------------------------------------------------------- // u32mulhi( EXP_0_BY_4_MINUS_1, h ) + EXP_0_BY_4_MINUS_1 + h // u32mulhi( 0 , h ) + 0 + h // h //---------------------------------------------------------------------- case 0: break; //---------------------------------------------------------------------- // i=1: multiply by exp( 1/4 ) //---------------------------------------------------------------------- // exp( 1/4 ) = 1 + EXP_1_BY_4_MINUS_1 //---------------------------------------------------------------------- case 1: h = u32mulhi( EXP_1_BY_4_MINUS_1, h ) + EXP_1_BY_4_MINUS_1 + h; break; //---------------------------------------------------------------------- // i=2: multiply by exp( 2/4 ) //---------------------------------------------------------------------- // exp( 2/4 ) = 1 + EXP_2_BY_4_MINUS_1 //---------------------------------------------------------------------- case 2: h = u32mulhi( EXP_2_BY_4_MINUS_1, h ) + EXP_2_BY_4_MINUS_1 + h; break; //---------------------------------------------------------------------- // i>=3: cannot happen //---------------------------------------------------------------------- default: break; } //------------------------------------------------------------------------------ // heuristic correction for truncation errors and series remainder //------------------------------------------------------------------------------ return h + ( u32mulhi( arg, 12 ) >> 1 ); }
Download: u32p32exp2m1.cc
Bei der Verwendung des Ganzzahl-Typs uint32_t
wird das Ergebnis der
Multiplikation nach unten gerundet. Dadurch sind Zwischenergebnisse und das
Endergebnis immer zu niedrig. Der Korrekturfaktor verschiebt das Endergebnis,
der Fehler ist dadurch symmetrisch um 0 verteilt und das Endergebnis hat einen
maximalen Fehler von ±4, was ±9.3·10-10 entspricht.
Durch kleine Justierungen an den Parametern dürfte die Genauigkeit noch verbessert werden
können, das war aber für meine Anwendung nicht nötig.
Sie können das Gesamtpaket als u32p32exp2m1.zip herunterladen.