我试图用它的c ++代码测试simulink块，simulink块包含一些algebratic，三角函数和积分器。在我的测试过程中，从simulink块输入使用随机数生成器，输入和输出都被记录到mat文件中（使用MatIO），它将由C ++代码读取，输出与C ++计算结果进行比较。对于仅包含代数函数的信号，结果是精确的，差值为零，对于包含三角函数的路径，差值约为10e-16。 matlab社区声称他们是正确的而glibc不是。
最近我发现在glibc中实现的三角函数的输出值不等于matlabs中产生的值，根据旧问题1 2 3和我的实验，差异相关1ulp> glibc的准确性。对于大多数块而言，这个10e-16误差感觉不大，但是在积分器的输出中，10e-16积累的越来越多，积分器的最终误差将是大约1e-3，这有点高，这种阻止是不可接受的。
经过对该问题的大量研究后，我决定使用其他方法来计算sin / cos函数，而不是glibc中提供的函数。
我实施了这些方法，
1-泰勒系列具有长双变量和-O2（强制使用x87 FPU及其80位浮点运算）
2-taylor系列与GNU quadmath库（128位精度）
3个MPFR库（128位）
4- CRLibm（正确舍入的libm）
5- Sun的LibMCR（就像CRLibm一样）
6- X86 FSIN / FCOS具有不同的舍入模式
7- Java.lang.math到JNI（我认为matlab使用）
8-fdlibm（根据我见过的一篇博文）
9开放的libm
10-通过mex / matlab引擎调用matlab函数
除了最后一个以外的实验不能生成等于matlab的值。我测试了所有这些库和方法的广泛输入，其中一些像libmcr和fdlibm将为一些输入产生NAN值（看起来他们没有良好的范围检查），其余的产生值错误10e-16及更高。与matlab相比，只有最后一个产生正确的值，但是调用matlab函数并不比本机实现快得多且慢得多。
另外我还说为什么MPFR和泰勒系列长双重和四重奏都出错了。
这是具有长双变量（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;
}
使用x87的所有四种舍入模式测试的泰勒系列方法，最好的一种方法具有10e-16的误差
这是X87的fpu之一
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下面的代码产生0x3fe1b071cef86fbe，但在这些apporoches中我得到了0x3fe1b071cef86fbf（最后一位的差异）
format hex;
sin(0.5857069572718263)
ans = 0x3fe1b071cef86fbe
要明确这个问题，正如我上面所描述的那样，当它输入积分器并且我正在寻找一种解决方案来获得与matlab完全相同的值时，这一点不准确是很重要的。有任何想法吗？
UPDATE1：
1 Ulp错误根本不会影响算法输出，但它会阻止使用matlab结果进行验证，特别是在积分器的输出中。
正如@John Bollinger所说，误差不会累积在具有多个算术块的直接路径中，而是在馈入离散积分器时不会累积
Update2：我计算了上述所有方法的不等结果数，显然openlibm与matlab相比会产生更少的不等值，但它不是零。
我的猜测是Matlab使用的是最初基于FDLIBM的代码。我能用朱莉娅（使用openlibm）获得相同的结果：你可以尝试使用它，或者musl，我相信它也使用相同的代码。
最接近的double / IEEE二进制64到0.5857069572718263是
0.5857069572718263117394599248655140399932861328125
（有位模式0x3fe2be1c8450b590）
这是sin
0.55278864311139114312806521962078480744570117018100444956428008387067038680572587...
这两个最接近的double / IEEE二进制64是
a）0.5527886431113910870038807843229733407497406005859375（0x3fe1b071cef86fbe），误差为0.5055 ulps
b）0.55278864311139119802618324683862738311290740966796875（0x3fe1b071cef86fbf），误差为0.4945 ulps
FDLIBM仅保证正确<1 ulp，因此要么可以接受，要么返回（a）。 crlibm是正确的圆形，glibc provides a tighter guarantee是0.55 ulps，所以两者都将返回（b）。