Algorithmus für Logarithmus-Bererechnung?

[Michael S. (michael8192)]

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.

[Karl H. (kbuchegg)]

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.

35 / 2 = 17.5          1 mal
17.5 / 2 = 8.75        2 mal
8.75 / 2 = 4.375       3 mal
4.375 / 2 = 2.1875     4 mal
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.

[Michael S. (michael8192)]

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

Gruß Michael

[Michael S. (michael8192)]

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?

[Johann L. (gjlayde)]

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:

#include <stdint.h>

uint32_t mul_Q18_14 (uint32_t a, uint32_t b)
{
    return (a * b) >> 14;
}

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:

uint32_t mul_Q18_14 (uint32_t a, uint32_t b)
{
    uint32_t ab = a * b + (1u << 13)
    return ab >> 14;
}

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

#include <stdfix.h>

unsigned accum mul_Q16_16 (unsigned accum a, unsigned accum b)
{
    return a * b;
}

[Michael S. (michael8192)]

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.

[Bitflüsterer (Gast)]

> 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.

[Michael S. (michael8192)]

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

[Michael S. (michael8192)]

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...

#include   <xc.h>
#include    <stdint.h>

int32_t fp(double value, int8_t Q)
{ return (int32_t) (value * (1 << Q)); }


int main(int argc, char** argv)
{
    float zahl1 = 1.09375;
    float zahl2 = 1.09375;
    int8_t g = 0;
    int8_t h = 0;
    int16_t y = 0;

    g = fp(1.094, 6);  // g (Q6) = 70
    h = fp(  2.5, 4);  // h (Q4) = 40
    y = g * h;         // y (Q10) = g (Q6) * h (Q4)

    g = fp(zahl1,6);           // g (Q6) = 70 (1.09375*2^6)
    h = fp(zahl2,6);           // h (Q6) = 70 (1.09375*2^6)
    y = g*h;                   // y (Q12) = g(Q6) * h(Q6) = 4900 (1.19629*2^12)

    while(1)
    {
        Nop();
    }
    return (0);
}

[Michael S. (michael8192)]

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?

#define    _XTAL_FREQ       16000000ULL 
#include   <xc.h>
#include    <stdint.h>
#include    <math.h>

int32_t fp(double value, int8_t Q)
{ 
    return (int32_t) (value * (1 << Q));
}

int main(int argc, char** argv)
{
    uint8_type i,j = 0;
    int32_t b = 0;
    uint32_type zerg = 0;
    uint16_type intval = 0; 
    float erg = 0;
    float value = 0;
 

    erg=0;
    intval = 20000;
    value=intval;
    i=0;
    while(intval>=2)
    {
        intval>>=1;
        i++;    /*i=14*/
    }
    zerg+=i;    /*Ganzzahl (Vorkommastellen) in Zwischenergebnis =14*2^15*/
    zerg<<=15;  /*Ganzzahl normieren auf 24bit float-Format (VZ+Exp(8)+M(15))*/
    value/=(1<<i); /*Normierung 1<value<=2, float, = 20000/2^14*/
    
    for(j=0;j<16;j++)  /*Berechnung Nachkommastellen*/
    {
        value = value*value;
        if(value>2)
        {
            value/=2;
            b|=(0x8000>>j);   /*Bit an Position j beginnend von MSB setzen, b=Nachkommastellen binär*/
        }
        else
        {
            Nop();
        }

    }
    b>>=1;          /*Nachkommastelle schieben fuer Mantisse M0..M14 (15bit)*/
    zerg+=b;        /*Nachkommastellen addieren*/
    erg=zerg/32768.;    /*Ergebnis Log2(20.000)*/
    
    while(1)
    {
        Nop();
    }
    return (0);
}

[Johann L. (gjlayde)]

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.

[Michael S. (michael8192)]

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

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

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

[Johann L. (gjlayde)]

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.

[Yalu X. (yalu) (Moderator)]

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.

#include <stdio.h>
#include <stdint.h>
#include <inttypes.h>

void print16_16(uint32_t x) {          // x: 16.16
  printf("%"PRIu32".%06"PRIu32"\n", x >> 16, (x & 0xffff) * 15625 >> 10);
}

uint32_t fixlog2(uint16_t x) {
  uint32_t q;                          // 2.30
  uint16_t y = 0;                      // 0.16
  uint8_t intpart = 15, i;             // 8.0

  while(!(x & 0x8000)) {
    x <<= 1;
    intpart--;
  }
  for(i=0; i<16; i++) {
    q = (uint32_t)x * x;               // x: 1.15, q: 2:30
    y <<= 1;
    if(q & (uint32_t)1 << 31)          // q >= 2.0?
      y |= 1;                          // 1-Bit -> y, q *= 2, q: 1.31
    else
      q <<= 1;                         // 0-Bit -> y, q: 1.31
    x = q >> 16;                       // x = q
  }
  return (uint32_t)intpart << 16 | y;  // return: 16.16
}

int main(void) {
  uint16_t x = 20000;                  // 16.0
  uint32_t y = fixlog2(x);             // 16.16
  print16_16(y);
  return 0;
}

Ergebnis (wie auch bei deinem Code):

14.287689

[Michael S. (michael8192)]

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:

    value/=(1<<i); /*Normierung 1<value<=2, float, = 20000/2^14*/
    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:

#include   <xc.h>
#include    <stdint.h>
#include    <math.h>


int32_t fp(double value, int8_t Q)
{
    return (int32_t) (value * (1 << Q));
}

int main(int argc, char** argv)
{
    int8_t i,j = 0;
    int16_t b = 0;
    int32_t erg = 0;
    uint16_t intval = 0;
    float value = 0;
    uint32_t fp_value = 0;
    
    erg=0;
    intval = 20000;
    value=intval;
    i=0;
    while(intval>=2)
    {
        intval>>=1;
        i++;    /*i=14*/
    }
    erg+=i;    /*Ganzzahl (Vorkommastellen) in Zwischenergebnis =14*2^15*/
    erg<<=15;  /*Ganzzahl normieren auf 24bit float-Format (VZ+Exp(8)+M(15))*/
    value/=(1<<i); /*Normierung 1<value<=2, float, = 20000/2^14*/
    fp_value = fp(value,8);

    for(j=0;j<16;j++)  /*Berechnung Nachkommastellen*/
    {
        fp_value = fp_value*fp_value;
        fp_value>>=8;
        if(fp_value>512)/*Wert >2? (*2^8)*/
        {
            fp_value>>=1;   /*Wert durch 2 teilen*/
            b|=(0x8000>>j);   /*Bit an Position j beginnend von MSB setzen, b=Nachkommastellen binär*/
        }
        else
        {
            Nop();
        }

    }
    b>>=1;          /*Nachkommastelle schieben fuer Mantisse M0..M14 (15bit)*/
    erg+=b;        /*Nachkommastellen addieren, Ergebnis = erg/2^15*/


    while(1)
    {
        Nop();
    }
    return (0);
}

Nochmals vielen Dank.

[Uwe (de0508)]

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.

[Michael S. (michael8192)]

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.

[micha (Gast)]

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)?!