快速定点pow，log，exp和sqrt
我有一个固定点类（10.22），我需要一个pow，一个sqrt，一个exp和一个日志函数。
唉，我不知道从哪里开始。 任何人都可以提供一些有用的文章的链接，或者，还没有，给我提供一些代码？
我认为，一旦我有了exp函数，那么实现pow和sqrt变得相对容易。
pow（x，y）=> exp（y * log（x））sqrt（x）=> pow（x，0.5）
它只是那些我发现很难的exp和日志函数（好像我记得我的一些日志规则，我记不起其他的很多了）。
据推测，顺便说一句，对于sqrt和pow也会有一个更快的方法，所以即使只是说使用我在上面概述的方法:)前面的任何指针都会受到赞赏:)
请注意：这是跨平台和纯C / C ++代码，所以我不能使用任何汇编程序优化。
一个非常简单的解决方案是使用适当的表驱动近似。 如果正确减少输入，实际上并不需要大量数据。 exp(a)==exp(a/2)*exp(a/2) ，这意味着你真的只需要为1 < x < 2计算exp(x) 。 在该范围内，runga-kutta近似将得到合理的结果，约16个条目IIRC。
类似地， sqrt(a) == 2 * sqrt(a/4) ，这意味着您只需要1 < a < 4表条目。 Log（a）有点难： log(a) == 1 + log(a/e) 。 这是一个相当慢的迭代，但log（1024）只有6.9，所以你不会有很多迭代。
你会对pow使用类似的“整数优先”算法： pow(x,y)==pow(x, floor(y)) * pow(x, frac(y)) 。 这是有效的，因为pow(double, int)是微不足道的（分而治之）。
[编辑]对于log(a)的积分组件，存储表1, e, e^2, e^3, e^4, e^5, e^6, e^7可能是有用的，所以你可以通过该表中a的简单硬编码二进制搜索来减少log(a) == n + log(a/e^n) 。 从7步到3步的改进并不是那么大，但这意味着你只需要用e^n除以e^n而不是n除以e 。
[编辑2]对于最后一个log(a/e^n)项，你可以使用log(a/e^n) = log((a/e^n)^8)/8 - 每次迭代产生3个以上通过表查找的位。 这使您的代码和表大小变小。 这通常是嵌入式系统的代码，并且它们没有大型缓存。
[编辑3]这不是我的聪明。 log(a) = log(2) + log(a/2) 。 您可以只存储定点值log2=0.30102999566 ，计算前导零的数量，将a移动到用于查找表的范围，并将该移位（整数）乘以定点常量log2 。 可以低至3条指令。
使用e进行还原步骤只会给你一个“漂亮”的log(e)=1.0常量但这是错误的优化。 0.30102999566和1.0一样好; 两者都是10.22固定点的32位常数。 使用2作为范围缩小的常量，可以使用位移来进行除法。
你仍然可以从编辑2获得技巧， log(a/2^n) = log((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号电源。）
[编辑4]仍然有一个开放的结局。 pow(x, frac(y)就是pow(sqrt(x), 2 * frac(y)) ，我们有一个不错的1/sqrt(x) 。这给了我们更有效的方法。说frac(y)=0.101二进制，即1/2加1/8。那意味着x^0.101是(x^1/2 * x^1/8) 。但x^1/2只是sqrt(x)和x^1/8是(sqrt(sqrt(sqrt(x))) 。再保存一次操作，Newton-Raphson NR(x)给出1/sqrt(x)所以我们计算1.0/(NR(x)*NR((NR(NR(x))) 。我们只反转最终结果，不要直接使用sqrt函数。
下面是Clay S. Turner的定点对数库2算法的示例C实现[ 1 ]。 该算法不需要任何类型的查找表。 这对于内存限制紧张且处理器缺少FPU的系统非常有用，例如许多微控制器就是这种情况。 然后还使用对数属性支持日志库e和日志库10，对于任何基数n ：
log (x) y log (x) = _______ n log (n) y
其中，对于此算法， y等于2。
这个实现的一个很好的特性是它支持变量精度：精度可以在运行时确定，代价是范围。 我实现它的方式，处理器（或编译器）必须能够进行64位数学运算以保持一些中间结果。 它可以很容易地适应不需要64位支持，但范围将会减少。
使用这些函数时， x是根据指定precision缩放的定点值。 例如，如果precision是16，那么x应该缩放2 ^ 16（65536）。 结果是一个定点值，其比例因子与输入相同。 INT32_MIN的返回值表示负无穷大。 返回值INT32_MAX表示错误， errno将设置为EINVAL ，表示输入精度无效。
#include #include #include "log2fix.h" #define INV_LOG2_E_Q1DOT31 UINT64_C(0x58b90bfc) // Inverse log base 2 of e #define INV_LOG2_10_Q1DOT31 UINT64_C(0x268826a1) // Inverse log base 2 of 10 int32_t log2fix (uint32_t x, size_t precision) { int32_t b = 1U << (precision - 1); int32_t y = 0; if (precision < 1 || precision > 31) { errno = EINVAL; return INT32_MAX; // indicates an error } if (x == 0) { return INT32_MIN; // represents negative infinity } while (x < 1U << precision) { x <<= 1; y -= 1U << precision; } while (x >= 2U << precision) { x >>= 1; y += 1U << precision; } uint64_t z = x; for (size_t i = 0; i < precision; i++) { z = z * z >> precision; if (z >= 2U << precision) { z >>= 1; y += b; } b >>= 1; } return y; } int32_t logfix (uint32_t x, size_t precision) { uint64_t t; t = log2fix(x, precision) * INV_LOG2_E_Q1DOT31; return t >> 31; } int32_t log10fix (uint32_t x, size_t precision) { uint64_t t; t = log2fix(x, precision) * INV_LOG2_10_Q1DOT31; return t >> 31; }
此实现的代码也存在于Github ，以及一个示例/测试程序，该程序说明如何使用此函数计算和显示从标准输入读取的数字的对数。
[ 1 ] CS Turner， “快速二进制对数算法” ， IEEE信号处理Mag。 ，pp.124,140，2010年9月。
一个很好的起点是Jack Crenshaw的书“实时编程的数学工具包” 。 它对各种超越函数的算法和实现进行了很好的讨论。
仅使用整数运算检查我的定点sqrt实现。 发明很有趣。 现在相当老。
否则，请检查CORDIC算法集。 这是实现列出的所有function和三角函数的方法。
编辑：我在这里发布了GitHub上的评论来源