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;
}

输出：

#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;
}

（与 GCC、Clang、MSVC 的结果相同，有和没有优化）

现场演示

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";
}

它输出：

#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";
}

我们可以看到它0.0001不能表示为 a doublesob实际上设置为0.00010000000000000000479217360238593。

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

（在 speedcrunch 中计算的数字可能与 IEEE 双精度值不完全匹配）

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

Answer 2

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

给定fmod(1001.0, 0.0001)使用 IEEE-754 binary64（最常用的格式double）的实现中的 C++ 源代码，源文本0.0001被转换为double值 0.00010000000000000000479217360238592959831294137984514295342

然后1001 = 10009999•0.000100000000000000004792173602385929598312941379845142364501953125 + 0.000099999999952030347032290447106817055100691504776477813720703125，所以fmod(1001, 0.0001)是完全0.000099999999952030347032290447106817055100691504776477813720703125。

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

Answer 3

这里的基本问题（ 的 IEEE-754 表示0.0001）已经很好地建立起来，但只是为了踢球，我从https://en.cppreference.com/w/cpp/numeric/math/fmod复制了fmodusing的实现并将其与.std::remainderstd::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;
}

结果是：

#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;
}

如最后两行输出所示，实际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) 结果等于 的 64 位 IEEE-754 值0.0001。也就是说，两个函数都返回与相同值 modulo 一致的结果0.0001000000000000000047921736，但std::fmod选择最小的正值，同时 std::remainder选择最接近零的值。