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

Question

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

根据旧问题1，最近我发现在 glibc 中实现的三角函数的输出值不等于在 matlabs 中生成的值。 2 3和我的实验与glibc的1ulp>准确性相关的差异。对于大部分块来说，这个 10e-16 误差没有多大意义，但是在积分器的输出中，10e-16 积累的越来越多，积分器的最终误差约为 1e-3，这有点高，并且对于这种块是 Not Acceptable 。

在对该问题进行了大量研究之后，我决定使用除 glibc 中提供的方法之外的其他方法来计算正弦/余弦函数。

我实现了这些方法，

1- 带有 long double 变量和 -O2 的 taylor 系列(强制使用 x87 FPU 及其 80 位浮点运算)

2- 带有 GNU quadmath 库的 taylor 系列(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 函数效率不高，而且比 native 实现慢得多。

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

这是带有 long double 变量(80 位精度)的 taylor 系列
并且应该使用 -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.55278864311139108700388078432229733407497406005859375 ( 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)。