快速定点战俘，对数，指数和开方 [英] Fast fixed point pow, log, exp and sqrt

问题描述

我有一个固定点类（10.22），我有需要的战俘，一个开方，一个exp和日志的功能。

I've got a fixed point class (10.22) and I have a need of a pow, a sqrt, an exp and a log function.

唉！我不知道在哪里，甚至在此启动。任何人都可以向我提供一些链接到有用的文章或bettter然而，我提供一些code？

Alas I have no idea where to even start on this. Can anyone provide me with some links to useful articles or, bettter yet, provide me with some code?

我assumign，一旦我有一个EXP函数则变得相对容易实现战俘和开方，因为他们刚刚成为

I'm assumign that once I have an exp function then it becomes relatively easy to implement pow and sqrt as they just become.

POW（X，Y）=> EXP（Y *日志（X））
   SQRT（x）=> POW（X，0.5）

pow( x, y ) => exp( y * log( x ) ) sqrt( x ) => pow( x, 0.5 )

它只是那些记录和日志我发现很难功能（好像我记得我的几个日志规则，我不记得关于他们的很多其他）。

Its just those exp and log functions that I'm finding difficult (as though I remember a few of my log rules, I can't remember much else about them).

presumably，顺便说一句，有也将是对开方和POW等等这一方面的任何指针更快的方法是将AP preciated即使它只是说使用的方法我概述以上：）

Presumably, btw, there would also be a faster method for sqrt and pow so any pointers on that front would be appreciated even if its just to say use the methods i outline above :)

请注意：这是跨平台，纯C / C ++ code，所以我不能使用任何汇编的优化

Please note: This HAS to be cross platform and in pure C/C++ code so I cannot use any assembler optimisations.

推荐答案

一个非常简单的解决方法是使用一个体面的表驱动逼近。你实际上并不需要大量的数据，如果你正确地降低您的投入。 EXP（一）== EXP（一/二）* EXP（一/ 2），这意味着你真的只需要计算 EXP（X ）为 1 LT; X - LT; 2 。在该范围，一个runga - 库塔逼近将给与合理的结果〜16 IIRC条目。

A very simple solution is to use a decent table-driven approximation. You don't actually need a lot of data if you reduce your inputs correctly. exp(a)==exp(a/2)*exp(a/2), which means you really only need to calculate exp(x) for 1 < x < 2. Over that range, a runga-kutta approximation would give reasonable results with ~16 entries IIRC.

同样，的sqrt（一）== 2 * SQRT（A / 4）这意味着你需要为 1只下表项; A＆LT; 4 。日志（a）是有点困难：日志（一）== 1 +日志（A / E）。这是一个相当缓慢的迭代，但日志（1024）只有6.9，所以你不会有很多反复。

Similarly, sqrt(a) == 2 * sqrt(a/4) which means you need only table entries for 1 < a < 4. Log(a) is a bit harder: log(a) == 1 + log(a/e). This is a rather slow iteration, but log(1024) is only 6.9 so you won't have many iterations.

您会使用POW类似整首算法： POW（X，Y）== POW（X，地板（Y））* POW（X，压裂（Y） ）。这工作，因为战俘（双，INT）是微不足道的（分而治之）。

You'd use a similar "integer-first" algorithm for pow: pow(x,y)==pow(x, floor(y)) * pow(x, frac(y)). This works because pow(double, int) is trivial (divide and conquer).

对于的组成部分日志（一），它可能是存储表有用1，E，E ^ 2 ，E ^ 3，E ^ 4，E ^ 5，E ^ 6，E ^ 7 这样就可以减少日志（一）== N +日志（A / E ^ N）通过简单的硬codeD二元该表中搜索的。从7到3个步骤的改善是没有那么大，但它意味着你只需要通过电子^ N 而不是 N 倍电子。

[edit] For the integral component of log(a), it may be useful to store a table 1, e, e^2, e^3, e^4, e^5, e^6, e^7 so you can reduce log(a) == n + log(a/e^n) by a simple hardcoded binary search of a in that table. The improvement from 7 to 3 steps isn't so big, but it means you only have to divide once by e^n instead of n times by e.

而对于去年日志（A / E ^ N）项，你可以使用日志（A / E ^ N）=日志（（A / E ^ n）的^ 8）/ 8 - 每一次迭代产生3个位<击>通过查表。 ，让您的code和表大小小。这通常是code为嵌入式系统，而他们没有大的缓存。

[edit 2] And for that last log(a/e^n) term, you can use log(a/e^n) = log((a/e^n)^8)/8 - each iteration produces 3 more bits by table lookup. That keeps your code and table size small. This is typically code for embedded systems, and they don't have large caches.

这STIL不聪明在我身边。 日志（一）=日志（2）+数（A / 2）。你可以只存储定点值 LOG2 = 0.30102999566 ，计数前导零，移 A 的数量成用于您的查找表范围，并乘以定点定 LOG2 这种转变（整数）。可低至3的说明。

[edit 3] That's stil not to smart on my side. log(a) = log(2) + log(a/2). You can just store the fixed-point value log2=0.30102999566, count the number of leading zeroes, shift a into the range used for your lookup table, and multiply that shift (integer) by the fixed-point constant log2. Can be as low as 3 instructions.

使用电子用于还原步骤只是给你一个好日志（E）= 1.0 恒定的，而是那是假的优化。 0.30102999566是一样好一个常数1.0;两者在10.22定点32位常数。使用2作为减少范围的不断让你用一个部门的移位。

Using e for the reduction step just gives you a "nice" log(e)=1.0 constant but that's false optimization. 0.30102999566 is just as good a constant as 1.0; both are 32 bits constants in 10.22 fixed point. Using 2 as the constant for range reduction allows you to use a bit shift for a division.

您仍然编辑2，日志（一/ 2 ^ N）=日志获得的伎俩（（A / 2 ^ n）的^ 8）/ 8 。基本上，这可以让你的结果（A + B / 8 + C / 64 + D / 512）* 0.30102999566 - 与B，C，D的范围在[0， 7。 a.bcd 真的是一个八进制数。不是因为我们用8为动力的惊喜。 （诀窍同样适用与电源2,4或16个）。

You still get the trick from edit 2, log(a/2^n) = log((a/2^n)^8)/8. Basically, this gets you a result (a + b/8 + c/64 + d/512) * 0.30102999566 - with b,c,d in the range [0,7]. a.bcd really is an octal number. Not a surprise since we used 8 as the power. (The trick works equally well with power 2, 4 or 16.)

仍然有一个开口端。 POW（X，压裂（Y）就是 POW（SQRT（X），2 *压裂​​（Y））我们有一个像样的 1 / SQRT（X）。这为我们提供了更为有效的方法。说压裂（Y）= 0.101 二进制，即1/2 1/8加。那意味着 X ^ 0.101 是（X ^ 1/2 * X ^ 1/8），但 X ^ 1 / 就是的sqrt（x）的和 X ^ 1月8日是（开方（开方（SQRT（X）））。保存一个更多的操作，牛顿迭代 NR（X）为我们提供了 1 / SQRT（X）让我们计算 1.0 /（NR（X）* NR（（NR（NR（X）））。我们只反转最终结果，不要直接使用sqrt函数。

[edit 4] Still had an open end. pow(x, frac(y) is just pow(sqrt(x), 2 * frac(y)) and we have a decent 1/sqrt(x). That gives us the far more efficient approach. Say frac(y)=0.101 binary, i.e. 1/2 plus 1/8. Then that means x^0.101 is (x^1/2 * x^1/8). But x^1/2 is just sqrt(x) and x^1/8 is (sqrt(sqrt(sqrt(x))). Saving one more operation, Newton-Raphson NR(x) gives us 1/sqrt(x) so we calculate 1.0/(NR(x)*NR((NR(NR(x))). We only invert the end result, don't use the sqrt function directly.

这篇关于快速定点战俘，对数，指数和开方的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！