matlab和c与cos函数不同

我有一个在matlab中实现的程序和c中的相同程序,结果不同。

我有点困惑的是cos函数没有返回完全相同的结果。

在这两种情况下,我使用相同的计算机,Intel Core 2 Duo和8字节双数据类型。

为什么结果不同?

这是测试:

c: double a = 2.89308776595231886830; double b = cos(a); printf("a = %.50f\n", a); printf("b = %.50f\n", b); printf("sizeof(a): %ld\n", sizeof(a)); printf("sizeof(b): %ld\n", sizeof(b)); a = 2.89308776595231886830106304842047393321990966796875 b = -0.96928123535654842068964853751822374761104583740234 sizeof(a): 8 sizeof(b): 8 matlab: a = 2.89308776595231886830 b = cos(a); fprintf('a = %.50f\n', a); fprintf('b = %.50f\n', b); whos('a') whos('b') a = 2.89308776595231886830106304842047393321990966796875 b = -0.96928123535654830966734607500256970524787902832031 Name Size Bytes Class Attributes a 1x1 8 double Name Size Bytes Class Attributes b 1x1 8 double So, b differ a bit (very slightly, but enough to make my debuging task difficult) b = -0.96928123535654842068964853751822374761104583740234 c b = -0.96928123535654830966734607500256970524787902832031 matlab 

我使用相同的计算机,Intel Core 2 Duo和8字节双数据类型。

为什么结果不同?

matlab不使用英特尔内置的cos函数硬件吗?

有没有一种简单的方法在matlab和c中使用相同的cos函数(具有精确的结果),即使有点慢,这样我可以安全地比较我的matlab和c程序的结果?


更新:

非常感谢你的回答!

因此,正如您所指出的,matlab和c的cos函数不同。 棒极了! 我以为他们正在使用英特尔微处理器内置的cos函数。

matlab的cos版本与matlab相同(至少对于此测试而言)。 你可以尝试从matlab:b = java.lang.Math.cos(a)

然后,我做了一个小的MEX函数来使用matlab中的cos c版本,它工作正常; 这允许我调试我的程序(在matlab和c中实现的相同程序),看看它们在哪些方面有所不同,这就是本文的目的。

唯一的事情是从matlab调用MEX c cos版本太慢了。

我现在试图从c调用Java cos函数(因为它与matlab相同),看看它是否更快。

使用http://www.mathworks.com/matlabcentral/fileexchange/1777-from-double-to-string上的脚本

两个数字之间的差异仅在最后一位:

 octave:1> bc = -0.96928123535654842068964853751822374761104583740234; octave:2> bm = -0.96928123535654830966734607500256970524787902832031; octave:3> num2bin(bc) ans = -.11111000001000101101000010100110011110111001110001011*2^+0 octave:4> num2bin(bm) ans = -.11111000001000101101000010100110011110111001110001010*2^+0 

其中一个必须更接近“正确”的答案,假设给出的值是准确的。

 >> be = vpa('cos(2.89308776595231886830)',50) be = -.96928123535654836529707365425580405084360377470583 >> bc = -0.96928123535654842068964853751822374761104583740234; >> bm = -0.96928123535654830966734607500256970524787902832031; >> abs(bc-be) ans = .5539257488326242e-16 >> abs(bm-be) ans = .5562972757925323e-16 

因此,C库结果更准确。

但是,出于您的问题的目的,您不应期望在matlab和您链接的C库中获得相同的答案。

浮点数以二进制forms存储,而不是十进制。 double精度浮点数具有52位精度,可转换为大约15位有效小数位。 换句话说,以十进制打印的double精度的前15个非零十进制数字足以唯一地确定打印了哪个double精度数。

作为一个二元理性double有一个十进制的精确表示,它需要比15代表更多的小数位(在你的情况下,52或53个位置,我相信)。 但是, printf和类似function的标准不要求超过15的数字是正确的; 他们可能完全是胡说八道。 我怀疑两种环境中的一种是打印精确值,另一种是打印不良的近似值,实际上两者都对应于完全相同的二进制double值。

结果是相同的最多15个小数位,我怀疑这对几乎所有的应用程序都足够了,如果你需要更多,你应该可能实现自己的余弦版本,这样你就可以控制细节并且你的代码是可移植的跨越不同的C编译器。

它们会有所不同,因为它们无疑使用不同的方法来计算结果的近似值或迭代不同的次数。 由于余弦定义为无穷大的一系列术语,因此必须使用近似值来实现其软件。 CORDIC算法是一种常见的实现方式。

不幸的是,在任何一种情况下我都不知道实现的细节,实际上C语言将取决于您正在使用的C标准库实现。

正如其他人所解释的那样,当您在源代码中直接输入该数字时,并不会使用所有小数位数,因为您只能得到15/16小数位数以获得精确度。 实际上,它们会以二进制forms转换为最接近的double值(超出固定数字限制的任何值都会被删除)。

更糟糕的是,根据@R,IEEE 754在使用余弦函数时容忍最后一位的错误。 我在使用不同的编译器时遇到了这个问题。

为了说明,我使用以下MEX文件进行了测试,一旦使用默认的LCC编译器编译,然后使用VS2010(我在WinXP 32位上)。

在一个函数中,我们直接调用C函数( mexPrintf只是一个宏#define作为printf )。 另一方面,我们调用mexEvalString来评估MATLAB引擎中的东西(相当于在MATLAB中使用命令提示符)。

prec.c

 #include  #include  #include  #include "mex.h" void c_test() { double a = 2.89308776595231886830L; double b = cos(a); mexPrintf("[C] a = %.25Lf (%16Lx)\n", a, a); mexPrintf("[C] b = %.25Lf (%16Lx)\n", b, b); } void matlab_test() { mexEvalString("a = 2.89308776595231886830;"); mexEvalString("b = cos(a);"); mexEvalString("fprintf('[M] a = %.25f (%bx)\\n', a, a)"); mexEvalString("fprintf('[M] b = %.25f (%bx)\\n', b, b)"); } void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { matlab_test(); c_test(); } 

与LCC合作

 >> prec [M] a = 2.8930877659523189000000000 (4007250b32d9c886) [M] b = -0.9692812353565483100000000 (bfef045a14cf738a) [C] a = 2.8930877659523189000000000 ( 32d9c886) [C] b = -0.9692812353565484200000000 ( 14cf738b) <--- 

用VS2010编译

 >> prec [M] a = 2.8930877659523189000000000 (4007250b32d9c886) [M] b = -0.9692812353565483100000000 (bfef045a14cf738a) [C] a = 2.8930877659523189000000000 ( 32d9c886) [C] b = -0.9692812353565483100000000 ( 14cf738a) <--- 

我使用以下命令编译以上内容: mex -v -largeArrayDims prec.c ,并使用以下命令在后端编译器之间切换: mex -setup

请注意,我还尝试打印数字的hex表示。 我只设法在C中显示二进制双数的下半部分(也许你可以使用某种位操作得到另一半,但我不确定如何!)

最后,如果在计算中需要更高的精度,请考虑使用库进行变量精度算术。 在MATLAB中,如果您可以访问Symbolic Math Toolbox ,请尝试:

 >> a = sym('2.89308776595231886830'); >> b = cos(a); >> vpa(b,25) ans = -0.9692812353565483652970737 

所以你可以看到实际值介于我上面得到的两个不同的近似值之间,实际上它们都等于第15个小数位:

 -0.96928123535654831.. # 0xbfef045a14cf738a -0.96928123535654836.. # <--- actual value (cannot be represented in 64-bit) -0.96928123535654842.. # 0xbfef045a14cf738b ^ 15th digit --/ 

更新:

如果要在C中正确显示浮点数的hex表示,请改用此辅助函数(类似于MATLAB中的NUM2HEX函数):

 /* you need to adjust for double/float datatypes, big/little endianness */ void num2hex(double x) { unsigned char *p = (unsigned char *) &x; int i; for(i=sizeof(double)-1; i>=0; i--) { printf("%02x", p[i]); } }