与 matlab 相比，获得三角函数的正确值

Question

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

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

在对该问题进行大量研究后，我决定使用其他方法来计算 sin/cos 函数，而不是 glibc 中提供的方法。

我实现了这些方法，

1- 带有 long double 变量和 -O2 的泰勒级数（强制使用带有 80 位浮点运算的 x87 FPU）

2- 带有 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-openlibm

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

除了最后一个实验之外，上面的所有实验都无法生成等于 matlab 的值。我对所有这些库和方法进行了广泛的输入测试，其中一些像 libmcr 和 fdlibm 会为一些输入产生 NAN 值（看起来它们没有很好的范围检查），其余的产生值10e-16 及更高的错误。 与预期的 matlab 相比，只有最后一个会产生正确的值，但调用 matlab 函数效率不高，而且比原生实现慢得多。

我也很惊讶为什么 MPFR 和 taylor 系列的 long double 和 quadmath 会出错。

这是具有长双变量（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 完全相同的值。有什么想法吗？

更新1：

1 Ulp 错误根本不影响算法输出，但它会阻止使用 matlab 验证结果，特别是在积分器的输出中。

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

更新 2： 我计算了上述所有方法的不相等结果的数量，显然与 matlab 相比，openlibm 会产生更少的不等值，但它不是零。

Answer 1

我的猜测是 Matlab 使用的代码最初基于 FDLIBM。我能够使用 Julia（使用 openlibm）获得相同的结果：您可以尝试使用它，或者 musl，我相信它也使用相同的代码。

double/IEEE binary64 与 0.5857069572718263 最接近的是

0.5857069572718263117394599248655140399932861328125

（具有位模式0x3fe2be1c8450b590）

这个sin是

0.55278864311139114312806521962078480744570117018100444956428008387067038680572587...

最接近的两个double/IEEE binary64 是

a) 0.5527886431113910870038807843229733407497406005859375 (0x3fe1b071cef86fbe)，误差为 0.5055 ulps

b) 0.55278864311139119802618324683862738311290740966796875 (0x3fe1b071cef86fbf)，误差为 0.4945 ulps

FDLIBM 只保证对 glibc provides a tighter guarantee 为 0.55 ulps，因此两者都将返回 (b)。