Code für Logarithmus gesucht

Jörg H. #704874
Lesenswert?
• ▲
▼ 
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!
Antwort Zitieren Bearbeiten Thread verschieben Thread sperren Anmeldepflicht aktivieren Anpinnen Thread löschen Thread mit anderem zusammenführen Beitrag melden
Stefan Gast #705045
Lesenswert?
• ▲
▼ 
> 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.
Antwort Zitieren Bearbeiten Löschen Beitrag melden
Stefan Gast #705055
Lesenswert?
• ▲
▼ 
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) :-)
Antwort Zitieren Bearbeiten Löschen Beitrag melden
yalu Gast #705080
Lesenswert?
• ▲
▼ 
@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).
Antwort Zitieren Bearbeiten Löschen Beitrag melden
Jörg H. #705215
Lesenswert?
• ▲
▼ 
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.
Antwort Zitieren Bearbeiten Löschen Beitrag melden
Stefan Gast #705282
Lesenswert?
• ▲
▼ 
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).
Antwort Zitieren Bearbeiten Löschen Beitrag melden
Jörg H. #705348
Lesenswert?
• ▲
▼ 
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
}
Antwort Zitieren Bearbeiten Löschen Beitrag melden
Thread beobachten

Antwort schreiben

Bitte melde dich an, um einen Beitrag zu schreiben.

Bestehender Account

oder

Mit Google-Account einloggen

Noch kein Account?

Die Registrierung ist kostenlos und dauert nur eine Minute.

Jetzt registrieren