Here's a collection of different ways to compute logarithms.
ln(x) - Natural logarithms:
x = e^t
ln(x) = t
x = e^ln(x)
log(x) - Common logarithms
if x = 10^t then log(x) = t
lg(x) - Binary Logarithms
if x = 2^t then lg(x) = t
Identities, where f () = ln(), log(), or lg()
f(a*b) = f(a) + f(b)
f(a/b) = f(a) - f(b)
f(a^b) = b*f(a)
We can use the identities to show how ln, log, and lg are related. Suppose we are given x and it is greater than zero. We can always find the three exponents that satisfy this equation:
x = e^t1 = 10^t2 = 2^t3
Using natural logarithms:
ln(x) = t1 = t2*ln(10) = t3*ln(2)
Using common logarithms:
log(x) = t1*log(e) = t2 = t3*log(2)
And finally binary logarithms:
lg(x) = t1*lg(e) = t2*lg(10) = t3
Combining these equations:
ln(x) = log(x) / log(e) = lg(x)/lg(e)
log(x) = ln(x) / ln(10) = lg(x)/lg(10)
lg(x) = ln(x)/ln(2) = log(x)/log(2)
Thus if we know the logarithm of x in any base, we can easily convert to another base by multiplying by a constant.
4. Table Look-up + Linear Interpolation
The easiest way to compute logarithms is to use the built in math libraries. Most contain code that will find the natural logarithm and the common logarithm. Converting to binary logarithms or an other base logarithm requires calling either one of the library functions and multiplication (or division) by a constant. This is too simple. If you want to complicate things then consider the following algorithms.
ln(x) ~ 6*(x-1)/ ( x + 1 + 4*(x^0.5))
The identity f(a*b) = f(a) + f(b) indicates that the logarithm of any number can be expressed as a sum of logarithms of the factors of that number. For example, the natural log of 72 can be expressed as:
ln(72) = ln(8) + ln(9) = 3*ln(2) + 2*ln(3)
In general,
______ _____
/ | | \ \
f(X) = f| | | x_i | = / f(x_i)
\ | | / -----
where x_i are the factors of N.
However this leaves two problems: 1) How can N be factored? and 2) Suppose N is not an integer? Well fortunately both of these stones can be killed by one bird. Finding the prime factors of an arbitrary integer will kill only one stone, not to mention the fact that this it is a very complicated bird. So let's assume that this approach is beyond the scope of the task at hand. Now consider non-integer factors. For example, 72 can be factor infinitely many ways. One that's conducive to finding the natural logarithm:
72 = (2.718281828459)^4 * 1.318725999989
ln(72) = 4 * ln(2.718281828459) + ln(1.318725999989)
= 4 + 0.276666119016
Another for binary logarithms:
72 = (2^6) * 1.125
lg(72) = 6*lg(2) + lg(1.125)
= 6 + 0.1699250014423
In both of these examples, 72 was reduced to two factors:
X = (b^n) * c
where b is the base of the logarithm, n is an integer, and c is a real number that is in the range 1.0 <= c <= b. And the logarithm is
log_b(X) = n + log_b(c)
Thus the problem has been reduced somewhat to just finding n and the logarithm of c. Here are two methods for doing this for binary logarithms:
Here are the steps one can go through to apply the factoring method to binary logarithms. (Actually, with a simple modification you can also apply this method to other based logarithms.)
Input: 16 bit unsigned integer x; 0 < x < 65536
(or 8 bit unsigned integer...)
Output: g, lg(x) the logarithm of x with respect to base 2.
log2arr[i] = lg(2^i / (2^i-1))
i = 1..M, M == desired size of the table.
This array will have to be scaled in such a way that the fixed point arithmetic that uses it below will work. (A convenient scale factor is 2^n, where n is the location of the decimal point in the fixed point arithmetic.)
The first few values of the array are
lg(2/1), lg(4/3), lg(8/7), lg(16/15),...
Note, if you wish to compute the logarithms in a different base, then simply use that base in computing the array. For example, use
log2arr[i] = log_b(2^i / (2^i-1))
where log_b() is the base b logarithm function.
In essence, you need to find the most significant bit of x. Let g, the output, equal lg(MSB(x)). For example, if bit 13 is the most significant bit of x then set g equal to lg(2^13) = 13. (Note, g must be scaled in such a way to make the fixed point arithmetic work.) Then, shift x left until the MSB occupies the bit 15. Here's some psuedo code:
g = 15;
while(g && !(x&0x8000) ) // loop until msb of x occupies bit 15
{
x <<= 1;
g--;
}
An efficiency comment:
Since the MSB of x is always set after the scaling process, we could shift x one more position to the left. This is analagous to how most normalized floating point number are stored in memory.
If you are not computing base 2 logarithms, then the scaling will have to be performed differently. One way is to change g to a floating point number and apply the following psuedo code:
g = log_b(2^15);
while(g && !(x&0x8000) ) // loop until msb of x occupies bit 15
{
x <<= 1;
g -= log_b(2);
}
Another way is to scale as though if binary logarithms were being computed. Then after scaling, convert from binary logarithms to base_b logarithms:
g = g * log_b(2);
Now, change your perspective on y. Consider it a fixed point number with the decimal point between bit positions 15 and 14. Since bit 15 is set, x is equal to one plus a fraction. The fraction is (x &0x7fff) / 0x8000.
The goal here is to find the set of numbers a,b,c,d,... such that:
x = a * b * c * d * ....
Then we can take the logarithm of x and use the identity:
lg(x) = lg(a) + lg(b) + lg(c) + ....
The difficulty is how do we find the factors? Well, suppose we consider the numbers:
2^i
--------- i = 2,3,4,...
2^i - 1
i=1: 2
i=2: 1.333333333333
i=3: 1.142857142857
i=4: 1.066666666667
etc.
If x is greater than one of these factors, then divide that factor out. As an example, suppose x = 1.9
a) 1.9 > 1.333,
x = x/1.333 = 1.425
b) 1.425 > 1.3333 so divide again
x = x/1.333 = 1.06875
c) 1.06875 < 1.1428, so don't divide
d) 1.06875 > 1.06666
x = 1.06875/1.06666 = 1.001953125
etc
So, x ~= 1.3333^2 * 1.0667 etc. and its base logarithm is
lg(1.9) ~= 2*lg(1.3333) + lg(1.0667)
or in terms of the log2arr[] created in step (I),
lg(1.9) ~= 2*log2arr[2] + log2arr[4]
Note that these are not perfect factors. But who cares? The reason these particular factors is that they lend themselves well to an efficient division algorithm.
Recall the multiplication or division by 2 is equivalent to shifting left or shifting right. So that,
x * 2^n = x << n
In the previous step, the division can be converted to shifting:
2^i
x / --------- = x * (2^i - 1) / (2^i)
2^i - 1
= (x * 2^i - x) / 2^i
= x - x / 2^i
= x - (x >> i)
The first 16 entries of log2arr[] created in step(i) are:
1 .4150375 .1926451 9.310941E-02 4.580369E-02 2.272008E-02 1.131531E-02 5.646563E-03 2.820519E-03 1.40957E-03 7.04613E-04 3.522635E-04 1.76121E-04 8.80578E-05 4.402823E-05 2.201395E-05 1.100693E-05
If all of the entries are added together, you get about 1.791906. In the process of iterating towards lg(x), the algorithm will add the entries of log2arr[] together. Now suppose the logarithm of x is greater than 1.791906. It appears the algorithm will under estimate lg(x) in this case. However, there are instances when a factor is used more that once. (Recall the example of factoring 1.9 in step 4).
The preceding steps may be expressed succinctly with the following algorithm. (Note that this comes from Knuth (2), section 1.2.3, exercise 25).
Given x, compute y=log_b(x), where log_b is the base b logarithm.
L1. [initialize.] Set y = 0, z = x >> 1, k = 1.
L2. [Test for end.] If x = 1, stop.
L3. [Compare.] If x-z < 1, go to L5.
L4. [Reduce values.] Set x = x - z, z = x >> k, y = y + log_b(2^k/(2^k-1)), go to L2.
L5. [Shift.] Set z = z >> 1, k = k + 1, go to L2.
This case is very similar to the previous one. The only difference is in how the input is factored. Like before we are given:
Input: 16 bit unsigned integer x; 0 < x < 65536
(or 8 bit unsigned integer...)
Output: g, lg(x) the logarithm of x with respect to base 2.
log2arr[i] = lg(1 + 2^(-(i+1)))
i = 0..M, M == desired size of the table.
The first few values of the array are
lg(3/2), lg(5/4), lg(9/8), lg(17/16),...
Recall that in the previous case the factors were
lg(2/1), lg(4/3), lg(8/7), lg(16/15),...
Again, if you wish to compute logarithms to a different base, then substitute the lg() function with the appropriate based logarithm function.
This is identical to 3(a).ii.
Again, this is identical to 3(a).iii.
This is very similar to step (iv) above. However, we now have different factors. Using the same example, x = 1.9, we can find the factors for this case.
a) 1.9 > 1.5,
x = x/1.5 ==> 1.266666
b) 1.26 > 1.25
x ==> 1.0133333
c) 1.0133 < 1.125 so don't divide
d) 1.0133 < 1.0625 " "
e) etc.
So, x ~= 1.5 * 1.25 * etc.
Like the previouse case, these factors are not perfect. Also, they're somewhat redundant in the sense that 1.5*1.25*1.125*1.0625*... spans a range that is larger than 2 ( ~2.38423 for i<=22). So unlike the previous factoring method, this one will not have repeated factors. Here's some psuedo code:
for(i=1,d=0.5; i<M; i++, d/=2)
if( x > 1+d)
{
x /= (1+d);
g += log2arr[i-1]; // log2arr[i-1] = log2(1+d);
}
Here, d takes on the values of 0.5, 0.25, 0.125, ... , 2^(-i). Then 1+d is the trial factor at each step. If x is greater than this trial factor, then we divide the trial factor out and add to g (ultimately the logarithm of x) the partial logarithm of the factor.
At this point, the answer is in the variable g.
Note that the division by 1+2^(-i) can be rearranged into a multiplication:
Using the expansion
x x / / e \ / e \2 / e \3 \ ------- = --- * | 1 - |---| + |---| - |---| + ... | v + e v \ \ v / \ v / \ v / /
let x = 1, v= 1, e = 2^(-i)
1
---------- = [ 1 - 2^(-i) + 2^(-2*i) - 2^(-3*i) + ...]
1 + 2^(-i)
So, if you want to divide x by 1+2^(-i) :
x
---------- = [ x - x>>i + x>>(2*i) - x>>(3*i) + ...]
1 + 2^(-i)
The number of terms need depends on i. For example, on the first pass, i =1, so you would need to keep 15 terms in the expansion, however for i=4 you would only need to keep the first four terms.
The logarithm function is a monotonic,slowly increasing function. This suggests that a straight line may be a good approximation to it. And in fact with a little work, it's an excellent approximation. The general linear interpolation formula is
b - x x - a
f(x) ~= ------ * f(a) + ------- * f(b), a <= x <= b
b - a b - a
Where in our case, f() is one of the logarithm functions. This approximates f(x) when x is in between a and b.
This can be arranged into the perhaps more familiar point-slope form:
/ f(b) - f(a) \ / a*f(b) - b*f(a) \
f(x) = x * |-------------| - |-----------------|
\ b - a / \ b - a /
Or into a significantly more useful form:
/ f(b) - f(a) \
f(x) = (x - a) * |-------------| + f(a) eq(4.3)
\ b - a /
There are a couple steps required before we can use this equation in a program. First, we need to map x (the range) into a sequence of integers: x -> 0..n. This is necessary to index the table. Second we need to choose the spacing between the tabulated function values so that the division by (b - a) is simplified. Fortunately, we can kill both of these stones with one bird.
For example, suppose that x is a 16-bit integer. We obviously don't want to have 65k look-up table (at least not with this algorithm). So a real simple way to map x into a smaller range is to shift it right. If we have decided to tabulate only 16 values of f(x), then we need only the 4 most significant bits of x (because 2^4 = 16). This is easily obtained by shifting x right 12 bit positions. Then all we need is to tabulate the 16 values: f(0x0000), f(0x1000), f(0x2000), etc. (BTW, the f(0x0000) value is undefined, but there's a way around that problem... see below).
This kills one stone, now for the other. Note that in the point slope equation that we divide by (b - a). However, this is always an integer power of 2. For the example of 16 tabulated values, a and b take on the values of 0, 0x1000, 0x2000, etc. Their difference (assuming we interpolate between consecutive table entries) is always 0x1000. So the division reduces to a shift right!
As an example, suppose we wish to find lg(0x3456) given that we have a table of lg(0x1000*i) for i=0..15. Using equation 4.3:
/ lg(0x4000) - lg(0x3000) \
lg(0x3456) = (0x3456 - 0x3000) * |-------------------------| + lg(0x3000)
\ 0x4000 - 0x3000 /
/ lg(4/3) \
= 0x456 * |---------| + lg(0x3000)
\ 0x1000 /
= (0x456 * lg(4/3) ) >> 12 + lg(0x3000)
= 13.697
This is to within 0.1% of the true answer (~13.710).
There is another trick we can perform. It's to reduce the range of x like we did above in case 3. While not necessary, it does have the benefit of reducing the dynamic range of the tabulated function values. This has two effects. First, the arithmetic can use smaller variables, e.g. 16-bit integers instead of 32-bit integers. Second, if we decide to interpolate with many line segments, then the memory to store the array of points that approximate the function is smaller. If we are working with 16-bit integers then x ranges from 1 to 2^16-1 = 65535. The range reduction method discussed above can scale x between a value of 1.0 and 2.0. (Recall that x is shifted left until its most significant bit occupies the most significant bit location of the memory word.)
So like the algorithms discussed in case 3, let's assume we are given:
Input: 16 bit unsigned integer x; 0 < x < 65536
(or 8 bit unsigned integer...)
Output: g, lg(x) the logarithm of x with respect to base 2.
The following (untested!) C-program implements table look-up plus first order interpolation. It's been designed to use integer arithmetic.
g = 15;
while(g && !(x&0x8000) ) // loop until msb of x occupies bit 15
{
x <<= 1;
g--;
}
// If g==0, then we're done.
if(!g) return(0);
// We have the integer portion of the log(x). Now get the fractional part.
x <<= 1; // "normalize" x, get rid of the MSB.
j = x >> 12; // Get the array index
x_minus_a = x & 0xfff; // The lower bits of x are all that's left
// after subtracting "a".
gf = log_table[j] + (x_minus_a * (log_table[j+1] - log_table[j])) >> 12;
Check out PIC Logarithms for an implementation of this method in PIC assembly.
Here are some handy series expansions.(3)
ln(1+x) = x - x^2/2 + x^3/3 - x^4/4 + ... - (-x)^n/n for -1 < x < 1
/ x-1 \ 1 / x-1 \2 1 / x-1 \3
ln(x)=|-----| + ---|-----| + ---|-----| + ..., for x>= 1/2
\ x / 2 \ x / 3 \ x /
/ x+1 \ / 1 1 1 \
ln|-----| = 2*| --- + ----- + ----- + ... |, for |x| >= 1
\ x-1 / \ x 3*x^3 5*x^5 /
Working software.
Other Theoretical stuff.
(1) Doerfler, Ronald W., "Dead Reconing: Calculating without instruments",Gulf Publishing Company, Houston. ISBN 0-88415-087-9
(2) Knuth, Donald E., "The Art of Computer Programming Vol 1", Addison-Wesley Publishing Company, ISBN 0-201-03822-6
(3) Abramowitz and Stegan, "Handbook of Mathematical Functions", Dover, ISBN 0-486-61272-4