Wie pow implementieren?

Hallo, ich wollte mal fragen, wie normalerweise die power-function für reelle Zahlen implementiert wird? Ich habe es mal über Taylorreihen implementiert, aber das konvergiert super langsam.
Wisst ihr, welchen Algorithmus dazu normalerweise verwendet?

hustbaer

Vermutlich einfach mit exp() und log()?

Ethon

double pow(double a, double b) {
    int tmp = (*(1 + (int *)&a));
    int tmp2 = (int)(b * (tmp - 1072632447) + 1072632447);
    double p = 0.0;
    *(1 + (int * )&p) = tmp2;
    return p;
}

Ich verstehe den Algo nicht, er ist aber schneller als alles andere. Laut Author 8x schneller als pow() aus cmath.
http://martin.ankerl.com/2007/10/04/optimized-pow-approximation-for-java-and-c-c/

TGGC

Mit diesem FPU-Befehl:
http://www.maartensz.org/computing/Assembly/RosaHelps/Mnemonics_Reference/FYL2X.htm

Oder meinst du wie dieser Befehl in der Hardware realisiert wird? Oder die Mathematik dahinter?

krümelkacker

Ethon schrieb:

double pow(double a, double b) {
    int tmp = (*(1 + (int *)&a));
    int tmp2 = (int)(b * (tmp - 1072632447) + 1072632447);
    double p = 0.0;
    *(1 + (int * )&p) = tmp2;
    return p;
}
Ich verstehe den Algo nicht, er ist aber schneller als alles andere. Laut Author 8x schneller als pow() aus cmath.

Blöd nur, dass es
- undefiniertes Verhalten hervorruft (nach §3.10/15).
- auf diverse Implementierungsdetails setzt (Intel-Byte-Ordnung und IEEE-754 floats)
- den Wert nur approximiert

Die Grundidee des Algorithmus basiert darauf, dass die oberen Bits der IEEE-754 Fießkommazahlen-Darstellung einer positiven, normalisierten Zahl -- wenn man das Bitmuster als Ganzzahl interpretiert -- nach Abzug eines Offsets (hier 1072632447) annähernd proportional zum Logarithmus der Fließkommazahl ist. Das, "1+", was hier zweimal auftaucht, ist die Zeigerarithmetik, um auf die "hinteren" 32Bits der double-Zahl zuzugreifen, was mit der Intel-Byte-Ordnung die 32 höherwertigen Bits der Fielßkommazahl sind. Hier wird also quasi a logarithmiert, mit b multipliziert und das wird dann wieder delogarithmiert.

TGGC schrieb:

Oder meinst du wie dieser Befehl in der Hardware realisiert wird? Oder die Mathematik dahinter?

Ich habe mich vielleicht ein wenig kurz gefasst Ich möchte eine eigene Fließkommazahlenklasse implementieren und dafür auch eine pow Funktion bereitstellen.

hustbaer schrieb:

Vermutlich einfach mit exp() und log()?

Ja, das sieht sehr gut aus! Hab mich gleich mal über die Implementierung von ln informiert. Was ich dabei nicht verstehe:
Mein Ansatz sähe jetzt so aus:

Verwende eine Reihe um Entwicklungspunkt 1. Normalisiere Zahl in der Form:
x umschreiben in 2^m * a
Daher: m*ln(2) + ln(a)

Das Beispiel auf Wikipedia verwendet diese Form:

2*m*ln(sqrt(2)) + ln(a)

Ich verschtehe schon, wie die Umformung zu stande kommt. Aber rechtfertigt die schnellere Konvergenz von 1.41... (im Gegensatz zu 2) wirklich die komplette Berechnung der Wurzel von 2?
Oder gibt es irgend ein "Standard-Verfahren", wie man ln ausrechnet?
Ich hoffe, ihr versteht, was ich meine

Michael E.

PowerRower schrieb:

Das Beispiel auf Wikipedia verwendet diese Form:
2*m*ln(sqrt(2)) + ln(a)
Ich verschtehe schon, wie die Umformung zu stande kommt. Aber rechtfertigt die schnellere Konvergenz von 1.41... (im Gegensatz zu 2) wirklich die komplette Berechnung der Wurzel von 2?

Den Wert darf man ruhig hartkodiert in den Quelltext schreiben

TGGC

PowerRower schrieb:

TGGC schrieb:

Oder meinst du wie dieser Befehl in der Hardware realisiert wird? Oder die Mathematik dahinter?

Ich habe mich vielleicht ein wenig kurz gefasst Ich möchte eine eigene Fließkommazahlenklasse implementieren und dafür auch eine pow Funktion bereitstellen.

Das erklaert trotzdem deine Frage nicht. Deine pow Funktion koennte einfach die pow Funktion aus der clib aufrufen. Du koenntest den Code, der in der clib steht, reinkopieren. Du koenntest eigenen Code schreiben, der die nuetzlichen FPU Befehle dafuer aufruft. Du koenntest eigenen Code schreiben, der das auf der CPU nachbildet mit einer Menge Bitgefummel.

Keiner weiss, was du brauchst.

zunächst kann man unterscheiden ob der Exponent ganzzahlig sein soll oder nicht.

Falls ja: Rückführung auf Multiplikation:

b = Basis
e = Exponent (in binärer Darstellung: e_0*2^0 + e_1*2^1 + ... + e_n*2^n)

b^e = b^(e_0*2^0 + e_1*2^1 + ... + e_n*2^n) 
    = b^e_0 * (b^2)^e_1 * (b^4)^e_2 * ... * (b^(2^n))^e_n

Anfang:

B := b  (Potenzen von b)
p := 1  (Zwischenprodukt)

der Exponent e wird nun bitweise nach rechts verschoben. Jedesmal, wenn rechts ein 1-Bit herausfällt, wird

p <- p*B und B <- B*B

gerechnet.

ist der Exponent nicht ganzzahlig, wird es komplizierter.

b^f = exp(f*log b)

führt b^f auf exp und log zurück, und für exp und log gibt es numerisch Standardverfahren.

krümelkacker

Nur um noch ein Stichwort für den Fall von ganzzahligen Exponenten zu nennen: Square-And-Multiply heißt der Algorithmus, auch bekannt unter binärer Exponientiation.

ist sogar ein Bug drin und keiner hat's gemerkt

B <- B*B

muß jedes Mal ausgeführt werden, wenn der Exponent um eine Bitstelle nach rechts verschoben wird, egal ob das herausfallende Bit 1 ist oder 0.