Forum: Compiler & IDEs Algorithmus für Logarithmus-Bererechnung?


Announcement: there is an English version of this forum on EmbDev.net. Posts you create there will be displayed on Mikrocontroller.net and EmbDev.net.
von Michael S. (michael8192)


Angehängte Dateien:

Lesenswert?

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
Als erstes soll x in das Gleitkommaformat gewandelt und anschliessend 
auf 1<=x<2 normalisiert werden? (Ich verwende den XC8-Compiler mit der 
modified 24-bit IEEE754 Darstellung).
Das bedeutet also 35d entspricht 1,000110000000000*2^5
Exponent = 132 = 2^(5+127) = 0x84, Mantisse = 0x0C00 = 3072d, Vorzeichen 
pos.

35d = 35.0 wird also vom Compiler als:
0 10000100 000110000000000 (VZ+Exp+M) = 0x420C00 dargestellt.

Jetzt müsste ich doch gemäss dem Algorithmus x quadrieren und Abfragen, 
ob die Zahl >2 ist. Wenn >2 dann die Zahl durch zwei teilen.
Irgendwie komme ich da nicht weiter. Vielleicht hat das schon mal jemand 
gemacht und kann mir weiterhelfen?

Besten Dank.

von Karl H. (kbuchegg)


Lesenswert?

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.

: Bearbeitet durch User
von Michael S. (michael8192)


Lesenswert?

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

Gruß Michael

von Michael S. (michael8192)


Lesenswert?

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?

: Bearbeitet durch User
von Johann L. (gjlayde) Benutzerseite


Lesenswert?

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
}

: Bearbeitet durch User
von Michael S. (michael8192)


Lesenswert?

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.

von Bitflüsterer (Gast)


Lesenswert?

> Wie soll das denn funktionieren ...

Schau' Dir doch mal die Dokumentation zu stdfix.h an und ein paar 
Beispiele lassen sich sicher auch dazu finden.

von Michael S. (michael8192)


Lesenswert?

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

von Michael S. (michael8192)


Lesenswert?

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
}

von Michael S. (michael8192)


Lesenswert?

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
}

von Johann L. (gjlayde) Benutzerseite


Lesenswert?

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.

von Michael S. (michael8192)


Lesenswert?

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.

von Johann L. (gjlayde) Benutzerseite


Lesenswert?

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.

von Yalu X. (yalu) (Moderator)


Lesenswert?

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

von Michael S. (michael8192)


Lesenswert?

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.

von Uwe (de0508)


Lesenswert?

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.

: Bearbeitet durch User
von Michael S. (michael8192)


Lesenswert?

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.

von micha (Gast)


Lesenswert?

Michael S. schrieb:
> 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.

Wäre für sowas nicht ein einfacher LUT besser (einfacher, schneller, 
simpler)?!

von Henrik H. (Firma: TU Chemnitz) (heha)


Lesenswert?

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;

: Bearbeitet durch User
von Helmut -. (dc3yc)


Lesenswert?

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?

von Christoph M. (mchris)


Lesenswert?

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 :-)

von Andi (chefdesigner)


Lesenswert?

Henrik H. schrieb:
> [i] = byte(256*log(1+float(i)/N)/log(2)+0.499);

Ziemlich mutwillige Form der Rundung!

von Peter D. (peda)


Lesenswert?

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.

von Frank M. (ukw) (Moderator) Benutzerseite


Lesenswert?

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

Bitte 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
Noch kein Account? Hier anmelden.