Sierpinski (<1K) 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[23(L[23x]+L[23y])]
    (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):

 These will contribute fractional errors of to xy respectively.
     We can prepare the following table of contributions to the net fractional error in xy of the four error sources:
xypLog
Sparse
Log
Thin
Alog
Sparse
Alog
Thin
Total
Frac
Total
Abs
<216Big2-122-9BigBig
2160232-122-92-123219
2182212-122-92-321219
220 42-1 2-122-92-5 2-1219
222 62-3 2-122-92-7 2-3219
224 82-5 2-122-92-9 2-5219
226102-7 2-122-92-112-7219
228122-9 2-122-92-132-8220
230142-112-122-92-152-9221
232162-122-122-92-172-9223

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 As 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]
xypLog
Sparse
Log
Thin
Alog
Sparse
Alog
Thin
Total
Frac
Total
Abs
216  022  2-122-92-1 22    218
218  220  2-122-92-3 20    218
220  42-2 2-122-92-5 2-2   218
222  62-4 2-122-92-7 2-4   218
224  82-6 2-122-92-9 2-6   218
226 102-8 2-122-92-113*2-9 3*217
228 122-102-122-92-133*2-103*218
230 142-122-122-92-152-9   221
232 162-132-122-92-172-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:
xypLog
Sparse
Log
Thin
Alog
Sparse
Alog
Thin
Total
Frac
Total
Abs
All2-112-122-802-82-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:
xyLog
Sparse
Log
Thin
Alog
Sparse
Alog
Thin
Total
Frac
Total
Abs
All2-113*2-132-1602-102-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(2k1)<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:

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.