为什么浮点数的整数表示提供对数的分段线性近似？
Ada*_*and 21 math floating-point logarithm approximation
如果您正在阅读有关 20 世纪 90 年代图形学发展的新闻，您可能会关注 Jim Blinn 在IEEE 计算机图形与应用中的专栏“Jim Blinn 的角落”。1997 年夏天，您可能会在名为“浮点技巧”的专栏中打开您的问题。所有的技巧都基于 Blinn 从微软的两位老手 Steve Gabriel 和 Gideon Yuval 那里传授的以下知识：
\n“如果只处理正数，则将浮点数的位模式解释为整数，给出对数函数的分段线性近似”
\n我的问题是，浮点数是什么导致了这个特性？
\n示范
\n\n// These will result in undefined behavior\n// and are not the \'right\' way to do this\n// but I am following Blinn 1997 here\nint AsInteger(float f) {\n return * ( int * ) &f;\n}\n\nfloat AsFloat(int i) {\n return * ( float * ) &i;\n}\n\n// FISR demonstrating the logaritmic nature of the floating-point numbers\n// The constant chosen is not "magic" but that which replaces the lost\n// exponents from the shift of the number as an integer.\nfloat BlinnISR(float x) {\n int i;\n float y;\n // 0x5F400000 = 1598029824\n // which is (AsInteger(1.0f) + (AsInteger(1.0f) >> 1))\n i = 0x5F400000 - ( AsInteger(x) >> 1);\n y = AsFloat(i);\n // Blinn uses 1.47 and 0.47 not 1.5 and 0.5 \n // which are what you would get from differentiating\n y = y * (1.47f - 0.47f * x * y * y);\n return y;\n}\n上面是Blinn 示例代码的重新创建，展示了它如何实际用于平方根倒数。
\n请记住 1 / sqrt(x) 与 x (-1/2)相同。因此，如果我们可以访问对数表，求平方根的一种方法是取 x (-1/2)的自然对数，得到 -1/2 * ln(x)。如果对其求幂，则得到 1 / sqrt(x) = e ((-1/2) * ln(x))。
\n右移浮点位的整数表示形式为我们提供了除以输入浮点数的对数的二的近似值（这是我们的 1/2 * ln(x)）。这是从恢复常数中减去的，这里只是0x5F400000，还没有“魔法”，所以我们基本上有 -1/2 * ln(x)。将其返回到浮点数作为我们的求幂步骤。然后我们就有了一个数字，它可以作为牛顿法的良好第一个输入。
它看起来像什么
\n从视觉上看，您可以看到这里的关系（Blinn 1997，第 81 页）：
\n
Mitchell 1961（第 513 页，链接如下）中出现了类似的理想化图表：
\n
Jean-Michel Muller，\xe2\x80\x9d基本函数和近似计算，\xe2\x80\x9d中显示了同样的关系，IEEE 会议录，卷。108，没有。12，p。2020 年 12 月 2137
\n
适用的情况
\n这个特性并不是 IEEE 754 浮点数所独有的，尽管这是 Quake III Arena 的快速平方根反比中最著名的利用，它使用了一个标准 float。然而，它可以像 Knuth 的MIX教学格式这样的抽象格式工作，或者与其他格式一起工作，包括下面的例子。
据我所知，这种关系首先在学术媒体上被记录，约翰·米切尔 (John Mitchell) 于 1961 年发表的《使用二进制对数的计算机乘法和除法》用抽象的二进制对数系统演示了该方法。William Kahan 于 1961 年或 1962 年左右在 IBM 7090 上实现了类似的方案。我们也知道，部分由于 Noah Hellman 的硕士论文《基于Mitchell 的浮点数近似运算》，它也在 UNIX 中实现从 1974 年开始一直作为系统平方根，直到 20 世纪 80 年代中期为 BSD 4.3 创建 libm。大约在同一时间，Kahan 和他的研究生 KC Ng 使用这种方法制作了一个“神奇”平方根，但没有发表，该方法旨在对 not 进行操作，floats但doubles会减半并作为 32 位部分进行操作（请参阅相关的 stackoverflow）问题）。
但是，它并不适用于所有浮点格式。DEC 的VAX 浮点格式交错符号、指数和分数部分（将其分成两部分），因此它不会这样做。
\n我不问什么
\n我不关心（对于这个问题）工程师如何了解这种关系，或者它在哪里被利用（参见 FISR）。我想知道它为何存在。我问过哪些与该主题相关的历史问题可以在逆计算堆栈交换中找到：早期在 C 中索引浮点位的示例，类型双关语和 C 标准/编译器
\n我也不是问如何做到这一点，尽管一个好的答案可能包括带有一些代码的说明性示例。然而，您应该首先看看codegolf上这个令人惊奇的 FISR 问题，看看一些真正有趣的问题。
\n一个很好的答案
\n一个糟糕的答案可能是“任何拥有计算尺的工作工程师都会理解这种关系”，尽管我可以想象一个好的答案可能会诉诸计算尺。另一个糟糕的答案可能是“因为 Knuth 是这么说的”，而一个非常好的答案可能会在解释中包含这种情绪，这取决于 TAOCP 第 2 卷第 4.2.4 节（我没有页码，我有 O\ '赖利的分页令人厌恶，它是印刷版第一版的第 240-247 页）。
\n一个好的答案应该清楚地理解这种关系的性质及其产生的原因。尽管我可以预见到一些挑战，但证明这种关系一定如此的证据将会很有趣。
\n我为什么要问
\n我正在写快速平方根反比的历史，相关信息收集在0x5f37642f.com上。作为历史学家调查这一问题的一部分意味着向社区提出这样的问题。我没有足够的专业知识来独自提炼其中的大部分内容，因为我不是数学家或计算机科学家。
\nFISR 基本上有三个有趣的部分：对数近似、“神奇”常数”和牛顿法。在这三者之间，当人们写 FISR 时，他们会写神奇常数，然后是牛顿法，然后花很少的时间谈论近似值。我从更广泛的“文献”中得到的感觉是，近似值已广为人知，但还没有被清楚地理解。我从答案中学到的知识将帮助我理解它。
\n出色的答案将在本文或演讲的致谢中注明。
\n我也问这个问题，因为拜托，还没有人问过这个问题，而且我们已经在这里呆了多久了？
\nstackoverflow 上的相关问题和解答
\n\nSte*_*mit 19
我将限制我对这个问题的回答，“为什么浮点数的\n整数表示提供\n对数的分段线性近似？”。\n快速平方根反比是一个令人着迷的话题，但我对它了解不够，无法谈论它。
\n但其实回答标题中的问题很容易。首先我将从“TL;DR”摘要开始：
\n\n\n这种近似之所以有效，是因为当一个数字是底数的精确幂时，当以指数表示法表示时，它的对数（简单地）就是它的指数。像 IEEE-754 这样的浮点表示法是以2 为基数的指数表示法，指数几乎占据了格式中最高精度的位，只有符号位除外，符号位对于正数为 0，因此可以忽略。所有较低精度位（即指数下方）形成完美的线性级数，理论上从 000 到 999，部分原因是 IEEE-754 格式使用带有“隐藏”前导 1 位的归一化有效数。
\n
较长的解释从回顾开始：对数是什么意思？\n嗯，它是指数或幂函数的反函数。\n指数表示诸如 10 的二次方\n(10 2 ) 是100. 因此，（以 10 为底）对数函数\n会提出这样的问题：\n我们需要对 10 进行多少次幂才能得到\n数字 100？答案当然是2。
\n您还可以问这样的问题：“我们需要加注 10 到多少才能得到\n数字 50？” 这显然更难，因为 50 不是 10 的偶数次方。这背后有一些非常酷的数学原理，我不会尝试解释，但将其提高到小数次方是有定义的这样 10 1.69897就是 50。\n因此 50 的对数（以 10 为底）是 1.69897。
\n然后你也不必将自己限制在 10 的幂。\n数学家喜欢使用特殊数字 e \ n(2.71828) 的幂，当然计算机程序员也喜欢（只是喜欢）\n使用2. 我们需要将 2 进行多少次方才能得到数字 32？每个计算机极客都知道它是 5。因此 32 的对数（以 2 为底）是 5。
\n然后是科学记数法。我们以数字 12345 为例。\n用科学计数法表示该数字的一种方法是 1.2345\n\xc3\x97 10 4。事实证明，这给了我们一个关于 12345 以 10 为底的对数是什么的提示。它\n将大于或等于 4，且小于 5。（事实上，它约为 4.09149。）\n只要我们选择以某种方式标准化的科学表示，\n这就有效，有效数字部分（通常称为“尾数”）选择在 1 到 10 之间。（也就是说，如果我们查看的话，我们将无法以这种方式得到对数的良好估计在其他等效\n符号 0.12345 \xc3\x97 10 5或 12.345 \xc3\x97 10 3。）
\n现在让我们转向计算机浮点表示。\n这些通常使用 \xe2\x80\x94 ，而 IEEE-754 肯定使用 \xe2\x80\x94 base-2\n（二进制）指数表示，即具有归一化\n有效数值乘以 2 的某个幂的值。例如，\n数字 1.125 表示为 1.125 \xc3\x97 2 0，\n数字 25.25 表示为 1.578125 \xc3\x97 2 4，\n数字 0.1875表示为 1.5 \xc3\x97 2 -3。
\n而且，就像以 10 为底一样，这些指数再次为我们提供了以 2 为底的对数 \xe2\x80\x94 的粗略的、仅整数的近似值，\只要我们使用标准化有效数，在本例中为 1 到 2 之间.\n因此 log 2 (25.25) 将介于 4 和 5 之间。
\n现在也是我预示\n我们将自己限制在标准化有效数字值这一事实所暗示的超级特殊秘密事实的时候了。\n我说过它们将“在 1 和 2 之间”，但是\n更准确地说，它们的范围是从 1.00000 到 1.99999。\n但请注意，无论有效数值是什么，\n对于基数 2，第一个数字始终是 1。\n稍后我们将看看我们能用它做什么。
\nIEEE-754 二进制浮点表示均由符号位、后跟一定数量的指数位、再后跟一定数量的有效数位组成。例如，对于我们的数字 25.25，我们将使用一个 0 位来表示正号，然后是指数 4 的某种表示，然后是有效数字 1.578125 的某种表示。
\n但现在我们已经了解到浮点数的原始位解释可能与对数有什么关系。对于正数，我们不必担心符号位（因为它是 0），接下来我们看到的是指数，而且我们已经看到指数是对数的粗略（仅整数）近似值。因此，剩下要问的是，我们如何\n在具有精确\n以 2 为底的对数\n(log 2 (8) = 3, log 2 (16)的二次方值之间进行插值= 4，\nlog 2 (32) = 5，等等）\n以及中间的小数？
\n如果我给你均匀间隔的数字，例如 100、200、300，并要求\n你填写之间的数字，那么你就不会遇到什么麻烦：\n你可以用 01、02、03 替换 00 部分， ...，最多 99。\n如果您了解有关二进制数的任何信息，\n我要求您对均匀间隔的二进制数执行相同的操作，例如
\n00100000\n01000000\n01100000\n10000000\n您可以做或多或少相同的事情，00000用连续的二进制数00001, 00010,00011等替换 \n 部分。
值得注意的是，当我们用 IEEE-754 浮点数对有效数字进行编码时，我们将看到与此完全相同的情况。
\n让我们开始看看实际的 IEEE-754\n浮点格式。我们将使用单精度 float32，以使数字保持合理。此格式使用 1 位表示符号 (S)，使用 8 位表示指数 (e)，使用 23 位表示有效数 (s)：
\nS e e e e e e e e s s s s s s s s s s s s s s s s s s s s s s s\n让我们开始弄清楚数字 16.000 = 2 4和 32.000 = 2 5的表示形式会是什么样子。
\n16 是 1.0 \xc3\x97 2 4。所以尾数是 1.0，指数是 4。这些将如何存储？我将首先关注有效数。还记得观察到的第一个数字（小数点左边的那个）总是 1 吗？IEEE-754 采取的立场是，这个始终为 1 的前导位根本不需要存储。为 1.0 有效数存储的\n实际位将为 00000。\n指数将是 4 的某种表示形式。因此，我\n将\n表示数字 16.000，如下所示：
\n0 4 4 4 4 4 4 4 4 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n这些“4”是近似值；我还没来得及解释指数是如何编码的。
\n数字 32.000 实际上类似于 \xe2\x80\x94，除了指数之外，它是相同的：
\n0 5 5 5 5 5 5 5 5 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n但我要构建的关键点是 \n\n\n\n\xe2\x80\x94 之间那些大于 16.0 但小于 32.0 的数字，都将是\n的变体
\n0 4 4 4 4 4 4 4 4 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n除了用所有连续的 23 位二进制数字 \xe2\x80\x94 替换 0\xe2\x80\x94 中的所有 2^23 或 8388608 \xe2\x80\x94 从 000\xe2\x80\xa6001 到 111\ xe2\x80\xa6111。
\n现在我们已经准备好了所有的部分，看看“浮点数的整数表示如何提供对数的分段线性近似”。我们已经证明，16 到 32 之间的数字的二进制表示形式可以从以 4 等开头的数字到以 5 等开头的数字形成一个良好的线性斜坡。
\n我一直在回避指数如何编码的问题。\n这是使用偏差方案完成的。对于binary32，存储的位模式\n是实际指数值加上127。因此我们的指数值4\n将存储为4+127=131的位模式，而5\n将存储为132。因此实际的指数值\n将存储为4+127=131的位模式。 16.0 和\n32.0 的编码将是
\n0 1 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n\n0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n其中我填写了编码指数 131 和 132 的二进制值，即\n10000011和10000100。现在我们可以看到实际的原始值
01000001100000000000000000000000和\n 01000010000000000000000000000000、或0x41800000和\n 0x42000000、或 1098907648 和 1107296256。
这些数字与 16 和 32（即 4 和 5）的以 2 为底的对数没有任何明显的关系，但它们实际上非常接近。它们被偏差的影响所抵消，并且它们的大小完全错误。但这并不奇怪：\n通过采用 32 位量并将它们视为整数，我们\n得到了很大的数字（它们确实是整数），但这些对数的实际\n值应该是 4.0 和 5.0，我们想要\n小数部分的可能性。因此，如果我们要将这些大整数视为一些较小数字的对数的某种再现，那么它将涉及一些非常类似于定点表示的东西。
\n推导该定点表示的细节非常简单：缩放因子显然是 2 23 = 8388608，并且还有一个偏移量，对应于 binary32 格式中内置的偏差。如果在缩放之前应用，该偏移量将为 1065353216\n(127 \xc3\x97 2 23 )\n，或者在缩放之后应用 127。
\n让我们看看这是否有效。我们得出了对应于 16.0 和 32.0 的整数值\n1098907648 和 1107296256。所以：
\n(1098907648 - 1065353216) / 8388608 = 4\n(1107296256 - 1065353216) / 8388608 = 5\n或者其他方式：
\n1098907648 / 8388608 - 127 = 4\n1107296256 / 8388608 - 127 = 5\n所以它起作用了！
\n据说一张图片胜过 1,000 个单词，而且早就该用这种方式进行解释了。这是一张图表，\xe2\x80\x94 并非巧合，\xe2\x80\x94 看起来非常像您在 1961 年 Mitchell 发布的图表：
\n
红色曲线是实际的以 2 为底的对数函数。黑点是 log 2 (X) 为整数的地方。在这些点上，binary32 浮点值的原始位（解释为整数并适当缩放）恰好等于 log 2。然后在蓝色中，我们有直接连接这些点的直线段。
\n回顾一下：如果我们采用与 2 的精确幂相对应的二进制 32 浮点值的原始位，我们会得到：
\n| X | 二进制32（十六进制） | 二进制 32（十进制） | 记录2 |
|---|---|---|---|
| 0.25 | 3e800000 | 1048576000 | -2 |
| 0.5 | 3f000000 | 1056964608 | -1 |
| 1 | 3f800000 | 1065353216 | 0 |
| 2 | 40000000 | 1073741824 | 1 |
| 4 | 40800000 | 1082130432 | 2 |
| 8 | 41000000 | 1090519040 | 3 |
| 16 | 41800000 | 1098907648 | 4 |
| 32 | 42000000 | 1107296256 | 5 |
2 的精确幂的有效数为 1.0（也就是说，这些数字都是 1.0 \xc3\x97 2 N），这意味着它们在 IEEE-754 中被编码为全 0（由于“隐藏 1 位”） " 属性）。\n并且由于它们是正数，因此它们的符号位也为 0，唯一不为 0 的部分是指数。由于对于精确的幂，指数是对数，因此对于这些数字，在移位和缩放之后，原始位模式的整数解释始终是对数。（我认为“移位和缩放”的某些方面就是神奇数字 0x5f37642f 的用武之地。）
\n然后，对于所有不是 2 的精确幂且对数不是整数的数字，有效数字值的布置方式意味着原始值（如中的蓝线所示）上图）在精确幂之间形成完美的线性插值。因此，在（几乎）所有情况下，原始位模式确实最终成为“对数的分段线性近似”。（我说“几乎”是因为这个结果对于次正规以及特殊值 Inf 和 NaN 不成立。）
\n对“为什么”问题的另一种解释是“为什么 IEEE-754 的设计者要这样设置？” 我不认为这是故意的，也就是说，我不认为原始位模式（解释为整数）对应于以 2 为底的对数的分段线性近似。然而，考虑到这三个相关事实：
\n- \n
- 指数符号的定义基本上是对数的分段近似 \n
- IEEE-754 确实有一个目标，即保持正浮点值按其原始位值单调递增 \n
- IEEE-754隐藏了有效数的高位“1”位；所有有效数均采用 1.xxx 形式，但仅 xxx 部分如此存储 \n
“分段线性近似”的结果相当整齐。
\n- 这是一个很好的答案，您为此付出的努力是显而易见的。我已经开始悬赏，但我希望你从文字中看到你是击败的答案。我希望本页能够成为我 2009 年遇到 FISR 时所希望的参考资料。 (3认同)
nju*_*ffa 11
问题中引用的文字实际上是不正确的。更准确的版本将是：正标准 IEEE-754 浮点数的位模式，解释为定点数并针对指数偏差进行调整，给出了二进制对数函数的分段线性近似。
\n为了使该技术发挥作用，浮点格式必须是二进制的，并且有效数的最高有效位必须是隐式的并且不实际存储。此外，它必须将有效数存储在最低有效位的连续块中，并将指数存储在直接相邻的较高有效位的连续块中。IEEE-754 之前的一些浮点格式满足这些要求，但许多却不能。
\n该技术之所以有效，是因为 IEEE-754 浮点格式使用半对数表示，其中指数提供对数标度，但有效数提供对数标度之间的线性标度细分。这与计算尺不同，计算尺是完全对数的，即主要分区和细分都遵循对数刻度。
\n如果我们有一个正态IEEE-754binary32数y，其存储的指数E是真指数e加上指数偏差 127。有效数s位于 [1,2) 中，仅f存储 [0, 1) 中的小数部分，即s = 1 + f，因为有效数的最高有效位是隐式的。该数字的值为y2 e * (1 + f)。其二进制对数为 log 2 (2 e * (1 + f)) = log 2 (2 e ) + log 2 (1 + f)。给定 的域f，一个非常简单粗略的近似是 log 2 (1 + f) = f （请注意，这与端点处的数学结果匹配），因此 log 2 (2 e * (1 + f)) \xe2 \x89\x88 e + f。
换句话说：减去指数偏差后，该数字表示 的二进制对数的粗略近似值的s8.23 定点y表示。要将其转换为浮点数，请除以 2 23或乘以 2 -23。为了使这个过程在大多数现代处理器上更加高效，可以在过程的最后而不是一开始就减去指数偏差。这导致以下示例性 ISO-C99 实现，最大绝对误差为 8.6e-2：y
#include <stdint.h>\n#include <string.h>\n\nuint32_t float_as_uint32 (float x) \n{ \n uint32_t r; \n memcpy (&r, &x, sizeof r); \n return r; \n} \n\nfloat fastest_log2 (float y)\n{\n const int FP32_EXPO_BIAS = 127;\n const float two_to_m23 = 1.192092896e-7f; // 0x1.0p-23\n float x = (float)(int) float_as_uint32 (y);\n return x * two_to_m23 - (float) FP32_EXPO_BIAS;\n}\n在支持融合乘加运算并具有适当编译器设置的硬件上，最后一行fastest_log2()应映射到某种 FMA 风格。
这个答案讨论了二进制浮点格式的正的、有限的正常数。（该前提不适用于零、负数、无穷大、次正规数和 NaN。）
\n对于固定基数b、整数指数e和尾数s ，浮点数y表示为 \xc2\xb1 b e s， 0 \xe2\x89\xa4 s < b。对于正常形式， 1 \xe2\x89\xa4 s < b。对于二进制格式，b是 2。对于正数，符号当然是 +。对于任何特定的浮点格式，e都有界限， s的大小和精度（ s中以b为基数的位数）也有界限。对于这个答案的一般讨论，我们并不关心具体的特征，尽管它们可能会在某些时候被介绍。
\n给定y = 2 e s，log 2 y = log 2 (2 e s ) = e + log 2 s。
\n在 IEEE 754 \xe2\x80\x9c 单精度\xe2\x80\x9d 格式中，binary32，通常用于float，普通数 2 e s被编码为 0 位作为符号，8 位为e +的二进制127 和 23 位是 ( s \xe2\x88\x921)2 23的二进制。将这 32 位重新解释为二进制数字（按照上面列出的顺序从最高有效位到最低有效位）生成值 ( e +127)2 23 +( s \xe2\x88\x921)2 23 = ( e +126+ s ) 2 23。
\n\n我的问题是，浮点数是什么导致了这个特性？
\n
正如我们在上面所看到的，重新解释浮点数的编码，特别是 IEEE-754 二进制32，并没有给出对数的近似值，因为它完全超出了比例范围并且存在偏移量。
\n然而，如果我们考虑调整 f( x ) = x/2 23 \xe2\x88\x92127，则 f( y重新解释为二进制) = f(( e +126+ s )2 23 ) = ( e +126+ s )2 23 /2 23 \xe2\x88\x92127 = e + s \xe2\x88\x921。现在我们看到y的指数e；该值e + s \xe2\x88\x921 至少接近 log 2 y = e + log 2 s。因此，它是 log 2 y的近似值。（另请注意，在s = 1 和 2 时s \xe2\x88\x921 分别为 0 和 1 ，与 log 2 s相同。因此它可以作为有效数覆盖的区间内log 2 s的近似值。）
\n如果我们将e的每个值作为分段线性近似的一个部分，则近似e + s \xe2\x88\x921 是线性的，因为它符合线性模型e + s \xe2\x88\x921 = ay + b，其中a = 2 \xe2\x88\x92 e和b = e \xe2\x88\x921，正如我们所见，因为ay + b = 2 \xe2\x88\x92 e y + e \xe2\x88\x92 1 = 2 \ xe2\x88\x92 e (2 e s ) + e \xe2\x88\x92 1 = e + s \xe2\x88\x921。
\n关于调整函数 f( x ) = x/2 23 \xe2\x88\x92127，您经常会在操作浮点数编码的代码中看到这一点\xe2\x80\x94，它会以各种方式添加或减去以考虑指数和有效数的前 1 位的偏差，并且它将以各种方式移动以将指数字段移动到所需的位置。
\n