5 ms·

From a practical point of view, in case you have Mathematica (which is a proprietary software unfortunately), you can do this easily for any function in a gener

by tagrun 11y ago

From a practical point of view, in case you have Mathematica (which is a proprietary software unfortunately), you can do this easily for any function in a general way: http://reference.wolfram.com/language/FunctionApproximations/tutorial/FunctionApproximations.html http://reference.wolfram.com/language/FunctionApproximations... (MiniMaxApproximation minimizes the relative maximum error) For the case discussed in write up, it gives MiniMaxApproximation[Log2[x + 1], {x, {1, 2}, 2, 0}] log2(x+1) ~ 0.177392 + 0.943626 x - 0.120232 x^2 with a maximum error of -0.000786282. In the 4th order (that is, (4,0)), it gives 0.0399284 + 1.27027 x - 0.391233 x^2 + 0.0909038 x^3 - 0.00986308 x^4 max error: -4.83031 10^-6. However, if the division isn't too expensive, a (2,2) approximant also works for an even better precision: (0.00383323 + 1.42412 x + 0.336161 x^2)/(1 + 0.704317 x + 0.0598015 x^2) max error: -1.33902 10^-7 As for [0.75, 1.5) interval, I get 0.108715 + 1.05621 x - 0.165315 x^2 which gives a max error -0.000651425 and the max relative error is around 0.0006, which is again way better than what he gives in the write up, although the error is non zero at points 0.75 and 1.5 (which is not really a necessary constraint, but an artifact of the particular polynomial he chose). The overall plot for relative error in a wider range is smooth. Bottom line: 1) Unless you are an expert in numerical methods, consult to a software/article/book by an expert (or to the expert him/herself if you can) in the field at some point. 2) While I appreciate his efforts and the write up, I wouldn't recommend using his polynomials in practice. (BTW, one can use frexp instead of bit-twiddling to make the C code friendlier since not everyone is familiar with the bit layout of IEEE 754 numbers.)

stagger87 11y ago
I am assuming a vectorized version of the function will appear in a later update, which is why frexp is not used. Range reduction -> polynomial -> scale is the process an expert/book will recommend for log functions.
- tagrun 11y ago
  That is the general procedure for any function of course. The problem usually lies in implementing these steps, which requires an expert.
tagrun 11y ago
I realized that I missed the part that he's talking about log2(x) rather than log(x+1), which changes the interval. Since MiniMaxApproximation gets quirky when the function is 0 in the given interval (it minimizes the relative error, so a division by 0 occurs, unless you cancel it), I'm rather using the interval [0.25, 0.5] with t=1/4: -3.64631 + 7.93796 x - 5.30406 x^2 This gives a maximum relative error of 0.00334047, which is still an order better than his. As for the (2,2) approximant: (-5.83261 - 16.7012 x + 22.394 x^2)/(1 + 11.2638 x + 7.81124 x^2) max rel error: 1.69763 10^-6 I checked several other values of t, and the max rel error for (2,0) approximant is roughly the same. This leads me to believe that the reason he is getting 10 times larger error is because he is imposing that errors should vanish at certain points. Not only this increases the max rel error, but it also distorts the smoothness of the function in a bad way (the kinks in the plots are a bad sign in general).
- gjm11 11y ago
  Your polynomial, since it doesn't get it exactly right at 1/4 & 1/2, will have an enormous relative error when you use it for computing logs of numbers very close to 1. Goldberg's initial function b ignores this issue, and a glance at the graph of its absolute errors will confirm that it's the minimax poly for absolute error or very close to it. His next candidate t half-fixes this, and again looking at the graph of the relative error makes it clear that it's right. (Note that avoiding intervals containing zero will indeed allow you to get better max relative errors -- keeping the relative error finite where the function is zero effectively reduces your polynomial degree by 1.) Then finally he fully fixes it with the polynomial g, solving the more interesting problem of minimizing the max relative error for the log function everywhere, and chooses a better functional form by shifting the interval he reduces to with the polynomial s. (The fact that his earlier polynomials do the wrong thing is clearly not a mistake but a paedagogical strategy: start with something simple and incrementally improve it.) My apologies for spelling all this out, but it really appears that you're accusing him of making elementary mistakes without actually understanding what he's doing. "Imposing that errors should vanish at certain points" is not a mistake, it's a necessity if you're implementing log by argument reduction plus local approximation and you care about relative error in the resulting approx_log function, and you are getting better numbers than Goldberg by solving the wrong problem. Your (2,2) rational approximant will be much more expensive to compute than the polynomial and won't yield the speed gains Goldberg is looking for. And, again, it will give you a log function with infinite relative error at 1. [EDITED to fix a typo and clarify something.]

by tagrun 11y ago

stagger87 11y ago

I am assuming a vectorized version of the function will appear in a later update, which is why frexp is not used. Range reduction -> polynomial -> scale is the process an expert/book will recommend for log functions.

tagrun 11y ago
That is the general procedure for any function of course. The problem usually lies in implementing these steps, which requires an expert.

tagrun 11y ago

I realized that I missed the part that he's talking about log2(x) rather than log(x+1), which changes the interval. Since MiniMaxApproximation gets quirky when the function is 0 in the given interval (it minimizes the relative error, so a division by 0 occurs, unless you cancel it), I'm rather using the interval [0.25, 0.5] with t=1/4: -3.64631 + 7.93796 x - 5.30406 x^2 This gives a maximum relative error of 0.00334047, which is still an order better than his. As for the (2,2) approximant: (-5.83261 - 16.7012 x + 22.394 x^2)/(1 + 11.2638 x + 7.81124 x^2) max rel error: 1.69763 10^-6 I checked several other values of t, and the max rel error for (2,0) approximant is roughly the same. This leads me to believe that the reason he is getting 10 times larger error is because he is imposing that errors should vanish at certain points. Not only this increases the max rel error, but it also distorts the smoothness of the function in a bad way (the kinks in the plots are a bad sign in general).

gjm11 11y ago
Your polynomial, since it doesn't get it exactly right at 1/4 & 1/2, will have an enormous relative error when you use it for computing logs of numbers very close to 1. Goldberg's initial function b ignores this issue, and a glance at the graph of its absolute errors will confirm that it's the minimax poly for absolute error or very close to it. His next candidate t half-fixes this, and again looking at the graph of the relative error makes it clear that it's right. (Note that avoiding intervals containing zero will indeed allow you to get better max relative errors -- keeping the relative error finite where the function is zero effectively reduces your polynomial degree by 1.) Then finally he fully fixes it with the polynomial g, solving the more interesting problem of minimizing the max relative error for the log function everywhere, and chooses a better functional form by shifting the interval he reduces to with the polynomial s. (The fact that his earlier polynomials do the wrong thing is clearly not a mistake but a paedagogical strategy: start with something simple and incrementally improve it.) My apologies for spelling all this out, but it really appears that you're accusing him of making elementary mistakes without actually understanding what he's doing. "Imposing that errors should vanish at certain points" is not a mistake, it's a necessity if you're implementing log by argument reduction plus local approximation and you care about relative error in the resulting approx_log function, and you are getting better numbers than Goldberg by solving the wrong problem. Your (2,2) rational approximant will be much more expensive to compute than the polynomial and won't yield the speed gains Goldberg is looking for. And, again, it will give you a log function with infinite relative error at 1. [EDITED to fix a typo and clarify something.]