Exponentialfunktion zur Basis 2
u32p32exp2m1()
Geschrieben von Wolf am .
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 Multiplikationsbefehle implementiert werden kann. 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.