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

Question

我正在将一些程序从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++中都是不同的)所以最终的差异是巨大的)(如果你很好奇,看看差异如何开始变大,这里是完整输出[很快链接],但这与问题无关)

Answer 1

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

IEEE754 2008标准不要求函数被正确舍入(尽管它"推荐"它).一些C libms确实提供了正确的舍入三角函数:我相信glibc libm(通常用于大多数Linux发行版),就像CRlibm一样.大多数其他现代libms将提供1 ulp以内的触发函数(即真值两侧的两个浮点值之一),通常称为忠实舍入,计算速度更快.

您打印的这些值中没有一个实际上可以作为IEEE 64位浮点值出现(即使是舍入的):最接近的3(打印到全精度)是:

0.477158760259608 405451814405751065351068973541259765625

0.477158760259608 46096296563700889237225055694580078125

0.477158760259608 516474116868266719393432140350341796875

您可能需要的值是:

1. 十进制的确切罪.497418836818383950,这是

0.477158760259608 433132061388630377105954125778369485736356219 ...

(这似乎是Mathematica给出的).

2. 最接近的64位浮点数.497418836818383950:

0.477158760259608 430531153841011107415427334794384396325832953 ...

在这两种情况下,上面列表中的第一个是最近的(尽管只有1的情况).