matlab和c的区别在于cos函数

Question

我有一个用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不使用Intel内置的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 中的函数相同），看看是否会更快。

Answer 1

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

作为二元有理，double 具有精确 的十进制表示形式，它需要比 15 多得多的小数位来表示（在您的情况下，我相信有 52 或 53 个地方）。然而，printf 和类似函数的标准并不要求第 15 位之后的数字是正确的；它们可能完全是无稽之谈。我怀疑这两个环境之一正在打印精确值，另一个环境正在打印较差的近似值，而实际上两者都对应于完全相同的二进制 double 值。

Answer 2

使用 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

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

>> 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 库中得到相同的答案。

Answer 3

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

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

不幸的是，我不知道这两种情况下的实现细节，实际上 C 实现将取决于您正在使用的 C 标准库实现。

Answer 4

正如其他人所解释的，当您直接在源代码中输入该数字时，不会使用所有小数位，因为您只能获得 15/16 位小数的精度。事实上，它们被转换为最接近的二进制双精度值（超出固定数字限制的任何内容都会被丢弃）。

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

为了说明这一点，我使用以下 MEX 文件进行了测试，使用默认的 LCC 编译器进行编译，然后使用 VS2010（我使用的是 WinXP 32 位）。

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

prec.c

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#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)    <---

copmiled 与使用 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

我还尝试打印数字的十六进制表示形式。我只设法在 C 中显示二进制双精度数的下半部分（也许您可以使用某种位操作获得另一半，但我不确定如何！）

最后，如果您需要更精确的计算，请考虑使用可变精度算术库。在 MATLAB 中，如果您可以访问符号数学工具箱，请尝试：

>> 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 中正确显示浮点数的十六进制表示形式，请改用此辅助函数（类似于 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]);
    }
}