Logarithms

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 log_a(x) for this number we have x = a^log_a(x).
    The inverse of log_a is written antilog_a and we have antilog_a(y) = a^y.
    Logarithms are useful because xy = a^log(x)a^log(y) = a^{log(x)+ log(y)} and so log_a(xy) = log_a(x) + log_a(y).
    Further x^y = (a^log(x))^y = a^ylog(x) and so log_a(x^y) = ylog_a(x).
    Combining these two results we also have log_a(x/y) =log_a(x)log_a(y).
    Clearly log_a(a)=1 and log_a(1)=0 and so antilog_a(1)=a and antilog_a(0)=1.
    Common values for a are two, ten, and the natural base e » 2.718. We can change between these using: log_b(x) = log_a(x)/log_a(b).
    We write ln(x) for log_e(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)=Klog_a(x) is d and that the maximum fractional error of A(x) approximating A(x)=a^x/K is e. Noting that the first two derivatives of L(x)=Klog_a(x) are L'(x) = K/( ln(a)x) ; L"(x) = -K/( ln(a)x²) and that the first two derivatives of A(x)=a^x/K are A'(x) = ln(a)a^x/K/K ; A"(x) = ln(a)2a^x/k/K², we can write:

A(L(x)+L(y))	£	A(L(x)+L(y))(1+e)
	£	A(L(x)+L(y)+2d>)(1+e)
	=	`a`^2d/K(1+e)A(L(x)+L(y))
	=	`a`^2d/K(1+e)xy.

Similarly A(L(x)+L(y)) ³ a^2d/K(1-e)xy.
    Now a^g » 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 N² 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 N_L,W_L,N_A,W_A 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 N_L and N_A powers of two.
    Let us tentatively assign 16 Kbytes = 2¹⁴ to the log table and consider W_L=2, N_L=2¹³.
    We will use X = 2^-3Z_2¹⁶ as a representative portion of Z_2¹⁶ and a contraction mapping l1(x)=[2^-3x] obtained by combining l11(x)=2^-3x with l₁₂(x)=[x].
    We tabulate L[i]=[Klog₂(i)| i Î Z_2¹⁶-{0}, L[0]=0.
    (This is better than setting X_L^*=X_L and l12(x)=[2^-3x] since then we must tabulate [Klog₂(2³i)| which is just the same table with 3K added to every entry).
    We choose K so that L[N_L] is as close as possible to 2¹⁵ while still below it (so that the sum of two table entries is always less than 2¹⁶). Since log₂(2¹³)=13 we set K=2¹⁵/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 N_A>2¹⁶ and we will consider W_A=2, N_A=2¹³.
    We need XA= 2*(dynamic range of L)= 6K+Z_2¹⁶_-1. We will take X_A^*=Z_2¹⁶-1 as a representative portion of this using a₁₁(x)=(x-6K). Combining this with a₁₂(x)=[2^-3x] yields a contraction function of a₁(x)=[2^-3(x-6K)].
    To tabulate A(x)=2^x/K over X_A* = Z_2¹⁶-1 using a₁₂(x)=[2^-3x] we would store A[i] = 2^8i/K. However, this has a dynamic range of 1 to 2^(2¹⁶)/K = 2²⁶ so to to avoid exceeding our 2-byte width we tabulate
A[i] = [2^-102^8i/K|.
    Our expansion function is the composition of a₂₁(y) = 2¹⁰y and a₂₂(y,x) = 2⁶y (compensating for a₁₁(x) = x-6K), ie. a₂(y,x) = 2¹⁶y.
    We thus have
    xy = 2¹⁶A[2³(L[2³x]+L[2³y])]
    (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 2^m, x=2¹⁶ y=2^m-16 gives the gretest fractional error.
    This will be due to (writing p=m-16):

An absolute sparseness error in L of Klog₂(1+1/2^P-3)+Klog₂(1+1/2^16-3)»Klog₂(1+2^3-p+2^-13);
An absolute thinness error in L of 2^-1+2^-1=1;
A fractional sparseness error in A of approx. 2³ln(2)/K;
A fractional thinness error in A of 2¹⁵/xy = 2^15-m = 2^-1-p.

These will contribute fractional errors of

2K( log(1+2^(3-p)+2^-13))/K-1 = 2^3-p+2^-13;
2^1/K-1 = 1.13*2^-12 » 2^-12;
2³ln(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	Log Sparse	Log Thin	A‘log Sparse	A‘log Thin	Total Frac	Total Abs
<2¹⁶	Big	2^-12	2^-9	Big	Big
2¹⁶	0	2³	2^-12	2^-9	2^-1	2³	2¹⁹
2¹⁸	2	2¹	2^-12	2^-9	2^-3	2¹	2¹⁹
2²⁰	4	2^-1	2^-12	2^-9	2^-5	2^-1	2¹⁹
2²²	6	2^-3	2^-12	2^-9	2^-7	2^-3	2¹⁹
2²⁴	8	2^-5	2^-12	2^-9	2^-9	2^-5	2¹⁹
2²⁶	10	2^-7	2^-12	2^-9	2^-11	2^-7	2¹⁹
2²⁸	12	2^-9	2^-12	2^-9	2^-13	2^-8	2²⁰
2³⁰	14	2^-11	2^-12	2^-9	2^-15	2^-9	2²¹
2³²	16	2^-12	2^-12	2^-9	2^-17	2^-9	2²³

Reducing N_A to 2¹² and doubling W_A (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 W_L, (ii) Double N_L, (iii) Double W_A, or (iv) Double N_A.
(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	Log Sparse	Log Thin	A‘log Sparse	A‘log Thin	Total Frac	Total Abs
2¹⁶	0	2²	2^-12	2^-9	2^-1	2²	2¹⁸
2¹⁸	2	2⁰	2^-12	2^-9	2^-3	2⁰	2¹⁸
2²⁰	4	2^-2	2^-12	2^-9	2^-5	2^-2	2¹⁸
2²²	6	2^-4	2^-12	2^-9	2^-7	2^-4	2¹⁸
2²⁴	8	2^-6	2^-12	2^-9	2^-9	2^-6	2¹⁸
2²⁶	10	2^-8	2^-12	2^-9	2^-11	3*2^-9	3*2¹⁷
2²⁸	12	2^-10	2^-12	2^-9	2^-13	3*2^-10	3*2¹⁸
2³⁰	14	2^-12	2^-12	2^-9	2^-15	2^-9	2²¹
2³²	16	2^-13	2^-12	2^-9	2^-17	2^-9	2²³

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 N_A=N_L=2¹³, W_A=W_L=2. We will set K=2¹⁵/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 l₁₁(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)ÎX_L*=[1,2¹³).
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 N_L/2 so our absolute sparseness error in L is at most Klog₂(1+2/N_L)»2K/(N_Lln(2)) contributing a fractional error in xy of 4/N_L, ie. 2^-11 in this case.
    Rather than incorporating an addition of Km(x) into l2 we multiply the final result by 2^Km(x)/K=2^m(x) thus:
    xy = 2^m(x)+m(y)+16A[2^-3(L[2^-m(x)x]+L[2^-m(y)y])].
where L[i]=[Klog₂(i)| iÎZ_N-{0}, A[i]=[2^-102^8i/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 N_L to 2¹⁰ (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 W_A at the expense of N_A    eliminates them entirely but doubles the fractional Antilog Thiness errors to 2^-8 for all products (including those less than 2¹⁶). At this stage we have:

`xy`	p	Log Sparse	Log Thin	A‘log Sparse	A‘log Thin	Total Frac	Total Abs
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=2⁴, p=4.
    We set y=[2^-4x], z=2^-4x-y, so that x=2⁴(y+z) and express z as R_z/2⁴, ie. R_z=bottom 4 bits of x.
    We set G=K(log₂(2^16/K-1)-4)»-11.83K»-&746E and define

g(x)	=	A(2⁴y) + A(2⁴y + Klog₂(R_z) + G)
	»	A(2⁴y) + A(24y + L[R_z] + G)
	»	2¹⁶A[y] + 2¹⁶A[y+[2^-4(L[R_z]+G)]]
	=	2¹⁶(A[2^-4x] + A[[2^-4x]+[2^-4(L[x AND 15]-&746E)]]

Our product function is then xy = 2^m(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=9e²/8 = 9*2^-16/8 » 2^-16. So our final errors are:

`xy`	Log Sparse	Log Thin	A‘log Sparse	A‘log Thin	Total Frac	Total Abs
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 N_A.

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=2^k where k=12 (the largest integer giving 2^klog₂(2^k1)<2¹⁶).
    We will use an N_L-entry 16-bit LUT L[i]=[2^klog₂(i)| 1<i<N_L ; and an N_A=2^a-entry b-bit LUT A[i]=[2^{b-2^(16-k)+(2^(16-a))i/K}| 0<i<N_A.
    We will use the contraction function l₁(x)=[2^-m(x)x]

where m(x) is some shift designed so that m(x)<N_A.
    Our composite contraction function for A is
    a₁(x)=[2^a-16(x MOD 2¹⁶)]
and our expansion function is
    a₂(y,x) = 2^{32-b+(x DIV 2¹⁶)2^(16-k)}y.
    Since L(N_L-1) is greater than 2¹⁵ it is possible that adding L(x) to L(y) will generate a result greater than 2¹⁶.
    Writing z for (L[2^-m(x)x]+L[2^-m(y)y])MOD 2¹⁶ ; C for (L[2^-m(x)x]+L[2^-m(y)y])DIV 2¹⁶ (=0 or 1, the "carry") we have: xy » 2^{m(x)+m(y)+32-b+2^(16-k)C}A[2^a-16z].
    Our errors in xy are:

4/N_L fractional error due to sparseness of L;
2^1/K-1 = 2^-kln(2) fractional error due to thinness of L;
2^{2^16-a-k-1} = 2^16-a-kln(2) fractional error due to sparseness of A;
2^31-b absolute error due to thinness of A.

If N_L=2¹³,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.