# tl;dr

No.

Well, perhaps, but it can be the lesser of two evils.

It depends.

Leave a reply

In this post I finally get to the punch line and construct a set of approximations to (actually ) that cover a range of speeds and accuracies. As execution time decreases, so does the accuracy. You can choose an approximation with a tradeoff that meets your needs.

When making a fast but approximate function, the design parameter is the form of the approximating function. Is it a polynomial of degree 3? Or the ratio of two linear functions? The table below has a systematic list of possibilities: the top half of the table uses only multiplication, the bottom half uses division.

Table 1: Approximating Functions

The *rewritten* column rewrites the expression in a more efficient form. The expressions are used in the approximation procedure as follows: to get the approximation to , first is reduced to the interval , then is substituted into the expression. As I explained in Part I and Part II, this procedure gives good accuracy and avoids floating-point cancellation problems.

For each form, the minimax values of are determined—that is, the values that minimize the maximum relative error. The *bits* column gives the bits of accuracy, computed as where is the maximum relative error. This value is computed using the evaluation program from my Part II post, invoked as .

The *“cost”* column is the execution time. In olden days, floating-point operations were a lot more expensive than integer instructions, and they were executed sequentially. So cost could be estimated by counting the number of floating-point additions and multiplications. This is no longer true, so I estimate cost using the evaluation program. I put cost in quotes since the numbers apply only to my MacBook. The numbers are normalized so that the cost of log2f in the standard C library is 1.

The table shows that using division increases accuracy about the same amount as does adding an additional free parameter. For example, line 4 has 11.3 bits of accuracy with four free parameters *a*, *b*, *c*, and *d*. Using division, you get similar accuracy (11.6) with line 8, which has only three free parameters. Similarly, line 6 has about the same accuracy as line 10, and again line 10 has one less free parameter than line 6.

It’s hard to grasp the cost and time numbers in a table. The scatter plot below is a visualization with lines 2–10 each represented by a dot.

You can see one anomaly—the two blue dots that are (roughly) vertically stacked one above the other. They have similar cost but very different accuracy. One is the blue dot from line 9 with cost .43 and 7.5 bits of accuracy. The other is from line 3 which has about the same cost (0.42) but 8.5 bits of accuracy. So clearly, line 9 is a lousy choice. I’m not sure I would have guessed that without doing these calculations.

Having seen the problem with line 9, it is easy to explain it using the identity . If is approximated by on , then in , and this works out to be

which means that the denominators must be roughly equal. Multiplying each by , this becomes . The only way this can be true is if and are very large, so that it’s as if the term doesn’t exist and the approximation reduces to the quotient of two linear functions. And that is what happens. The optimal coefficients are on the order of , which makes (now writing in terms of ) . Plugging in the optimal values , , gives , , which are very similar to the coefficients in line 7 shown in Table 2 later in this post. In other words, the optimal rational function in line 9 is almost identical to the one in line 7. Which explains why the bit accuracy is the same.

In the next plot I remove line 9, and add lines showing the trend.

The lines show that formulas using division outperform the multiplication-only formulas, and that the gain gets greater as the formulas become more complex (more costly to execute).

You might wonder: if one division is good, are two divisions even better? Division adds new power because a formula using division can’t be rewritten using just multiplication and addition. But a formula with two divisions can be rewritten to have only a single division, for example

Two divisions add no new functionality, but could be more efficient. In the example above, a division is traded for two multiplications. In fact, using two divisions gives an alternative way to write line 10 of Table 1. On my Mac, that’s a bad tradeoff: the execution time of the form with 2 divisions increases by 40%.

Up till now I’ve been completely silent about how I computed the minimax coefficients of the functions. Or to put it another way, how I computed the values of , , , etc. in Table 1. This computation used to be done using the Remez algorithm, but now there is a simpler way that reduces to solving a convex optimization problem. That in turn can then be solved using (for example) CVX, a Matlab-based modeling system.

Here’s how it works for line 8. I want to find the minimax approximation to . As discussed in the first post of this series, it’s the relative error that should be minimized. This means solving

This is equivalent to finding the smallest for which there are , , and satisfying

For the last equation, pick . Of course this is not exactly equivalent to being true for all but it is an excellent approximation. The notation may be a little confusing, because are constants, and , and are the variables. Now all I need is a package that will report if

has a solution in , , . Because then binary search can be used to find the minimal . Start with that you know has no solution and that is large enough to guarantee a solution. Then ask if the above has a solution for . If it does, replace ; otherwise, . Continue until has the desired precision.

The set of satisfying the equation above is convex (this is easy to check), and so you can use a package like CVX for Matlab to quickly tell you if it has a solution. Below is code for computing the coefficients of line 8 in Table 1. This matlab/CVX code is modified from http://see.stanford.edu/materials/lsocoee364a/hw6sol.pdf.

It is a peculiarity of CVX that it can report for a value of , but then report for a smaller value of . So when seeing I presume there is a solution, then decrease the upper bound and also record the corresponding values of , , and . I do this for each step of the binary search. I decide which to use by making an independent calculation of the minimax error for each set .

You might think that using a finer grid (that is, increasing so that there are more ) would give a better answer, but it is another peculiarity of CVX that this is not always the case. So in the independent calculation, I evaluate the minimax error on a very fine grid that is independent of the grid size given to CVX. This gives a better estimate of the error, and also lets me compare the answers I get using different values of . Here is the CVX code:

format long format compact verbose=true bisection_tol = 1e-6; m=500; lo=0.70; % check values a little bit below 0.75 hi=1.5; xi = linspace(lo, hi, m)'; yi = log2(xi); Xi = linspace(lo, hi, 10000); % pick large number so you can compare different m Xi = Xi(Xi ~= 1); Yi = log2(Xi); xip = xi(xi >= 1); % those xi for which y = x-1 is positive xin = xi(xi = 0.75); xinn = xi(xi = 1); yin = yi(xi = 0.75); yinn = yi(xi = 0.75); Xin = Xi(Xi = bisection_tol gamma = (l+u)/2; cvx_begin % solve the feasibility problem cvx_quiet(true); variable A; variable B; variable C; subject to abs(A*(xip - 1).^2 + B*(xip - 1) - yip .* (xip - 1 + C)) <= ... gamma * yip .* (xip-1 + C) abs(A*(xin - 1).^2 + B*(xin - 1)- yin .* (xin - 1 + C)) <= ... -gamma * yin .* (xin-1 + C) abs(A*(2*xinn - 1).^2 + B*(2*xinn - 1) - (1 + yinn) .* (2*xinn - 1 + C)) <= ... -gamma * yinn .* (2*xinn - 1 + C) cvx_end if verbose fprintf('l=%7.5f u=%7.5f cvx_status=%s\n', l, u, cvx_status) end if strcmp(cvx_status,'Solved') | strcmp(cvx_status, 'Inaccurate/Solved') u = gamma; A_opt(k) = A; B_opt(k) = B; C_opt(k) = C; lo = (A*(2*Xin - 1).^2 + B*(2*Xin - 1)) ./ (2*Xin - 1 + C) - 1; hi = (A*(Xip - 1).^2 + B*(Xip -1)) ./ (Xip - 1 + C); fx = [lo, hi]; [maxRelErr(k), maxInd(k)] = max(abs( (fx - Yi)./Yi )); k = k + 1; else l = gamma; end end [lambda_opt, k] = min(maxRelErr); A = A_opt(k) B = B_opt(k) C = C_opt(k) lambda_opt -log2(lambda_opt)

Here are the results of running the above code for the expressions in the first table. I don’t bother giving all the digits for line 9, since it is outperformed by line 7.

Table 2: Coefficients for Table 1

So, what’s the bottom line? If you don’t have specific speed or accuracy requirements, I recommend choosing either line 3 or line 7. Run both through the evaluation program to get the cost for your machine and choose the one with the lowest cost. On the other hand, if you have specific accuracy/speed tradeoffs, recompute the cost column of Table 1 for your machine, and pick the appropriate line. The bits column is machine independent as long as the machine uses IEEE arithmetic.

If you want a rational function with more accuracy than line 10, the next choice is cubic/quadratic which gives 20.7 bits of accuracy. That would be with coefficients *A* = 0.1501692, *B* = 3.4226132, *C* = 5.0225057, *D* = 4.1130283, *E* = 3.4813372.

Finally, I’ll close by giving the C code for line 8 (almost a repeat of code from the first posting). This is bare code with no sanity checking on the input parameter . I’ve marked the lines that need to be modified if you want to use it for a different approximating expression.

float fastlog2(float x) // compute log2(x) by reducing x to [0.75, 1.5) { /** MODIFY THIS SECTION **/ // (x-1)*(a*(x-1) + b)/((x-1) + c) (line 8 of table 2) const float a = 0.338953; const float b = 2.198599; const float c = 1.523692; #define FN fexp + signif*(a*signif + b)/(signif + c) /** END SECTION **/ float signif, fexp; int exp; float lg2; union { float f; unsigned int i; } ux1, ux2; int greater; // really a boolean /* * Assume IEEE representation, which is sgn(1):exp(8):frac(23) * representing (1+frac)*2^(exp-127). Call 1+frac the significand */ // get exponent ux1.f = x; exp = (ux1.i & 0x7F800000) >> 23; // actual exponent is exp-127, will subtract 127 later greater = ux1.i & 0x00400000; // true if signif > 1.5 if (greater) { // signif >= 1.5 so need to divide by 2. Accomplish this by // stuffing exp = 126 which corresponds to an exponent of -1 ux2.i = (ux1.i & 0x007FFFFF) | 0x3f000000; signif = ux2.f; fexp = exp - 126; // 126 instead of 127 compensates for division by 2 signif = signif - 1.0; lg2 = FN; } else { // get signif by stuffing exp = 127 which corresponds to an exponent of 0 ux2.i = (ux1.i & 0x007FFFFF) | 0x3f800000; signif = ux2.f; fexp = exp - 127; signif = signif - 1.0; lg2 = FN; } // last two lines of each branch are common code, but optimize better // when duplicated, at least when using gcc return(lg2); }

In response to comments, here is some sample code to compute logarithms via table lookup. A single precision fraction has only 23 bits, so if you are willing to have a table of 2^{23} floats (2^{25} bytes) you can write a logarithm that is very accurate and very fast. The one thing to watch out for is floating-point cancellation, so you need to split the table into two parts (see the `log_biglookup()`

code block below).

The `log_lookup()`

sample uses a smaller table. It uses linear interpolation, because ordinary table lookup results in a step function. When *x* is near 0, log(1+*x*) ≈ *x* is linear, and any step-function approximation will have a very large relative error. But linear interpolation has a small relative error. In yet another variation, `log_lookup_2()`

uses a second lookup table to speed up the linear interpolation.

float log_lookup(float x) // lookup using NDIGIT_LOOKUP bits of the fraction part { int exp; float lg2, interp; union { float f; unsigned int i; } ux1, ux2; unsigned int frac, frac_rnd; /* * Assume IEEE representation, which is sgn(1):exp(8):frac(23) * representing (1+frac)*2ˆ(exp-127) Call 1+frac the significand */ // get exponent ux1.f = x; exp = ((ux1.i & 0x7F800000) >> 23); // -127 done later // top NDIGITS frac = (ux1.i & 0x007FFFFF); frac_rnd = frac >> (23 - NDIGITS_LOOKUP); // for interpolating between two table values ux2.i = (frac & REMAIN_MASK) << NDIGITS_LOOKUP; ux2.i = ux2.i | 0x3f800000; interp = ux2.f - 1.0f; if (frac_rnd < LOOKUP_TBL_LN/2) { lg2 = tbl_lookup_lo[frac_rnd] + interp*(tbl_lookup_lo[frac_rnd+1] - tbl_lookup_lo[frac_rnd]); return(lg2 + (exp - 127)); } else { lg2 = tbl_lookup_hi[frac_rnd] + interp*(tbl_lookup_hi[frac_rnd+1] - tbl_lookup_hi[frac_rnd]); return(-lg2 + (exp - 126)); } } static float log_lookup_2(float x) // use a second table, tbl_interp[] { int exp; float lg2; union { float f; unsigned int i; } ux1; unsigned int frac, frac_rnd, ind; /* * Assume IEEE representation, which is sgn(1):exp(8):frac(23) * representing (1+frac)*2ˆ(exp-127) Call 1+frac the significand */ // get exponent ux1.f = x; exp = ((ux1.i & 0x7F800000) >> 23); // -127 done later // top NDIGITS frac = (ux1.i & 0x007FFFFF); frac_rnd = frac >> (23 - NDIGITS_LOOKUP_2); // for interpolating between two table values ind = frac & REMAIN_MASK_2; // interp = tbl_inter[ind] if (frac_rnd < LOOKUP_TBL_LN_2/2) { lg2 = tbl_lookup_lo[frac_rnd] + tbl_interp[ind]*(tbl_lookup_lo[frac_rnd+1] - tbl_lookup_lo[frac_rnd]); return(lg2 + (exp - 127)); } else { lg2 = tbl_lookup_hi[frac_rnd] + tbl_interp[ind]*(tbl_lookup_hi[frac_rnd+1] - tbl_lookup_hi[frac_rnd]); return(-lg2 + (exp - 126)); } } static float log_biglookup(float x) // full lookup table with 2^23 entries { int exp; float lg2; union { float f; unsigned int i; } ux1; unsigned int frac; /* * Assume IEEE representation, which is sgn(1):exp(8):frac(23) * representing (1+frac)*2ˆ(exp-127) Call 1+frac the significand */ // get exponent ux1.f = x; exp = ((ux1.i & 0x7F800000) >> 23); // -127 done later frac = (ux1.i & 0x007FFFFF); if (frac < TWO_23/2) { lg2 = tbl_lookup_big_lo[frac]; return(lg2 + (exp - 127)); } else { lg2 = tbl_lookup_big_hi[frac - TWO_23/2]; return(-lg2 + (exp - 126)); } } #define NDIGITS_LOOKUP 14 #define LOOKUP_TBL_LN 16384 // 2^NDIGITS #define REMAIN_MASK 0x1FF // mask with REMAIN bits where REMAIN = 23 - NDIGITS static float *tbl_lookup_lo; static float *tbl_lookup_hi; static void init_lookup() { int i; tbl_lookup_lo = (float *)malloc(( LOOKUP_TBL_LN/2 + 1)*sizeof(float)); tbl_lookup_hi = (float *)malloc(( LOOKUP_TBL_LN + 1)*sizeof(float)); tbl_lookup_hi = tbl_lookup_hi - LOOKUP_TBL_LN/2; for (i = 0; i <= LOOKUP_TBL_LN/2; i++) // <= not < tbl_lookup_lo[i] = log2f(1 + i/(float)LOOKUP_TBL_LN); for (i = LOOKUP_TBL_LN/2; i < LOOKUP_TBL_LN; i++) // log2, not log2f tbl_lookup_hi[i] = 1.0 - log2(1 + i/(float)LOOKUP_TBL_LN); tbl_lookup_hi[LOOKUP_TBL_LN] = 0.0f; } /* two tables */ #define NDIGITS_LOOKUP_2 12 #define LOOKUP_TBL_LN_2 4096 // 2^NDIGITS #define TWO_REMAIN_2 2048 // 2^REMAIN, where REMAIN = 23 - NDIGITS #define REMAIN_MASK_2 0x7FF // mask with REMAIN bits static float *tbl_interp; static void init_lookup_2() { int i; tbl_lookup_lo = (float *)malloc(( LOOKUP_TBL_LN_2/2 + 1)*sizeof(float)); tbl_lookup_hi = (float *)malloc(( LOOKUP_TBL_LN_2/2 + 1)*sizeof(float)); tbl_lookup_hi = tbl_lookup_hi - LOOKUP_TBL_LN_2/2; tbl_interp = (float *)malloc(TWO_REMAIN_2*sizeof(float)); // lookup for (i = 0; i <= LOOKUP_TBL_LN_2; i++) tbl_lookup_lo[i] = log2f(1 + i/(float)LOOKUP_TBL_LN_2); for (i = LOOKUP_TBL_LN_2/2; i < LOOKUP_TBL_LN_2; i++) // log2, not log2f tbl_lookup_hi[i] = 1.0 - log2(1 + i/(float)LOOKUP_TBL_LN_2); tbl_lookup_hi[LOOKUP_TBL_LN_2] = 0.0f; for (i = 0; i < TWO_REMAIN_2; i++) tbl_interp[i] = i/(float)TWO_REMAIN_2; } #define TWO_23 8388608 // 2^23 static float *tbl_lookup_big_lo; static float *tbl_lookup_big_hi; static void init_biglookup() { int i; tbl_lookup_big_lo = (float *)malloc((TWO_23/2) * sizeof(float)); tbl_lookup_big_hi = (float *)malloc((TWO_23/2) * sizeof(float)); for (i = 0; i < TWO_23/2; i++) { tbl_lookup_big_lo[i] = log2f(1 + i/(float)TWO_23); // log2, not log2f tbl_lookup_big_hi[i] = 1.0 - log2(1 + (i+ TWO_23/2)/(float)TWO_23); } }

Powered by QuickLaTeX

In Part I, I assumed all computations were done exactly. It’s time to bring in floating-point rounding error. There’s a key fact about rounding error that might seem surprising: If , then the computation of in floating-point has no rounding error. The surprise comes for those who’ve been taught to beware of rounding error when subtracting two similar quantities. For example, the computation in floating-point of might be far from the true value of when . That’s because can have rounding error, so the subtraction might cancel most of the good bits leaving mainly the incorrect (due to rounding) bits. But in , neither nor 1 have rounding error, so the computation of is exactly equal to when .

Here’s the explanation. Rounding error in *x+y* can happen in two ways. First if *x > y* and *x* has a different exponent from *y*, then *x* will have its fractional part shifted right to make the exponents match, and so *x* might drop some bits. Second, even if the exponents are the same, there may be rounding error if the addition of the fractional parts has a carry-out from the high order bit. In the case of *x − 1*, the exponents are the same if 1 ≤ *x* < 2. And if 1/2 ≤ *x* < 1, the larger number is 1 and it does not drop bits when shifted. So there is no rounding error in *x − 1* if 1/2 < *x* < 2.

The rule of thumb is that when approximating a function with , severe rounding error can be reduced if the approximation is in terms of . For us, so , and the rule suggests polynomials written in terms of rather than . By the key fact above, there is no rounding error in when .

Let me apply that to two different forms of the quadratic polynomials used in Part I: the polynomial can be written in terms of or .

If they are to be used on the interval and I want to minimize relative error, it is crucial that the polynomial be 0 when , so they become

The second equation has no constant term, so they both cost the same amount to evaluate, in that they involve the same number of additions and multiplications.

But one is much more accurate. You can see that empirically using an evaluation program (code below) that I will be using throughout to compare different approximations. I invoked the program as and and got the following:

SPACING IS 1/1024 using x bits 5.5 at x=0.750000 2.06 nsecs nsec/bit=0.372 bits/nsec=2.69 using x-1 bits 5.5 at x=0.750000 2.10 nsecs nsec/bit=0.380 bits/nsec=2.63 SPACING IS 1/4194304 = 2^-22 using x bits 1.7 at x=1-2.4e-07 2.08 nsecs nsec/bit=1.222 bits/nsec=0.82 using x-1 bits 5.5 at x=0.750000 2.12 nsecs nsec/bit=0.384 bits/nsec=2.61

When the approximating polynomials are evaluated at points spaced 1/1024 apart, they have similar performance. The accuracy of both is 5.5 bits, and the one using is slightly slower. But when they are evaluated at points spaced apart, the polynomial using has poor accuracy when is slightly below 1. Specifically, the accuracy is only 1.7 bits when .

To see why, note that when , is summing two numbers that have rounding error, but are of different sizes, since . But is summing two numbers of similar size, since and the sum of the first two terms is about . This is the bad case of subtracting two nearby numbers (cancellation), because they both have rounding error.

I suppose it is an arguable point whether full accuracy for all is worth a time performance hit of about 2%. I will offer this argument: you can reason about your program if you know it has (in this case) 5.5 bits of accuracy on *every* input. You don’t want to spend a lot of time tracking down unexpectedly low accuracy in your code that came about because you used a log library function with poor precision on a small set of inputs.

Here’s some more information on the output of the evalution program displayed above. The first number is accuracy in bits measured in the usual way as where is the maximum relative error. Following is the value of where the max error occured. The execution time (e.g. 2.06 nsecs for the first line) is an estimate of the time it takes to do a single approximation to , including reducing the argument to the interval . The last two numbers are self explanatory.

Estimating execution time is tricky. For example on my MacBook, if the argument must be brought into the cache, it will significantly affect the timings. That’s why the evaluation program brings and into the cache before beginning the timing runs.

For polynomials, using has almost the same cost and better accuracy, so there is a good argument that it is superior to using . Things are not so clear when the approximation is a rational function rather than a polynomial. For example, . Because , the numerator is actually . And because you can multiply numerator and denominator by anything (as long as it’s the same anything), it further simplifies to . This will have no floating-point cancellation, and will have good accuracy even when . But there’s a rewrite of this expression that is faster:

Unfortunately, this brings back cancellation, because when there will be cancellation between and the fraction. Because there’s cancellation anyway, you might as well make a further performance improvement eliminating the need to compute , namely

Both sides have a division. In addition, the left hand side has a multiplication and 2 additions. The right hand side has no multiplications and 2 additions ( is a constant and doesn’t involve a run-time multiplication, similarly for ). So there is one less multiplication, which should be faster. But at the cost of a rounding error problem when .

SPACING IS 1/1024 using x-1 bits 7.5 at x=0.750000 2.17 nsecs nsec/bit=0.289 bits/nsec=3.46 using x bits 7.5 at x=0.750000 2.07 nsecs nsec/bit=0.275 bits/nsec=3.64 SPACING IS 1/4194304 = 2^-22 using x-1 bits 7.5 at x=0.750000 2.24 nsecs nsec/bit=0.298 bits/nsec=3.36 using x bits 1.4 at x=1+2.4e-07 2.09 nsecs nsec/bit=1.522 bits/nsec=0.66

As expected, the rational function that has one less multiplication (the line marked *using x)* is faster, but has poor accuracy when is near 1. There’s a simple idea for a fix. When is small, use the Taylor series, . Using is a subtraction and a multiplication, which is most likely cheaper than a division and two additions, What is the size cutoff? The error in the Taylor series is easy to compute: it is the next term in the series, , so the relative error is about . And I want to maintain an accuracy of 7.5 bits, or . So the cutoff is , where or .

On my MacBook, the most efficient way to implement appears to be . In the evaluation program, most of the are greater than 1, so only the first of the inequalities is executed. Despite this, adding the check still has a high cost, but no more accuracy than using .

SPACING IS 1/1024 using x-1 bits 7.5 at x=0.750000 2.17 nsecs nsec/bit=0.289 bits/nsec=3.46 using x bits 7.5 at x=0.750000 2.07 nsecs nsec/bit=0.275 bits/nsec=3.64 cutoff bits 7.5 at x=0.750000 2.58 nsecs nsec/bit=0.343 bits/nsec=2.91 SPACING IS 1/4194304 = 2^-22 using x-1 bits 7.5 at x=0.750000 2.24 nsecs nsec/bit=0.298 bits/nsec=3.36 using x bits 1.4 at x=1+2.4e-07 2.09 nsecs nsec/bit=1.522 bits/nsec=0.66 cutoff bits 7.5 at x=0.989000 2.60 nsecs nsec/bit=0.347 bits/nsec=2.88

In Part I of this series, I noted that testing whether was in the range can be done with a bit operation rather than a floating-point one. The same idea could be used here. Instead of using the Taylor series when or , use it in a slightly smaller interval

The latter can be converted to bit operations on , the fraction part of x, as follows:

As bit operations, this is

(exp == 0 && (f & 111111100...) == 0) OR (exp = -1 && (f & 11111100...) == 111111000...)

When I tested this improvement ( in the table below). it was faster, but still slower than using , at least on my MacBook.

SPACING IS 1/1024 using x-1 bits 7.5 at x=0.750000 2.17 nsecs nsec/bit=0.289 bits/nsec=3.46 using x bits 7.5 at x=0.750000 2.07 nsecs nsec/bit=0.275 bits/nsec=3.64 cutoff bits 7.5 at x=0.750000 2.58 nsecs nsec/bit=0.343 bits/nsec=2.91 fcutoff bits 7.5 at x=0.750000 2.46 nsecs nsec/bit=0.327 bits/nsec=3.06 SPACING IS 1/4194304 = 2^-22 using x-1 bits 7.5 at x=0.750000 2.24 nsecs nsec/bit=0.298 bits/nsec=3.36 using x bits 1.4 at x=1+2.4e-07 2.09 nsecs nsec/bit=1.522 bits/nsec=0.66 cutoff bits 7.5 at x=0.989000 2.60 nsecs nsec/bit=0.347 bits/nsec=2.88 fcutoff bits 7.5 at x=0.750001 2.44 nsecs nsec/bit=0.325 bits/nsec=3.08

Bottom line: having special case code when appears to significantly underperform computing in terms of .

In the first post, I recommended reducing to instead of because you get one extra degree of freedom, which in turns gives greater accuracy. Rounding error gives another reason for preferring . When , reduction to will have cancellation problems. Recall the function that was optimal for the interval , . When , must be multiplied by two to move into , and then to compensate, the result is . When , , and so you get cancellation. Below are the results of running the evaluation program on . If there was no rounding error, would be accurate to 3.7 bits. As you get closer to 1 () the accuracy drops.

SPACING IS 1/2^19 g bits 3.5 at x=1-3.8e-06 2.05 nsecs nsec/bit=0.592 bits/nsec=1.69 SPACING IS 1/2^20 g bits 2.9 at x=1-9.5e-07 2.05 nsecs nsec/bit=0.706 bits/nsec=1.42 SPACING IS 1/2^21 g bits 2.9 at x=1-9.5e-07 2.05 nsecs nsec/bit=0.706 bits/nsec=1.42 SPACING IS 1/2^22 g bits 1.7 at x=1-2.4e-07 2.05 nsecs nsec/bit=1.203 bits/nsec=0.83

The goal of this series of posts is to show that you can create logarithm routines that are much faster than the library versions and have a minimum guaranteed accuracy for all . To do this requires paying attention to rounding error. Summarizing what I’ve said so far, my method for minimizing rounding error problems is to reduce to the interval and write the approximating expression using , for example ). More generally, the approximating expression would be a polynomial

or rational function

I close by giving the code for the evaluation program that was used to compare the time and accuracy of the different approximations:

#include <stdio.h> #include <stdlib.h> #include <sys/time.h> #include <math.h> /* * Usage: eval [hi reps spacing] * Evaluates an approximation to log2 in the interval [0.125, hi] * For timing purposes, repeats the evaluation reps times. * The evaluation is done on points spaced 1/spacing apart. */ int main(argc, argv) char **argv; { float x; struct timeval start, stop; float lo, hi, delta; int i, j, n, repetitions, one_over_delta; double xd; float *xarr, *lg2arr, *yarr; // parameters lo = 0.125; hi = 10.0; one_over_delta = 4194304.0; // 2^22 repetitions = 1; if (argc > 1) { hi = atof(argv[1]); repetitions = atoi(argv[2]); one_over_delta = atoi(argv[3]); } delta = 1.0/one_over_delta; // setup n = ceil((hi - lo)/delta) + 1; xarr = (float *)malloc(n*sizeof(float)); yarr = (float *)malloc(n*sizeof(float)); lg2arr = (float *)malloc(n*sizeof(float)); i = 0; for (xd = lo; xd <= hi; xd += delta) { x = xd; if (x == 1.0) // relative error would be infinity continue; xarr[i] = x; lg2arr[i++] = log2(x); } if (i >= n) // assert (i < n) fprintf(stderr, "Help!!!\n"); n = i; /* cache-in xarr[i], yarr[i] */ yarr[0] = 0.0; for (i = 1; i < n; i++) { yarr[i] = xarr[i] + yarr[i-1]; } fprintf(stderr, "cache-in: %f\n\n", yarr[n-1]); // to foil optimizer gettimeofday(&start, 0); for (j = 0; j < repetitions; j++) { for (i = 0; i < n; i++) { yarr[i] = approx_fn(xarr[i]); } } gettimeofday(&stop, 0); finish(&start, &stop, "name ", n, repetitions, xarr, yarr, lg2arr); exit(0); } // convert x to string, with special attention when x is near 1 char *format(float x) { static char buf[64]; float y; if (fabs(x - 1) > 0.0001) sprintf(buf, "%f", x); else { y = x-1; if (y < 0) sprintf(buf, "1%.1e", y); else sprintf(buf, "1+%.1e", y); } return(buf); } void finish(struct timeval *start, struct timeval *stop, char *str, int n, int repetitions, float *xarr, float *yarr, float *lg2arr) { double elapsed; // nanosecs float max, rel; int maxi, i; double bits; elapsed = 1e9*(stop->tv_sec - start->tv_sec) + 1000.0*(stop->tv_usec - start->tv_usec); max = 0.0; for (i = 0; i < n; i++ ) { rel = fabs( (yarr[i] - lg2arr[i])/lg2arr[i]); if (rel > max) { max = rel; maxi = i; } } bits = -log2(max); elapsed = elapsed/(n*repetitions); printf("%s bits %4.1f at x=%s %.2f nsecs nsec/bit=%.3f bits/nsec=%.2f\n", str, bits, format(xarr[maxi]), elapsed, elapsed/bits, bits/elapsed); }

Powered by QuickLaTeX