5

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

测试我的 fdlibm 的 sqrt 端口，我在 64 位双精度上使用 sqrt 得到以下计算：

sqrt(1977061516825203605555216616167125005658976571589721139027150498657494589171970335387417823661417383745964289845929120708819092392090053015474001800648403714048.0) = 44464159913633855548904943164666890000299422761159637702558734139742800916250624.0

（这个案例在我关于精度的测试中打破了一个简单的后置条件；我不确定这个后置条件是否可以通过 IEEE-754 实现）

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

sqrt(1977061516825203605555216616167125005658976571589721139027150498657494589171970335387417823661417383745964289845929120708819092392090053015474001800648403714048.0) =44464159913633852501611468455197640079591886932526256694498106717014555047373210.truncated

可以看到，左边的第 17 个数字是不同的，这意味着如下错误：

3047293474709469249920707535828633381008060627422728245868877413.0

问题 1：允许这么大的错误吗？

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

问题 2：这是否意味着，每个基本操作都应该有一个错误 < 2.220446e-16 和 64 位双精度数（机器 epsilon）？

我确实用 x86-32 linux 系统（glibc / eglibc）计算了相同的结果，并得到了与 fdlibm 相同的结果，这让我认为：

  • a：我做错了什么（但是如何：printf会成为候选人，但我不知道这是否可能是原因）
  • b：错误/精度在这些库中很常见
4

2 回答 2

16

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

在双精度中，数字具有 53 位二进制精度，因此正确答案是四舍五入到 53 位有效数字的确切答案。正如Rick Regan在他的回答中所表明的那样，这正是你得到的结果。

您的问题的答案是：

问题 1：允许这么大的错误吗？

是的，但是将此错误称为“巨大”是一种误导。事实是，没有可以返回的双精度值具有更小的误差。

问题 2：这是否意味着，每个基本操作都应该有一个错误 < 2.220446e-16 和 64 位双精度数（机器 epsilon）？

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

问题 3：您使用 x86 硬件和 gcc + libc 获得了哪个结果？

您所做的相同答案，因为sqrt在任何合理的平台上都正确四舍五入。

于 2010-11-30T22:04:20.473 回答
8

在二进制中，任意精度答案的前 58 位是 1011111111111111111111110101010101111111111111111011010001...

双精度值的 53 位有效位是

10111111111111111111111101010101011111111111111110111

这意味着双精度值正确舍入到 53 个有效位，并且在 1/2 ULP 之内。（错误是“大”只是因为数字本身很大）。

于 2010-11-30T21:44:36.403 回答