为何C++的sin()函数对同一参数返回不同值？

问题描述 #

C++的sin()函数有时对同一参数返回略有差异的结果（差异出现在第8位小数），原因是什么？

测试代码 #

#include <fstream>
#include <iostream>
#include <thread>
#include <vector>
#include <cmath>
#include <immintrin.h>

std::ofstream toFile("file.txt");
double xValue = 1.3808129283228576;

void func(int nValues) {
    std::vector<double> results(nValues);
    for (int i=0; i<nValues; i++) {
        results[i] = sin(xValue);
        uint64_t *ptrX     = (uint64_t *)&xValue;
        uint64_t *ptrValue = (uint64_t *)&results[i];
        toFile << xValue << " " << results[i] << "   " << std::hex << *ptrX << " " << *ptrValue << " " << 
                  _mm_getcsr() << "\n";
        std::this_thread::sleep_for(std::chrono::milliseconds(30));
    }
    for (int i=0; i+1<nValues; i++)
        if (fabs(results[i] - results[i+1]) > 1e-10) {
            toFile << "ERROR:\n" << results[i] << "\n" << results[i+1] << "\n";
            std::cout <<  "ERROR:\n" << results[i] << "\n" << results[i+1] << "\n";
            return;
        }
    toFile << "OK\n";
    std::cout << "OK\n";
}

int main() {
    toFile.precision(16);
    std::cout.precision(16);
    func(1000);
    return 0;
}

输出示例（file.txt） #

1.380812928322858 0.9820073703281703   3ff617cf4c10a691 3fef6c9ab87fb1c9 1fa0
1.380812928322858 0.9820073703281703   3ff617cf4c10a691 3fef6c9ab87fb1c9 1fa0
1.380812928322858 0.9820073703281703   3ff617cf4c10a691 3fef6c9ab87fb1c9 1fa0
1.380812928322858 0.982007359072273    3ff617cf4c10a691 3fef6c9ab274b21b 1fa0
1.380812928322858 0.982007359072273    3ff617cf4c10a691 3fef6c9ab274b21b 1fa0
1.380812928322858 0.982007359072273    3ff617cf4c10a691 3fef6c9ab274b21b 1fa0
1.380812928322858 0.9820073703281703   3ff617cf4c10a691 3fef6c9ab87fb1c9 1fa0
1.380812928322858 0.9820073703281703   3ff617cf4c10a691 3fef6c9ab87fb1c9 1fa0
1.380812928322858 0.9820073703281703   3ff617cf4c10a691 3fef6c9ab87fb1c9 1fa0

编译环境 #

Windows平台下Mingw的g++ 9.2及11.2版本，编译参数为-O3。

更新说明 #

已添加xValue和sin(xValue)的十六进制表示输出，以及_mm_getcsr()函数的返回结果。

原因分析 #

核心原因：指令集实现差异 #

从输出的十六进制结果和_mm_getcsr()返回值（始终为1fa0，对应MXCSR寄存器默认设置：舍入模式为最近偶数，精度控制为53位双精度）可知，输入参数没有变化，问题根源在于编译器在优化模式下对sin()函数的实现选择：

• 在-O3优化级别下，GCC会根据上下文自动切换不同指令集的实现：
  • 一种是基于x87 FPU的三角函数指令，另一种是基于SSE/AVX的SIMD指令
  • 这两种指令集的底层近似算法（如CORDIC迭代、多项式拟合）存在细微精度差异，导致最终结果在第8位小数出现偏差
• 即使浮点控制寄存器（MXCSR）的精度配置一致，不同指令集的硬件实现或编译器内置函数的近似精度也会有差别，非特殊角度的浮点数计算中，累积误差会放大这种差异

解决方案 #

• 固定浮点指令集：编译时添加-mfpmath=sse或-mfpmath=387参数，强制编译器使用指定的浮点运算指令集，确保sin()的实现一致
• 指定精确实现：使用编译器提供的高精度数学库，或明确调用__builtin_sin等内置函数，锁定三角函数的计算实现

内容的提问来源于stack exchange，提问作者Borys L.