Building a logarithm function in C without using float type

Question

I need to rewrite the log function (base 2 or base 10 doesn't matter which) without using float type, but I need to get the precision of few decimal digits after the decimal point. ( like a float * 100 to get 2 decimals inside integer type eg: if the 1.4352 would be the result, my function should return something like 143 (int type) and I know that the last 2 numbers are the decimals.

I found over the stackoverflow some methods like:

• How can I compute a base 2 logarithm without using the built-in math functions in C#?

but all of them returns int precision ( avoiding the decimals ).

I have no idea how to approach this so the question is:

How to encode (and or change) integer log implementation to support decimal result?

Cecrops

Is that a homework? – Tessatessellate

nope, it's not. i want to use it in a personal project. I done a function witch checks for the highest bit position, but using it i just return the integ part of the resoult, ignoring the zecimals, and i need at least 2 of them ... it will help me to a exponential function. P.S: I dont have to use any libraries – Cecrops

Well log2(1024^100) gives you log2(1024)*100 but that will not be much help (since you cannot effectivaly count that in C)... – Tessatessellate

You need to ask a specific question. What is it specifically that is preventing you from writing your code? Show our attempt, explain what is wrong with it and ask a specific question that will help you progress that code. – Haggadist

i used the ideas from here #12023507 but as i said, it generated me a result at int precision .. – Cecrops

i need any idea to generate a result with some decimal precision – Cecrops

Why not to use float as a result type and when you need just multiply it say by 100 and get integer part ? – Interferon

my bet is you are on MCU without FPU. anyway use fixed precision and binary search for this – Colemancolemanite

i would do that, but i can't get a float result by my method, it only returns the position of the most significant 1 bit from the number, and this is exactly the value of the integer part of the result. I don't know how to find the decimal values... thats why i asked you. – Cecrops

you are right @Spektre. – Cecrops

Do you want C++ example for log2 for example 24.8 bit ? – Colemancolemanite

what do you mean with "24.8 bit" ? Yes, show me senpai. – Cecrops

it means you got 32 bit uint and 24 bits are integer part and 8 decimal ... wil wrote one and post as answer – Colemancolemanite

that sounds great. I hope it will help me. – Cecrops

If you only need 2 decimals the easiest way is to precompute a table. – Darnelldarner

Fast log2(float x) implementation C++, graphics.stanford.edu/~seander/bithacks.html#IntegerLogObvious – Amid

log approximation in embedded systems, Fast Approximation of the Tangent, Hyperbolic Tangent, Exponential and Logarithmic Functions, A simple and fast look-up table method to compute the exp(x) and ln(x) functions, Fast fixed point pow, log, exp and sqrt – Amid

@Cecrops sorry chat does not notify I did not see your message until now. For LUT you need just DWORD LUT[_fx32_fract_bits] = { 0.5^0.5, 0.5^0.25, 0.5^0.125,...} and use that instead of b=fx32_sqrt(b); Do you need real code? – Colemancolemanite

Your MCU surely has a floating point emulation library and probably using that would be the easiest. If not, what is the MCU and the compiler? I could probably port something to it to get it going. How time constrained is that log2, ie. how many clock cycles you got for it? If you can spare a couple microseconds then full blown fp emulation is probably simplest, and you could just use floats. I use floats all the time on quite tiny MCUs - just 1k of RAM, and it’s mostly IEEE-754 compliant and dependable. – Dentelle

Answer 1

You need to use fixed point precision/arithmetics/math for this. It means you use integer type variables but some of the bits are after the decimal point.

for example let assume 8 decimal bits so operations are done like this:

a = number1*256
b = number2*256
c=a+b // +
c=a-b // -
c=(a*b)>>8 // *
c=(a/b)<<8 // /

Here simple fixed point log2 example via binary search in C++:

//---------------------------------------------------------------------------
const DWORD _fx32_bits      =32;                            // all bits count
const DWORD _fx32_fract_bits= 8;                            // fractional bits count
const DWORD _fx32_integ_bits=_fx32_bits-_fx32_fract_bits;   // integer bits count
//---------------------------------------------------------------------------
const DWORD _fx32_one       =1<<_fx32_fract_bits;           // constant=1.0 (fixed point)
const DWORD _fx32_fract_mask=_fx32_one-1;                   // fractional bits mask
const DWORD _fx32_integ_mask=0xFFFFFFFF-_fx32_fract_mask;   // integer bits mask
const DWORD _fx32_MSB_mask=1<<(_fx32_bits-1);               // max unsigned bit mask
//---------------------------------------------------------------------------
DWORD bits(DWORD p)             // count how many bits is p
    {
    DWORD m=0x80000000; DWORD b=32;
    for (;m;m>>=1,b--)
     if (p>=m) break;
    return b;
    }
//---------------------------------------------------------------------------
DWORD fx32_mul(DWORD x,DWORD y)
    {
    // this should be done in asm with 64 bit result !!!
    DWORD a=x,b=y;              // asm has access only to local variables
    asm {                       // compute (a*b)>>_fx32_fract
        mov eax,a               // eax=a
        mov ebx,b               // ebx=b
        mul eax,ebx             // (edx,eax)=eax*ebx
        mov ebx,_fx32_one
        div ebx                 // eax=(edx,eax)>>_fx32_fract
        mov a,eax;
        }
    return a;
    // you can also do this instead but unless done on 64bit variable will overflow
    return (x*y)>>_fx32_fract_bits;
    }
//---------------------------------------------------------------------------
DWORD fx32_sqrt(const DWORD &x) // unsigned fixed point sqrt
    {
    DWORD m,a;
    if (!x) return 0;
    m=bits(x);                  // integer bits
    if (m>_fx32_fract_bits) m-=_fx32_fract_bits; else m=0;
    m>>=1;                      // sqrt integer result is half of x integer bits
    m=_fx32_one<<m;             // MSB of result mask
    for (a=0;m;m>>=1)           // test bits from MSB to 0
        {
        a|=m;                   // bit set
        if (fx32_mul(a,a)>x)    // if result is too big
         a^=m;                  // bit clear
        }
    return a;
    }
//---------------------------------------------------------------------------
DWORD fx32_exp2(DWORD y)       // 2^y
    {
    // handle special cases
    if (!y) return _fx32_one;                    // 2^0 = 1
    if (y==_fx32_one) return 2;                  // 2^1 = 2
    DWORD m,a,b,_y;
    // handle the signs
    _y=y&_fx32_fract_mask;      // _y fractional part of exponent
     y=y&_fx32_integ_mask;      //  y integer part of exponent
    a=_fx32_one;                // ini result
    // powering by squaring x^y
    if (y)
        {
        for (m=_fx32_MSB_mask;(m>_fx32_one)&&(m>y);m>>=1);     // find mask of highest bit of exponent
        for (;m>=_fx32_one;m>>=1)
            {
            a=fx32_mul(a,a);
            if (DWORD(y&m)) a<<=1;  // a*=2
            }
        }
    // powering by rooting x^_y
    if (_y)
        {
        for (b=2<<_fx32_fract_bits,m=_fx32_one>>1;m;m>>=1)      // use only fractional part
            {
            b=fx32_sqrt(b);
            if (DWORD(_y&m)) a=fx32_mul(a,b);
            }
        }
    return a;
    }
//---------------------------------------------------------------------------
DWORD fx32_log2(DWORD x)    // = log2(x)
    {
    DWORD y,m;
    // binary search from highest possible integer power of 2 to avoid overflows (log2(integer bits)-1)
    for (y=0,m=_fx32_one<<(bits(_fx32_integ_bits)-1);m;m>>=1)
        {
        y|=m;   // set bit
        if (fx32_exp2(y)>x) y^=m; // clear bit if result too big
        }
    return y;
    }
//---------------------------------------------------------------------------

Here simple test (using floats just for loading and printing you can handle booth on integers too, or by compiler evaluated constants):

float(fx32_log2(float(125.67*float(_fx32_one)))) / float(_fx32_one)

This evaluates: log2(125.67) = 6.98828125 my win calc returns 6.97349648 which is pretty close. More precise result you need more fractional bits you need to use. Int and compile time evaluation float example:

(100*fx32_log2(125.67*_fx32_one))>>_fx32_fract_bits

returns 698 which means 6.98 as we multiplied by 100. You can also write your own load and print function to convert between fixed point and string directly.

To change precision just play with _fx32_fract_bits constant. Anyway if your C++ does not know DWORD it is just 32bit unsigned int. If you are using different type (like 16 or 64 bit) then just change the constants accordingly.

For more info take a look at:

• Power by squaring for negative exponents

[Edit2] fx32_mul on 32bit arithmetics without asm base 2^16 O(n^2)

DWORD fx32_mul(DWORD x,DWORD y)
    {
    const int _h=1; // this is MSW,LSW order platform dependent So swap 0,1 if your platform is different
    const int _l=0;
    union _u
        {
        DWORD u32;
        WORD u16[2];
        }u;
    DWORD al,ah,bl,bh;
    DWORD c0,c1,c2,c3;
    // separate 2^16 base digits
    u.u32=x; al=u.u16[_l]; ah=u.u16[_h];
    u.u32=y; bl=u.u16[_l]; bh=u.u16[_h];
    // multiplication (al+ah<<1)*(bl+bh<<1) = al*bl + al*bh<<1 + ah*bl<<1 + ah*bh<<2
    c0=(al*bl);
    c1=(al*bh)+(ah*bl);
    c2=(ah*bh);
    c3= 0;
    // propagate 2^16 overflows (backward to avoid overflow)
    c3+=c2>>16; c2&=0x0000FFFF;
    c2+=c1>>16; c1&=0x0000FFFF;
    c1+=c0>>16; c0&=0x0000FFFF;
    // propagate 2^16 overflows (normaly to recover from secondary overflow)
    c2+=c1>>16; c1&=0x0000FFFF;
    c3+=c2>>16; c2&=0x0000FFFF;
    // (c3,c2,c1,c0) >> _fx32_fract_bits
    u.u16[_l]=c0; u.u16[_h]=c1; c0=u.u32;
    u.u16[_l]=c2; u.u16[_h]=c3; c1=u.u32;
    c0 =(c0&_fx32_integ_mask)>>_fx32_fract_bits;
    c0|=(c1&_fx32_fract_mask)<<_fx32_integ_bits;
    return c0;
    }

In case you do not have WORD,DWORD add this to start of code

typedef unsigned __int32 DWORD;
typedef unsigned __int16 WORD;

or this:

typedef uint32_t DWORD;
typedef uint16_t WORD;

[Edit3] fx32_mul debug info

let call and trace/breakpoint this (15 fractional bits):

fx32_mul(0x00123400,0x00230056);

Which is:

0x00123400/32768 * 0x00230056/32768 =
36 * 70.00262451171875 = 2520.094482421875

So:

DWORD fx32_mul(DWORD x,DWORD y) // x=0x00123400 y=0x00230056
    {
    const int _h=1;
    const int _l=0;
    union _u
        {
        DWORD u32;
        WORD u16[2];
        }u;
    DWORD al,ah,bl,bh;
    DWORD c0,c1,c2,c3;
    // separate 2^16 base digits
    u.u32=x; al=u.u16[_l]; ah=u.u16[_h]; // al=0x3400 ah=0x0012
    u.u32=y; bl=u.u16[_l]; bh=u.u16[_h]; // bl=0x0056 bh=0x0023
    // multiplication (al+ah<<1)*(bl+bh<<1) = al*bl + al*bh<<1 + ah*bl<<1 + ah*bh<<2
    c0=(al*bl);        // c0=0x00117800
    c1=(al*bh)+(ah*bl);// c1=0x0007220C
    c2=(ah*bh);        // c2=0x00000276
    c3= 0;             // c3=0x00000000
    // propagate 2^16 overflows (backward to avoid overflow)
    c3+=c2>>16; c2&=0x0000FFFF; // c3=0x00000000 c2=0x00000276
    c2+=c1>>16; c1&=0x0000FFFF; // c2=0x0000027D c1=0x0000220C
    c1+=c0>>16; c0&=0x0000FFFF; // c1=0x0000221D c0=0x00007800
    // propagate 2^16 overflows (normaly to recover from secondary overflow)
    c2+=c1>>16; c1&=0x0000FFFF; // c2=0x0000027D c1=0x0000221D
    c3+=c2>>16; c2&=0x0000FFFF; // c3=0x00000000 c2=0x0000027D
    // (c3,c2,c1,c0) >> _fx32_fract_bits
    u.u16[_l]=c0; u.u16[_h]=c1; c0=u.u32; // c0=0x221D7800
    u.u16[_l]=c2; u.u16[_h]=c3; c1=u.u32; // c1=0x0000027D
    c0 =(c0&_fx32_integ_mask)>>_fx32_fract_bits; // c0=0x0000443A
    c0|=(c1&_fx32_fract_mask)<<_fx32_integ_bits; // c0=0x04FA443A
    return c0; // 0x04FA443A -> 83510330/32768 = 2548.53302001953125
    }

Colemancolemanite

Comments are not for extended discussion; this conversation has been moved to chat. – Liddie

During log2, there’s no need to run a full exponentiation: you already have a bunch of bits previously computed - all you need at each step is a trial fixed point multiplication by a constant from a table. In fact, you don’t need either exponentiation nor square roots. Just multiplies, adds and compares. But at that point you may see that you’ll get many less multiplies if you use the typical ways floating point is emulated, so taking verbatim code from a float software library and cutting down its accuracy would yield even faster code! – Dentelle