where : ibrtses embedded
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
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
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) + 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
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 ...
last updated: 12.may.07
Copyright (99,2007) Ing.Büro R.Tschaggelar