c++ - 快速定点 pow、log、exp 和 sqrt

Question

我有一个定点类 (10.22)，我需要一个 pow、一个 sqrt、一个 exp 和一个 log 函数。

唉，我什至不知道从哪里开始。谁能给我一些有用文章的链接，或者更好的是，给我一些代码？

我假设一旦我有了一个 exp 函数，那么实现 pow 和 sqrt 就变得相对容易了。

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

我发现困难的只是那些 exp 和 log 函数(好像我记得我的一些日志规则，但我记不起关于它们的更多内容)。

据推测，sqrt 和 pow 也会有一种更快的方法，因此即使它只是说使用我上面概述的方法，也将不胜感激这方面的任何指针。

请注意:这必须是跨平台和纯 C/C++ 代码，所以我不能使用任何汇编程序优化。

Answer 1

一个非常简单的解决方案是使用合适的表驱动近似值。如果您正确减少输入，您实际上并不需要大量数据。 exp(a)==exp(a/2)*exp(a/2) ，这意味着你真的只需要计算exp(x)对于 1 < x < 2 .在该范围内，runga-kutta 近似值将给出合理的结果，大约 16 个 IIRC 条目。

同样，sqrt(a) == 2 * sqrt(a/4) == sqrt(4*a) / 2这意味着您只需要 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 可能很有用所以你可以减少log(a) == n + log(a/e^n)通过对该表中的 a 进行简单的硬编码二进制搜索。从 7 步到 3 步的改进不是很大，但这意味着您只需除以一次 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.6931471805599 ，计算前导零的个数，移位a进入用于查找表的范围，并将该移位(整数)乘以定点常数 log2 .可以低至 3 条指令。

使用 e因为减少步骤只是给你一个“好”log(e)=1.0不变，但这是错误的优化。 0.6931471805599 和 1.0 一样好；两者都是 10.22 定点的 32 位常量。使用 2 作为常数来缩小范围允许您使用位移位进行除法。

[编辑 5]而且由于您将其存储在 Q10.22 中，因此您可以更好地存储 log(65536)=11.09035488。 (16 x 对数 (2))。 “x16”意味着我们还有 4 位可用的精度。

您仍然可以从编辑 2 中获得诀窍，log(a/2^n) = log((a/2^n)^8)/8 .基本上，这会给你一个结果(a + b/8 + c/64 + d/512) * 0.6931471805599 - 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 函数。