Ich suche einen AVR-tauglichen (Festkomma-)Algorithmus, der den Logarithmus einer 32-Bit Zahl bildet. Zu welcher Basis ist egal, das ist ja nur ein Faktor im Ergebnis. Vermutlich ist die Basis 2 am geeignetsten für Bitschiebereien. Meine Eingangswerte liegen im Bereich von ca. 1000 bis 100.000.000. Es muß nicht so fürchterlich genau sein, aber nur das höchste gesetzte Bit auszählen reicht mir nicht. Hat da jemand vielleicht was Nettes parat? Vielen Dank!
Der Logarithmus ist relativ einfach. Es beginnt damit zu zaehlen auf welcher Position das MSB ist. Fuer die folgenden bits nimmt man ein Tabelle. Eine Potenzreihe ist unpassend.
> Es muß nicht so fürchterlich genau sein, aber nur das höchste gesetzte > Bit auszählen reicht mir nicht. Ein genaueres Resultat kann man aber mit Festkomma-Zahlen gar nicht speichern. Das höchste Bit zu suchen ist also der einfachste, schnellste und genaueste Algorithmus der möglich ist mit Fix-Point-Arithmetik. Beweisskizze: x sei ein unsigned 32-bit Integer (ganz normal, kein Komma). Der Wertebereich liegt zwischen 0 und 2^32-1. log2(x) wird sich entsprechend zwischen -unendlich und 32 bewegen. Soll der log2 wiederum in dem selben Datentyp abgelegt werden, muss man auf ganze Zahlen runden (z.B. hier z.B. abrunden). Somit sind nur 32 Werte als Resultat möglich. Wenn du die wieder zurück-exponenzierst wirst du feststellen das dies genau jene 32 Fälle sind, die du per MSB Bit-Zählen bekommst. Nun zur Fix-Point-Arithmetik. Sei b=Komma-Position dann kann man eine 32-bit-Fix-Punkt Zahl y Schreiben als y = x / 2^b. Wobei x ein 32-bit-Int ist. Es folgt: log2(y) = log2(x / 2^b) = log2(x) - log2(2^b) = log2(x) - b. Sprich der fix-point-log unterscheidet sich nur duch eine ganzzahlige Subtraktion vom integer-log, eine höhere Genauigkeit ist also nicht möglich. Wenn dir die Genauigkeit nicht ausreicht, müssen input und output Datentyp des Logarithmus eine unterschiedliche Komma-Position haben.
Hmm Moment, irgenwas stimmt da doch ganz und gar nicht. Ich glaub was ich da geschrieben hab ist falsch, zumindest der fix-point Teil. Ich schau mir das Morgen nochmals an (brauch wohl etwas Schlaf) :-)
@Stefan: Endlich mal einer, der seine Aussagen ordentlich begründet. Aber ich glaube, in dem Beweis steckt ein kleiner Fehler: > log2(y) = log2(x / 2^b) = log2(x) - log2(2^b) = log2(x) - b = (log2(x) - b) * 2^b / 2^b (log2(x) - b) * 2^b (gerundet) ist also der Wert, der als Ergebnis in die 32-Bit-Variable geschrieben wird und repräsentiert eine Zahl mit b Nachkommastellen. Anders als in deiner ersten Beweishälfte ist log2(x) jetzt natürlich keine Ganzzahl mehr. Folglich ist es kein Problem, aus einer Festkommazahl mit b Nachkomma- stellen den Logarithmus mit ebenfalls b Nachkommastellen zu berechnen. Oder denke ich da falsch? > Wenn dir die Genauigkeit nicht ausreicht, müssen input und output > Datentyp des Logarithmus eine unterschiedliche Komma-Position haben. Was sicher auch kein Problem wäre. Mein Favorit wäre übrigens der erste Vorschlag von sechseinssechs (Schieben und Tabelle).
Ich habe mich an den Vorschlag Schieben+Tabelle gehalten. Per Tabelle nähere ich aus 16 Geradenstücken an. Meine Annäherung liegt wegen der Krümmung unter der Originalkurve, für kleineren durchschnittlichen Fehler könnte man die etwas höher schieben, ich weiß aber nicht wie man diesen Ausgleich berechnen täte. Die Genauigkeit ist schon bemerkenswert hoch, nahe 10e-5. Hier ist mein erster Code:
1 | uint32_t log2fix(uint32_t n) |
2 | {
|
3 | static const uint16_t lut[17] = { // log2(1.0)...log2(2.0) in 1/16 interval |
4 | 0, 5732, 11136, 16248, 21098, 25711, 30109, 34312, 38336, 42196, 45904, 49472, 52911, 56229, 59434, 62534, 65535 |
5 | };
|
6 | uint8_t ri = 31; // result integer part |
7 | uint16_t rf; // result fractional part |
8 | |
9 | uint8_t j; // table index |
10 | uint16_t nf; // input value, fractional part |
11 | |
12 | while (!(n & 0x80000000)) // count bits to MSB |
13 | {
|
14 | n <<= 1; |
15 | ri--; |
16 | }
|
17 | |
18 | n <<= 1; // shift out the mantissa '1', now only fractional part left |
19 | |
20 | j = (uint8_t)(n >> (32-4)); // use 4 bit for table entry selection |
21 | nf = (uint16_t)(n >> (16-4)); // use the bits below for linear interpolation |
22 | |
23 | rf = lut[j] + ((nf * (lut[j+1] - lut[j])) >> 16); |
24 | |
25 | return ((uint32_t)ri << 16) + rf; // return as 8.16 fixpoint |
26 | }
|
Die Hi/Lo Teile kann man mit Unions wohl noch optimieren.
Merci yalu Genau da lag mein Fehler. Offensichtlich kann man mit Fix-Point Arithmetik den Logarithmus genauer darstellen. Das ist mir beim erneuten durchlesen dann auch aufgefallen, nur um mit dem Finger darauf zu zeigen war ich dann doch zu müde. Sorry wenn ich für Verwirrung gesorgt habe :) Schliesse mich 616 an. Die Zahl ins 1..2 Intervall schieben und dann eine der üblichen Approximationsmethoden verwenden. Der Logartihmus ist in dem Bereich ja ziemlich "zahm". Für eine erste grobe Approximation könnte man einfach den Wert übernehmen (Approx durch Gerade mit Steigung 1).
Mit Unions zum partiellen Zugriff auf die Teile eines 32-bit Wertes kriege ich es effizenter, im Hinblick auf 8 Bit Controller ohne Barrel-Shifter. Außer der MSB-Suche sind keine Shifts mehr nötig, wenn ich damit passend für eine Bytegrenze aufhöre. Mit diesem Code kann ich die oberen 3 Bit des Eingangswertes nicht nutzen, aber so große Werte treten bei mir nicht auf. (Bei Verdopplung der Tabellengröße würde man ein Bit davon zurückgewinnen)
1 | #pragma pack(1) // for MSVC
|
2 | // gcc needs __attribute__((packed))
|
3 | |
4 | // this struct allows efficient access to the parts of a 32 bit value, without shifting/masking
|
5 | typedef union |
6 | { // union for little endian CPU, make sure this gets packed |
7 | struct
|
8 | {
|
9 | uint8_t ll8; |
10 | uint8_t lh8; |
11 | uint8_t hl8; |
12 | uint8_t hh8; |
13 | };
|
14 | struct
|
15 | {
|
16 | uint16_t l16; |
17 | uint16_t h16; |
18 | };
|
19 | struct
|
20 | {
|
21 | uint8_t l8; |
22 | uint16_t m16; // this needs unaligned read capability |
23 | uint8_t h8; |
24 | };
|
25 | uint32_t val32; |
26 | } variant_t; |
27 | |
28 | |
29 | // logarithm to base 2, returns a 16.16 fixed point value
|
30 | // restrictions: input value must not be 0 and not exceed 0x1FFFFFFF (because of an implementation shortcut)
|
31 | uint32_t log2fix(uint32_t n) |
32 | {
|
33 | static const uint16_t lut[17] = { // log2(1.0)...log2(2.0) in 1/16 interval |
34 | 0, 5732, 11136, 16248, 21098, 25711, 30109, 34312, 38336, 42196, 45904, 49472, 52911, 56229, 59434, 62534, 65535 |
35 | };
|
36 | variant_t v; // alias for the input value |
37 | variant_t r; // result |
38 | uint8_t j; // table index |
39 | |
40 | r.hh8 = 0; // upper integer part is always zero |
41 | r.hl8 = 28; // integer part is a bit countdown |
42 | |
43 | v.val32 = n; |
44 | |
45 | while (!(v.hh8 & 0xF0)) // count bits until MSB reaches place 28 |
46 | {
|
47 | v.val32 <<= 1; |
48 | r.hl8--; |
49 | }
|
50 | |
51 | j = v.h8 & 0x0F; // use 4 bit in places 27...24 for table entry selection |
52 | |
53 | // use the bits below (in v.m16) for linear interpolation
|
54 | r.l16 = lut[j] + ((v.m16 * (lut[j+1] - lut[j])) >> 16); |
55 | |
56 | return r.val32; // return as 16.16 fixpoint |
57 | }
|
Thread beobachten
Seitenaufteilung abschaltenBitte melde dich an um einen Beitrag zu schreiben. Anmeldung ist kostenlos und dauert nur eine Minute.
Bestehender Account
Schon ein Account bei Google/GoogleMail? Keine Anmeldung erforderlich!
Mit Google-Account einloggen
Mit Google-Account einloggen
Noch kein Account? Hier anmelden.