std::fmod 糟糕的双精度

Question

fmod(1001.0, 0.0001) 给出了0.00009999999995，考虑到0 的预期结果，这似乎是一个非常低的精度（10^-5）。

根据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 的结果相同，有和没有优化）

Live demo

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 不能表示为双精度所以它在传递到地板之前四舍五入到10010000。

Answer 2

这里的基本问题（0.0001 的 IEEE-754 表示）已经很成熟，但只是为了好玩，我使用来自 https://en.cppreference.com/w/cpp/numeric/math/fmod 的 std::remainder 复制了 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 页面上建议的实现相匹配，至少在此示例中如此。

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

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

Answer 3

fmod 产生准确的结果，没有错误。

给出了C ++源代码@ 987654322使用IEEE-754 Binary64（最常用的格式）（double），源文本0.0001 double value 0.0001000000000000 0003839320385929531238592953123859295329364501953123941379531230385953123941379531239413795312323645019531238595312323645953123236439531232930323 >

然后1001 = 10009999•0.000100000000000000004792173602385929598312941379845142364501953125 + 0.000099999999952030347032290447106817055100691504776477813720703125，所以fmod(1001, 0.0001)正是0.000099999999952030347032290447106817055100691504776477813720703125。 P>

将源文本中的十进制数字转换为基于二进制的double 格式时会出现唯一错误。 fmod操作没有错误。