使用固定点而非浮点数实现快速反平方根

Question

我正在尝试实现快速反平方根算法，但对于固定点数却一无所获。

我的实现原理与文章中完全相同，只是将数字以浮点格式 x = (-1) ^ s * (1 + M) * 2 ^ (E-127) 改为以32位固定点数格式 x = M * 2 ^ -16，具有16个小数位和16个分数位。

问题在于，我找不到“魔数”的值。根据我的计算，它不存在，但我不是数学家，我认为我做错了什么。

为了解决Y = 1 / sqrt (x)，我使用了以下推理（我不知道是否正确）。

在原始代码中，我们得出了用于牛顿近似的Y0：

i = 0x5f3759df - (i >> 1);

这意味着我们将得到一个浮点数，其值为：

y0 = (1 + R2 - M / 2) * 2 ^ (R1 - E / 2);

这是因为操作符>>将指数和尾数除以2，然后我们将这两个数作为整数进行减法运算。

按照文章中所展示的步骤，我将x的格式设置为：

x = M * 2 ^ -16

为了执行相同的逻辑，我尝试为以下内容定义Y0：

Y0 = (R2 - M / 2) * 2 ^ (R1 - (-16/2));

我正在尝试找到一个能够最小化下列误差的数字：

 

error = (Y - Y0) / Y

无论 R1 的值是多少，我都可以进行移位操作来纠正最终结果的指数值，从而在定点处获得正确的结果。

我错在哪里？

- vicente cesar

---

也许这个问题应该在数学交换平台上提出，如果问题是关于数学的。 - Christian Gibbons

当然不行。"Fast Inverse Square Root"算法使用IEEE-754二进制32位浮点表示的具体特征，以及浮点运算的特性，将该算法适应于固定点算术是一个软件工程问题。 - Eric Postpischil

Answer 1

无法做到。

快速反平方根是由于浮点表示，已经将数字分为2的幂（指数）和显著数字。

可以做到。

通过与浮点数相同的技巧，可以将固定点转换为2^exp * x。给定uint32_t a，uint8_t exp = bias- builtin_count_leading_zeros(a)；uint32_t b = a << exp，常数（以及a的域）要非常小心地选择，这样就不会有下溢或上溢。

因此，您实际上将拥有自定义浮点表示，该表示适用于此特定目的，至少省略符号位，并具有指数的最佳可能位数，可能是8位。

- Aki Suihkonen

---

我看到过使用相同原理的技术，但是将来我会在FPGA中实现这个算法，“builtin_count_leading_zeros”可能不是一个好选择。这将涉及到“优先编码器”，据我所知，它会限制电路上可用的时钟速度。这个算法是错误的选择吗？ - vicente cesar

@vincente cesar 最好在问题中提到这一点。我们不知道您的FPGA的具体情况。您可以尝试一种逐位算法，恢复或非恢复变体，并可能对其进行流水线处理。如果您有一个计数前导零电路、乘法器和一些ROM空间，您可以尝试基于牛顿-拉弗森的算法，就像我在这个答案中展示的那样。 - njuffa

在约翰·卡马克 （John Carmack）的“快速求平方根倒数”方法中，使用精心选择的常量来讨论。其美妙之处在于该解决方案通常在3次迭代内收敛，在所选择的准确度范围内非常高效。 - David C. Rankin

在FPGA中，我可能会通过每次向左移动16、8、4、2、1位来归一化值，在每次迭代中得到一个指数的位。这可能不会限制时钟速度，但会增加延迟。 - Aki Suihkonen

@AkiSuihkonen 我正在尝试学习链接中的算法。是否可能将其调整为计算数字的倒数？我仍在努力理解。我正在阅读尽可能多的解决方案，它的算法看起来是计算倒数和平方根倒数的最佳选择。我正在使用FPGA Spartan 6 xc6slx16。 - vicente cesar

是的，一个2^n *(1+m)的数的倒数大约是2^-n * (1-m)。 - Aki Suihkonen