Logarithmic Arithmetic
Logarithms can be used to quickly multiply two numbers or to raise a number to a (possibly non-integral) power.
Logarithms are defined with respect to a positive base a. The  logarithm base a of a positive number x is that number which a must be raised to the power of to give x. Writing loga(x) for this number we have x = aloga(x).
The inverse of loga is written antiloga and we have antiloga(y) = ay.
Logarithms are useful because xy = a log(x)a log(y) = a log(x)+ log(y) and so loga(xy) = loga(x) + loga(y).
Further xy = (a log(x))y = aylog(x) and so loga(xy)  = yloga(x).
Combining these two results we also have loga(x/y) =loga(x)­loga(y).
Clearly loga(a)=1 and loga(1)=0 and so antiloga(1)=a and antiloga(0)=1.
Common values for a are two, ten, and the  natural base e » 2.718. We can change between these using: logb(x) = loga(x)/loga(b).
We write ln(x) for loge(x).

How errors in Log and Antilog tables effect final results
We are primarily interested in how the errors in our tabulated log and antilog functions effect the accuracy of operations performed using them. We will assume for the moment that the maximum absolute error of L(x) approximating L(x)=Kloga(x) is d and that the maximum fractional error of A(x) approximating A(x)=ax/K is e. Noting that the first two derivatives of L(x)=Kloga(x) are L'(x) = K/( ln(a)x) ; L"(x) = -K/( ln(a)x2) and that the first two derivatives of A(x)=ax/K are  A'(x) =  ln(a)ax/K/K ;  A"(x) =  ln(a)2ax/k/K2, we can write:
 A(L(x)+L(y)) £ A(L(x)+L(y))(1+e) £ A(L(x)+L(y)+2d>)(1+e) = a2d/K(1+e)A(L(x)+L(y)) = a2d/K(1+e)xy.

Similarly A(L(x)+L(y)) ³ a2d/K(1-e)xy.
Now ag » 1 + g ln(a) for small g so, provided d is small compared to K, the fractional error in A(L(x)+L(y)) as an approximation to xy is approximately ln(a)2d/K+e.
For a=2 this is just under 1.39d/K+e.
If the A LUT has an index range from 0 to 2L(N) and if it is the direct inverse of the log table so that A[L[x]]=x, then the product of two positive integers in the range 1 to N is approximated by A[L[x]+L[y]]. However, a table of 2L(N) cells each wide enough to hold N2 may require more memory than is available.

Examples

To illustrate the above we will consider multiplying two positive nonzero 16-bit numbers using logs to give a 32-bit result. We will assume that a total of 32Kbytes is available for log and antilog tables. Our first method will attempt to optimise speed, our second accuracy.
We will use binary logarithms (a=2) and assume that our options for widths of LUT entries are 16 or 32 bits.
We will write NL,WL,NA,WA for the length and width (in bytes) of the log and antilog tables respectively.

Case 1 - Speed
We take zero base functions and insist that all contraction and expansion functions be simple binary shifts (ie. multiplication by integral powers of 2). This effectively forces us to make NL and NA powers of two.
Let us tentatively assign 16 Kbytes = 214 to the log table and consider WL=2, NL=213.
We will use X = 2-3Z216 as a representative portion of Z216 and a contraction mapping l1(x)=[2-3x] obtained by combining l11(x)=2-3x with l12(x)=[x].
We tabulate L[i]=[Klog2(i)| i Î Z216-{0},  L[0]=0.
(This is better than setting XL*=XL and l12(x)=[2-3x] since then we must tabulate [Klog2(23i)| which is just the same table with 3K added to every entry).
We choose K so that L[NL] is as close as possible to 215 while still below it (so that the sum of two table entries is always less than 216). Since log2(213)=13 we set K=215/13.
Our expansion function is l2(y,x)=y+3K but we will never actually compute this.
Since only 16 bits of data are coming out of the log table there is no point in having NA>216 and we will consider WA=2, NA=213.
We need XA= 2*(dynamic range of L)= 6K+Z216-1. We will take XA*=Z216-1 as a representative portion of this using a11(x)=(x-6K). Combining this with a12(x)=[2-3x] yields a contraction function of a1(x)=[2-3(x-6K)].
To tabulate A(x)=2x/K over XA* = Z216-1 using a12(x)=[2-3x] we would store A[i] = 28i/K. However, this has a dynamic range of 1 to 2(216)/K = 226 so to to avoid exceeding our 2-byte width we tabulate
A[i] = [2-1028i/K|.
Our expansion function is the composition of a21(y) = 210y and a22(y,x) = 26y (compensating for a11(x) = x-6K), ie. a2(y,x) = 216y.
We thus have
xy = 216A[2­3(L[2­3x]+L[2­3y])]
(Where T[x] is taken to mean T[[x]]. T[[x|] will give greater accuracy but is not a simple binary shift.). However the errors are high:
Of all the pairs x,y whose product gives 2m, x=216 y=2m-16 gives the gretest fractional error.
This will be due to (writing p=m-16):

• An absolute sparseness error in L of  Klog2(1+1/2P-3)+Klog2(1+1/216-3)»Klog2(1+23-p+2-13);
• An absolute thinness error in L of 2-1+2-1=1;
• A fractional sparseness error in A of approx. 23ln(2)/K;
• A fractional thinness error in A of 215/xy = 215-m = 2-1-p.
These will contribute fractional errors of
• 2K( log(1+2(3-p)+2-13))/K-1 = 23-p+2-13;
• 21/K-1 = 1.13*2-12        » 2-12;
• 23ln(2)/K = 1.13*2-9         » 2-9;
• and  2-1-p
to xy respectively.
We can prepare the following table of contributions to the net fractional error in xy of the four error sources:
 xy p LogSparse LogThin A‘logSparse A‘logThin TotalFrac TotalAbs <216 Big 2-12 2-9 Big Big 216 0 23 2-12 2-9 2-1 23 219 218 2 21 2-12 2-9 2-3 21 219 220 4 2-1 2-12 2-9 2-5 2-1 219 222 6 2-3 2-12 2-9 2-7 2-3 219 224 8 2-5 2-12 2-9 2-9 2-5 219 226 10 2-7 2-12 2-9 2-11 2-7 219 228 12 2-9 2-12 2-9 2-13 2-8 220 230 14 2-11 2-12 2-9 2-15 2-9 221 232 16 2-12 2-12 2-9 2-17 2-9 223

Reducing NA to 212 and doubling WA (ie. 32 bit entries) eliminates the thinness errors in A at the expense of doubling the sparseness ones but since the sparseness errors in L dominate A‘s thinness errors this does not improove things.
If we allocate another 16 Kbytes to the LUTs we can either:
(i) Double WL, (ii) Double NL, (iii) Double WA, or (iv) Double NA.
(ii) is the most sensible since other than for large xy, the sparseness of L is the major source of error and this would be halved giving: bs]
 xy p LogSparse LogThin A‘logSparse A‘logThin TotalFrac TotalAbs 216 0 22 2-12 2-9 2-1 22 218 218 2 20 2-12 2-9 2-3 20 218 220 4 2-2 2-12 2-9 2-5 2-2 218 222 6 2-4 2-12 2-9 2-7 2-4 218 224 8 2-6 2-12 2-9 2-9 2-6 218 226 10 2-8 2-12 2-9 2-11 3*2-9 3*217 228 12 2-10 2-12 2-9 2-13 3*2-10 3*218 230 14 2-12 2-12 2-9 2-15 2-9 221 232 16 2-13 2-12 2-9 2-17 2-9 223

Case Two - Accuracy
We now sacrifice a little evaluation speed and do not restrict our expansion and contraction functions to binary shifts. We will assume 32 Kbytes and again consider NA=NL=213, WA=WL=2. We will set K=215/13 as before.
We know that the greatest fractional error in xy derives from sparseness errors in L so we will tackle these first. The problem lies in our choice of l11(x)=2-3x. When this is followed by l12(x)=[x] the bottom three bits of x are lost which for low x is disastrous. A better contraction mapping is based on l11(x)=2-m(x)x where m(x) is the minimal positive ineteger such that l11(x)ÎXL*=[1,213).
This is the contraction function discussed in the theoretical section on Log Tables with q=0 and t=13. Application of l12(z)=[z] where z=l11(x) will then only cause data loss if m(x)>0. But if this is so z will be at least NL/2 so our absolute sparseness error in L is at most Klog2(1+2/NL)»2K/(NLln(2)) contributing a fractional error in xy of 4/NL, ie. 2-11 in this case.
Rather than incorporating an addition of Km(x) into l2 we multiply the final result by 2Km(x)/K=2m(x) thus:
xy = 2m(x)+m(y)+16A[2-3(L[2-m(x)x]+L[2-m(y)y])].
where L[i]=[Klog2(i)| iÎZN-{0}, A[i]=[2-1028i/K|
The disadvantage of this contraction function is that the operation "shift x right until x is first less than y and remember the number of shifts performed" is usually slow to impliment with time proportional to the number of shifts required. This problem aside, we could reduce NL to 210 (a 4 Kbyte LUT) before the L errors would again become significant compared to th A ones.
Having reduced the Log Sparseness errors to 2-11 our next priority is to tackle the Antilog Thinness ones. Doubling WA at the expense of NA    eliminates them entirely but doubles the fractional Antilog Thiness errors to 2-8 for all products (including those less than 216). At this stage we have:
 xy p LogSparse LogThin A‘logSparse A‘logThin TotalFrac TotalAbs All 2-11 2-12 2-8 0 2-8 2-8xy

To reduce the sparseness errors in A requires linear interpolation. This can be done relatively painlessly as in the Sparseness Errors section of Antilog Tables above with d=0, b=24, p=4.
We set y=[2-4x], z=2-4x-y, so that x=24(y+z) and express z as Rz/24, ie. Rz=bottom 4 bits of x.
We set G=K(log2(216/K-1)-4)»-11.83K»-&746E and define
 g(x) = A(24y) + A(24y + Klog2(Rz) + G) » A(24y) + A(24y + L[Rz] + G) » 216A[y] + 216A[y+[2-4(L[Rz]+G)]] = 216(A[2-4x] + A[[2-4x]+[2-4(L[x AND 15]-&746E)]]

Our product function is then xy = 2m(x)+m(y)g(L[2-m(x)x]+L[2-m(y)y]).
As described above, this doubles our thinness errors (zero in this case) and reduces the sparseness ones to e=9e2/8 = 9*2-16/8 » 2-16. So our final errors are:
 xy LogSparse LogThin A‘logSparse A‘logThin TotalFrac TotalAbs All 2-11 3*2-13 2-16 0 2-10 2-10xy

a fractional error of less than 0.1%. This  may be reduced further by introducing linear interpolation for L or increasing NA.

Case Three - Compromise
The following represents only one possible balance between speed and accurracy. The preceeding rather lengthy discussion should enable the reader to find for herself one more suited to any given application. A sensible approach is often to have several different multiplication routines using the same tables in different ways.
We will take K=2k where k=12 (the largest integer giving 2klog2(2k­1)<216).
We will use an NL-entry 16-bit LUT L[i]=[2klog2(i)| 1<i<NL ; and an NA=2a-entry b-bit LUT A[i]=[2b-2^(16-k)+(2^(16-a))i/K| 0<i<NA.
We will use the contraction function l1(x)=[2-m(x)x]

where m(x) is some shift designed so that m(x)<NA.
Our composite contraction function for A is
a1(x)=[2a-16(x MOD 216)]
and our expansion function is
a2(y,x) = 232-b+(x DIV 216)2(16-k)y.
Since L(NL-1) is greater than 215 it is possible that adding L(x) to L(y) will generate a result greater than 216.
Writing z for (L[2-m(x)x]+L[2-m(y)y])MOD 216 ; C for (L[2-m(x)x]+L[2-m(y)y])DIV 216  (=0 or 1, the "carry") we have:  xy » 2m(x)+m(y)+32-b+2^(16-k)CA[2a-16z].
Our errors in xy are:

• 4/NL fractional error due to sparseness of L;
• 21/K-1 = 2-kln(2) fractional error due to thinness of L;
• 22^16-a-k-1 = 216-a-kln(2) fractional error due to sparseness of A;
• 231-b absolute error due to thinness of A.
If NL=213,k=11,a=13,b=32 these are 0.05%,0.03%,0.27%, and 0.5% respectively.

Glossary   Contents   Author
Copyright (c) Ian C G Bell 1998
Web Source: www.iancgbell.clara.net/maths or www.bigfoot.com/~iancgbell/maths
18 Nov 2006.