std::fmod 函数的双精度计算精度极低

Question

fmod(1001.0, 0.0001)返回的是0.00009999999995，精度非常低（10^-5），与期望结果0相比。

根据cppreference所述，fmod()可以使用remainder()实现，但remainder(1001.0, 0.0001)返回的是-4.796965775988316e-14（仍然远低于double精度，但比10^-5好得多）。

fmod为什么会这么大程度地依赖于输入参数的精度？这正常吗？

MCVE：

#include <cmath>
#include <iomanip>
#include <iostream>
using namespace std;

int main() {
    double a = 1001.0, b = 0.0001;
    cout << setprecision(16);
    cout << "fmod:      " << fmod(a, b) << endl;
    cout << "remainder: " << remainder(a, b) << endl;
    cout << "actual:    " << a-floor(a/b)*b << endl;
    cout << "a/b:       " << a / b << endl;
}

输出：

fmod:      9.999999995203035e-05
remainder: -4.796965775988316e-14
actual:    0
a/b:       10010000

（使用GCC、Clang和MSVC，无论是否进行优化，结果相同）

演示链接

- rustyx

---

问题在于 0.0001 无法在二进制浮点数中精确地表示。 - Barmar

我已经重新打开了，但我认为如果你只是打印b并带有许多位精度，你会茅塞顿开。 - Barmar

@Barmar 我要毁掉惊喜了：假设几乎普遍采用的IEEE-754表示法，0.0001实际上是0.000100000000000000004792173602385929598312941379845142364501953125。与此同时，1001.0是精确的。 - Mark Ransom

另一种方法是使用％a来发布所有FP值以确保它们的精度，而不是使用数十个数字。 - chux - Reinstate Monica

故事的寓意似乎是不要认为你懂浮点数 :P - Passer By

如果想继续使用std::cout，你可以使用std::hexfloat修饰符（需要C++11），这样可以将输出格式转换为十六进制浮点数。请注意，不要改变原有的意思。 - chtz

Answer 1

如果我们修改你的程序为：

#include <cmath>
#include <iomanip>
#include <iostream>

int main() {
    double a = 1001.0, b = 0.0001;
    std::cout << std::setprecision(32) << std::left;
    std::cout << std::setw(16) << "a:" << a << "\n"; 
    std::cout << std::setw(16) << "b:" << b << "\n"; 
    std::cout << std::setw(16) << "fmod:" << fmod(a, b) << "\n";
    std::cout << std::setw(16) << "remainder:" << remainder(a, b) << "\n";
    std::cout << std::setw(16) << "floor a/b:" << floor(a/b) << "\n";
    std::cout << std::setw(16) << "actual:" << a-floor(a/b)*b << "\n";
    std::cout << std::setw(16) << "a/b:" << a / b << "\n";
    std::cout << std::setw(16) << "floor 10009999:" << floor(10009999.99999999952) << "\n";
}

它的输出为：

a:              1001
b:              0.00010000000000000000479217360238593
fmod:           9.9999999952030347032290447106817e-05
remainder:      -4.796965775988315527911254321225e-14
floor a/b:      10010000
actual:         0
a/b:            10010000
floor 10009999: 10010000

我们可以看到，0.0001无法表示为一个double类型，因此b实际上被设置为0.00010000000000000000479217360238593。

这导致a/b为10009999.9999999995203034224，因此意味着fmod应该返回1001 - 10009999*0.00010000000000000000479217360238593，即9.99999999520303470323e-5。

(在SpeedCrunch中计算的数字可能不完全匹配IEEE双精度值)

你的“实际”值之所以不同，是因为floor(a/b)返回的是10010000而不是fmod使用的确切值10009999，这本身是由于10009999.99999999952不能表示为双精度，因此在传递给floor之前会四舍五入为10010000。

- Alan Birtles

Answer 2

fmod函数能够准确计算，没有误差。

对于给定的C++源代码fmod(1001.0, 0.0001)，在使用IEEE-754二进制64位格式（双精度浮点型最常用的格式）的实现中，源文本中的0.0001会被转换成双精度浮点数值0.000100000000000000004792173602385929598312941379845142364501953125。

然后可以得到等式1001 = 10009999 • 0.000100000000000000004792173602385929598312941379845142364501953125 + 0.000099999999952030347032290447106817055100691504776477813720703125。因此，fmod(1001, 0.0001)的结果是精确等于0.000099999999952030347032290447106817055100691504776477813720703125。

唯一的误差发生在将源文本中的十进制数转换为基于二进制的双精度浮点数格式时，而fmod运算本身没有误差。

- Eric Postpischil

Answer 3

这里的基本问题（IEEE-754表示的0.0001）已经被充分解决了，但为了好玩，我复制了使用std::remainder实现fmod的代码，来自https://en.cppreference.com/w/cpp/numeric/math/fmod并将其与std::fmod进行了比较。

#include <iostream>
#include <iomanip>
#include <cmath>

// Possible implementation of std::fmod according to cppreference.com
double fmod2(double x, double y)
{
#pragma STDC FENV_ACCESS ON
    double result = std::remainder(std::fabs(x), (y = std::fabs(y)));
    if (std::signbit(result)) result += y;
    return std::copysign(result, x);
}

int main() {
    // your code goes here
    double b = 0.0001;
    std::cout << std::setprecision(25);
    std::cout << "              b:" << std::setw(35) << b << "\n"; 
    
    double m = 10010000.0;
    double c = m * b;
    double d = 1001.0 - m * b;
    std::cout << std::setprecision(32);
    std::cout << "     10010000*b:" << std::setw(6) << c << "\n"; 
    std::cout << std::setprecision(25);
    std::cout << "1001-10010000*b:" << std::setw(6) << d << "\n";
    
    long double m2 = 10010000.0;
    long double c2 = m2 * b;
    long double d2 = 1001.0 - m2 * b;
    std::cout << std::setprecision(32);
    std::cout << "     10010000*b:" << std::setw(35) << c2 << "\n"; 
    std::cout << std::setprecision(25);
    std::cout << "1001-10010000*b:" << std::setw(35) << d2 << "\n";
    
    std::cout << "      remainder:" << std::setw(35) << std::remainder(1001.0, b) << "\n"; 
    std::cout << "           fmod:" << std::setw(35) << std::fmod(1001.0, b) << "\n"; 
    std::cout << "          fmod2:" << std::setw(35) << fmod2(1001.0, b) << "\n"; 
    std::cout << " fmod-remainder:" << std::setw(35) <<
                 std::fmod(1001.0, b) - std::remainder(1001.0, b) << "\n"; 
    return 0;
}

结果为：

              b:     0.0001000000000000000047921736
     10010000*b:  1001
1001-10010000*b:     0
     10010000*b:  1001.0000000000000479616346638068
1001-10010000*b:    -4.796163466380676254630089e-14
      remainder:    -4.796965775988315527911254e-14
           fmod:     9.999999995203034703229045e-05
          fmod2:     9.999999995203034703229045e-05
 fmod-remainder:     0.0001000000000000000047921736

正如输出的最后两行所示，在这个实现中，实际的 std::fmod 与 cppreference 页面上建议的实现匹配，至少对于这个例子是如此。

我们还看到，64位 IEEE-754 不足以精确表示 10010000 * 0.0001 与整数之间的差异。但如果我们使用128位，则小数部分可以明确表示，当我们从 1001.0 中减去它时，我们发现余数与 std::remainder 的返回值大致相同。（差异可能是由于 std::remainder 使用少于128位的计算；它可能使用80位算术。）

最后，需要注意的是std::fmod(1001.0, b) - std::remainder(1001.0, b)的结果等于64位IEEE-754值0.0001。 也就是说，这两个函数返回的结果在模0.0001000000000000000047921736下同余， 但std::fmod选择最小正值，而std::remainder选择最接近零的值。

- David K