`std::sin` 最后一位是错误的答案

【问题标题】：`std::sin` is wrong in the last bit`std::sin` 最后一位是错误的
【发布时间】：2015-08-12 08:15:47
【问题描述】：

为了提高效率，我正在将一些程序从 Matlab 移植到 C++。两个程序的输出必须完全相同 (**)。

此操作我面临不同的结果：

std::sin(0.497418836818383950)   = 0.477158760259608410 (C++)
sin(0.497418836818383950)        = 0.47715876025960846000 (Matlab)
N[Sin[0.497418836818383950], 20] = 0.477158760259608433 (Mathematica)

所以，据我所知，C++ 和 Matlab 都使用 IEEE754 定义的双精度算法。我想我在某处读过 IEEE754 允许在最后一位产生不同的结果。使用mathematica 来决定，似乎C++ 更接近结果。 如何强制 Matlab 精确计算包含最后一位的 sin，以使结果相同？

在我的程序中，这种行为会导致很大的错误，因为数值微分方程求解器会在最后一位中不断增加这个错误。但是我不确定 C++ 移植版本是否正确。我猜 即使 IEEE754 允许最后一位不同，以某种方式保证在更多 IEEE754 定义的双重操作中使用结果时此错误不会变大（因为否则，两个不同的程序正确根据 IEEE754 标准可以产生完全不同的输出）。所以另一个问题是我说得对吗？

我想得到两个粗体问题的答案。 编辑：第一个问题颇有争议，但不太重要，有人可以评论第二个问题吗？

注意：这不是打印错误，以防万一你想检查，这是我得到这些结果的方式：

http://i.imgur.com/cy5ToYy.png

注意（**）：我的意思是最终输出，即一些计算的结果，显示一些小数点后 4 位的实数，需要完全相同。我在问题中谈到的错误变得更大（因为更多的操作，在 Matlab 和 C++ 中的每一个都是不同的）所以最终的差异是巨大的）（如果你很好奇看到差异是如何开始变大的，这里是完整的输出[链接很快]，但这与问题无关）

【问题讨论】：

IEEE-754 doubles 的精度仅为 15–17 significant decimal digits。
@Trollkemada It is important for the output of both programs to be exactly the same 你会等着戈多。甚至不能保证用不同的编译器编写的 C++ 程序的行为与浮点数完全相同。
Re 重要的是两个程序的输出完全相同。 相反，重要的是你应该容忍一些最小的差异。您期望完全匹配但没有得到它是您的问题，而不是 C++、Matlab 或 Mathematica。你需要改变你的期望。最好不要期望跨平台、产品或什至完全相同产品的不同编译的浮点数完全匹配。
基于“在我的程序中，这种行为会导致很大的错误，因为数值微分方程求解器会在最后一位中不断增加这个错误。”你的问题比正弦结果的最后一点要大得多。请记住，原始输入数据可能来自测量精度远低于 10^15 中的一部分。
您需要考虑您尝试求解混沌方程组的可能性。微小的差异迅速变得显着是混乱的典型。

标签： c++ matlab floating-point ieee-754

【解决方案1】：

首先，如果您的数值方法依赖于 sin 到最后一位的精度，那么您可能需要使用任意精度库，例如 MPFR。

IEEE754 2008 标准不要求对函数进行正确舍入（尽管它确实“推荐”了它）。一些 C libm 确实提供了正确舍入的三角函数：我相信 glibc libm 和 CRlibm 一样（通常用于大多数 linux 发行版）。大多数其他现代 libms 将提供 1 ulp 以内的三角函数（即真值两侧的两个浮点值之一），通常称为 忠实舍入，计算速度要快得多。

您打印的这些值实际上都不是 IEEE 64 位浮点值（即使四舍五入）：最接近的 3 个（打印到全精度）是：

0.477158760259608 405451814405751065351068973541259765625

0.477158760259608 46096296563700889237225055694580078125

0.477158760259608 516474116868266719393432140350341796875

您可能想要的值是：

小数点 0.497418836818383950 的正弦值，即

0.477158760259608 433132061388630377105954125778369485736356219...

（这似乎是 Mathematica 提供的）。

最接近 0.497418836818383950 的 64 位浮点数的确切 sin：

0.477158760259608 430531153841011107415427334794384396325832953...

在这两种情况下，上面列表中的第一个是最近的（尽管仅在 1 的情况下几乎没有）。

【讨论】：

【解决方案2】：

你写的double常数的正弦值大约是0x1.e89c4e59427b173a8753edbcb95p-2，它最近的double是0x1.e89c4e59427b1p-2。到小数点后 20 位，最接近的两个 doubles 是 0.47715876025960840545 和 0.47715876025960846096。

也许 Matlab 显示的是截断值？（编辑：我现在看到倒数第四个数字是 6，而不是 0。Matlab 给你的结果仍然是四舍五入的，但它是最接近所需结果的两个 doubles 中较远的一个。而且它仍然打印出错误的数字。

我还应该指出，Mathematica 可能正在尝试解决一个不同的问题——将十进制数 0.497418836818383950 的正弦计算到小数点后 20 位。您不应该期望这与 C++ 代码的结果或 Matlab 的结果相匹配。

【讨论】：

从 Matlab 和 C++ 输出的两个双精度常量以及与实际结果最接近的可表示实数值（检查 i.imgur.com/ZSUyPmM.png）。（如问题中所述，计算并未精确到小数点后 20 位，而是使用 IEEE 754 标准）
@Trollkemada：不，Matlab 结果的最后几位数字是伪造的。没有double 有一个以这样开头的小数扩展。