matlab and c differ with cos function

Question

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 function from c (as it is the same from matlab), see if that goes faster.

Answer 1

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.

Answer 2

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.

Answer 3

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.