working with logarithms

notation : we consider x^y to be the y-th power of x, x^y = x*x*x*x*... (y times)

Logarithms measure orders of magnitudes. A useful property is that on logarithmic scales, multiplications become additions and power-of becomes a multiplication. It is the inverse of the power-of function in general. And power of e in particular. What is it that makes the power of e so special ? It is the differential calculus that distinguishes e as base from all others.

diff(a^b, b) = b * log(a)
diff(e^b, b) = b        // with log (e) = 1

We thus use exp() and log() for the base e only. All others are marked. Nevertheless, there are other bases which can be consideres useful, mainly 10 and 2. The power-of are related as

x^y = exp(y*log(x)), thus 10^a = exp(a*log(10)) and 2^a = exp(a*log(2))

likewise the logarithms are related as

log(x) = log(2)*log(base2)(x) = log(10)*log(base10)(x)

working with the power of 2 and log2

Apparently working with power of 2 and log2 are suboptimal compared with the e based stuff. But with computers, they hold more advantages. Power of 2 means, how many binary positions it the number shifted in either direction. Consider a generic binary number, msb to the left, lsb to the right, each position to the right is half to the one on the left :

 ...00000001001101001,0100100000000...   forget leading and trailing zeros
 ...0000000xxxxxxxxxxxxxxxx00000000...   now the x holding the number which is 
                     
  1,00110100101001 shifted to the left by 9 positions, it is bigger than 
  1,00000000000000 and smaller than
 10,00000000000000 

 thus we know the logarithm(2) being between 9 and 10, short 9,xxxxx . The fracional ",xxxxx" 
 is the log of the useful digits.

calculating the log2 of a number

This approach was posted by Peter Dickerson, email unknown, on the comp.arch.embedded newsgroup on 11.may.07.

First we shift the input left until a one is shifted out while subtracting 
the fixed point representation of ln(2) from the answer - initialised to 
32*ln(2). The typical normization step. Obviously this stage can be speeded 
up with binary search approaches.

Start with a 32bit floating point register initialized to 32*log(2) = 32 * 0.69314718.. = 22.180710 . Now in a loop, shift the input to the left and subtract log(2) from the register until the input overflows to the left with a one. If started with a 1.00, the register holds now a 0.00 . For each power of 2 greater 1, once log2 is in the register.

So now we have the value 1+x+y in the form of a 32-bit fixed point number 
x+y with the binary point at the left. x, here, is the top eight bits and y 
the next 24. So 0 <= x < 1 is in steps of  1/256, and 0 <= y < 1/256.

input = 1.xxxxxxxx yyyyyyyy yyyyyyyy yyyyyyyy

We want ln(1+x+y)
  = ln((1+x)*(1+y/(1+x))
  = ln(1+x) + ln(1+y/(1+x)). (a)

The first term of (a) has only 256 values, and so can be looked up from a 
table. The second term is small, which we'll now estimate.

Define z = y/(1+x), then

ln(1+z) is approximately z - z*z/2 +...

note that the z cubed and higher terms are smaller than 1 lsb and to a good 
approximation we can correct ln(1+z) from the linear with a small table 
because no more than eight bits are significant.

OK, to make this a bit more explicit, lets use upper case letters for the 
binary (whole) numbers representing the various fractions. X = x*(2**8); Y = 
y*(2**32). Then the original normalized number is 1 + X/(2**8) + Y/(2**32).

Look up ln(1+X/(2**8)) in a table indexed by X using a 24-bit fraction 
format - call this L.

Set Z = Y/((2**8)+X) = z * (2**24), so its already in the 24-bit fraction 
format and less than 2**16.

Add Z to L.

Subtract a correction corresponding to the z*z/2 term by using (Z>>8) as an 
index into a byte-wide table. In the 24-fraction format the correction is 
roughly Z*Z/2/(2**24) or (Z/(2**8))*(Z/(2**8))/(2**9). So max correction is 
< 128.

So the longest step is the divide, which is actually 24-bit/9-bit. This 
could be speeded up on ARM7 by using the 32x32->64 multiply, looking up the 
multipliers corresponding to 1/(1+x).

That sounds doable, doesn't it.
He had x**y as the y-th power of x, by the way.

To be continued ...

Questions ?
Suggestions?
Feedback ?

sponsored links

embedded
home

last updated: 12.may.07

Copyright (99,2007) Ing.Büro R.Tschaggelar