|
John DeVos presents a compact logarithm algorithm that trades off accuracy for size and performance.
March 01, 1996
URL:http://www.drdobbs.com/database/algorithm-alley/184409852
John is an engineer at Ircon, a manufacturer of optical thermometers. He can be contacted at [email protected].
When you become accustomed to high-level languages, it is easy to forget that the whole point of using them is to free you from worrying about the bare-metal details of processing. In a low-level embedded environment, however, a language like C can be a reminder of exactly how free from this worry we've become. In these environments, questions such as "Is this causing the right kind of divide?" and "Is the For or the While smaller?" crop up regularly. It is also likely that before too long, you will miss a library routine or two.
This was the situation I faced while calculating the natural logarithm, ln(x), in a program for an embedded controller; see Figure 1. My resources were limited to internal RAM and PROM, and timing limitations were a significant consideration. An "infinite" series (the kind a high-level-language library routine probably uses) seemed like a good starting point, but unfortunately would have required the calculation of a large number of terms to provide accuracy over a useful range of x. This would have resulted in a large, time-consuming routine.
I suspected, however, that if I could tolerate less-than-stellar accuracy, there might be a compact way to calculate logarithms. Converting logs to different bases is trivial, and log-base-2, log2(x), is easier to get at in a binary environment than ln(x) or log10(x).
In my application, it was necessary to evaluate an expression of the form in Figure 1. This expression can be rewritten assuming the use of log2(x), as in Figure 2, in which the constant ln(2) can be included in the constants A and B, making run-time base conversion unnecessary.
In another part of the same program, I needed an exponential function only requiring rough accuracy. The function was implemented using the equivalence of the expressions 2x and 1<<x (one left-shifted x bits). This equivalence suggests that a rough estimate of log2(x) could be obtained by determining the number of binary digits in x. Refining this idea slightly reveals that the integer part of log2(x) is equal to the position of the most-significant bit in x; the least-significant bit (LSB) = 0. I'll refer to this initial estimate as "L0(x)."
Figure 3 shows that the error between L0(x) and log2(x) depends on x's relationship to the surrounding powers of 2 and not its absolute magnitude, suggesting that an additive correction based on this relative value of x is in order. I would like to report that I found a correction after some methodical analysis, but the truth is that I guessed a workable solution on the second try. Trials with linear interpolation showed that the residual error after such a correction was not only large, but difficult to approximate and therefore difficult to cancel.
An interesting, and in this case useful, property of logarithms is that the slope of a logarithmic curve is easy to calculate. Instead of using a straight line to interpolate, why not use a slope based on that of the function you hope to match? Examination of the plot of log2(x) and L0(x) shows that a linear addition to L0(x) would be about right if the slope used was that of the curve log2(x) at an "effective" x half-way between x and the power of 2 below it. The slope of log2(x) is given by Figure 4(a) and the effective x is simply the average of x and 2L0(x); see Figure 4(b). Combining these expressions with the relative x value of (x-2L0(x)) gives the first correction term C0(x); see Figure 4(c).
After this correction, C0(x) has about the right shape, but is increasingly low as x increases from one power of 2 to the next. There is also a discontinuity as x approaches each power of 2 from below. The limit of C0(x) as x->2n from below evaluates to 2/(3*ln(2)), so dividing C0(x) by this factor would make the discontinuity vanish. The resulting expression for the correction term C1(x) resembles a textbook problem in which a complicated looking expression evaluates to an integer, as in Figure 5(a), and the first corrected estimate in Figure 5(b).
A plot of the difference L1(x)-log2(x) shows that the residual error is (conveniently) very nearly parabolic between powers of 2 with the apex of each curve lying very close to half-way between the powers of 2; see Figure 6. Also, the height of each apex, like the error in L0(x), is independent of x.
The expression in Figure 7(a) describes a parabola that closely approximates this residual error, while Figure 7(b) is the twice-corrected approximation. Since this parabola does not exactly match the actual error, selecting PA to equal the maximum value of the residual error results in L2(x) having error with a positive average. In my application, the constant PA was chosen numerically to make the remaining error symmetric in amplitude. Although the average error over the range of x is nonzero (albeit small), this approach minimizes the obligatory +/- error spec.
Listing One is the log2(x) routine implemented in Microsoft Visual C/C++ 1.0. Although this source code works with the Intel iC-96 compiler I use for embedded development, the embedded version of the routine uses a smaller set of local variables (to conserve stack space) and is therefore less readable. This routine works for signed integers (2-bytes) although long integers could be accommodated without much additional code. Table 1 summarizes the routine's performance.
Although I haven't tested it, I suspect that an algorithm similar to this one could be developed for floating-point numbers that would probably be faster and smaller (if less accurate) than a library routine or an algorithm based on a series expansion. Another possible use for such an algorithm is calculating the square root of x, since Figure 8 is true for an arbitrary base B.
Argument range tested: 3 to 32767 Highest relative error: @ log2(3) is 0.0153 percent Lowest relative error: @ log2(7) is -0.0102 percent Highest absolute error: @ log2(10815) is 0.000438 Lowest absolute error: @ log2(15199) is -0.000514 Average absolute error: 0.000002 RMS absolute error: 0.000272
A
f(x) = ----------- for constants A, B, or C
B+C*ln(x)
A*ln(2)
f(x) = ------------------
B*ln(2)+C*log2(x)
(a)
d 1
m(x) = ---(log2(x)) = -------
dx ln(2)*x
(b)
x+(2L0(x))
xeff = ----------
2
(c)
C0(x) = m(xeff)*(x-2L0(x))
2*(x-2L0(x))
C0(x) = ----------------
ln(2)*(x+2L0(x))
(a)
(x-2L0(x))
C1(x) = 3*----------
(x+2L0(x))
(b)
L1(x) = L0(x)+C1(x)
/* log2.h */
/* Prototypes */
long log2K (int);
/* Defines & Macros */
#define LOG2_ERROR 666L /* Impossible result */
#define BASE_FACTOR 16384
#define BASE_SHIFT 14
#define ERR_C1 242 /* PA * BASE_FACTOR */
#define ERR_C2 726 /* 3 * ERR_C1 */
/* log2K -- Returns the log-base-2 of the argument relative to
BASE_FACTOR. If the parameter is not positive, it returns a defined
impossible result. */
long log2K (int x) {
int j, hiBitPos, C1, C2;
long L0;
if (x <= 0) return LOG2_ERROR;
/* Zero-order estimate */
j = 1;
while (x >> j++);
hiBitPos = j - 2;
j = 1 << hiBitPos; /* j = 2 ^ hiBitPos(x) */
L0 = (long)hiBitPos << BASE_SHIFT;
if (x == j) return L0; /* Catch powers of 2 */
/* First correction */
C1 = (int)(((((long)x - j) * 3) << BASE_SHIFT) / ((log)x + j));
/* Second correction */
j >>= 1; /* j = 2 ^ (hiBit(x) - 1) */
C2 = (int)(((long)x * ERR_C1) / j) - ERR_C2;
C2 = ERR_C1 - (int)(((long)C2 * C2) / ERR_C1);
return (L0 + C1 - C2);
}
Terms of Service | Privacy Statement | Copyright © 2024 UBM Tech, All rights reserved.