floating-point - IEEE-754 浮点精度 : How much error is allowed?

Question

我正在移植 sqrt函数(用于 64 位 double )来自 fdlibm到我目前使用的模型检查器工具( cbmc )。
作为我工作的一部分，我阅读了很多关于 ieee-754 标准的内容，但我认为我不了解基本操作(包括 sqrt)的精度保证。

测试我的 fdlibm 的 sqrt 端口，我在 64 位 double 上使用 sqrt 进行了以下计算:

sqrt(1977061516825203605555216616167125005658976571589721139027150498657494589171970335387417823661417383745964289845929120708819092392090053015474001800648403714048.0) = 44464159913633855548904943164666890000299422761159637702558734139742800916250624.0

(这个案例在我的关于精度的测试中打破了一个简单的后置条件；我不确定这个后置条件是否适用于 IEEE-754)

为了进行比较，几个多精度工具计算如下:

sqrt(1977061516825203605555216616167125005658976571589721139027150498657494589171970335387417823661417383745964289845929120708819092392090053015474001800648403714048.0) =44464159913633852501611468455197640079591886932526256694498106717014555047373210.truncated

可以看到，从左边数第 17 个数字是不同的，意思是这样的错误:

3047293474709469249920707535828633381008060627422728245868877413.0

问题 1:允许如此巨大的误差吗？

标准是说每个基本操作 (+,-,*,/,sqrt) 应该在 0.5 ulps 之内，这意味着它应该等于一个数学上精确的结果四舍五入到最接近的 fp 表示(维基说一些库仅保证 1 ulp，但目前这并不重要)。

问题 2:这是否意味着每个基本操作都应该有一个错误 < 2.220446e-16 和 64 位 double (machine-epsilon)？

我确实用 x86-32 linux 系统 (glibc/eglibc) 计算了相同的结果，并得到了与 fdlibm 相同的结果，这让我想到:

a:我做错了什么(但如何:printf 将成为候选人，但我不知道这是否可能是原因)

b:错误/精度在这些库中很常见

Answer 1

IEEE-754 标准要求正确舍入所谓的“基本运算”(包括加法、乘法、除法和平方根)。这意味着有一个唯一的允许答案，并且它是与所谓的“无限精确”运算结果最接近的可表示浮点数。

在 double 中，数字具有 53 位二进制精度，因此正确答案是四舍五入为 53 位有效数字的精确答案。如 里克·里根 在他的回答中显示，这正是您得到的结果。

你的问题的答案是:

问题 1:允许如此巨大的误差吗？

是的，但称此错误为“巨大”是非常误导的。事实是，没有可以返回具有较小错误的 double 值。

问题 2:这是否意味着每个基本操作都应该有一个错误 < 2.220446e-16 和 64 位 double (machine-epsilon)？

不。这意味着每个基本操作都应该根据当前的舍入模式舍入到(唯一的)最接近的可表示浮点数。这与说相对误差受机器 epsilon 限制并不完全相同。

问题 3:您使用 x86 硬件和 gcc + libc 得到哪个结果？

同样的答案，因为 sqrt在任何合理的平台上正确四舍五入。