matlab and c differ with cos function

I have a program implemented in matlab and the same program in c, and the results differ.

I am bit puzzled that the cos function does not return the exact same result.

I use the same computer, Intel Core 2 Duo, and 8 bytes double data type in both cases.

Why does the result differ?

Here is the test:

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

I use the same computer, Intel Core 2 Duo, and 8 bytes double data type.

Why does the result differ?

does matlab do not use the cos function hardware built-in in Intel?

Is there a simple way to use the same cos function in matlab and c (with exact results), even if a bit slower, so that I can safely compare the results of my matlab and c program?

Update:

thanks a lot for your answers!

So, as you have pointed out, the cos function for matlab and c differ. That's amazing! I thought they were using the cos function built-in in the Intel microprocessor.

The cos version of matlab is equal (at least for this test) to the one of matlab. you can try from matlab also: b=java.lang.Math.cos(a)

Then, I did a small MEX function to use the cos c version from within matlab, and it works fine; This allows me to debug the my program (the same one implemented in matlab and c) and see at what point they differ, which was the purpose of this post.

The only thing is that calling the MEX c cos version from matlab is way too slow.

I am now trying to call the Java cos functi开发者_StackOverflow社区on from c (as it is the same from matlab), see if that goes faster.

Floating point numbers are stored in binary, not decimal. A double precision float has 52 bits of precision, which translates to roughly 15 significant decimal places. In other words, the first 15 nonzero decimal digits of a double printed in decimal are enough to uniquely determine which double was printed.

As a diadic rational, a double has an exact representation in decimal, which takes many more decimal places than 15 to represent (in your case, 52 or 53 places, I believe). However, the standards for printf and similar functions do not require the digits past the 15th to be correct; they could be complete nonsense. I suspect one of the two environments is printing the exact value, and the other is printing a poor approximation, and that in reality both correspond to the exact same binary double value.

Using the script at http://www.mathworks.com/matlabcentral/fileexchange/1777-from-double-to-string

the difference between the two numbers is only in the last bit:

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

One of them must be closer to the "correct" answer, assuming the value given for a is exact.

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

So, the C library result is more accurate.

For the purposes of your question, however, you should not expect to get the same answer in matlab and whichever C library you linked with.

The result is the same up to 15 decimal places, I suspect that is sufficient for almost all applications and if you require more you should probably be implementing your own version of cosine anyway such that you are in control of the specifics and your code is portable across different C compilers.

They will differ because they undoubtedly use different methods to calculate the approximation to the result or iterate a different number of times. As cosine is defined as an infinite series of terms an approximation must be used for its software implementation. The CORDIC algorithm is one common implementation.

Unfortunately, I don't know the specifics of the implementation in either case, indeed the C one will depend on which C standard library implementation you are using.

As others have explained, when you enter that number directly in your source code, not all the fraction digits will be used, as you only get 15/16 decimal places for precision. In fact, they get converted to the nearest double value in binary (anything beyond the fixed limit of digits is dropped).

To make things worse, and according to @R, IEEE 754 tolerates error in the last bit when using the cosine function. I actually ran into this when using different compilers.

To illustrate, I tested with the following MEX file, once compiled with the default LCC compiler, and then using VS2010 (I am on WinXP 32-bit).

In one function we directly call the C functions (mexPrintf is simply a macro #define as printf). In the other, we call mexEvalString to evaulate stuff in the MATLAB engine (equivalent to using the command prompt in 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();
}

copmiled with LCC

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

compiled with VS2010

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

I compile the above using: mex -v -largeArrayDims prec.c, and switch between the backend compilers using: mex -setup

Note that I also tried to print the hexadecimal representation of the numbers. I only managed to show the lower half of binary double numbers in C (perhaps you can get the other half using some sort of bit manipulations, but I'm not sure how!)

Finally, if you need more precision in you calculations, consider using a library for variable precision arithmetic. In MATLAB, if you have access to the Symbolic Math Toolbox, try:

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

So you can see that the actual value is somewhere between the two different approximations I got above, and in fact they are all equal up to the 15th decimal place:

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

UPDATE:

If you want to correctly display the hexadecimal representation of floating point numbers in C, use this helper function instead (similar to NUM2HEX function in MATLAB):

/* 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]);
    }
}