将三角函数的正确值与matlab进行比较

我正在尝试用它的c++代码测试一个simulink块，该simulink区块包含一些代数、三角函数和积分器。在我的测试过程中，一个随机数生成器被用于simulink块输入，输入和输出都被记录到mat文件中(使用MatIO)，该文件将由C++代码读取，并将输出与C++计算结果进行比较。对于只包含代数函数的信号，结果是精确的，并且差为零；对于包含三角函数的路径，差约为10e-16。matlab社区声称它们是正确的，而glibc则不然。

最近，我发现在glibc中实现的三角函数的输出值不等于在matlabs中产生的值，根据旧问题1 2 3和我的实验，这些差异与glibc的1ulf>精度有关。对于块的大部分，这个10e-16误差没有太大意义，但在积分器的输出中，10e-16累积得越来越多，积分器的最终误差将是大约1e-3，这有点高，对于那种块来说是不可接受的。

经过对这个问题的大量研究，我决定使用其他方法来计算sin/cos函数，而不是glibc中提供的方法。

我实施了这些方法，

具有长双变量和-O2的1-taylor级数(强制使用x87 FPU及其80位浮点运算)

带有GNU四进制库的2-泰勒级数(128位精度)

3-MPFR库(128位)

4-CRLibm(正确取整的libm)

5-Sun的LibMCR(就像CRLibm)

具有不同舍入模式的6-X86 FSIN/FCOS

7-Java.lang.math通过JNI(我认为matlab使用)

8-fdlibm(根据我看到的一篇博客文章)

9-openlibm

10-通过mex/matlab引擎调用matlab函数

除了最后一个实验外，上面的所有实验都不能产生与matlab相等的值。我针对各种输入测试了所有这些库和方法，其中一些库和方法(如libmcr和fdlibm)将为一些输入生成NAN值(看起来它们没有很好的范围检查)，其余库和方法生成的值的误差为10e-16或更高。与matlab相比，只有最后一个函数产生了正确的值，但调用matlab函数并不有效，而且比原生实现慢得多。

并对长二重和四重数学的MPFR和taylor级数为什么会出错表示支持。

这是具有长双变量的泰勒级数(精度为80位)并且应该使用-O2进行编译，这可以防止将FPU堆栈中的值存储到寄存器中(80位到64位=精度损失)，在进行任何计算之前，x87的舍入模式将设置为最接近的

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

taylor级数方法在x87的所有四种舍入模式下进行了测试，其中最好的一种具有10e-16的误差

这是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;
}

而且x87 fpu代码并不比以前的更准确

这是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;
}

我不明白为什么MPFR版本没有像预期的那样工作

此外，我测试过的所有其他方法的代码都是相同的，与matlab相比，它们都有错误。

所有的代码都经过了广泛的数字测试，我发现了一些简单的失败案例。例如：

在matlab中，以下代码生成0x3fe1b0071cef86fbe，但在这些近似中，我得到了0x3fe1 b0071cep86fbf(最后一位的差异)

format hex;
sin(0.5857069572718263)
ans = 0x3fe1b071cef86fbe

为了明确这个问题，正如我上面所描述的，当它被输入积分器时，这个一位的不准确度是很重要的，我正在寻找一个解决方案，以获得与matlab完全相同的值。有什么想法吗？

更新1:

1 Ulp误差完全不影响算法的输出，但它阻止了与matlab结果的验证，特别是在积分器的输出中。

正如@John Bollinger所说，错误不会在多个算术块的直接路径中累积，但不会在输入离散积分器时累积

Update2:我统计了上述所有方法的不等结果的数量，很明显，与matlab相比，openlibm将产生更少的不等值，但它不是零。