binary logarithm

calculating binary logarithm of 32-bit unsigned word;

log(2,x) = lb(x)

for more precision will calc 8*lb(x)

basic idea:

> #include <stdio.h>
> #include <math.h>
> #include <unistd.h>
>
> double ln2f(double x)
> {
>        const   double eps = 1.0/256.0;
>        int     ival;
>        double  part, fval;
>
>        ival=0;
>        while(x < 1.0) { ival--; x *= 2.0; }
>        while(x >= 2.0){ ival++; x /= 2.0; }
>
>        fval = 0;
>        part = 1.0/2.0;
>        x *= x;
>        while (part > eps) {
>                if (x >= 2.0) { fval += part; x /= 2.0; }
>                part /= 2.0; x *= x;
>        }
>        return ival + fval;
> }

3-bit fraction version:

> uint8_t ln2fx8(double x)
> {
>        uint8_t ival;
>
>        ival=0;
>        while(x < 1.0) { ival-=8; x *= 2.0; }
>        while(x >= 2.0){ ival+=8; x /= 2.0; }
>        // here x E [1;2); i.e. x^2 => [1;4)
>        x *= x; if (x >= 2.0) { ival += 4; x /= 2.0; }
>        x *= x; if (x >= 2.0) { ival += 2; x /= 2.0; }
>        x *= x; if (x >= 2.0) { ival += 1; x /= 2.0; }
>
>        return ival;
> }

now, basic idea two is 8*lb(x) == lb(x^8); i.e. number of significant bits in x^8. x^8 is represented as 256-bit number. o-ho-ho. but in this can be used approimated multiplication, to get number of bits*

now, multiplication improved basic idea. for 3-bit fraction, use estimation of F2.6 multiplication.

> uint8_t ln2x8(uint32_t x)
> {
>        uint8_t ival;
>        if (!x) return 0;
>
>        ival=30*8; // F2.30(1) == 0x40000000; F2.30(2) == 0x80000000; 
>        while(x < 0x40000000L) { ival-=8; x <<= 1; }
>        while(x >=0x80000000L) { ival+=8; x >>= 1; }
>        // here x E [1;2); i.e. x^2 => [1;4)
>
>        //now, switch to F2.6; F2.6(2) == 0x80
>        x >>=24;
>
>        // squaring of F2.6 move to F4.12; but x^2 is in [1;4), cut'em
>        x *= x; x>>=8-2; if (x >= 0x80) { ival += 4; x >>=1; }
>        x *= x; x>>=8-2; if (x >= 0x80) { ival += 2; x >>=1; }
>        x *= x; x>>=8-2; if (x >= 0x80) { ival += 1; x >>=1; }
>
>        return ival;
> }

new, improved version, more precision(:

> uint8_t ln2x8(uint32_t x)
> {
>        uint8_t ival;
>        if (!x) return 0;
>
>        // F1.31(1) == 0x80000000; F1.31(2) == 0x100000000 overflow 
>        ival=31*8;
>        while(x < 0x80000000L) { ival-=8; x <<= 1; }
>
>        // here x E [1;2); i.e. x^2 => [1;4)
>        //now, switch to F1.7; F1.7(2) == 0x100 byte overflow
>        x >>=24;
>
>        // squaring of F1.7 move to F2.14; but x^2 is in [1;4), cut'em
>        x *= x; x>>=8-1; if (x >= 0x100) { ival += 4; x >>=1; }
>        x *= x; x>>=8-1; if (x >= 0x100) { ival += 2; x >>=1; }
>        x *= x; x>>=8-1; if (x >= 0x100) { ival += 1; x >>=1; }
>
>        return ival;
> }

newest, re-improved version, faster. no 32 shifts:

> uint8_t ln2x8(uint32_t x)
> {
>        uint8_t ival;
>        if (!x) return 0;
>
>        ival=0;
>        if (x >= 0x00010000L) ival+=16*8; else x<<=16;
>        if (x >= 0x01000000L) ival+= 8*8; else x<<=8;
>        if (x >= 0x10000000L) ival+= 4*8; else x<<=4;
>        if (x >= 0x40000000L) ival+= 2*8; else x<<=2;
>        if (x >= 0x80000000L) ival+= 1*8; else x<<=1;
>
>        // F1.31(1) == 0x80000000; F1.31(2) == 0x100000000 overflow 
>        // here x E [1;2); i.e. x^2 => [1;4)
>        //now, switch to F1.7; F1.7(2) == 0x100 byte overflow
>        x >>=24;
>
>        // squaring of F1.7 move to F2.14; but x^2 is in [1;4), cut'em
>        x *= x; x>>=8-1; if (x >= 0x100) { ival += 4; x >>=1; }
>        x *= x; x>>=8-1; if (x >= 0x100) { ival += 2; x >>=1; }
>        x *= x; x>>=8-1; if (x >= 0x100) { ival += 1; x >>=1; }
>
>        return ival;
> }