Atan2 Sin Cos Matlab vs C

Hello, I have an identical implementation of a code, both in Matlab and C. My problem is that the values that C and Matlab returns are sometimes different, after the 10th digit or so, but I would like the two programs to return identical values. For instance:

    x=0.29695007709080329
    y= -0.96764419257510348
    Matlab atan2= 2.8438364863237604
    C atan2= 2.8438364863237608

I have tried to even write mex files that implement the cos, sin and atan2, and I even put the whole computation that I am doing in a mex file and the result is the same. I have set the compiler of the mex file to be the same as the one from C and still, the difference is there. Do you have any idea how can I solve it? Does it have to do with operation precision, or the way that the two programs interpret numbers? Thank you

Answer 1

[Jan]

The documentation of ACOS is determined by the FDLIBM library, as explained in the documentation. But there is no description for the ATAN2 function.

I tried to get equal results for a MEX and a MATLAB function based on ACOS, but this was impossible even when I include the identical FDLIBM function. This is caused by the fact, that the result depends on the compilers optimization flags, the internal floating point precision (the Intel and BCC compilers can store temporary values in a 80 bit register, while MSVC rounds them to 64 bit), the rounding direction, and the usage of SSE2/3/4 commands. The order of execution can be modified automatically by the compiler to increase the speed or accuracy.

If you find a compiler and corresponding flags on your computer, such that the C and MATLAB results are identical, they can be different, if you run the program on another processor.

Another cause for the difference can be the definition of the inputs: If the decimal value has no exact representation as binary number, even this line can create different results in C and MATLAB (just an example, I did not check the number):

x = 0.29695007709080329;

Conclusion: It is usually not possible to reproduce results exactly when different engines (computers, processors, compilers, languages, floating point settings) are used. All floating point calculations suffer from this limited accuracy. Even calculating with quadrupel or octupel precision will not change this.

However, this is very helpful to cross-check an algorithm! If the results are very different on different engines, the method has either a bug or is numerically instable. Or, fortunately not very frequent, the optimizer of the compiler fails to create proper code. E.g. there have been several bugs in MATLAB, which could be solved by "feature accel off", what disables the JIT optimization.

[Walter Roberson]

I do not recall who referred me to this article, but it is relevant to the discussion:

http://scientificcomputing.com/article-hpc-Numerical-Precision-How-Much-is-Enough-063009.aspx

I also saw a more detailed exploration of this topic; if I recall correctly, that more detailed version was by the person who invented the "ULP" (unit in the last place).

Answer 2

[Walter Roberson]

That isn't "after the 10th digit or so": the difference is in the 16th digit, and represents a difference of one unit in the last place (ULP).

Why this difference of 1 ULP? Probably difference in the rounding settings.

The rounding setting that MATLAB uses is not, as far as I can recall, documented. On some systems, on some versions, it can be controlled using an undocumented call: see this CSSM discussion and thank Yair for his detailed work in keeping track of such matters.

[Walter Roberson]

Jan's comments about compiler flags are good ones.

[Dragos Catalin]

But having implemented the function in a mex file, shouldn't it give the same results as in C? I am also using the same compilers for both C and mex files.

[Jan]

@Dragos: Mex files written in C are actually C files, such that the results should be (and are usually) equal. But the floating point environment is set by MATLAB when a mex function is entered. E.g. the floating point precision is set to 53 bit mantissa, such that intermediate results are not stored in the 80 bit registers of the Intel processors. So if you want to guarantee identical results, use _control87 or similiar platform dependent commands to control the floating point directives exactly. In addition you have to check, if the mexopts.bat file and the C-compiler use identical flags for the compilation: E.g. changing /O2 to /Ox or using /fp:fast influence the results for the MSVC compilers.