I am trying to test a simulink block with its c++ code, the simulink block contains some algebratic, trigonometric functions and integrators. In my test procedure, a random number generator is used form the simulink block input and both input and outpus are recorded into mat file (using MatIO) which will be read by C++ code and the outputs are compared with C++ calculate ones. for signals containing only algebraic functions the results are exact and the difference is zero, for paths which contains the trigonometric functions the difference is about 10e-16. The matlab community claim they are correct and glibc isn't.
Recently i discovered the output values of trigonometric functions implemented in glibc isn't equal to values produced in matlabs, according to old questions 1 2 3 and my experiments the differences related 1ulp> accuracy of glibc. for most of the block this 10e-16 error doesn't sense much, but in the output of the integrator the 10e-16 accumulated more and more and the final error of the integrator will be about 1e-3 which is a bit high and isn't acceptable for that kind of block.
After lot of research about the problem, i decided to use other approaches for calculating sin/cos functions than those provided in glibc.
I implemented these apporaches,
1- taylor series with long double variables and -O2 (which force using x87 FPU with its 80bit floating point arithmetic)
2- taylor series with GNU quadmath library (128bit precision)
3- MPFR library (128 bit)
4- CRLibm (correctly rounded libm)
5- Sun's LibMCR ( just like CRLibm )
6- X86 FSIN/FCOS with different rounding modes
7- Java.lang.math through JNI ( as i think matlab uses )
8- fdlibm ( according to one the blogpost i have seen )
9- openlibm
10- calling matlab function through mex/matlab engine
Non of the experiments above except last one couldn't generate values equal to matlab. i tested all of those libraries and approaches for wide range of inputs, some of them like libmcr and fdlibm will produce NAN value for some of the inputs (looks like they doesn't have good range checkings) and the rest of them producing values with error of 10e-16 and higher. Only the last one produces correct values compared to matlab as was expeceted, but calling matlab function isn't efficent and much slower than native implementations.
also i supprised why the MPFR and taylor series with long double and quadmath are getting in error.
This is taylor series with long double variables (80bit precision) and should be compiled with -O2 which prevents storing values from FPU stack into registers (80bit to 64bit = precision loss), also before doing any calculations the rounding mode of x87 will be set to nearest
typedef long double dt_double;
inline void setFPUModes(){
unsigned int mode = 0b0000111111111111;
asm(
"fldcw %0;"
: : "m"(mode));
}
inline dt_double factorial(int x) //calculates the factorial
{
dt_double fact = 1;
for (; x >= 1 ; x--)
fact = x * fact;
return fact;
}
inline dt_double power(dt_double x, dt_double n) //calculates the power of x
{
dt_double output = 1;
while (n > 0)
{
output = (x * output);
n--;
}
return output;
}
inline double sin(double x) noexcept //value of sine by Taylors series
{
setFPUModes();
dt_double result = x;
for (int y = 1 ; y != 44; y++)
{
int k = (2 * y) + 1;
dt_double a = (y%2) ? -1.0 : 1.0;
dt_double c = factorial(k);
dt_double b = power(x, k);
result = result + (a * b) / c;
}
return result;
}
the taylor series approach tested with all four rounding modes of x87, the best one has error of 10e-16
This is the X87 fpu one
double sin(double x) noexcept
{
double d;
unsigned int mode = 0b0000111111111111;
asm(
"finit;"
"fldcw %2;"
"fldl %1;"
"fsin;"
"fstpl %0" :
"+m"(d) : "m"(x), "m"(mode)
);
return d;
}
also the x87 fpu code isn't more accurate than previous one
Here is the code for MPFR
double sin(double x) noexcept{
mpfr_set_default_prec(128);
mpfr_set_default_rounding_mode(MPFR_RNDN);
mpfr_t t;
mpfr_init2(t, 128);
mpfr_set_d(t, x, MPFR_RNDN);
mpfr_t y;
mpfr_init2(y, 128);
mpfr_sin(y, t, MPFR_RNDN);
double d = mpfr_get_d(y, MPFR_RNDN);
mpfr_clear(t);
mpfr_clear(y);
return d;
}
I can't understand why the MPFR version didn't work as expected
also the codes for all other approaches i have tested is same and all of them have errors compared to matlab.
all the codes are tested for wide range of numbers, and i found simple cases which they fail. for example :
in matlab following code produces 0x3fe1b071cef86fbe but in these apporoches i got 0x3fe1b071cef86fbf ( difference in last bit )
format hex;
sin(0.5857069572718263)
ans = 0x3fe1b071cef86fbe
To be clear about the question, As i described above, this one bit inaccuracy is important when it fed into the integrator and i am looking for a solution to get values exactly same as the matlab one. any ideas?
Update1:
1 Ulp error doesn't affect the alogorithm output at all, but it prevents verfication with matlab results specially in the output of integrators.
As @John Bollinger said, the errors doesn't accumulate in direct path with multiple arithmetic blocks, but not when feed into the discrete integrator
Update2: I counted the number of inequal results for all of the approaches above, clearly openlibm will produce less inequal values compared to matlab ones, but it isn't zero.
