`std :: sin`在最后一点是错误的 [英] `std::sin` is wrong in the last bit

问题描述

为了提高效率，我正在将一些程序从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 ++中每个操作都不相同)，所以最终的差异很大)(如果您很好奇地看到差异如何开始变大，请在这里是完整的输出[很快就会链接]，但这与问题无关)

解决方案

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

IEEE754 2008标准不要求对函数进行正确的四舍五入(尽管确实推荐"了它).一些C libms确实提供了正确的舍入三角函数:我相信glibc libm确实(通常在大多数linux发行版中使用)，CRlibm也是如此.大多数其他现代libms都将提供1 ulp以内的trig函数(即，真值两侧的两个浮点值之一)，通常称为忠实舍入，它的计算速度更快.

您打印的所有值实际上都不会出现为IEEE 64位浮点值(即使是四舍五入的):最接近的3个(打印为全精度)是:

0.477158760259608 405451814405751065351068973541259765625

0.477158760259608 46096296563700889237225055694580078125

0.477158760259608 516474116868266719393432140350341796875

您可能想要的值是:

1. 十进制.497418836818383950的确切正弦，即

0.477158760259608 433132061388630377105954125778369485736356219 ...

(这似乎是Mathematica所提供的).

2. 最接近.497418836818383950的64位浮点数的确切正弦:

0.477158760259608 430531153841011107415427334794384396325832953 ...

在这两种情况下，上面列表中的第一个都是最近的(尽管在1的情况下才勉强是这样的.)

I am porting some program from Matlab to C++ for efficiency. It is important for the output of both programs to be exactly the same (**).

I am facing different results for this operation:

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

So, as far as I know both C++ and Matlab are using IEEE754 defined double arithmetic. I think I have read somewhere that IEEE754 allows differents results in the last bit. Using mathematica to decide, seems like C++ is more close to the result. How can I force Matlab to compute the sin with precision to the last bit included, so that the results are the same?

In my program this behaviour leads to big errors because the numerical differential equation solver keeps increasing this error in the last bit. However I am not sure that C++ ported version is correct. I am guessing that even if the IEEE754 allows the last bit to be different, somehow guarantees that this error does not get bigger when using the result in more IEEE754 defined double operations (because otherwise, two different programs correct according to the IEEE754 standard could produce completely different outputs). So the other question is Am I right about this?

I would like get an answer to both bolded questions. Edit: The first question is being quite controversial, but is the less important, can someone comment about the second one?

Note: This is not an error in the printing, just in case you want to check, this is how I obtained these results:

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

Note (**): What I mean by this is that the final output, which are the results of some calculations showing some real numbers with 4 decimal places, need to be exactly the same. The error I talk about in the question gets bigger (because of more operations, each of one is different in Matlab and in C++) so the final differences are huge) (If you are curious enough to see how the difference start getting bigger, here is the full output [link soon], but this has nothing to do with the question)

解决方案

Firstly, if your numerical method depends on the accuracy of sin to the last bit, then you probably need to use an arbitrary precision library, such as MPFR.

The IEEE754 2008 standard doesn't require that the functions be correctly rounded (it does "recommend" it though). Some C libms do provide correctly rounded trigonometric functions: I believe that the glibc libm does (typically used on most linux distributions), as does CRlibm. Most other modern libms will provide trig functions that are within 1 ulp (i.e. one of the two floating point values either side of the true value), often termed faithfully rounded, which is much quicker to compute.

None of those values you printed could actually arise as IEEE 64bit floating point values (even if rounded): the 3 nearest (printed to full precision) are:

0.477158760259608 405451814405751065351068973541259765625

0.477158760259608 46096296563700889237225055694580078125

0.477158760259608 516474116868266719393432140350341796875

The possible values you could want are:

1. The exact sin of the decimal .497418836818383950, which is

0.477158760259608 433132061388630377105954125778369485736356219...

(this appears to be what Mathematica gives).

2. The exact sin of the 64-bit float nearest .497418836818383950:

0.477158760259608 430531153841011107415427334794384396325832953...

In both cases, the first of the above list is the nearest (though only barely in the case of 1).

这篇关于`std :: sin`在最后一点是错误的的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！