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.

Rather than relying on both programs doing precisely the same steps in their calculations (which can be broken easily, at least on the C code, with a different compiler version or different compiler flags), how about comparing the results like one usually compares floats - `abs(x - y) < epsilon`? — , Oct 01 '11 at 18:10
In the age of SSE instructions, the x87 FPU isn't used much anymore. There's no SSE instruction that does `cos` or other trig functions. So they are implemented in software more often than not now. — Mysticial, Oct 01 '11 at 18:15
@Mysticial, say *what*?! The SSE et al hacks are a way to use the *floating point registers* (i.e., the registers used to do, see, *floating point*) to futz around with short vectors of integers (whch is useful for multimedia). The point is that even if your operating system comes from the '80s (and knows nothing of such newfangled stuff) it will save and restore said registers. But for exactly that same reason, the functions are still implemented in silicon. In any case, it isn't like the CPUs are getting smaller... — vonbrand, Jan 23 '13 at 00:43

score 6 · Answer 1 · answered Oct 01 '11 at 18:44

6

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.

answered Oct 01 '11 at 18:44

R.. GitHub STOP HELPING ICE

208,859
35
376
711

1

It turns out that two doubles are not the same; the last bit is different. – stardt Oct 01 '11 at 19:44
3

IEEE 754 does allow 1ulp ("unit in the last place") of error for transcendental functions... – R.. GitHub STOP HELPING ICE Oct 01 '11 at 19:59
By "the same" I mean having the exact same binary value. I agree that both are sufficiently accurate, but only one is closest to the exact answer. – stardt Oct 01 '11 at 20:04
Right. I was just pointing out that the standard does not require "closest to exact answer" for transcendental functions, just that the error be at most 1ulp (or something pretty similar; I forget the exact wording). – R.. GitHub STOP HELPING ICE Oct 01 '11 at 20:24
1

+1 very concise take on this issue. The difference is still probably in the calc like "R.." says, not printf, because Matlab's printf seems to call down to C's printf, so should be comparable. (I infer this because I've seen Matlab `printf` do conversions that the underlying C lib supports but Matlab doco doesn't mention.) I think you can use the `%x` conversion in printf in both languages to get exact hex output of the double values, eliminating possible printf differences. – Andrew Janke Oct 02 '11 at 05:45
@AndrewJanke: `%x` in C only gives you the lower 32-bits, while in MATLAB we have `%bx` and `%tx` which can be applied to 64bit `double` and 32bit `single` respectively (see the example I posted) – Amro Oct 02 '11 at 17:41
`%x` in C has nothing whatsoever to do with printing floating point. If you want to see the exact value of a float, use `%a` or else read it byte-by-byte with the overlaid `unsigned char [sizeof(double)]` and print each byte with `%x`... – R.. GitHub STOP HELPING ICE Oct 02 '11 at 22:07
@R.., `%x` will give you the underlying byte pattern so you can see if the two values are identical, confirming that the difference originated in the `cos` call. If you don't care about endianness, you can do a single `%x` call instead of looping over the bytes. You're right, `%x` output here isn't meaningful as a number. I didn't know about `%a`, that looks better. – Andrew Janke Oct 03 '11 at 02:24
@Amro: The `%lx` modifier in C will show 64 bits, like `%bx` in Matlab. – Andrew Janke Oct 03 '11 at 02:28
3

@Andrew: No, it will not. `printf` has **undefined behavior** if the type required by the conversion does not match the type of the actual argument passed. This is not just theoretical language lawyering. On real-world systems like x86_64 and basically anything but i386, variadic arguments are passed, like non-variadic ones, in registers, and the integer/pointer registers are separate from the floating point registers. Thus trying to use `%lx` to print a `double` argument will just print the old contents of some integer register. – R.. GitHub STOP HELPING ICE Oct 03 '11 at 02:37
@R..: Ok so `%x` should only be used with integers. But is `%a` part of the standard? because while it seems to work in GCC (I am using Cygwin), it does not give the correct output in VS2010 and completely doesn't work in MinGW-GCC... So, and I guess this should be a posted as a separate question, *how do you correctly convert floats and doubles to IEEE-754 hexadecimal representation in C?* (in a way that works across different compilers/architectures) – Amro Oct 03 '11 at 13:07
@R..: Oh, I see. You're right, `%x` is broken for `double` args. Amro, printing each byte with `%x` as in R..'s earlier comment will work portably and can be used to produce the same output as Matlab's `%bx`. – Andrew Janke Oct 03 '11 at 14:42
2

`%a` is part of the C standard. MSVC does not support the C standard (neither the current one, C99, which introduced `%a`, nor the original one). The solution is to use a C compiler (and standard library) if you want to program in C. – R.. GitHub STOP HELPING ICE Oct 03 '11 at 14:48

stardt · Accepted Answer · 2011-10-01T19:34:54.560

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.

that helped a lot to debug! do you know if there is an equivalent num2bin function for matlab and/or c? — user974735, Oct 02 '11 at 16:06
in matlab there is num2hex(b) which already helps to see the difference, but the num2bin from octave is much more cool for debugging. — user974735, Oct 02 '11 at 16:32
@user974735 I just happened to use num2bin in octave. The same code can be used in matlab as well. Writing an equivalent for C shouldn't be too hard. Let me know if you want such code. — stardt, Oct 09 '11 at 19:40

Charles Keepax · Answer 3 · 2011-10-01T18:27:37.657

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.

score 2 · Answer 4 · edited May 23 '17 at 12:07

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

Try lower case `%lx` instead of `%Lx` to show full width of doubles like Matlab's `%bx`; the `L`/`l` modifiers are different case for integer and floating point conversions. — Andrew Janke, Oct 03 '11 at 02:14
@AndrewJanke: it seems that unlike MATLAB, `%x` in C is to be used with integers only, otherwise it's undefined behavior... (for what it's worth, `%lx` and `%Lx` gave the same thing in the example above) — Amro, Oct 03 '11 at 13:13
My mistake - I ran the test code on my x64 OS X as `printf("%lx",*((unsigned long *)&a))`, where %lx happens to correspond to an 8 byte long. Not portable. `%llx` and a cast to unsigned long long will probably work for your 32-bit Win; but printing individual bytes like discussed in the other answer would be better. And yeah, `%x` without a conversion is broken in C. I'd gotten lucky with that writing against x86; will have to be more careful. — Andrew Janke, Oct 03 '11 at 14:59
@AndrewJanke: I added a C function to correctly print the hex values (print byte-by-byte as suggested) — Amro, Oct 05 '11 at 23:51

matlab and c differ with cos function

4 Answers4

prec.c

copmiled with LCC

compiled with VS2010

UPDATE:

Linked