计算机如何进行浮点运算？

Posted 2023-03-06

【中文标题】计算机如何进行浮点运算？

【英文标题】：How computer does floating point arithmetic?

【发布时间】：2011-05-17 15:24:40

【问题描述】：

我看过很长的文章解释如何存储浮点数以及这些数字的算术是如何完成的，但是请在我写的时候简要解释一下为什么

cout << 1.0 / 3.0 <<endl;

我看到 0.333333，但是当我写的时候

cout << 1.0 / 3.0 + 1.0 / 3.0 + 1.0 / 3.0 << endl;

我看到 1。

计算机是如何做到这一点的？请解释这个简单的例子。这对我来说已经足够了。

【问题讨论】：

这听起来是个好消息，不是吗？

@close voters：我在这个问题中没有看到任何离题或过于宽泛的内容。如果你这样做，请在 cmets 中解释自己。

我很惊讶是什么让这四位亲近的选民也生气了:)！

顺便说一句，几乎每个答案都错了。大多数答案都假设 FPU 得到了错误的位，但输出转换四舍五入到正确的值。一点也不真实。在十进制和二进制中，有理但重复的分数 1/3 不能精确表示，但是当添加两次时，它确实会四舍五入到完全正确的答案。许多答案出错的地方：这个舍入发生在最后一个操作中，第二个加法，它发生在最低有效位（2**-23 或 2**-53）位。无论输出转换如何完成，它都会产生精确的 1.0 (0x3f800000)。

【参考方案1】：

查看"What every computer scientist should know about floating point arithmetic"上的文章

【讨论】：

我看过这篇文章，里面有一些定理和证明。似乎不是太简短，因此我问了这个问题！ :)

@Narek：浮点运算不简单，如果你打算做任何不平凡的事情，你绝对应该阅读这篇文章。我建议你问一些具体的问题，而不是问这么宽泛的问题。

@Alexandre C，现在我不需要整个理论。我只需要理解这个简单的例子。

@Narek：你需要它来理解这个简单的例子。

好的，那我需要推迟到不确定的时间，直到我有空闲时间！ :(

【参考方案2】：

问题是浮点格式表示以 2 为底的分数。

第一个小数位是 ½，第二个 ¼，它继续为 1 / 2ⁿ。

那个的问题在于，并非每个有理数（可以表示为两个整数之比的数字）实际上都具有这种以 2 为底的格式的有限表示。

（这使得浮点格式难以用于货币值。虽然这些值总是有理数 (n/100)，但实际上只有 .00、.25、.50 和 .75具有以二为底的分数的任意位数的精确表示。 )

无论如何，当您将它们添加回来时，系统最终有机会将结果四舍五入为一个可以准确表示的数字。

在某些时候，它发现自己将 .666... 数字添加到 .333... 数字中，如下所示：

  00111110 1  .o10101010 10101010 10101011
+ 00111111 0  .10101010 10101010 10101011o
------------------------------------------
  00111111 1 (1).0000000 00000000 0000000x  # the x isn't in the final result

最左边的位是符号，接下来的八位是指数，其余位是小数。在指数和分数之间是一个假定的“1”，它始终存在，因此实际上并未存储为归一化的最左边的分数位。我写的零实际上并不像o那样作为单独的位出现。

这里发生了很多事情，在每一步，FPU 都采取了相当英勇的措施来完善结果。已经保留了两个额外的精度数字（超出了结果中的精度），并且在许多情况下，FPU 知道剩余最右边的位中是否有或至少 1 个是 1。如果是这样，则分数的那部分大于 0.5（按比例缩放），因此向上取整。中间四舍五入的值允许 FPU 将最右边的位一直传送到整数部分，最后四舍五入到正确的答案。

这并没有发生，因为有人添加了 0.5； FPU 只是在格式的限制范围内尽力而为。实际上，浮点并不是不准确的。它非常准确，但是我们期望在以 10 为底的有理数世界观中看到的大多数数字都不能用格式的以 2 为底的分数来表示。事实上，很少。

【讨论】：

我想强调的是，1.0 的答案绝对是不是在输出转换例程中最终四舍五入的结果。在这种情况下，所有当前的 FPU 单元实际上都会产生一个精确的 1.0。

【参考方案3】：

让我们算一下。为简洁起见，我们假设您只有四个有效（base-2）数字。

当然，由于gcd(2,3)=1，1/3 在以 base-2 表示时是周期性的。尤其是它不能准确表示，所以我们需要满足于近似值

A := 1×1/4 + 0×1/8 + 1×1/16 + 1*1/32

这比1/3的实际值更接近

A' := 1×1/4 + 0×1/8 + 1×1/16 + 0×1/32

因此，以十进制打印 A 会得到 0.34375（您在示例中看到 0.33333 的事实只是证明了 @ 中有效数字的数量较多987654327@).

将这些相加三次，我们得到

A + A + A
= ( A + A ) + A
= ( (1/4 + 1/16 + 1/32) + (1/4 + 1/16 + 1/32) ) + (1/4 + 1/16 + 1/32)
= (   1/4 + 1/4 + 1/16 + 1/16 + 1/32 + 1/32   ) + (1/4 + 1/16 + 1/32)
= (      1/2    +     1/8         + 1/16      ) + (1/4 + 1/16 + 1/32)
=        1/2 + 1/4 +  1/8 + 1/16  + 1/16 + O(1/32)

O(1/32) 术语不能在结果中表示，所以它被丢弃了，我们得到

A + A + A = 1/2 + 1/4 + 1/8 + 1/16 + 1/16 = 1

QED :)

【讨论】：

但是结果正好是 1.0 只是运气。它们很容易偏离 1 或 2 LSB，而他的程序仍会显示 1，因为默认情况下，浮点转换的精度仅为 6。

是的，这个答案是不正确的，但正确结果的原因不是“运气”。事实上，正如我已经解释的那样，所有三个分数的加法都会产生精确的 1.0，无论您继续打印多少小数位。

@DigitalRoss：想解释一下“这个答案不正确”的地方吗？

@UpAndAdam 再次阅读。核心点是结果的 O(1/32) 部分由于有效 base-2 数字的数量有限，无法用 FP 类型表示，因此产生了精确的 1.0。

【参考方案4】：

至于这个具体的例子：我认为现在的编译器太聪明了，如果可能的话，会自动确保原始类型的const 结果是准确的。我没能骗过 g++ 做一个像这样错误的简单计算。

但是，使用非常量变量很容易绕过这些事情。不过，

int d = 3;
float a = 1./d;
std::cout << d*a;

将准确地产生 1，尽管这不应该是真正的预期。如前所述，原因是 operator<< 将错误舍入。

至于为什么它可以这样做：当您添加相似大小的数字或将 float 乘以 int 时，您几乎可以获得浮点类型可以最大程度地为您提供的所有精度 - 这意味着，比率错误/结果非常小（换句话说，假设您有一个正错误，错误发生在小数点后）。

所以3*(1./3)，尽管作为浮点数，不完全是==1，具有很大的正确偏差，这会阻止operator<< 处理小错误。但是，如果您随后通过仅减去 1 来消除此偏差，则浮点数将滑落到错误处，并且突然之间它不再是可以忽略的了。正如我所说，如果您只输入3*(1./3)-1，则不会发生这种情况，因为编译器太聪明了，但请尝试

int d = 3;
float a = 1./d;
std::cout << d*a << " - 1 = " <<  d*a - 1 << " ???\n";

我得到的（g++，32 位 Linux）是

1 - 1 = 2.98023e-08 ???

【讨论】：

我在 Linux 64 位和 g++ 中得到 0 :)。但我真的很感兴趣如何找到错误？？？

这很有趣。我没有期望它在 32 位和 64 位上的工作方式不同，但我不确定...

你试过用比 3 更大的素数吗？嵌套分数怎么样？

您正在使用不同版本的 g++ 和/或不同的优化级别。尝试将其设为extern int d; extern float a，并仅使用链接到可执行文件的int d = 3; float a = 1./d; 编译一个单独的源文件。这应该会挫败其所做的任何精度优化优化。

【参考方案5】：

这是可行的，因为默认精度为 6 位，四舍五入为 6 位结果为 1。请参阅 C++ draft standard (n3092) 中的 27.5.4.1 basic_ios 构造函数。