Algorithmus für Logarithmus-Bererechnung?

Michael S. schrieb:
> Hallo,
>
> ich möchte den Logarithmus einer Zahl mit dem Prinzip berechnen wie er
> in Wikipedia beschrieben ist (siehe Anhang). Leider komme ich da nicht
> weiter.
> Konkretes Beispiel:
> Log2(x) mit x=35

Die Binärdarstellung von 35 lautet (8Bit) 00100011

Die linkste 1 ist an der Bitposition 6, also lautet der Vorkommaanteil 
des log2 von 35 schon mal 5, weil ja 2 hoch 5 gleich 32 ist. 32 ist ein 
bischen weniger als 35, also gibt es noch Nachkommastellen. Welche sind 
das?

Und genau da kommt jetzt das Verfahren ins Spiel.

praktisch gesehen kannst du dir das Abzählen der Bits sparen, wenn du 
einfach sukzessive durch 2 teilst, bis das Ergebnis im Bereich 1 bis 2 
liegt und mitzählst, wie oft du dividieren musstest.

1

35 / 2 = 17.5          1 mal

2

17.5 / 2 = 8.75        2 mal

3

8.75 / 2 = 4.375       3 mal

4

4.375 / 2 = 2.1875     4 mal

5

2.1875 / 2 = 1.09375   5 mal

Das Ergebnis, während der Normierung erhalten, ist dasselbe wie vorher: 
Der ganzzahlige Anteil des Logarithmus von 35 lautet 5.
Um jetzt noch den Nachkommaanteil (binär) zu erhalten, lässt du den 
Algorithmus wie in deinem Link beschrieben, weiterlaufen.

Super Erklärung, vielen Dank Karl Heinz. Ich glaube ich habs jetzt 
verstanden. Ich mache gleich mal noch ein paar Tests.

Gruß Michael

Aber die Quadratur- und Divisionsoperationen für die Nachkommastellen 
muss ich dann doch alle in Gleitkommaformat durchführen, d.h. der 
Rechenaufwand ist auch sehr hoch? Oder muss man diese Operationen mit 
Festkomma-Arithmetik durchführen?

Die Division ist doch nur durch 2, also in Fixed-Point trivial zu 
implementieren.

Ansonsten brauch man nur Vergleiche (easy) und Quadrieren, also 
Multiplikation. Die in Fixed-Point ist auch kein Hexenwerk, z.b. für 2 
32-Bit Zahlen, die den Wert in Q18.14 enthalten:

1

#include <stdint.h>

2

3

uint32_t mul_Q18_14 (uint32_t a, uint32_t b)

4

{

5

    return (a * b) >> 14;

6

}

Weil a und b bereits < 2 sind, sind sie als uint32 < 2^16.  Die 
Multiplikation läuft also nicht über.  Wenn man zusätzlich noch runden 
will:

1

uint32_t mul_Q18_14 (uint32_t a, uint32_t b)

2

{

3

    uint32_t ab = a * b + (1u << 13)

4

    return ab >> 14;

5

}

Und noch einfacher wird's mit Fixed-Point Unterstützung im Compiler, 
z.B. mit avr-gcc 4.8+:

1

#include <stdfix.h>

2

3

unsigned accum mul_Q16_16 (unsigned accum a, unsigned accum b)

4

{

5

    return a * b;

6

}

Wie soll das denn funktionieren z.B. mit a=b=1.09375?
Daraus wird ja beim Datentyp uint_32t oder überhaut uint eine "1".
Oder muss ich gedanklich das Komma schieben? Z.Bsp. rechnen 1093*1093? 
Dann gibts allerdings irgendwann einen Überlauf.

Ich verwende den XC8-Compiler, die stdfix.h wird glaub ich erst beim 
XC16-Compiler unterstützt.

Aha, ich glaube jetzt ist der "Cent" (ich glaube früher hieß der 
Groschen :-) gefallen:
Ich habe mal das hier probiert, das sollte glaub ich gemeint sein mit 
Fixpoint-Arithmetik:
Das habe ich bisher noch nicht verwendet, man lernt doch jeden Tag was 
dazu...

1

#include   <xc.h>

2

#include    <stdint.h>

3

4

int32_t fp(double value, int8_t Q)

5

{ return (int32_t) (value * (1 << Q)); }

6

7

8

int main(int argc, char** argv)

9

{

10

    float zahl1 = 1.09375;

11

    float zahl2 = 1.09375;

12

    int8_t g = 0;

13

    int8_t h = 0;

14

    int16_t y = 0;

15

16

    g = fp(1.094, 6);  // g (Q6) = 70

17

    h = fp(  2.5, 4);  // h (Q4) = 40

18

    y = g * h;         // y (Q10) = g (Q6) * h (Q4)

19

20

    g = fp(zahl1,6);           // g (Q6) = 70 (1.09375*2^6)

21

    h = fp(zahl2,6);           // h (Q6) = 70 (1.09375*2^6)

22

    y = g*h;                   // y (Q12) = g(Q6) * h(Q6) = 4900 (1.19629*2^12)

23

24

    while(1)

25

    {

26

        Nop();

27

    }

28

    return (0);

29

}

Ich habe jetzt folgendes Programm, dass mir den Log2 einer Zahl (20000) 
berechnet und das Ergebnis stimmt bis auf kleine Abweichungen mit dem 
des Taschenrechners überein:
Log2(20000) = Log10(20000)/Log10(2) = 14.28771
Aber die Berechnungszeit bis zum Ergebnis ist alles Andere als 
brauchbar, das geht schneller mit
erg = Log10(20000)*3.322
Aber so in etwa sollte das dem Prinzip entsprechen?

1

#define    _XTAL_FREQ       16000000ULL 

2

#include   <xc.h>

3

#include    <stdint.h>

4

#include    <math.h>

5

6

int32_t fp(double value, int8_t Q)

7

{ 

8

    return (int32_t) (value * (1 << Q));

9

}

10

11

int main(int argc, char** argv)

12

{

13

    uint8_type i,j = 0;

14

    int32_t b = 0;

15

    uint32_type zerg = 0;

16

    uint16_type intval = 0; 

17

    float erg = 0;

18

    float value = 0;

19

20

21

    erg=0;

22

    intval = 20000;

23

    value=intval;

24

    i=0;

25

    while(intval>=2)

26

    {

27

        intval>>=1;

28

        i++;    /*i=14*/

29

    }

30

    zerg+=i;    /*Ganzzahl (Vorkommastellen) in Zwischenergebnis =14*2^15*/

31

    zerg<<=15;  /*Ganzzahl normieren auf 24bit float-Format (VZ+Exp(8)+M(15))*/

32

    value/=(1<<i); /*Normierung 1<value<=2, float, = 20000/2^14*/

33

34

    for(j=0;j<16;j++)  /*Berechnung Nachkommastellen*/

35

    {

36

        value = value*value;

37

        if(value>2)

38

        {

39

            value/=2;

40

            b|=(0x8000>>j);   /*Bit an Position j beginnend von MSB setzen, b=Nachkommastellen binär*/

41

        }

42

        else

43

        {

44

            Nop();

45

        }

46

47

    }

48

    b>>=1;          /*Nachkommastelle schieben fuer Mantisse M0..M14 (15bit)*/

49

    zerg+=b;        /*Nachkommastellen addieren*/

50

    erg=zerg/32768.;    /*Ergebnis Log2(20.000)*/

51

52

    while(1)

53

    {

54

        Nop();

55

    }

56

    return (0);

57

}

Sinn und Zweck des Algorithmus' ist float zu vermeiden.  Ansonsten 
bist du mit der Potenzreihe von logf besser bedient da sie bestimmt 
weniger (float-)Multiplikationen braucht.

Was und wie müsste ich im Code noch ändern, damit ich ganz auf float 
verzichten kann? Auf die letzte Berechnung

1

erg=zerg/32768.;    /*Ergebnis Log2(20.000)*/

kann ich verzichten, die habe ich nur zum Testen eingebaut.

Du überlegst dir welche Genauigkeit du brauchst (sowohl für's Ergebnis 
als auch während der Berechnung) und ersetzt float durch einen adäquaten 
Fixed-Point Typ nebst passender Arithmetik.

Hier ist eine für 8-Bit-Mikrocontroller optimierte Version, die keine
FP-Berechnungen und keine Shift-Operationen mit variabler Weite
verwendet.

Die Funktion fixlog2 berechnet den Zweierlogarithmus für ganzahlige x
von 1 bis 65535. Das Ergebnis ist eine 16.16-Zahl (jeweils 16 Bits von
und nach dem Komma), die mit der Funktion print16_16 ausgegeben wird.
Der maximale absolute Fehler beträgt 0.000051.

Je nachdem, wie gut der Compiler die einzelnen Berechnungen optimiert,
kann man den Code evtl. durch andere Formulierungen noch etwas schneller
machen. Man kann auch durch eine bessere Rundung der Zwischenergebnisse
die Genauigkeit noch etwas verbessern, allerdings braucht das dann auch
wieder etwas mehr Rechenzeit.

Da in der For-Schleife x immer in [1,2) liegt, d.h. das höchstwertige
Bit immer gesetzt istund damit keine Information trägt, könnte man evtl.
die Formeln so umformen, dass stattdessen mit x'=x-1 gerechnet wird.
x' hat dadurch eine zusätzliche Nachkommastelle, was möglicherweise auch
das Ergebnis im 1 Bit genauer macht. Man muss dabei allerdings mögliche
Überläufe von q=x² geschickt verhindern, was den Algorithmus wieder
etwas aufwendiger macht.

Ich bin mir nicht sicher, ob der Code tatsächlich schneller als die
(wesentlich genauere) log2-Funktion aus math.h ist. Am besten probierst
du es aus. Immerhin werden hier 16 32-Bit-Multiplikationen ausgeführt,
was auf einem 8-Bit-µC ohne Hardwaremultiplizierer auch eine ganze Weile
dauert.

Damit du siehst, in welchem Festkommaformat jeweils gerechnet wird, habe
ich dieses als Kommentar neben die einzelnen Zeilen geschrieben.

1

#include <stdio.h>

2

#include <stdint.h>

3

#include <inttypes.h>

4

5

void print16_16(uint32_t x) {          // x: 16.16

6

  printf("%"PRIu32".%06"PRIu32"\n", x >> 16, (x & 0xffff) * 15625 >> 10);

7

}

8

9

uint32_t fixlog2(uint16_t x) {

10

  uint32_t q;                          // 2.30

11

  uint16_t y = 0;                      // 0.16

12

  uint8_t intpart = 15, i;             // 8.0

13

14

  while(!(x & 0x8000)) {

15

    x <<= 1;

16

    intpart--;

17

  }

18

  for(i=0; i<16; i++) {

19

    q = (uint32_t)x * x;               // x: 1.15, q: 2:30

20

    y <<= 1;

21

    if(q & (uint32_t)1 << 31)          // q >= 2.0?

22

      y |= 1;                          // 1-Bit -> y, q *= 2, q: 1.31

23

    else

24

      q <<= 1;                         // 0-Bit -> y, q: 1.31

25

    x = q >> 16;                       // x = q

26

  }

27

  return (uint32_t)intpart << 16 | y;  // return: 16.16

28

}

29

30

int main(void) {

31

  uint16_t x = 20000;                  // 16.0

32

  uint32_t y = fixlog2(x);             // 16.16

33

  print16_16(y);

34

  return 0;

35

}

Ergebnis (wie auch bei deinem Code):

1

14.287689

Vielen Dank Yalu für Deine Hilfe. Das Programm ist sehr kompakt und auch 
schnell. Ich habe gestern abend auch noch mal einen Versuch gestartet 
und bin auf u.a. Lösung gekommen. Wenn ich beide Berechnungen 
vergleiche, braucht meine ca. 500 Cycles länger als Deine mit ca. 7000.
Es sind noch zwei Berechnungen drin, die jeweils ca. 1000 Cycles 
brauchen:

1

    value/=(1<<i); /*Normierung 1<value<=2, float, = 20000/2^14*/

2

    fp_value = fp(value,8);

Vielleicht kann man die noch optimieren.
Wenn man die for-Schleife auf 8 begrenzt, wird es nochmal ein Stück 
schneller, allerdings auch etwas ungenauer.

Hier z.Vgl. mein Code:

1

#include   <xc.h>

2

#include    <stdint.h>

3

#include    <math.h>

4

5

6

int32_t fp(double value, int8_t Q)

7

{

8

    return (int32_t) (value * (1 << Q));

9

}

10

11

int main(int argc, char** argv)

12

{

13

    int8_t i,j = 0;

14

    int16_t b = 0;

15

    int32_t erg = 0;

16

    uint16_t intval = 0;

17

    float value = 0;

18

    uint32_t fp_value = 0;

19

20

    erg=0;

21

    intval = 20000;

22

    value=intval;

23

    i=0;

24

    while(intval>=2)

25

    {

26

        intval>>=1;

27

        i++;    /*i=14*/

28

    }

29

    erg+=i;    /*Ganzzahl (Vorkommastellen) in Zwischenergebnis =14*2^15*/

30

    erg<<=15;  /*Ganzzahl normieren auf 24bit float-Format (VZ+Exp(8)+M(15))*/

31

    value/=(1<<i); /*Normierung 1<value<=2, float, = 20000/2^14*/

32

    fp_value = fp(value,8);

33

34

    for(j=0;j<16;j++)  /*Berechnung Nachkommastellen*/

35

    {

36

        fp_value = fp_value*fp_value;

37

        fp_value>>=8;

38

        if(fp_value>512)/*Wert >2? (*2^8)*/

39

        {

40

            fp_value>>=1;   /*Wert durch 2 teilen*/

41

            b|=(0x8000>>j);   /*Bit an Position j beginnend von MSB setzen, b=Nachkommastellen binär*/

42

        }

43

        else

44

        {

45

            Nop();

46

        }

47

48

    }

49

    b>>=1;          /*Nachkommastelle schieben fuer Mantisse M0..M14 (15bit)*/

50

    erg+=b;        /*Nachkommastellen addieren, Ergebnis = erg/2^15*/

51

52

53

    while(1)

54

    {

55

        Nop();

56

    }

57

    return (0);

58

}

Nochmals vielen Dank.

Hallo Michael S.,

ich habe diesen Thread mit Interesse verfolgt und möchte gerne von Dir 
noch erfahren für was benötigst Du den log zur Basis 2 ?

Besonderen dank an Yalu X. für den Fixpunkt Algorithmus.

Den Logarithmus brauche ich für einen Lichtsensor, mit dem ich die 
Helligkeit eines Displays ändere. Bei den verschiedenen 
Beleuchtungsstärken (Nacht, Zimmerbeleuchtung, Sonne) ist für mich ein 
logarithmisches Verhalten besser geeignet.

Eine Lookup-Table lohnt sich nur für AVRs ohne MUL-Befehl, also ATtinys sowie ATmega16U2. Hierbei wird nur der Bereich [1..2> in die Tabelle gesteckt, also für die Nachkommastellen. Der ganzzahlige Teil wird wie oben mittels clz (count leading zeros) ermittelt.

Statt die Tabelle mit einem externen Skript zu generieren lässt sich diese mit C++14 vom Compiler erstellen:

1

template<unsigned N> struct Logtab{

2

  byte a[N] {};  // lb(x) für x = [1..2>

3

  constexpr Logtab() {

4

    for (byte i=0; i<N; i++)

5

    a[i] = byte(256*log(1+float(i)/N)/log(2)+0.499);

6

  }

7

}

8

9

// Hier: 64 Byte lange Tabelle

10

static constexpr PROGMEM Logtab<64> logtab;

Henrik H. schrieb:

Eine Lookup-Table lohnt sich nur für AVRs ohne MUL-Befehl, also ATtinys sowie ATmega16U2.

Dauert dein Studium schon so lange, dass du erst nach 10 Jahren auf diese Erkenntnis kommst?

Helmut -. (dc3yc)

Dauert dein Studium schon so lange, dass du erst nach 10 Jahren auf diese Erkenntnis kommst? Die Zeit hättest du ja verwenden können, um sein Ergebnis früher zu posten.
Ich finde es äußerst nützlich, wenn funktionierende Codebeispiele ausgeschrieben werden. @Henrik Danke :-)

Michael S. schrieb:

Den Logarithmus brauche ich für einen Lichtsensor, mit dem ich die Helligkeit eines Displays ändere.

Ein schönes Beispiel, warum die konkrete Anwendung immer an den Anfang einer Frage gehört. Dann hätte man sich die ganze Arbeit sparen können und einfach eine Tabelle mit max 16 Stützwerten hingeschrieben.

Dieser Artikel LED-Fading berücksichtigt die Kennlinie des Auges und benutzt letztendlich auch Stützpunkte in Form einer Tabelle.

Andi schrieb:

Peter D. schrieb:

Ein schönes Beispiel, warum die konkrete Anwendung immer an den Anfang einer Frage gehört.

... am Besten auch gleich mitsamt dem Definitionsbereich der benötigten Formel. Kein Mensch kommt auf die Idee, für so weniger Werte einen LOG einzubauen. Das tun nur Programmierer. Physiker nehmen dafür dann sogar noch MATLAB und einen C-Coder.

Ja, das geht sogar noch weiter: In nahezu allen Fällen liefert auch ein Log2-basierender Algorithmus hinreichend gute Ergebnisse und das sehr schnell und ohne speicherfressende Tabellen. Der reale Gamma-Wert wird dadurch zwar nicht korrekt abgebildet, aber es ist eben für die meisten Fälle doch hinreichend nah dran.

Das liegt daran, dass schon der jeweilige "Aktor" (meist ja PWM) sowieso oft viel zu wenig Auflösung hat, um der tatsächlichen Gamma-Kurve korrekt folgen zu können.

Michael S. schrieb:

Den Logarithmus brauche ich für einen Lichtsensor, mit dem ich die Helligkeit eines Displays ändere.

Dann ist dein Problem nicht die Berechnung des Logarithmus, sondern die Auflösung deines Stellgliedes im unteren Helligkeitsbereich.

Welchen Dynamikbereich willst du abdecken und was schafft deine Helligkeitssteuerung?